Non-equilibrium quantum dynamics with collisional correlations

Our formal starting point is Stochastic Time-Dependent Hartree–Fock (STDHF) [28], which provides a description of dynamical correlations in terms of a stochastic ensemble of TDHF (mean-field) states. Although originally formulated in nuclear dynamics, it can be used in any system where a quantum mean field provides a sound description of ground state and low-energy properties, such as nuclei, molecules, or clusters. In the electronic case the underlying practical framework will be DFT [29] with local density approximation (LDA) as associated effective mean field. Time dependent LDA (TDLDA) is meanwhile widely used for simulating various dynamical scenarios [30, 31]. STDHF as proposed here is equally applicable to DFT thus delivering a Stochastic TDLDA (STDLDA).

Originally, STDHF was formulated in terms of density matrices with the von Neumann equation as formal starting point. It was also shown that STDHF can be reduced to a quantum Boltzmann–Langevin equation, namely a kinetic equation for the one-body density matrix containing a stochastic term beyond standard kinetic picture [28]. STDHF thus formally contains all ingredients for describing the path to thermalization and the thermal fluctuations around it [2, 28]. But it has never been explored in practice at full quantum level because of its complexity. Our aim here is to reformulate the theory in simpler terms, restricting correlations to their leading two-particle-two-hole ($2ph$) contribution.

Means of the description is an ensemble of $\mathcal{N}$ Slater states $\left\{ \left| {{\Phi }^{\alpha }} \right\rangle ,\alpha =1,...\;\mathcal{N} \right\}$. The $\left| {{\Phi }^{\alpha }} \right\rangle$ are built from single-particle (s.p.) states $\left\{ \varphi _{i}^{\alpha },i=1...\Omega \right\}$, where $i=1,...N$ stands for the occupied (hole) states. We also carry a sufficient amount of unoccupied (particle) states $i=N+1,...\Omega$ to supply space for the stochastic jumps to come. The enlarged space now allows to unfold the hierarchy of n-particle-n-hole (nph) excitations. The $1ph$ excitations need not to be taken care of, as they are already accounted for in the mean-field ground state and propagation. The first true excitations beyond mean field come along with $2ph$ states

$$\begin{eqnarray}&&\left| \Phi _{pp^{\prime} hh^{\prime} }^{\alpha } \right\rangle =\hat{a}_{p}^{\dagger }\hat{a}_{p^{\prime} }^{\dagger }\hat{a}_{h^{\prime} }^{{}}\hat{a}_{h}^{{}}\left| {{\Phi }^{\alpha }} \right\rangle ,\end{eqnarray} \tag{ 1 }$$

which as such are, again, Slater states. The idea is then to propagate for some time a correlated state $\left| {{\Psi }^{\alpha }}\left( t \right) \right\rangle =\left| {{\Phi }^{\alpha }}\left( t \right) \right\rangle +{{\sum }_{pp^{\prime} hh^{\prime} }}c_{pp^{\prime} hh^{\prime} }^{\alpha }\left( t \right)\left| \Phi _{pp^{\prime} hh^{\prime} }^{\alpha }\left( t \right) \right\rangle$ starting from an uncorrelated situation, $c_{pp^{\prime} hh^{\prime} }^{\alpha }\left( {{t}_{0}} \right)=0$. After a certain propagation time τ this coherent state $\left| {{\Psi }^{\alpha }}\left( t \right) \right\rangle$ is sampled in terms of an ensemble of pure Slater states $\left\{ \left| \Phi _{\kappa }^{\alpha } \right\rangle ,w_{\kappa }^{\alpha } \right\}$, where $\kappa \in \left\{ 0,pp^{\prime} hh^{\prime} \right\}$ and $w_{\kappa }^{\alpha }=|c_{pp^{\prime} hh^{\prime} }^{\alpha }{{|}^{2}}$ is the probability with which $\left| \Phi _{\kappa }^{\alpha } \right\rangle$ appear. The weight of the original state $\left| {{\Phi }^{\alpha }} \right\rangle$ is the complement $w_{0}^{\alpha }=1-\sum w_{pp^{\prime} hh^{\prime} }^{\alpha }$ (we attribute $w_{0}^{\alpha }$ to the 'no transition' case).

We evaluate $c_{pp^{\prime} hh^{\prime} }^{\alpha }$ in first order time-dependent many-body perturbation theory, propagated over a time span τ and delivering eventually the transition probability according to Fermiʼs golden rule [28]. In the $2ph$ picture it reads in detail

$$\begin{eqnarray}&&w_{pp^{\prime} hh^{\prime} }^{\alpha }=\tau {{\left| \left\langle \Phi _{pp^{\prime} ,hh^{\prime} }^{\alpha }|\hat{W}|{{\Phi }^{\alpha }} \right\rangle \right|}^{2}}{{\delta }_{\Gamma }}\left( \varepsilon _{p}^{\alpha }+\varepsilon _{p^{\prime} }^{\alpha }-\varepsilon _{h}^{\alpha }-\varepsilon _{h^{\prime} }^{\alpha }+E_{pp^{\prime} hh^{\prime} }^{\left( {\rm rearr} \right)} \right)\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\delta }_{\Gamma }}\left( \epsilon \right)=\vartheta \left( \epsilon -\Gamma /2 \right)\vartheta \left( \Gamma /2-\epsilon \right),\end{eqnarray} \tag{ 3 }$$

where $\hat{W}$ is the residual interaction complementing the mean-field Hamiltonian $\hat{h}$. The rearrangement energy $E_{pp^{\prime} hh^{\prime} }^{\left( {\rm rearr} \right)}$ corrects for the change in the mean field due to the actual $2ph$ transition (the actual transition preserving the total energy by construction [28]). It is negligible in large systems because 2 changing states out of hundredth of particles will not matter much. However, it cannot be ignored in small systems we want to address. The ϑ in equation (2) is the Heaviside function. The finite width is realistic as it can be related to a finite sampling time and, even for long sampling times, to the fluctuations of the s.p. energies in a mean-field (TDHF/TDLDA) evolution. The Heaviside function is used because it is technically the simplest finite-width function. Moreover, we are not trying to adjust the width to propagation time τ. In this first (proof of principle) paper, we simply employ the width $\Gamma$ as a numerical parameter to deal with collisions in a discrete spectrum. Of course, the width spoils to some extent energy conservation. One has thus to find a good compromise for $\Gamma$. It has to be larger than the average distance of excitation levels to provide an appropriate coverage of reaction rates. On the other hand, a small $\Gamma$ is preferable for minimal violation of energy conservation.

Note that the scheme is based on time-dependent perturbation theory. It thus stays in the weak coupling limit typical of standard kinetic theory. This implies restrictions on the strength of the residual interaction $\hat{W}$ and on the choice of the sampling interval τ [28, 32]. The restriction on the amplitude of $\hat{W}$ is governed by the nature of the considered system. As long as mean field provides a good starting point, $\hat{W}$ is expected to be small enough. The restriction on τ is more of practical nature. Indeed τ has to be long enough to randomize the quantum phase oscillations in the $2ph$ correlations, thus justifying stochastic reductions and loss of coherence, but short enough to maintain $\sum w_{pp^{\prime} hh^{\prime} }^{\alpha }\ll 1$. This means that in practice it has to be checked that it fulfills the above conditions. Once fixed it then allows to preserve perturbation theory on short time (of order τ) while allowing the progressive building up of (possibly large) perturbations (well beyond the pertubative regime) on long times ($\gg \tau$).

The stepping scheme thus proceeds as follows. An initial state $\left| {{\Phi }_{0}} \right\rangle$ is defined and each member of the ensemble is initially set to $\left| {{\Phi }^{\alpha }}\left( 0 \right) \right\rangle =\left| {{\Phi }_{0}} \right\rangle$. Each $\left| {{\Phi }^{\alpha }}\left( t \right) \right\rangle$ is then propagated individually. First, we propagate from 0 to τ by TDHF/TDLDA. At time τ one scans all $2ph$ states about $\left| {{\Phi }^{\alpha }}\left( \tau \right) \right\rangle$ and evaluates the jump probabilities $w_{\kappa }^{\alpha }$ according to equation (3). One state $\left| \Phi _{\kappa }^{\alpha } \right\rangle$ out of the possible κ is then chosen randomly according to its weight $w_{\kappa }^{\alpha }$. Subsequently, we propagate $\left| \Phi _{\kappa }^{\alpha } \right\rangle$ by TDHF/TDLDA from τ to $2\tau$ and perform another stochastic sampling at $2\tau$ and so on. The above procedure is restarted from initial time for each member of the ensemble separately delivering eventually the ensemble $\left\{ \left| {{\Phi }^{\alpha }}\left( t \right) \right\rangle ,\alpha =1,...\;\mathcal{N} \right\}$:

The s.p. energies entering crucially the ${{\delta }_{\Gamma }}$ function in equation (3) require a word of caution. A Slater state $\left| {{\Phi }^{\alpha }} \right\rangle$ is invariant under an arbitrary unitary transformation amongst the occupied s.p. states $\left| \varphi _{i}^{\alpha } \right\rangle ,i=1...N$. Thus the definition of s.p. energies is ambiguous. But the δ function in equation (3) is, in fact, an approximation. The full expression involves an operator δ function of the mean-field Liouvillean [28]. The replacement by s.p. energies is justified only if the uncertainties on the s.p. energies are small. To minimize these uncertainties, we diagonalize the actual mean-field Hamiltonian amongst the set of occupied states and also amongst the unoccupied states (as we have the freedom to do so) and use the emerging s.p. states and eigenvalues in equation (3). This renders the choice of s.p. energies unique and achieves the best possible approximation to the mean-field Liouvillean.

Once the ensemble $\left\{ \left| {{\Phi }^{\alpha }} \right\rangle ,\alpha =1,...,\mathcal{N} \right\}$ constructed, we can compute any observable by standard statistical averages. One-body observables, in particular, are deduced from the (correlated) one-body density operator built from the total ensemble as

$$\begin{eqnarray}&&\hat{\rho }=\frac{1}{\mathcal{N}}\sum \limits_{\alpha =1}^{\mathcal{N}}\sum \limits_{i=1}^{N}|\varphi _{i}^{\alpha }\rangle \langle \varphi _{i}^{\alpha }|\equiv \sum \limits_{\nu =1}^{\Omega }|\varphi _{\nu }^{\left( {\rm nat} \right)}\rangle {{n}_{\nu }}\langle \varphi _{i}^{\left( {\rm nat} \right)}|.\end{eqnarray} \tag{ 4 }$$

The second representation therein employs the natural s.p. orbitals $|\varphi _{\nu }^{\left( {\rm nat} \right)}\rangle$. These are the orbitals which diagonalize $\hat{\rho }$. The eigenvalues are the (fractional) occupation numbers ${{n}_{\nu }}$ of the natural orbitals and provide a natural tool for analyzing thermal effects. A useful quantity deduced from that is the distribution of occupation number $n\left( \varepsilon \right)$ built from identifying $n\left( {{\varepsilon }_{\nu }} \right)={{n}_{\nu }}$ (see figures 1 and 2 below). We also use the ${{n}_{\nu }}$ to evaluate the one-body entropy

$$\begin{eqnarray}&&S=-{{k}_{B}}\sum \limits_{\nu }\left[ {{n}_{\nu }}{\rm log} {{n}_{\nu }}+\left( 1-{{n}_{\nu }} \right){\rm log} \left( 1-{{n}_{\nu }} \right) \right]\end{eqnarray} \tag{ 5 }$$

which characterizes the approach to thermal equilibrium (${{k}_{B}}$ is Boltzmann constant).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A word is in order about the level of analysis. Equation (5) defines what is called the one-body entropy, which is the standard entropy associated to a (one-body) mean field description. In a pure TDHF time evolution the occupation numbers ${{n}_{\nu }}$ remain strictly 0 or 1, as fixed at initial time, and the entropy strictly 0. Coherent correlations (on top of a mean field description) would lead to some fractional ${{n}_{\nu }}$, even in the ground state. They would produce pattern which would stay far away from a thermal Fermi distribution and possibly pollute thermal observables. This problem is nevertheless negligible in the present case as we consider high excitations where thermal effects dominate anyway. Thus we think that the definition of entropy and temperature at one-body level suffices for the present purposes and clearly shows how the present theory allows to accommodate the building up of thermal effects, especially as compared to a pure mean field approach.

For a first exploration and proof of principle, we consider a simple 1D model mimicking a moleculer case with the mean field Hamiltonian (in x representation, $\hbar =1$)

$$\begin{eqnarray}&&{{\hat{h}}^{\alpha }}=-\frac{\Delta }{2m}+{{V}_{{\rm ext}}}\left( x \right)+\lambda {{\left( {{\varrho }^{\alpha }}\left( x \right) \right)}^{2}}\end{eqnarray} \tag{ 6 }$$

where ${{\varrho }^{\alpha }}\left( x \right)=\sum _{1}^{N}|\varphi _{i}^{\alpha }\left( x \right){{|}^{2}}$ is the local one-body density. For the external potential ${{V}_{{\rm ext}}}\left( x \right)$ we use a Woods-Saxon profile ${{V}_{{\rm ext}}}\left( x \right)={{V}_{0}}/\left( 1+{\rm exp} \left( \left( x-{{x}_{0}} \right)/a \right) \right)$ with ${{V}_{0}}=-5$ Ry, ${{x}_{0}}=15$ a₀, a = 2 a₀, complemented by a confining harmonic oscillator $\frac{m}{2}{{\omega }^{2}}{{\left( |x|-{{x}_{0}}-a \right)}^{2}}\vartheta \left( |x|-{{x}_{0}}-a \right)$ with $\hbar \omega =4$ Ry, acting outside the potential well ($|x|>{{x}_{0}}+a$) and ensuring soft reflecting boundary conditions. This allows to focus the analysis on thermalization by avoiding competition with ionization. For the self-consistent term we use $\lambda =2$ Ry ${\rm a}_{0}^{2}$. This model resembles a typical situation in molecules or clusters, where the ionic background provides a spatially fixed external potential. With the chosen parameters we obtain a sufficiently dense spectrum of s.p. energies in spite of the 1D restriction, while producing typical molecular/cluster values of energies and time scales. As residual interaction we choose a simple zero-range force $W\left( x,x^{\prime} \right)={{W}_{0}}\delta \left( x-x^{\prime} \right)$ with ${{W}_{0}}=3$ Ry, which leads to realistic relaxation times for the entropy (see below). This choice differs from the mean-field (TDLDA) residual interaction according to the mean field used (which should be $\delta {{\hat{h}}^{\alpha }}/\delta {{\rho }^{\alpha }}$). A separate modeling of the residual interaction is legitimate, even necessary, to calibrate the collision rates properly in realistic cases. This is common practice in BUU and it is justified by the observation that the collision term should employ the screened interaction and the mean-field part the interaction before screening [33–35].

The test case has nine physical particles with the same spin. Our working space thus has 9 hole (occupied) levels and an arbitrary number of particle (unoccupied) levels. In practice we used 8 particle states. No significant differences were observed when using 16 or 24 particle states. The distance between s.p. level energies typically varies between 10⁻⁴ and 2.5 × 10⁻³ Ry in the range of excitation energies we have explored. In order to fix the optimal width $\Gamma$ of ${{\delta }_{\Gamma }}$ we have studied its impact on observables of interest, such as the entropy (equation (5)) as well as the a priori highly $\Gamma$-sensitive number of collisions and find that they are rather stable up to $\Gamma =0.1$ Ry turning to moderate growth above that value. The actual choice of $\Gamma$ is thus not so critical. All values up to $\Gamma =0.1$ Ry are acceptable. Very small values of $\Gamma$ require larger samples. Thus we choose $\Gamma =0.1$ Ry as the compromise for our present studies.

The period associated to the dominant optical transition lies in the fs range (1.15 fs). The time for propagating two-body correlations in perturbation theory is then chosen as $\tau =1\;{\rm fs}$, following [28]. We found identical results for $\tau =0.5$ and $\tau =1.5$ fs. We follow the dynamics on rather long times, typically 100 fs, which is long enough for thermal relaxation studies and much shorter than Poincaré recurrence time [36]. We propagate ensembles with 100 members.

There are several ways to excite the system. We use here an instantaneous initial random (multi)ph excitation. Looking at s.p. energies as a function of time one observes many level crossings between occupied and unoccupied levels. This gives rise to a sufficient number of $2ph$ energies matchings according to equation (2) and sufficiently large total probabilities (up some 10ʼs of %), to allow enough transitions into new configurations, while remaining in the weak coupling regime.

The STDHF dynamics is then generated for each member of the ensemble. The energy stays constant during TDHF/TDLDA evolution, but may change a bit in stochastic jumps due to the finite width of the ${{\delta }_{\Gamma }}$ functions in equation (2). The ensemble thus develops a small energy uncertainty and a negligible trend to growing energy. The energy spread in the ensemble only grows within 100 fs to about 0.05 % (${{E}^{*}}$ = 1.9 Ry) to about 0.3% (${{E}^{*}}$ = 3.1 Ry) of the total energy. A small center-of-mass fluctuation correspondingly appears with an amplitude of order 0.05 to 0.1 a₀ around average value 0 which is in the average perfectly preserved (no drift).

The relevant one-body observables are extracted from the ensemble averaged one-body density matrix $\hat{\rho }$, see equation (4). Figure 1 shows the distribution of occupation numbers $n\left( \varepsilon \right)$ at a few selected times for the case of an initial $1ph$ excitation delivering ${{E}^{*}}=2.3$ Ry. The s.p. energies ${{\varepsilon }_{\nu }}=\langle \varphi _{\nu }^{\left( {\rm nat} \right)}|\hat{h}|\varphi _{\nu }^{\left( {\rm nat} \right)}\rangle$ are evaluated as the energy expectation values of the natural orbitals (equation (4)) from the (ensemble) Hamiltonian $\hat{h}={{\sum }_{\alpha }}{{\hat{h}}^{\alpha }}/\mathcal{N}$ (${{\hat{h}}^{\alpha }}$ is the mean field associated to state $\left| {{\Phi }^{\alpha }} \right\rangle$, (equation (6)). The first snapshot at time 5 fs carries still much of the initial $1ph$ pattern. The subsequent snapshots show nicely the steady evolution towards a thermal equilibrium distribution. The finite sampling adds, of course, some fluctuations around the leading monotonic trend. The final snapshot looks already very smooth and resembles the typical Fermi distribution of a thermalized state. We can associate a temperature T to this state by fitting a thermal Fermi distribution to this final $n\left( \varepsilon \right)$ profile. This yields here $T\sim 0.48$ Ry. The corresponding Fermi function is also indicated on the plot for completeness. Note the gradual convergence of the actual distribution towards the profile, in particular on long times.

Figure 2 shows the final distributions of occupation numbers for various excitation energies ${{E}^{*}}$. Increasingly longer tails develop with increasing ${{E}^{*}}$. All profiles resemble thermal equilibrium distributions, particularly for higher energies. Fits to Fermi distributions deliver the following temperatures: $T\left( {{E}^{*}}=1.3\;{\rm Ry} \right)\sim 0.34$ Ry, $T\left( 1.9\;{\rm Ry} \right)\sim 0.46$ Ry, $T\left( 3.1\;{\rm Ry} \right)\sim 0.56$ Ry, $T\left( 3.5\;{\rm Ry} \right)\sim 0.66$ Ry and $T\left( 4.3\;{\rm Ry} \right)\sim 0.78$ Ry.

As we are dealing with reflecting boundary conditions the initial excitation energy can be fully converted into thermal energy. It thus makes sense to use it in a purely thermal picture. Figure 3 shows the temperatures deduced from fits of the $n\left( \varepsilon \right)$ to Fermi factors and the final values of the entropies S as a function of the excitation energy. The Fermi gas model provides an analytical estimate for both quantities [37] in the low temperature limit in terms of the level density parameter a:

$$\begin{eqnarray}&&{{E}^{*}}=a{{T}^{2}},\quad S=2aT.\end{eqnarray} \tag{ 7 }$$

We have evaluated a for the present model by a direct fit of final S versus initial ${{E}^{*}}$ (lower panel of figure 3). The corresponding estimates of temperatures are also shown in figure 3 (upper panel). The results agree perfectly for the entropy S. There is some acceptable deviation for the temperature. Mind that the temperature was determined from fit of a Fermi distribution to the occupations $n\left( \varepsilon \right)$. The fluctuations in $n\left( \varepsilon \right)$ due to the finite sampling cause some uncertainties which somewhat degrade the computation of the temperature. The entropy is probably more robust in that respect because it is obtained by summing up over occupation numbers.

Zoom In Zoom Out Reset image size

Finally, we extract thermalization time scales from the one-body entropy S, see equation (5). Figure 4 shows the time evolution of entropy for various ${{E}^{*}}$. We see the typical pattern of a thermal relaxation process with an early rapid growth turning later to steadily slower increase. Except for the lowest excitation energy probably suffering from too large fluctuations or too low statistics the pattern resemble an exponential relaxation towards the maximum entropy for a given ${{E}^{*}}$. We have fitted exponential relaxation profiles $S={{S}_{0}}\left( 1-{\rm exp} \left( t/{{\tau }_{{\rm rel}}} \right) \right)$ to the S(t)ʼs. This delivers relaxation times ${{\tau }_{{\rm rel}}}$ decreasing with increasing ${{E}^{*}}$, from ∼22 fs at ${{E}^{*}}=1.5$ Ry down to ∼5 fs at ${{E}^{*}}=3.5$ Ry, representing realistic values for such a system [38].

Zoom In Zoom Out Reset image size