Quantum microscopic dynamical approaches

The quantum state of a single-particle is noted $|\varphi\rangle$. It belongs to the Hilbert space of one particle, noted $\mathcal{H}_{1}$. In the framework of second quantisation, such a state is written $|\varphi\rangle=\hat{a}^{\dagger}_{\varphi}|-\rangle$, where $|-\rangle$ is the state of the vacuum and $\hat{a}^{\dagger}_{\varphi}$ creates a particle in the state $|\varphi\rangle$.

When two particles are distinguishable, such as a proton and an electron in the Hydrogen atom, or a proton and a neutron in a deuteron, the state of each particle can be labelled, e.g., with a number. Two cases are considered: one where the particles are independent and one where they are correlated.

Now consider the case of two identical fermions which are indistinguishable. In this case, the states $|1:\varphi_{1},2:\varphi_{2}\rangle$ and $|1:\varphi_{2},2:\varphi_{1}\rangle$ are equivalent and must both appear with the same probability in the two-body state. For two independent fermions,

$$|\Phi\rangle=|\varphi_{1}\varphi_{2}\rangle=(|1:\varphi_{1},2:\varphi_{2}\rangle-|1:\varphi_{2},2:\varphi_{1}\rangle)/\sqrt{2}.$$

The minus sign accounts for the fact that the state must be antisymmetric for fermions (there would be a plus sign for bosons). The antisymmetry is obvious as $|\varphi_{1}\varphi_{2}\rangle=-|\varphi_{2}\varphi_{1}\rangle$. In particular, $|\varphi\varphi\rangle=0$ which is a manifestation of the Pauli principle.

The state $|\Phi\rangle$ can be written in the form of a Slater determinant

$$|\Phi\rangle=\frac{1}{\sqrt{2}}\,\begin{vmatrix}|1:\varphi_{1}\rangle&|1:\varphi_{2}\rangle\\ |2:\varphi_{1}\rangle&|2:\varphi_{2}\rangle\\ \end{vmatrix}=\mathcal{A}|1:\varphi_{1},2:\varphi_{2}\rangle,$$

where $\mathcal{A}$ is the antisymmetrisation operator.

In second quantisation, $|\Phi\rangle=\hat{a}^{\dagger}_{2}\hat{a}^{\dagger}_{1}|-\rangle$, with the antisymmetry ensured by the anticommutation relationship $\{{\hat{a}^{\dagger}}_{1},{\hat{a}^{\dagger}}_{2}\}=0$. The generalisation to correlated two-particle states becomes $|\Psi\rangle=\sum_{\alpha}C_{\alpha}\hat{a}^{\dagger}_{\alpha_{2}}\hat{a}^{\dagger}_{\alpha_{1}}|-\rangle$. Both $|\Phi\rangle$ and $|\Psi\rangle$ belong to the Hilbert space of two identical fermions $\mathcal{H}_{2}$.

The generalisation to the case of $N$ identical fermions is now straightforward. For independent particles, the many-body state is described by a Slater determinant

$$|\Phi\rangle=\frac{1}{\sqrt{N!}}\,\begin{vmatrix}|1:\varphi_{1}\rangle&\cdots&|1:\varphi_{N}\rangle\\ \vdots&&\vdots\\ |N:\varphi_{1}\rangle&\cdots&|N:\varphi_{N}\rangle\\ \end{vmatrix}=\mathcal{A}|1:\varphi_{1},\cdots,N:\varphi_{N}\rangle=\hat{a}^{\dagger}_{N}\cdots\hat{a}^{\dagger}_{1}|-\rangle.$$

As before, a correlated state is written as a sum of independent particle states $|\Psi\rangle=\sum_{\alpha}C_{\alpha}|\Phi_{\alpha}\rangle$. As a result, a complete set of orthonormal Slater determinants constitute a possible basis of the Hilbert space of $N$ particles $\mathcal{H}_{N}$.

There are several variational principles that could be used. For the moment, consider a variational principle based on Dirac’s action

$$S=\int_{0}^{T}dt\,\langle\Psi|\left(i\hbar\frac{d}{dt}-\hat{H}\right)|\Psi\rangle.$$

Requiring the stationarity of the action, $\delta S=0$, leads to the Schrödinger equation.

Remember that, for a complex variable $z={\mathrm{Re}}(z)+i\,{\mathrm{Im}}(z)$, there are two independent variables, ${\mathrm{Re}}(z)$ and ${\mathrm{Im}}(z)$, or, equivalently, $z={\mathrm{Re}}(z)+i\,{\mathrm{Im}}(z)$ and $z^{*}={\mathrm{Re}}(z)-i\,{\mathrm{Im}}(z)$. Similarly, $|\Psi\rangle$ and $\langle\Psi|$ can be treated as independent variational quantities and the property $\langle\Psi|=(|\Psi\rangle)^{\dagger}$ can be restored at the end. The stationarity condition is then expressed as

$$\,\,\,\frac{\delta S}{\delta\langle\Psi|}=0\,\,\,\mbox{ and }\frac{\delta S}{\delta|\Psi\rangle}=0\,\,\,\mbox{ for }\,\,\,0\leq t\leq T.$$

The first condition gives

	$\displaystyle\frac{\delta}{\delta\langle\Psi(t)|}\int_{0}^{T}dt^{\prime}\,\langle\Psi(t^{\prime})|\left(i\hbar\frac{d}{dt}-\hat{H}\right)|\Psi(t^{\prime})\rangle$	$\displaystyle=$	$\displaystyle 0$
	$\displaystyle\int_{0}^{T}dt^{\prime}\,\frac{\delta\langle\Psi(t^{\prime})|}{\delta\langle\Psi(t)|}\left(i\hbar\frac{d}{dt}-\hat{H}\right)|\Psi(t^{\prime})\rangle$	$\displaystyle=$	$\displaystyle 0$
	$\displaystyle\left(i\hbar\frac{d}{dt}-\hat{H}\right)|\Psi\rangle$	$\displaystyle=$	$\displaystyle 0$

where $\frac{\delta\langle\Psi(t^{\prime})|}{\delta\langle\Psi(t)|}=\delta(t-t^{\prime})$ was used. The last line is the Schrödinger equation.

The second condition can be expressed as

$$\frac{\delta}{\delta|\Psi(t)\rangle}\int_{0}^{T}dt^{\prime}\,\langle\Psi(t^{\prime})|\left(i\hbar\frac{d}{dt}-\hat{H}\right)|\Psi(t^{\prime})\rangle=0.$$

One has to be careful with the $\frac{d}{dt}$ term. Using integration by part,

	$\displaystyle\frac{\delta}{\delta|\Psi(t)\rangle}\int_{0}^{T}dt^{\prime}\,\langle\Psi(t^{\prime})|\left(\frac{d}{dt}|\Psi(t^{\prime})\rangle\right)$	$\displaystyle=$
	$\displaystyle\frac{\delta}{\delta|\Psi(t)\rangle}\left[\left[\langle\Psi(t^{\prime})|\Psi(t^{\prime})\rangle\right]_{0}^{T}-\int_{0}^{T}dt^{\prime}\,\left(\frac{d}{dt}\langle\Psi(t^{\prime})|\right)|\Psi(t^{\prime})\rangle\right]$	$\displaystyle=$
	$\displaystyle\left[\langle\Psi(t^{\prime})|\frac{\delta|\Psi(t^{\prime})\rangle}{\delta|\Psi(t)\rangle}\right]_{0}^{T}-\int_{0}^{T}dt^{\prime}\,\left(\frac{d}{dt}\langle\Psi(t^{\prime})|\right)\frac{\delta|\Psi(t^{\prime})\rangle}{\delta|\Psi(t)\rangle}$	$\displaystyle=$
	$\displaystyle\left[\langle\Psi(t^{\prime})|\frac{\delta|\Psi(t^{\prime})\rangle}{\delta|\Psi(t)\rangle}\right]_{0}^{T}-\frac{d}{dt}\langle\Psi(t)|.$				(2)

(Note that $\frac{\delta\langle\Psi(t^{\prime})|}{\delta\langle\Psi(t)|}=\delta(t-t^{\prime})$ in the $[\cdots]_{0}^{T}$ term was not used intentionally.)

The term with $\hat{H}$ is trivial if the latter does not depend on $|\Psi\rangle$:

$$\frac{\delta}{\delta|\Psi(t)\rangle}\int_{0}^{T}dt^{\prime}\,\langle\Psi(t^{\prime})|\hat{H}|\Psi(t^{\prime})\rangle=\langle\Psi(t)|\hat{H}.$$

(3)

Combining Eqs. (2) and (3),

$$i\hbar\left(\langle\Psi(T)|\frac{\delta|\Psi(T)\rangle}{\delta|\Psi(t)\rangle}-\langle\Psi(0)|\frac{\delta|\Psi(0)\rangle}{\delta|\Psi(t)\rangle}-\frac{d}{dt}\langle\Psi(t)|\right)-\langle\Psi(t)|\hat{H}=0.$$

(4)

Now “restore” the property $\langle\Psi|=(|\Psi\rangle)^{\dagger}$. As shown earlier, imposing $\frac{\delta S}{\delta\langle\Psi|}=0$ leads to the Schrödinger equation. Taking its Hermitian conjugate, $-i\hbar\frac{d}{dt}\langle\Psi|=\langle\Psi|\hat{H}$. Comparing with Eq. (4), it is seen that, to have $\frac{\delta S}{\delta|\Psi\rangle}=0$, one needs to have

$$\langle\Psi(T)|\frac{\delta|\Psi(T)\rangle}{\delta|\Psi(t)\rangle}-\langle\Psi(0)|\frac{\delta|\Psi(0)\rangle}{\delta|\Psi(t)\rangle}=0.$$

This is obtained by forbidding variations at the initial and final times, i.e., $\delta|\Psi(0)\rangle=\delta|\Psi(T)\rangle=0$.

To summarise, defining the Dirac action $S=\int_{0}^{T}dt\,\langle\Psi|\left(i\hbar\frac{d}{dt}-\hat{H}\right)|\Psi\rangle$ and requiring $\frac{\delta S}{\delta\langle\Psi|}=0$ leads to Schrödinger equation, while $\frac{\delta S}{\delta|\Psi\rangle}=0$ gives its Hermitian conjugate if it is required that $\Psi$ cannot vary at $t=0$ and $t=T$.

To get Schrödinger from the variational principle $\delta S=0$, the variational space for $\Psi$ was not restricted, i.e., it spans the entire Hilbert space $\mathcal{H}_{N}$. If the variational space was restricted to a sub-space $\mathcal{F}_{N}$ of $\mathcal{H}_{N}$, one would expect to get a solution to the variational principle $|\Psi(t)\rangle$ that is not solution to the Schrödinger equation, but which is an optimised approximation of the exact evolution within $\mathcal{F}_{N}$. The role of a theorist is then to choose a sub-space $\mathcal{F}_{N}$ which (i) contains enough physics to give a good approximation to the evolution of the system of interest, and (ii) is simple enough so that the resulting equations of motion can be solved analytically or numerically.

The Dirac action is

$$S_{t_{0},t_{1}}[\Phi]=\int_{t_{0}}^{t_{1}}dt\,\langle\Phi(t)|\left(i\hbar\frac{d}{dt}-\hat{H}\right)|\Phi(t)\rangle$$

where $|\Phi\rangle\in\mathcal{F}_{N}$ is a Slater determinant built from the $A$ occupied single-particle states $|\varphi_{i}\rangle$, $1\leq i\leq A$.

The expectation value of the time derivative is computed first. From

$$\frac{d}{dt}|\Phi\rangle=\frac{d}{dt}\mathcal{A}|\varphi_{A}\rangle\cdots|\varphi_{1}\rangle=\mathcal{A}\sum_{j=1}^{A}|\varphi_{A}\rangle\cdots\left(\frac{d}{dt}|\varphi_{j}\rangle\right)\cdots|\varphi_{1}\rangle$$

one gets

	$\displaystyle\langle\Phi|\frac{d}{dt}|\Phi\rangle$	$\displaystyle=$	$\displaystyle\left[\langle\varphi_{1}|\cdots\langle\varphi_{A}|\mathcal{A}\right]\mathcal{A}\sum_{j=1}^{A}|\varphi_{A}\rangle\cdots\left(\frac{d}{dt}|\varphi_{j}\rangle\right)\cdots|\varphi_{1}\rangle$
		$\displaystyle=$	$\displaystyle\sum_{j=1}^{A}\langle\varphi_{1}|\varphi_{1}\rangle\cdots\left(\langle\varphi_{j}|\frac{d}{dt}|\varphi_{j}\rangle\right)\cdots\langle\varphi_{A}|\varphi_{A}\rangle$
		$\displaystyle=$	$\displaystyle\sum_{j=1}^{A}\langle\varphi_{j}|\frac{d}{dt}|\varphi_{j}\rangle.$

Now compute the expectation value of the Hamiltonian and write it as an energy density functional (EDF) $E[\rho]$:

$$\langle\Phi|\hat{H}|\Phi\rangle=\langle\varphi_{1}\cdots\varphi_{A}|\hat{H}|\varphi_{A}\cdots\varphi_{1}\rangle=E[\varphi_{1}\cdots\varphi_{A}]\equiv E[\rho],$$

where $\rho$ is the one-body density matrix. Its elements are defined as

$$\rho_{\alpha\beta}=\langle\hat{a}_{\beta}^{\dagger}\hat{a}_{\alpha}\rangle.$$

(5)

For a Slater,

$$\rho_{\alpha\beta}=\langle\varphi_{1}\cdots\varphi_{A}|\hat{a}^{\dagger}_{\beta}\hat{a}_{\alpha}|\varphi_{A}\cdots\varphi_{1}\rangle=\sum_{j=1}^{A}\langle\varphi_{j}|\beta\rangle\langle\alpha|\varphi_{j}\rangle.$$

(6)

For instance, in the position-spin-isospin basis of $\mathcal{H}_{1}$, $\{|\mathbf{r},s,q\rangle\}$, one has

$$\rho(\mathbf{r}sq,\mathbf{r}^{\prime}s^{\prime}q^{\prime})=\sum_{j=1}^{A}\varphi_{j}^{*}(\mathbf{r}^{\prime}s^{\prime}q^{\prime})\varphi_{j}(\mathbf{r}sq).$$

The one-body density matrix is expressed as a function of all the occupied single-particle states. It contains the same information on the system as the Slater. This is why the functional $E[\varphi_{1}\cdots\varphi_{A}]$ can be replaced by the energy density functional $E[\rho]$. The EDF accounts for the kinetic energy and the interaction between the particles. The most popular EDF are the Skyrme and Gogny functionals [see (Bender et al 2003) for a review].

The action is then written as

	$\displaystyle S$	$\displaystyle=$	$\displaystyle\int_{t_{0}}^{t_{1}}dt\,\left(i\hbar\sum_{i=1}^{N}\langle\varphi_{i}|\frac{d}{dt}|{\varphi}_{i}\rangle-E[\rho(t)]\right)$
		$\displaystyle=$	$\displaystyle\int_{t_{0}}^{t_{1}}dt\,\left(i\hbar\sum_{i=1}^{N}\,\int dx\,\varphi_{i}^{*}(x,t)\,\frac{d}{dt}{\varphi}_{i}(x,t)-E[\rho(t)]\right)$

where the notations $x\equiv(\mathbf{r}sq)$ and $\int dx=\sum_{qs}\int d\mathbf{r}$ are introduced.

The variational principle $\delta S=0$ is solved by considering the variation of all the independent variables $\varphi_{i}$ and $\varphi^{*}_{j}$, $1\leq j\leq A$:

$\displaystyle\frac{\delta S}{\delta\varphi_{j}(x,t)}=0\mbox{\hskip 28.45274pt and \hskip 28.45274pt}\frac{\delta S}{\delta\varphi^{*}_{j}(x,t)}=0.$

(8)

The variation over $\varphi^{*}$ gives

$$\frac{\delta S}{\delta\varphi^{*}_{j}(x,t)}=i\hbar\,\frac{d}{dt}\varphi_{j}(x,t)-\int_{t_{0}}^{t_{1}}dt^{\prime}\,\frac{\delta E[\rho(t^{\prime})]}{\delta\varphi^{*}_{j}(x,t)}.$$

(9)

The functional derivative of $E$ can be re-written thanks to a change of variable

$$\frac{\delta E[\rho(t^{\prime})]}{\delta\varphi^{*}_{j}(x,t)}=\int dy\,dy^{\prime}\,\frac{\delta E[\rho(t^{\prime})]}{\delta\rho(y,y^{\prime};t^{\prime})}\,\frac{\delta\rho(y,y^{\prime};t^{\prime})}{\delta\varphi^{*}_{j}(x,t)}.$$

(10)

Using

$$\frac{\delta\rho(y,y^{\prime};t^{\prime})}{\delta\varphi^{*}_{j}(x,t)}=\varphi_{j}(y,t^{\prime})\,\delta(y^{\prime}-x)\,\delta(t-t^{\prime})$$

(11)

and noting the single-particle Hartree-Fock Hamiltonian $h$ with matrix elements

$$h(x,y;t)=\frac{\delta E[\rho(t)]}{\delta\rho(y,x;t)},$$

(12)

one gets the TDHF equation for the set of occupied states

$$\framebox{$\displaystyle i\hbar\,\frac{d}{dt}\varphi_{j}(x,t)=\int dy\,h(x,y;t)\,\varphi_{j}(y,t)$}.$$

(13)

As in the Schrödinger case, the variation over $\varphi$ leads to the complex conjugate of the TDHF equation.

The TDHF equation can also be written in a matrix form

$$i\hbar\frac{d}{dt}\varphi_{j}(t)=h(t)\varphi_{j}(t)\,\,,\,\,\,\,\,1\leq j\leq A,$$

or, using Dirac notation,

$$i\hbar\frac{d}{dt}|\varphi_{j}(t)\rangle=\hat{h}(t)|\varphi_{j}(t)\rangle.$$

(14)

Note that the single-particle Hamiltonian $h$ is self-consistent, i.e., it depends on the density matrix $\rho$. It is then also a function of time. As a result, the TDHF equation is equivalent to a set of non-linear single-particle Schrödinger equations. The non-linearity comes from the self-consistency of $h[\rho]$: The Hamiltonian is a function of the states on which it acts.

The one-body density matrix can be written as an operator of $\mathcal{H}_{1}$:

$$\hat{\rho}=\sum_{j=1}^{A}|\varphi_{j}\rangle\langle\varphi_{j}|.$$

Indeed, its matrix elements

$$\rho_{\alpha\beta}=\langle\alpha|\hat{\rho}|\beta\rangle=\sum_{j=1}^{A}\langle\varphi_{j}|\beta\rangle\langle\alpha|\varphi_{j}\rangle$$

are the same as for a Slater determinant in Eq. (6).

Multiplying Eq. (14) by $\langle\varphi_{j}|$ on the right and summing over $j$,

$$\sum_{j=1}^{A}i\hbar\left(\frac{d}{dt}|\varphi_{j}\rangle\right)\langle\varphi_{j}|=\sum_{j=1}^{A}\hat{h}|\varphi_{j}\rangle\langle\varphi_{j}|.$$

(15)

Similarly, multiplying the Hermitian conjugate of Eq. (14) by $|\varphi_{j}\rangle$ on the left and summing over $j$ leads to

$$\sum_{j=1}^{A}-i\hbar|\varphi_{j}\rangle\left(\frac{d}{dt}\langle\varphi_{j}|\right)=\sum_{j=1}^{A}|\varphi_{j}\rangle\langle\varphi_{j}|\hat{h}.$$

(16)

Taking the difference between Eqs. (15) and (16) gives the Liouville form of the TDHF equation

	$\displaystyle\sum_{j=1}^{A}i\hbar\frac{d}{dt}\left(|\varphi_{j}\rangle\langle\varphi_{j}|\right)$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{A}\left(\hat{h}|\varphi_{j}\rangle\langle\varphi_{j}|-|\varphi_{j}\rangle\langle\varphi_{j}|\hat{h}\right)$
	$\displaystyle i\hbar\frac{d}{dt}\hat{\rho}$	$\displaystyle=$	$\displaystyle\left[\hat{h},\hat{\rho}\right].$		(17)

This can be written also in the matrix form

$$i\hbar\frac{d}{dt}{\rho}=\left[{h},{\rho}\right].$$

(18)

The easiest way to solve the TDHF equation numerically is to start from the non-linear Schrödinger equation for occupied single-particles wave-functions

$$i\hbar\frac{d}{dt}\varphi_{j}(\mathbf{r}sq,t)=\sum_{s^{\prime}q^{\prime}}\int d^{3}r^{\prime}\,h(\mathbf{r}sq,\mathbf{r}^{\prime}s^{\prime}q^{\prime};t)\,\varphi_{j}(\mathbf{r}^{\prime}s^{\prime}q^{\prime},t)\,,\,\,\,\,\,1\leq j\leq A.$$

Because $h\equiv h[\rho(t)]$ is time-dependent, one cannot use $e^{-ih(t_{1}-t_{0})/\hbar}$ as an evolution operator, unless one considers an evolution over a time interval that is short enough so that the variation of $h(t)$ during this time interval can be neglected, i.e., $\varphi_{j}(t+\Delta t)\simeq e^{-ih\Delta t/\hbar}\varphi_{j}(t)$.

To ensure that the evolution is unitary (required for energy and particle number conservations), one needs the same operator (up to Hermitian conjugation) to go from $t\rightarrow t+\Delta t$ as the one to do the time-reversed operation $t+\Delta t\rightarrow t$. Therefore, $h$ has to be estimated at $t+\frac{\Delta t}{2}$.

A possible algorithm is

$$\begin{array}[]{ccccc}\{\varphi_{j}(t)\}&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Rightarrow&\rho(t)&\Rightarrow&h(t)\\ \Uparrow&&&&\Downarrow\\ {\varphi}_{i}(t+\Delta t)=e^{-i\frac{\Delta t}{\hbar}h(t+\frac{\Delta t}{2})}{\varphi}_{i}(t)&&&&\!\!\!\!\!\!\!\!\!\!\tilde{\varphi}_{i}(t+\Delta t)=e^{-i\frac{\Delta t}{\hbar}h(t)}{\varphi}_{i}(t)\\ \Uparrow&&&&\Downarrow\\ h(t+\frac{\Delta t}{2})&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Leftarrow&{\rho}(t+\frac{\Delta t}{2})=\frac{\tilde{\rho}(t+\Delta t)+{\rho}(t)}{2}&\Leftarrow&\tilde{\rho}(t+\Delta t)\end{array}$$

(19)

The initial state $\{\varphi_{j}(t_{0})\}$ could be a Hartree-Fock (HF) ground-stale which is put into motion thanks to an external one-body time dependent potential added to the mean-field. For collisions, one usually starts with two non-overlapping HF ground-states at some finite distance in a large box and a Galilean boost $e^{i\mathbf{k}_{n}\cdot\mathbf{r}}$ ($n=1,2$ denotes the two nuclei) is applied to induce a momentum $\hbar\mathbf{k}_{n}$ to the nucleons at the initial time $t_{0}$.

Typical TDHF solvers use cartesian grids with mesh spacing $\Delta x=0.6-1.0$ fm. Several types of boundary conditions can be used:

• •

  Hard (particles are reflected by an infinite potential on the edges of the box)

• •

  Periodic (particles reaching one edge reappear on the other side)

• •

  Absorbing (particles are absorbed on the edge)

Absorbing boundary conditions can be obtained with an additional layer of imaginary potential outside the box, or by “twisting” the phase on the boundaries (Schuetrumpf et al 2016).

Static HF ground-states are required to build the TDHF initial condition. The static HF equation reads

$$\left[h[\rho],\rho\right]=0.$$

As they commute, one can find a single-particle states basis that diagonalises both $\rho$ and $h$ simultaneously. (Note that to solve TDHF, a basis needs to be chosen. TDHF equations are usually solved in the “canonical” basis in which $\rho$ is diagonal. In this basis, $h$ is not diagonal if $h$ and $\rho$ do not commute, i.e., when the state evolves in time: $[h,\rho]=i\hbar\dot{\rho}$.)

In principle the HF ground-states could be determined from any HF code as long as the same energy density functional is used as in the TDHF evolution. However, to minimise numerical errors, it is best to use a HF code that uses the same numerical approximations as in TDHF, in particular in term of mesh grid and spatial derivatives.

One way to do this is by solving the HF equation with the imaginary time method, i.e., replacing $t\rightarrow-i\tau$ in the TDHF equation. This makes an initial state (usually built from Nilsson or harmonic oscillator wave-functions) converge towards the ground-state.

Of course, in the Hartree-Fock theory, the Hamiltonian is self-consistent so an iterative process with imaginary time step $\Delta\tau$ is needed. In addition, $e^{-h\tau/\hbar}$ is not unitary so an orthonormalisation of the occupied single-particle wave-functions is required, e.g., with a Gram-Schmidt procedure.

Here, some applications are only briefly mentioned. See Reviews (Simenel 2012; Lacroix and Ayik 2014; Nakatsukasa et al 2016, Simenel and Umar 2018; Sekizawa 2019; Stevenson and Barton 2019) for more detailed examples.

TDHF calculations are now standard to investigate various aspects of nuclear dynamics. Collective vibrations can be studied by computing the time evolution of multipole moments of a nucleus following a collective boost excitation (see next section on RPA). Although most applications have been dedicated to giant resonances, low-lying collective vibrations have been also studied.

Fusion barriers can be obtained in a straightforward manner by searching for fusion thresholds in head-on collisions. Below this threshold, the two fragments re-separate after contact. Once this threshold is overcome, the fragments merge in a single one. Compared to barrier heights of the “bare” nucleus-nucleus potential, this fusion threshold is often found at a lower energy due to dynamical effects such as vibration and transfer in the entrance channel. Applying the same method for searching for fusion thresholds at finite impact parameters allow the computation of above barrier fusion cross-sections.

Even if the nuclei do not fuse after contact, part of the wave-function can be transferred from one fragment to the other in quasi-elastic collisions. The final fragments being a coherent superposition of eigenstates of the particle number operator, one obtains a distribution of probability to find a given number of nucleons in the final fragments. Naturally, both fragments are entangled so that “measuring” the number of protons or neutrons of one fragment projects the other fragment into a state with the corresponding number of nucleons as the total number of particles is conserved in TDHF.

Significant efforts have been dedicated in the past decade to study heavy-ion collision at energies well above the barrier, leading to deep inelastic collisions (DIC), as well as in heavy systems leading to quasi-fission where fission-like fragments are produced without the intermediate formation of an equilibrated compound nucleus. These reactions offer a wide variety of observables to compare with experiment, such as fragment mass and charge distributions, scattering angle and final kinetic energy. Comparison with theory offers valuable information on expected contact times between the fragments, as well as equilibration and dissipation mechanisms.

The TDHF theory is based on the mean-field approximation. It describes the evolution of independent particles in a mean-field potential produced by the ensemble of particles. As a result, some correlations are neglected. This is the case of pairing correlations that lead to the formation of Cooper pairs and induce a superfluid behaviour to some nucleons in the nucleus. These correlations are generated by the pairing residual interaction that is neglected in the Hartree-Fock approximation.

Another limitation is the lack of quantum tunnelling of the many-body wave-function. This is due to the restriction of the variational space to a single Slater determinant. In reality, colliding nuclei have a non-zero probability to fuse even at sub-barrier energies thanks to tunnelling. However, to describe a coherent superposition of a fused system and a system of two outgoing fragments, one would need at least two Slaters, one for the compound nucleus and for the outgoing fragments. Standard applications of TDHF are then unable to describe sub-barrier fusion or spontaneous fission. There are however, techniques such as Density-Constrained TDHF that allow to extract a nucleus-nucleus potential from TDHF simulations of two colliding nuclei (Umar and Oberacker 2016). Such potential can then be used to evaluate tunnelling probabilities with standard techniques, and then predict fusion cross-sections even below the barrier. Alternatively, tunnelling can be studied at the mean-field level through a Wick rotation changing real time into imaginary time (Levit 1980). However, this approach has not reached the level of realistic applications yet.

It is also well known that TDHF lacks of quantum fluctuations. This is evident from comparisons between theoretical predictions and experimental measurements of fragment particle number distributions in deep-inelastic collisions in which TDHF underpredicts widths of these distributions. This is again a limitation due to the restriction to a single Slater determinant. Indeed, in TDHF, one Slater is used to describe all exit channels. While the mean-field dynamics of this Slater is expected to give a reasonable description of the most likely outcome, it is unlikely to be realistic for exit channels that deviate significantly from the average one.

It may seem strange that an independent particle approximation works for the nucleus which is a very dense object of nucleons very close to each other and interacting via the strong and Coulomb interactions. What makes it work is the Pauli exclusion principle. The latter prevents the collision of two nucleons to happen inside the nucleus if the final state of the collision is already occupied.

Consider two nucleons with an energy $E_{1}$ and $E_{2}$ before they interact and an energy $E_{1}^{\prime}$ and $E_{2}^{\prime}$ after the collision. Because of energy conservation, and assuming the state of the other nucleons is not changed, one has $E_{1}+E_{2}=E_{1}^{\prime}+E_{2}^{\prime}$. In the ground state, or a low excitation energy state, the nucleons have initially an energy lower than the Fermi energy $E_{F}$ below which single-particle states are essentially occupied and above which they are mostly empty. This means that $E_{1}<E_{F}$ and $E_{2}<E_{F}$. Thus, even if one final state has $E_{i}^{\prime}>E_{F}$, the other must have an energy lower than $E_{F}$ and then one nucleon ends up in a state which is likely to be already occupied, which is forbidden by Pauli. This “Pauli blocking” increases the mean-free path of a nucleon to about the size of the nucleus, reducing the effect of collisions and thus allowing the mean-field approximation to work.

The RPA is used to study small amplitude oscillations of many-body systems. It is a “harmonic” approximation as in the small amplitude limit (i.e. the system is only allowed small deformations around the ground-state) the potential energy varies quadratically, $V\propto(Q-Q_{g.s.})^{2}$, with respect to the deformation $Q$ (e.g., multipole moment). In the harmonic picture, the relationship between the frequency $\frac{\omega}{2\pi}$ of the oscillation and the energy spectrum of the harmonic oscillator $E_{n}=(n+\frac{1}{2})\hbar\omega$ can then be used.

In addition to the energy $\hbar\omega$ of the vibration, one is also interested in its collectivity quantified by the transition amplitude $q_{\nu}=\langle\nu|\hat{Q}|0\rangle$ between the ground-state $|0\rangle$ and the first phonon $|\nu\rangle$ of the vibrational spectrum. $\hat{Q}=\sum_{i=1}^{A}\hat{q}(i)$ is a one-body operator, often chosen as a multipole moment

$$\hat{Q}_{LM}=\sum_{i=1}^{A}\hat{r}^{L+2\delta_{L,0}}Y_{LM}(\hat{\theta},\hat{\phi}).$$

In the expression for the transition amplitude $q_{\nu}=\langle\nu|\sum_{i=1}^{A}\hat{q}(i)|0\rangle$, the contributions of each nucleon to the many-body state $|\nu\rangle$ are summed over. Then, the larger the magnitude of $q_{\nu}$, the larger the collectivity. Both $\hbar\omega_{\nu}$ and $q_{\nu}$ can be extracted from TDHF using the linear response theory.

The time evolution of an observable $Q(t)$ is computed from the time-dependent state $|\Psi(t)\rangle$ after an excitation induced by a small boost on the ground-state $|0\rangle$:

$$|\Psi(0)\rangle=e^{-i\epsilon\hat{Q}}|0\rangle=|0\rangle-i\epsilon\hat{Q}|0\rangle+O(\epsilon^{2}).$$

(20)

Inserting $\hat{1}=\sum_{\nu}|\nu\rangle\langle\nu|$ where $\{|\nu\rangle\}$ are the eigenstates of $\hat{H}$ with eigenenergies $\{E_{\nu}\}$,

$$|\Psi(0)\rangle=|0\rangle-i\epsilon\sum_{\nu}q_{\nu}|\nu\rangle+O(\epsilon^{2}).$$

The time-evolution of the state reads

	$\displaystyle|\Psi(t)\rangle$	$\displaystyle=$	$\displaystyle e^{-i\hat{H}t/\hbar}|\Psi(0)\rangle$
		$\displaystyle=$	$\displaystyle e^{-iE_{0}t/\hbar}\left(|0\rangle-i\epsilon\sum_{\nu}q_{\nu}e^{-i\omega_{\nu}t}|\nu\rangle\right)+O(\epsilon^{2}),$

with $\hbar\omega_{\nu}=E_{\nu}-E_{0}$. The response to this excitation can be written

	$\displaystyle Q(t)$	$\displaystyle=$	$\displaystyle\langle\Psi(t)|\hat{Q}|\Psi(t)\rangle-\langle 0|\hat{Q}|0\rangle$		(22)
		$\displaystyle=$	$\displaystyle\left(\langle 0|+i\epsilon\sum_{\nu}q_{\nu}^{*}e^{i\omega_{\nu}t}\langle\nu|\right)\hat{Q}\left(|0\rangle-i\epsilon\sum_{\nu}q_{\nu}e^{-i\omega_{\nu}t}|\nu\rangle\right)-\langle 0|\hat{Q}|0\rangle+O(\epsilon^{2})$
		$\displaystyle=$	$\displaystyle-2\epsilon\sum_{\nu}|q_{\nu}|^{2}\sin\omega t+O(\epsilon^{2}).$

The energies $\hbar\omega_{\nu}$ and transition amplitudes $q_{\nu}$ can be extracted from the strength function

$$R_{Q}(\omega)=\lim_{\epsilon\rightarrow 0}\frac{-1}{\pi\epsilon}\,\int_{0}^{\infty}dt\,{Q}(t)\,\sin(\omega t).$$

(23)

Indeed, using the expression (22) for $Q(t)$,

	$\displaystyle R_{Q}(\omega)$	$\displaystyle=$	$\displaystyle\frac{2}{\pi}\sum_{\nu}\,|q_{\nu}|^{2}\int_{0}^{\infty}dt\,\sin(\omega t)\sin(\omega_{\nu}t)$		(24)
		$\displaystyle=$	$\displaystyle\sum_{\nu}\,|q_{\nu}|^{2}\delta(\omega-\omega_{\nu}).$		(25)

The evolution of $Q(t)$ can be computed directly from TDHF codes. The boost of Eq. (20) can be applied directly to the single-particle wave-functions according to

$$|\Psi(0)\rangle=e^{-i\epsilon\hat{Q}}|0\rangle=e^{-i\epsilon\sum_{i=1}^{A}\hat{q}(i)}|\varphi_{1}\cdots\varphi_{A}\rangle=|\left(e^{-i\epsilon{q}}\varphi_{1}\right)\cdots\left(e^{-i\epsilon{q}}\varphi_{A}\right)\rangle.$$

To ensure that the response is in the linear regime, various values are usually used to check that $\langle\hat{Q}(t)\rangle\propto\varepsilon$.

An example of octupole response $\langle\hat{Q}(t)\rangle$ in the linear regime is shown in Fig. 1 in ⁴⁰Ca and ⁵⁶Ni. The computation of the strength function can be done with Eq. (23). However, the computational time being finite, the integral cannot be performed up to $t\rightarrow\infty$. Performing the integration up to $T$ without modifying Eq. (23) would induce spurious oscillations in the spectrum as can be expected from the Fourier transform of a signal convoluted with a step function. To avoid such spurious oscillations, one usually multiply $Q(t)$ by a damping function $f(t)$ (e.g., a $\cos^{2}$ or a Gaussian function) decreasing slowly from $f(0)=1$ to $f(T)=0$. The damping function induces a small width proportional to $1/T$ in the strength function.

The strength functions associated with the octupole responses in ⁴⁰Ca and ⁵⁶Ni are shown in Fig. 1. The rapid oscillation of the octupole moment in ⁵⁶Ni is associated with a higher energy $\hbar\omega\sim 9$ MeV than in ⁴⁰Ca where the main low-lying octupole vibration is found at $\sim 4$ MeV. Other states with smaller strengths are also observed at higher energy. As these nuclei have ground-state spin and parity $0^{+}$, and as the octupole operator allows for a change of angular momentum $\Delta L=3\hbar$, these vibrational states have a spin-parity $3^{-}$. The area of the peaks correspond to the transition probabilities $|q_{\nu}|^{2}$.

The transition amplitudes can be used to evaluate other useful quantities such as the reduced electric transition probability and the deformation parameters. The former can be evaluated according to

$$B(E\lambda;0^{+}_{gs}\rightarrow\nu)=\frac{Z^{2}}{A^{2}}e^{2}|q_{\nu}|^{2},$$

where the proton density is assumed to be proportional to the neutron density. For instance, one can evaluate $B(E2;0^{+}_{g.s.}\rightarrow 2^{+})$ from the TDHF evolution of the quarupole moment $Q_{20}(t)$ following a quadrupole boost.

The deformation parameter can be computed from

$$\beta_{\lambda}^{(\nu)}=\frac{4\pi|q_{\nu}|}{3AR_{0}^{\lambda}},$$

where $R_{0}\simeq r_{0}A^{1/3}$ is the nuclear radius with $r_{0}\simeq 1.1-1.2$ fm. It is often useful to compare deformations of different nuclei, as well as to account for couplings to low-lying vibrations in coupled-channel calculations of heavy-ion collisions.

Giant resonances are highly collective nuclear vibrations that can usually decay by emitting nucleons. Their first phonons are often found between 10 and 30 MeV, depending on their multipolarity as well as their scalar/vector and isoscalar/isovector characteristics. Although most TDHF applications to nuclear vibrations are dedicated to giant resonances, the fact that they decay via particle emission means that such calculations has to be treated with care to avoid spurious finite size effects of the numerical box. To avoid using very large boxes (which often limit the calculations to spherical symmetric systems), one can use absorbing boundary conditions. In this case, the resulting strength functions correspond to continuum-RPA calculations.

An example of TDHF calculation of a monopole giant resonance in ²⁰⁸Pb is shown in Fig. 2. The direct decay of the oscillation amplitude is due to evaporation of nucleon wave-functions that are leaving the nucleus. This decay is exponential and contributes to the width of the peak in the strength function.

TDHF can be used in the small amplitude regime to extract energies and transition amplitudes of vibrational states at the RPA level. This is not, however, how a standard RPA code works. The goal is now to show the connection between linearised TDHF and RPA in more details.

When a small boost $e^{-i\varepsilon\hat{Q}}$ is applied to the wave function, it induces a perturbation to the one-body density matrix

$$\rho=\rho^{(0)}+\varepsilon\rho^{(1)}+O(\varepsilon^{2})$$

(26)

where $\rho^{(0)}=\sum_{h=1}^{A}|\varphi_{h}^{(HF)}\rangle\langle\varphi_{h}^{(HF)}|$ is the one-body density matrix of the Hartree-Fock ground-state. In TDHF, the many-body state remains a Slater with the corresponding one-body density matrix $\rho=\sum_{i=1}^{A}|\varphi_{i}\rangle\langle\varphi_{i}|$. One can see that $\rho$ is a projector onto the sub-space of occupied single-particle states by computing $\rho^{2}=\sum_{i,j=1}^{A}|\varphi_{i}\rangle\langle\varphi_{i}|\varphi_{j}\rangle\langle\varphi_{j}|=\rho$ as $\langle\varphi_{i}|\varphi_{j}\rangle=\delta_{ij}$. Using this property with the expression (26) gives

	$\displaystyle\rho^{2}$	$\displaystyle=$	$\displaystyle{\rho^{(0)}}^{2}+\varepsilon\left(\rho^{(0)}\rho^{(1)}+\rho^{(1)}\rho^{(0)}\right)+O(\varepsilon^{2})$
	$\displaystyle=\rho$	$\displaystyle=$	$\displaystyle\rho^{(0)}+\varepsilon\rho^{(1)}+O(\varepsilon^{2}).$

As a result, ${\rho^{(0)}}^{2}=\rho^{(0)}$ and

$$\rho^{(1)}=\rho^{(0)}\rho^{(1)}+\rho^{(1)}\rho^{(0)}.$$

(27)