Impact of the initial fluctuations on the dissipative dynamics of interacting Fermi systems: A model case study

TL;DR

This paper investigates how initial fluctuations influence the dissipative dynamics in interacting Fermi systems, demonstrating that non-Gaussian initial distributions improve long-term evolution predictions beyond traditional Gaussian assumptions.

Contribution

It introduces a modification to the stochastic mean-field approach by relaxing the Gaussian assumption, enhancing the accuracy of long-time dissipative dynamics modeling.

Findings

Non-Gaussian initial distributions reduce overdamping effects.

Distributions with smaller kurtosis better reproduce long-term evolution.

Gaussian assumptions tend to overestimate damping in the dynamics.

Abstract

Standard methods used for computing the dynamics of a quantum many-body system are the mean-field (MF) approximations such as the time-dependent Hartree-Fock (TDHF) approach. Even though MF approaches are quite successful, they suffer some well-known shortcomings, one of which is insufficient dissipation of collective motion. The stochastic mean-field approach (SMF), where a set of MF trajectories with random initial conditions are considered, is a good candidate to include dissipative effects beyond mean field. In this approach, the one-body density matrix elements are treated initially as a set of stochastic Gaussian c numbers that are adjusted to reproduce first and second moments of collective one-body observables. It is shown that the predictive power of the SMF approach can be further improved by relaxing the Gaussian assumption for the initial probabilities. More precisely, using Gaussian or uniform distributions for the matrix elements generally leads to overdamping for long times, whereas distributions with smaller kurtosis lead to much better reproduction of the long time evolution.

Equations

$$i\hslash \frac{\partial }{\partial t}\rho =[h[\rho ],\rho ],$$

$$i\hslash \frac{\partial }{\partial t}{\rho }^{\lambda }=[h[{\rho }^{\lambda }],{\rho }^{\lambda }].$$

$$⟨A{⟩}^{\lambda }=\mathrm{Tr}({\rho }^{\lambda }A),$$

$$\overline{⟨A{⟩}^{\lambda }}=$$

$$=$$

$$\overline{{\rho }_{ij}^{\lambda }}=\underset{\mathcal{N}\to \infty }{\lim }\frac{1}{\mathcal{N}}\sum _{\lambda =1}^{\mathcal{N}}{\rho }_{ij}^{\lambda }.$$

$${\sigma }_A^{(2)}=$$

$$=$$

$$⟨A(0)⟩=\sum _i{n}_i{A}_{ii},$$

$$⟨A^2(0)⟩-⟨A(0){⟩}^2=\frac{1}{2}\sum _{ij}\left[n_i(1-n_j)+n_j(1-n_i)\right]A_{ji}{A}_{ij},$$

$$\overline{{\rho }_{ij}^{\lambda }(0)}=$$

$$\overline{\delta {\rho }_{ij}^{\lambda }(0)\delta {\rho }_{kl}^{\lambda }(0)}=$$

$${\sigma }_A^{(m)}=$$

$$=$$

$${\sigma }_A^{(3)}=$$

$${\sigma }_A^{(4)}=$$

$$⟨(\Delta A{)}^3⟩=$$

$$⟨(\Delta A{)}^4⟩=$$

$$+3\sum _{ijkl}{\Lambda }_{ijkl}^{(4b)}{A}_{ij}{A}_{ji}{A}_{kl}{A}_{lk},$$

$${\Lambda }_{ijk}^{(3)}=$$

$$+n_j(1-3n_k)(1-n_k{n}_i)],$$

$${\Lambda }_{ijkl}^{(4a)}=$$

$$+n_l(1-4n_i)(1-3n_j)(1-n_i{n}_j{n}_k)$$

$$+n_k(1-4n_l)(1-3n_i)(1-n_l{n}_i{n}_j)$$

$$+n_j(1-4n_k)(1-3n_l)(1-n_k{n}_l{n}_i)],$$

$${\Lambda }_{ijkl}^{(4b)}=$$

$$+n_j{n}_k\left[(1-2n_l)(1-n_l{n}_i)+(1-2n_i)(1-n_i{n}_l)\right]$$

$$+n_i{n}_l\left[(1-2n_k)(1-n_k{n}_j)+(1-2n_j)(1-n_j{n}_k)\right]$$

$$+n_j{n}_l\left[(1-2n_k)(1-n_k{n}_i)+(1-2n_i)(1-n_i{n}_k)\right]\}.$$

$$\overline{\delta {\rho }_{ij}^{\lambda }(0)\delta {\rho }_{kl}^{\lambda }(0)\delta {\rho }_{mn}^{\lambda }(0)}=$$

$$\overline{\delta {\rho }_{ij}^{\lambda }(0)\delta {\rho }_{kl}^{\lambda }(0)\delta {\rho }_{mn}^{\lambda }(0)\delta {\rho }_{pr}^{\lambda }(0)}=$$

$$+{\delta }_{il}{\delta }_{jk}{\delta }_{rm}{\delta }_{pn}3{\Lambda }_{jipr}^{(4b)}.$$

$$\delta {\rho }_{ii}=0$$

$$\displaystyle\overline{r_{ij}^{2}}+\overline{s_{ij}^{2}}=\left\{\begin{array}[]{cc}\frac{1}{2}&\qquad\text{if}~{}(i,j)=(p,h),\\ \\ 0&\text{otherwise},\end{array}\right.$$

$${\Lambda }_{pph}^{(3)}$$

$${\Lambda }_{phh}^{(3)}$$

$${\Lambda }_{ppph}^{(4a)}$$

$${\Lambda }_{hhhp}^{(4a)}$$

$${\Lambda }_{pphh}^{(4a)}$$

$${\Lambda }_{phph}^{(4a)}$$

$${\Lambda }_{phph}^{(4b)}$$

$$\overline{\delta {\rho }_{ij}\delta {\rho }_{ki}\delta {\rho }_{jk}}=$$

Full text

I Introduction

In many situations, the evolution of a quantum system can be replaced by a set of classical evolutions with random initial conditions optimized to reproduce at best the initial quantum zero-point motion Davis and Heller (1984). This quantum-to-classical mapping is particularly suitable in the absence of interference or tunneling and can be exact in some cases, like the free wave expansion or the quantum harmonic oscillator. This idea is employed in many fields of physics to describe the out-of-equilibrium motion of complex systems, like in quantum optics Gardiner and Zoller (2004). This also includes the possibility to describe many-body interacting systems. An illustration in bosonic systems is the truncated Wigner approximation (TWA) Sinatra et al. (2002) (see also Polkovnikov (2010)). For Fermi systems, as underlined in ref. Davidson et al. (2017), such mapping is more tricky due to the absence of natural classical representation, in contrast to bosonic systems. Despite this inherent difficulty for fermionic systems, one can mention two attempts to map the complex many-body problem of interacting systems as a set of "classical trajectories". The first one is the stochastic mean-field (SMF) approach Ayik (2008); Lacroix and Ayik (2014), where the stochastic one-body density matrices are treated as classical objects evolving through the time-dependent mean-field (TDMF) equation of motion. A second formulation was made more recently in ref. Davidson et al. (2017) based on the mapping between fermionic and bosonic operators, leading to the fermionic TWA (f-TWA). We actually realized recently that the equation of motion obtained in the f-TWA often coincides with the TDMF evolution, and therefore these two approaches are relatively close to each other. In recent decades, the SMF approach has been successfully applied to model cases Lacroix et al. (2012, 2014); Yilmaz et al. (2014a) as well as realistic simulations of dynamical phenomena in nuclear physics Ayik et al. (2009, 2010); Yilmaz et al. (2011, 2014b); Tanimura et al. (2017); Ayik et al. (2017, 2018, 2019). It has also been extended to describe superfluid systems in ref. Lacroix et al. (2013). The SMF approach and the f-TWA approach have also in common that they both assume Gaussian probabilities for the initial fluctuations. However, this assumption turns out to be guided more by practical arguments than by first principles. In the present work, we further explore the possibility of using an alternative probability distribution for the initial conditions in the SMF approach. We show that, while the second moment of a one-body observable can be generally interpreted classically, the fourth moment of a one-body observable is more problematic, since it can lead to negative values for the fourth moments of the stochastic matrix elements of one-body densities when considering many-body Fermi systems. Such behavior cannot be reproduced by a classical mapping. However, probability distributions with smaller kurtosis compared to the Gaussian distribution turn out to be more efficient to describe the time evolution in the SMF approach. Such finding is illustrated on a modified version of the Lipkin-Meshkov-Glick model Lacombe et al. (2016).

One of the interesting aspects of SMF is that it provides a fully microscopic description of fluctuation and dissipation. It also can be used to make connection with the more phenomenological Langevin approach in collective space where few relevant collective degrees of freedom are preselected Ayik (2008). The Langevin technique has been successfully applied to dissipative nuclear process like fission and/or heavy-ion collisions Abe et al. (1996); Karpov et al. (2001); Zagrebaev and Greiner (2005); Aritomo (2009); Wen et al. (2013); Karpov and Saiko (2017); Sekizawa and Hagino (2019); Lin and Danielewicz (2019). Although we are still trying to quantify and/or improve the predictive power of the SMF approach, in the long term it might provide a fully microscopic framework for dissipative process without needing assumptions inherent to more phenomenological approaches, like the selection of collective degrees of freedom or the assumption on the nature of the noise.

The paper is organized as follows. In Sec. II, the third and fourth moments of the matrix elements of the stochastic one-body densities are derived within the SMF approach and it is shown that initial probability distribution functions with small kurtosis can provide a better approximation to the fourth central moment. In Sec. III, the SMF dynamics is applied to a modified version of the Lipkin-Meshkov-Glick model. Finally, the conclusions are given in Sec. IV.

II Stochastic mean-field approach

The Schrödinger equation of a many-body fermionic system can be cast into a hierarchy of differential equations for reduced densities (one-body, two-body, etc.) known as BBGKY (Bogoliubov-Born-Green-Kirkwood-Yvon) hierarchy. The truncation of BBGKY equations at lowest order gives the mean-field equation for the dynamics of one-body density,

[TABLE]

where $h[\rho]$ is the mean-field Hamiltonian. In the mean-field approximation, one-body density contains all the information on the system, hence the many-body state is restricted to be a Slater determinant during the dynamical evolution. MF dynamics (or TDHF) is known to provide good approximation to one-body observables but severely underestimate their quantum fluctuations. Beyond -mean-field approaches are necessary to overcome these shortcomings and provide a more accurate description of the dynamical evolution.

In the SMF approach Ayik (2008), the initial one-body density $\rho(0)$ is replaced by an ensemble of stochastic initial one-body densities $\rho^{\lambda}(0)$, where $\lambda$ stands for event label. Each of these densities evolves with its self-consistent mean-field equation,

[TABLE]

For each event $\lambda$, the ”event“ expectation value of a one-body observable $A$ is given by

[TABLE]

where ”event“ means one trajectory associated with a specific distribution for the stochastic initial one-body densities.

In the SMF approach, the expectation values are obtained by ensemble averages. Hence, the expectation value of the one-body observable is defined as,

[TABLE]

where the overline stands for ensemble averaging explicitly given by

[TABLE]

Here, $\mathcal{N}$ is the number of events in the ensemble. The quantum variance of the one-body operator is defined by

[TABLE]

where $\delta\rho^{\lambda}=\rho^{\lambda}-\overline{\rho^{\lambda}}$.

In the SMF approach, the distribution of the stochastic matrix elements, $\rho_{ij}^{\lambda}(0)$, is chosen such that the initial mean and variances of the observables are the same with those obtained with the initial density $\rho(0)$. If $\{|i\rangle\}$ stands for the natural basis which satisfies $\langle i|\rho(0)|j\rangle=n_{i}\delta_{ij}$, the mean and variance of $A$ are the expectation values given by

[TABLE]

and

[TABLE]

respectively. Comparing Eqs. (II) and (II) with Eqs. (7) and (8), we see that the quantum averages match the ensemble averages if we have,

[TABLE]

In the original formulation of the SMF approach Ayik (2008), the stochastic matrix elements $\delta\rho^{\lambda}_{ij}(0)$ are assumed to be Gaussian random numbers satisfying Eqs. (9) and (10).

Higher order moments of the stochastic matrix elements of the one-body density

Here, we derive higher order moments for the stochastic matrix elements and test the Gaussian assumption of the original formulation of the SMF approach when the initial state is a Slater determinant. For this purpose, similarly to the second central moment Eq. (II), we define the expectation value of the mth central moment of a one-body operator $A$ as

[TABLE]

where $\delta\rho^{\lambda}=\rho^{\lambda}-\overline{\rho^{\lambda}}$. Hence, the third and fourth central moments can be written as,

[TABLE]

The corresponding quantum central moments, for an initial Slater determinant, are given by (see the Appendix for details)

[TABLE]

where $\Delta A=A-\langle A\rangle$ and

[TABLE]

Comparing Eqs. (12) and (13) with Eqs. (14) and (II), we see that the quantum and ensemble averages coincide if we have

[TABLE]

By considering that all initial densities $\rho^{\lambda}(0)$ in the ensemble are Hermitian [in the following, we drop the labels for brevity, i.e., $\delta\rho_{ij}\equiv\delta\rho_{ij}^{\lambda}(0)$.], it immediately follows from Eq. (10) that $\overline{(\delta\rho_{ii})^{2}}=0$ which means that

[TABLE]

and hence the matrix elements $\rho_{ii}$ are not fluctuating. Furthermore, for $i\neq j$, we have

[TABLE]

where we introduced the real and imaginary parts of the matrix element dispersions, $\delta\rho_{ij}=r_{ij}+i\,s_{ij}$. Here, $(i,j)=(p,h)$ means that one of the states is a particle state and the other a hole state. The variance cannot be negative and hence, from Eq. (25), we have $\overline{r_{ij}^{2}}=\overline{s_{ij}^{2}}=0$ for $(i,j)=(p,p)$ or $(i,j)=(h,h)$ which means that these matrix elements are zero. Note that the stochastic matrix elements satisfy the equations $r_{ji}=r_{ij}$ and $s_{ji}=-s_{ij}$ due to Hermiticity.

In order to deduce the properties required for $\delta\rho$ to fulfill the moments Eqs. (19) and (II), we need to analyze the possible values of the $\Lambda$ terms. In general, from Eqs. (II)-(II), the $\Lambda$ terms can assume the following nonzero values:

[TABLE]

where p stands for a particle and h stands for a hole state. The nonzero terms of Eq. (19) are given by

[TABLE]

Similarly, the nonzero terms of the first and second terms of Eq. (II), by taking into account the Kronecker delta functions, read as

[TABLE]

and

[TABLE]

respectively.

The form of the last three equations suggests that these equations can be satisfied if one considers that the stochastic matrix elements are correlated with each other. However, averages that contain terms with the same indices $\delta\rho_{jj}$ in Eqs. (27) and (28) lead to, for instance, $\overline{\delta\rho_{ij}\delta\rho_{ji}\delta\rho_{jj}}=\Lambda^{(3)}_{jij}$ and $\overline{\delta\rho_{ij}\delta\rho_{ii}\delta\rho_{ji}\delta\rho_{jj}}=\Lambda^{(4a)}_{jiij}$, where left-hand side of these equations is zero due to Eq. (21) and the right-hand side of the equations is non-zero as seen from Eq. (II). Hence, Eqs. (27) and (28) cannot be satisfied completely even with correlated matrix elements. The reason for such a discrepancy is related to the fact that probability distributions cannot simulate quantum mechanical systems. In conclusion, probability distributions can only approximately match higher quantum mechanical moments. Here, we explore probability distributions that can at best approximate the third and fourth quantum mechanical moments.

In this study, we investigate the consequences of the simple case of uncorrelated matrix elements. Hence, we assume that each stochastic matrix element is statistically independent of the others as well as that the real and imaginary parts of a matrix element are also statistically independent. Then, for simplicity, we set $\overline{\delta\rho_{ij}\delta\rho_{ji}\delta\rho_{ij}}=\overline{r_{ij}^{3}}+i\overline{s_{ij}^{3}}=0$ for the third central moments, which means

[TABLE]

This is in accordance with the original formulation of the SMF approach since the third central moment of a mean-zero Gaussian random number is zero. Note that Eq. (27) does not provide any restrictions for the $\overline{\delta\rho_{ij}\delta\rho_{ji}\delta\rho_{ij}}$ moments. However, a condition for the fourth central moments can be obtained by assuming uncorrelated matrix elements. It is seen from Eqs. (28) and (29) that one has

[TABLE]

which leads to the result,

[TABLE]

where Eq. (II) was used.

No probability distribution can have a negative second and/or fourth central moment. Hence, the left-hand side of Eq. (35) can never be negative. However, the distributions that we want to obtain are quantum mechanical in nature. Actually the noncommutativity of operators within quantum mechanics is the main reason for inadequacy of probability distributions to simulate the quantum systems. In the phase-space formulation of quantum mechanics, we already know such behaviors coming from the Wigner distribution which can assume negative values and therefore is called a quasiprobability distribution Bednorz and Belzig (2011). It is important to note that positive-definite phase-space distributions such as the Husimi distribution exist; however, they are still quasiprobability distributions since the averages are taken by the Weyl symbol of the operators rather than the operators themselves Yilmaz et al. (2014a).

Based on the discussion above, since a probability distribution cannot lead to a negative fourth central moment, we anticipate that an efficient approximation would be to consider a distribution which minimizes the left-hand side of Eq. (35) while satisfying Eq. (25). Equation (25) sets a condition on the sum of the variances and hence it does not impose any conditions on the weights to the variances of the real and imaginary parts. In order to investigate the dependence to the weights, we define a parameter $\chi=\overline{r^{2}_{ij}}$. From Eq. (25), we have $\overline{s^{2}_{ij}}=1/2-\chi$, where $0\leq\chi\leq 1/2$. Substituting these results into the left-hand side of Eq. (35) we obtain the function

[TABLE]

where we used the kurtosis of the distribution, which is defined as

[TABLE]

The kurtosis is an invariant form of the fourth central moment that only depends on the distribution and does not depend on the particular values of the variance. Figure 1 shows the plot of the function (36) versus the weight parameter $\chi$ for the kurtosis values $3$ and $1.8$, which are the kurtoses of the Gaussian and the uniform distributions, respectively. The kurtosis value of $1$ corresponds to the minimum possible value. The probability distribution function with the minimum kurtosis is the so-called two-point distribution function which can take only two values, $\sigma$ and $-\sigma$, with equal probabilities Moors (1988). This probability function can be written as

[TABLE]

where $\delta(x)$ is the Dirac delta function and $\sigma$ is the standard deviation of $x$. It is seen from Fig. 1 that the function $F$ decreases as the kurtosis $\gamma$ decreases. We show below that the SMF dynamics with a small kurtosis probability distribution provides a better approximation to the exact dynamics than the SMF dynamics with a large kurtosis distribution. In particular, for the three distributions we considered here, the two-point distribution provides the best approximation and the Gaussian distribution provides the worst approximation. The choice of the weight parameter $\chi$ also has an impact on the values of the function $F$, as seen from Fig. 1. The equal weight case, $\chi=\overline{r^{2}_{ij}}=\overline{s^{2}_{ij}}=1/4$, provides the minimal value of $F$ for a fixed value of the kurtosis $\gamma$ except for the two-point distribution which is constant. Hence, the case with equal weights for the variances of the real and imaginary parts of the stochastic matrix elements is anticipated to provide a better approximation than any unequal weights case.

III Applications

In this section, we apply the SMF approach with the two-point distribution as well as the Gaussian and uniform distributions on an exactly solvable model and show that indeed the SMF dynamics with the two-point distribution for the matrix elements gives a better agreement with the exact dynamics than the SMF dynamics with the uniform or Gaussian distributions.

III.1 A modified Lipkin-Meshkov-Glick model

We consider a modified Lipkin-Meshkov-Glick model (mLMG) that was introduced in ref. Lacombe et al. (2016). It has a pairing Hamiltonian structure and the single particle energies are stochastically distributed over the levels allowing for a description of dissipation. The Hamiltonian of this model reads

[TABLE]

with

[TABLE]

where $\alpha=(s_{\alpha},m_{\alpha})$, $s_{\alpha}\in\{-1,+1\}$, $m_{\alpha}\in\{-j,-j+1,...,j-1,j\}$, $S_{+}=\sum_{\alpha>0}a^{\dagger}_{\alpha}a^{\dagger}_{\underline{\alpha}}$, and $S_{-}=S_{+}^{\dagger}$. Here, $\underline{\alpha}=(s_{\alpha},-m_{\alpha})$ and $\alpha>0$ means $m_{\alpha}>0$. $v_{0}$ is the interaction strength and $\varepsilon_{\alpha}/2$ are the single-particle energies which are random numbers with some variance $\sigma^{2}_{\varepsilon}$. The mean values of the single-particle energies in the upper level $s_{\alpha}=1$ and in the lower level $s_{\alpha}=-1$ are $\Delta/2$ and $-\Delta/2$, respectively. In order to avoid any ambiguity, we should state that the stochastic single-particle energies are only set once and these values are used for the dynamics. Figure 2 shows a schematic illustration of the mLMG model. Two transitions are indicated on the figure. Note that, contrary to the original LMG model Lipkin et al. (1965), there are transitions within the same level, which is due to the fact that the single-particle energies are different.

III.2 Exact, mean-field, and stochastic mean-field dynamics

The exact dynamics, for an initial state $|\Psi(0)\rangle$, is formally given by

[TABLE]

In practice, there are several numerical methods to solve this equation of motion. We use the iterated Crank-Nicolson method for computing the exact dynamics Teukolsky (2000).

The mean-field equation of motion (or TDHF equation) can be derived by using the Ehrenfest theorem,

[TABLE]

where the expectation value is taken with respect to a Slater determinant. By using Eq. (39), the last equation gives

[TABLE]

The SMF equation is directly obtained from Eq. (III.2) by replacing all one-body densities with the stochastic ones, $\rho\rightarrow\rho^{\lambda}$.

There are two different single-particle bases that are used in this study. One of the bases is the fixed basis of the model, which consists of the single-particle states associated with the creation and annihilation operators $a^{\dagger}_{\alpha}$, $a_{\alpha}$ appearing in the Hamiltonian Eq. (39). We use the labels $\alpha$, $\beta$, etc. for the fixed basis. The other basis is the natural basis which is obtained by diagonalizing the initial one-body density and we use the labels $i$, $j$, etc. for the states in this basis. The fixed basis do not change over time, whereas the natural basis change over time. Note that the MF equation (III.2) is written in the fixed basis, whereas the properties of the matrix elements of the initial stochastic one-body densities in the SMF approach are derived in the natural basis in Sec. II. We solve the SMF equations of motion in the fixed basis. Hence, for each event we sample an initial stochastic one-body density in the natural basis, we write it in the fixed basis, and then evolve it with the SMF equation which is just the MF equation.

III.3 The initial state

We consider the half-filled case for which there are $N$ particles for $2N$ single-particle states with initial states of the form

[TABLE]

where $\mu$ is a real parameter and $|\Phi^{(\eta)}\rangle$ is a Slater determinant given by

[TABLE]

Here, $n_{\alpha}^{(\eta)}$ are occupation numbers with values 0 or 1 corresponding to the single-particle states indicated in Fig. 2. The dipole operator $D$ reads

[TABLE]

Since $D$ is a one-body operator the initial state $|\Psi(0)\rangle$ is also a Slater determinant due to the Thouless theorem Thouless (1960). The effect of the exponential term in Eq. (45) is to introduce some excitation by an instant dipole boost Lacombe et al. (2016).

We consider two different initial states for $N=6$ particles. The first state corresponds to the ground state, $|\Phi^{(0)}\rangle$, with all the particles in the lower level ($s_{\alpha}=-1$) excited by the dipole boost in Eq. (45) with $\mu=0.8$. The second state is found by applying again the dipole boost with $\mu=0.8$ to the state $|\Phi^{(\eta)}\rangle$ for which

[TABLE]

single-particle states are occupied (see Fig. 2). The first and second states are denoted by $|\Psi(0)\rangle_{1}$ and $|\Psi(0)\rangle_{2}$, respectively.

Note that the schematic model in Fig. 2 is shown for 10 particles in 20 states. Here, we consider 6 particles in 12 states, therefore $m_{\alpha}$ can take the values $\{-5/2,-3/2,-1/2,1/2,3/2,5/2\}$.

III.4 Results

In the following computations, the units of energy will be written in terms of the level spacing $\Delta$ and the time units in terms of $\Delta^{-1}$. The single-particle energies are chosen as

[TABLE]

The mean of the single-particle energies in the upper level $\varepsilon_{(+1,m_{\alpha})}/2$ and in the lower level $\varepsilon_{(-1,m_{\alpha})}/2$ are $\Delta/2$ and $-\Delta/2$, respectively. The standard deviation of the $\varepsilon_{\alpha}/2$ values in each level is $\sigma_{\varepsilon}=0.2\,\Delta$. All the SMF results are obtained by considering $\mathcal{N}=10^{6}$ events.

We follow the dynamics of the dipole operator $D$ and the fermionic one-body entropy given by

[TABLE]

where $\rho$ is the one-body density operator. The one-body entropy is a measure of departure from an uncorrelated (Slater) state and a measure of thermalization Lacombe et al. (2016). For a Slater state, the entropy gets the value $S=0$. Hence, during the MF dynamics the entropy remains zero. The maximum value of the entropy for the half-filled case is obtained when all single particle states are half-filled, leading to the value

[TABLE]

In the SMF approach, the one-body entropy is computed by the expression

[TABLE]

Figure 3 shows the expectation value of the dipole operator $D$, given by Eq. (47), and the one-body entropy $S$, given by Eq. (50), per particle versus time for the weak coupling $v_{0}=0.05\,\Delta$. The initial state is $|\Psi(0)\rangle_{1}$ which is explained in the Sec. III.3. In Fig. 3(a), it is observed that MF evolution starts to deviate from the exact one at around $t=6\,\Delta^{-1}$, whereas in Fig. 3(b) the SMF evolutions start to deviate later at around $t=20\,\Delta^{-1}$. Comparing the three SMF dynamics with the exact one in Fig. 3(b), we see that the SMF approach with the two-point distribution is much better than the other two distributions at all times. A similar behavior is seen for the one-body entropy $S$ in Fig. 3(c). The SMF evolution of the entropy with the two-point distribution follows the exact solution very closely, whereas SMF solutions with the Gaussian and uniform initial distributions are underestimating the exact entropy until $t=30\,\Delta^{-1}$. After that time, the exact entropy decreases, which is followed by the SMF solution with two-point distribution, whereas the solutions with the Gaussian and uniform distributions continue to increase and approach the maximum value $2\ln 2=1.38$. The solution with the uniform distribution has a small decrease from $t=40\,\Delta^{-1}$ to $t=70\,\Delta^{-1}$, hence it is slightly better than the solution with the Gaussian distribution. Note that as the kurtosis, Eq. (37), of the distribution decreases the SMF dynamics becomes better. The best approximation is attained for the two-point distribution, Eq. (38), with the minimum kurtosis value of 1. We clearly see in Fig. 3 that the use of the two-point distribution strongly increases the timescale over which the SMF approach is predictive.

Figure 4 is the same as Fig. 3 except that the long time dynamics is shown. In Fig. 4(a), the exact solution exhibits oscillations with decreasing amplitude for long times; on the other hand the MF solution has beating-like oscillations with almost constant amplitude. This result and the fact that the one-body entropy in MF approximation remains zero at all times clearly show the well-known shortcoming of MF dynamics, that is the underestimation of dissipation and hence thermalization in quantum many-body systems. In Fig. 4(b), the SMF solution with two-point distribution follows the amplitudes of the oscillation of the exact dynamics for a longer time than the SMF solutions with Gaussian and uniform distributions, which become almost completely damped around $t=50\,\Delta^{-1}$. Figure 4(c) shows the similar behavior for the entropy $S$. The SMF solutions with Gaussian and uniform distributions attain the maximum entropy value around $t=50\,\Delta^{-1}$ whereas the SMF solution with two-point distribution is in better agreement with the exact entropy for longer times.

The dynamics in the strong coupling case, $v_{0}=0.5\,\Delta$, is indicated in Fig. 5. It is seen that the MF solution of $\langle D\rangle$ deviates from the exact one at $t=0.7\,\Delta^{-1}$ and SMF solutions start to deviate later at around $t=2\,\Delta^{-1}$. In Fig. 5(a) and 5(b), the difference between the SMF results with the three distributions is almost negligible. However, when the plots at $t=2.1\,\Delta^{-1}$ are zoomed in, as seen from the inset figures, the SMF solution with the two-point distribution turns out to be the best approximation to the exact solution and the SMF solution with the Gaussian distribution is the worst approximation. It is known that the validity time of MF as well as beyond-MF approximations decreases in inverse proportion to the coupling strength, $\Delta t_{\text{val}}\propto v_{0}^{-1}$, Lacroix et al. (2014); Polkovnikov (2003). Here, we observe the same behavior by comparing Figs. 3 and 5. Furthermore, as the coupling strength increases the difference between SMF evolutions with different initial distributions decreases.

In Fig. 6, we investigate the dependence of the SMF dynamics on the parameter $\chi$ that controls the weights of the real and imaginary parts of the stochastic matrix elements of the one-body densities $\rho^{\lambda}$. Three cases are considered: the first one corresponds to the equal weights case, labeled as R+I, for which $\chi=\overline{r^{2}_{ij}}=\overline{s^{2}_{ij}}=1/4$; the second one corresponds to full weight to the real parts, labeled as R, for which $\chi=\overline{r^{2}_{ij}}=1/2$ and $\overline{s^{2}_{ij}}=0$; and the third case is the opposite, full weight to the imaginary parts, labeled as I, for which $\overline{s^{2}_{ij}}=1/2$ and $\chi=\overline{r^{2}_{ij}}=0$. It is clearly seen that the equal weight case gives the best dynamics with all three distributions. When the imaginary matrix elements are set to zero, case R, the SMF entropy is overestimating the exact one and when the real matrix elements are set to zero, case I, the SMF entropy is underestimating the exact entropy. As discussed in Sec. II, the quality of the SMF approximation increases when the function $F$, given by Eq. (36), assumes smaller values. In Figs. 6(a) and 6(b), the SMF solutions with the Gaussian and uniform distributions are consistent with this result. However, the weight parameter $\chi$ should not affect the quality of the SMF evolution with the two-point distribution since $F$ is constant for this distribution, as seen from Fig. 1. In Fig. 6(c), it is observed that the equal weight case is still much better than unequal weight cases. This can be explained by observing the distributions of real and imaginary parts of the stochastic matrix elements. As an illustration, in Fig. 7, we show the distributions of the real $r_{\alpha\beta}$ and imaginary $s_{\alpha\beta}$ parts of a matrix element $\rho^{\lambda}_{\alpha\beta}$ of the stochastic one-body densities $\rho^{\lambda}$ at the times $t=0$, $t=2\,\Delta^{-1}$, and $t=4\,\Delta^{-1}$. In this figure, the single particle states are arbitrarily chosen as $\alpha=(+1,+1/2)$ and $\beta=(-1,+3/2)$. The initial distribution of the matrix elements is the two-point distribution. The same cases as in Fig. 6 namely R+I, R, and I, are indicated. It is observed that when the imaginary (real) parts of the matrix elements are set to zero, case R (I), the real (imaginary) parts in the figure assume three values with equal probability whereas, the imaginary (real) parts are zero at the initial time $t=0$. After a very short time interval, at $t=2\,\Delta^{-1}$, both the real and imaginary parts develop similar three-peak structures. However, after the same amount of time, at $t=4\,\Delta^{-1}$, opposite distribution structures for the R and I cases are observed with respect to the distributions at time $t=0$. At time $t=4\,\Delta^{-1}$, the real (imaginary) parts develop a narrow distribution around the value zero compared to the imaginary (real) parts that have a wide three-peak distribution for the case R (I). The case R+I exhibits almost symmetrical distributions for the real and imaginary parts at all times. These results clearly show how the dynamical correlations between the real and imaginary parts of the matrix elements build up in time. In particular, the real and imaginary parts are balancing each other for the equal weight case R+I, whereas the distributions are changing between narrow single-peak around zero and a wide three-peak distribution for the unequal weights cases R and I. In Fig. 6, it is seen that the entropies for the cases R and I start to deviate from the exact and R+I case solutions at very short times such as $t=2\,\Delta^{-1}$ for all the distributions. At this time, even the MF solution agrees well with the exact solution as seen from Fig. 3(a). Based on these observations and the fact that the correlations between the real and imaginary parts of the matrix elements of the one-body density are governed by the MF equation, we think that the quality of the SMF approach for different weights cases, such as the cases R and I, is unrelated to the choice of the initial distributions in the SMF approach. Hence, the equal weight case R+I should always be considered as a better approximation over unequal cases due to dynamical correlations, governed by the MF equation, between the real and imaginary parts of the matrix elements.

Figure 8 shows the distribution of the ”event“ expectation value, defined by Eq. (3), of the collective observable $D$ for the three probability distributions at three different times. The distributions of $\langle D\rangle^{\lambda}$ are continuous and Gaussian for the initial Gaussian and uniform distributions at all times. However, for the initial two-point distribution, the resulting distribution of $\langle D\rangle^{\lambda}$ is discrete. This result is directly related to the fact that the two-point distribution is discrete. It is seen that in a very short time interval $4\,\Delta^{-1}$ the result with the two-point distribution becomes also Gaussian since the correlations between the stochastic matrix elements of the one-body density, governed by the stochastic version of the mean-field equation (III.2), build up very fast. Figure 9 is the same as Fig. 8 except that we consider a system with $N=12$ particles in 24 single-particle states. Here, it is observed that the number of discete peaks increases linearly with particle number and also that the time interval after which a Gaussian distribution is obtained is even shorter than that for $N=6$ particles. The behaviors here are closely related to the central limit theorem, which states that the sum of many independent random variables tends toward a Gaussian distribution. Note that Figs. 8 and 9 show how fast the discrete initial distribution of $\langle D\rangle^{\lambda}$ transforms to a Gaussian distribution and hence how fast the correlations between the matrix elements of the one-body density build up. If the dynamical evolution of the distribution of $\langle D\rangle^{\lambda}$ in Figs. 8 or 9 is followed for longer times (not shown for brevity), it is observed that the Gaussian shape is preserved; however, centroids and widths of the distributions of $\langle D\rangle^{\lambda}$ for the two-point, Gaussian, and uniform distributions deviate from each other since the corresponding SMF means and variances of $D$ deviate from each other as well.

Similar results are reached for the SMF evolutions with different initial states. Here, we also demonstrate a result with a different initial state $|\Psi(0)\rangle_{2}$, which is explained in Sec. III.3. The corresponding solutions are illustrated in Figs. 10 and 11. Note that these figures are the same as Figs. 3 and 4 except that the initial state is different. The MF evolution with this initial state results in a sinusoidal oscillation with constant amplitude for the dipole moment $\langle D\rangle$, and the exact evolution demonstrates a beating oscillation with almost equal amplitudes for the beatings, as seen in Figs. 10(a) and 11(a). We see that taking a different initial state leads to the same conclusion that the SMF evolution with the two-point distribution leads to a better agreement with the exact evolution compared to the SMF evolution with the Gaussian or uniform distributions. In Figs. 10(c) and 11(c) the system does not reach the maximum entropy value $2\ln 2$ due to the symmetry of the initial state. This is also in accordance with observation of the beating oscillations of the exact evolution. SMF evolution with the two-point distribution is able to reproduce the beatings with decreasing amplitudes in time whereas SMF evolutions with Gaussian and uniform distributions get damped in a very short time.

IV Conclusions

In the SMF approach, expectation values of observables are obtained by statistically averaging over an ensemble of stochastic one-body densities. It is through matching these expectation values to the corresponding quantum expectation values of collective observables at initial time that provides information about the properties of the initial distribution of the matrix elements of the stochastic one-body densities. In the original formulation of the SMF approach, the matching of the expectation values is performed for the quantum means and variances of collective observables, which provides relations for the means and (co)variances (first and second moments) of the stochastic one-body density matrix elements that are then assumed to be Gaussian random numbers Ayik (2008). In the present work, we have investigated the properties of the stochastic one-body densities within the SMF approach by considering higher order moments. An expression for the fourth central moments of the matrix elements is derived which suggest that a probability distribution with minimal kurtosis for the initial density matrix elements provides the best approximation to the exact dynamics when the stochastic matrix elements are assumed to be statistically independent random numbers. The distribution with minimal kurtosis is the so-called two-point distribution, which can take only two values and hence is discrete. The SMF approach with the two-point distribution as well as the Gaussian and uniform distributions are applied to an exactly solvable model that is inspired by the Lipkin-Meshkov-Glick model that is familiar from nuclear physics Lacombe et al. (2016). It is shown that indeed the two-point distribution is the best approximation for the matrix elements of the stochastic one-body densities in the weak as well as strong coupling strengths. However, in the strong coupling case, the difference among the SMF evolutions with different distributions is almost negligible.

MF dynamics in its various forms such as TDHF and density functional theory is widely used in many fields of physics. However, MF dynamics is known to contain some drawbacks. Quantum correlations in collective observables are mostly missing hampering the use of the MF dynamics for long time evolutions. Furthermore, the dissipation of the collective motion is severly underestimated. The SMF approach cures these drawbacks to some extent by taking into account the initial quantum fluctuations. It can be shown by following ref. Lacroix et al. (2016) that SMF is equivalent to solving a BBGKY-like hierarchy of equations describing the coupling between the different moments of the one-body density. As a consequence, the kth moment is coupled to the (k+1)th moment. The dissipation in one-body space emerges from the coupling with the second moment that is also coupled to higher order moments. This situation is rather similar to the case with the standard BBGKY hierarchy of equations. The equations of motion of the different moments are independent of the assumption for the initial probability. However, obviously, the quality of the result itself is anticipated to depend on how accurately not only the second moments but also higher order moments of the collective observables are described by the initial probability. This is actually what we show here: that a better account of moments higher than the second moment improves the SMF evolution. Hence, it is observed that the long time correlations are much better reproduced with the two-point distribution than the Gaussian or uniform distributions, which result in overdamping and hence faster thermalization.

The discrete form of the two-point distribution results in a discrete initial distributions of the collective observables with the tips of the discrete values having a Gaussian shape. Due to correlations between the matrix elements supplied by the self-consistent SMF equation, the distribution of the collective observables reaches a Gaussian distribution very fast, which is in accordance with the central limit theorem. The discrete nature of the physical observables and the matrix elements of the initial stochastic one-body densities opens up the possibility to perform SMF dynamics with small numbers of events.

Acknowledgements.

B.Y. gratefully acknowledges useful discussions with Phuong Mai Dinh. This work was supported in part by the Scientific and Technological Research Council of Turkey (TUBITAK).

Appendix A Quantum moments of a one-body operator

Let us consider that a many-body system is described by a Slater determinant and that the corresponding one-body density is $\rho$. In terms of the natural basis $\{|i\rangle\}$, which satisfies $\langle i|\rho|j\rangle=n_{i}\delta_{ij}$ where the occupation numbers $n_{i}$ can take values 0 or 1, a one-body operator can be written as $A=\sum_{ij}A_{ij}a_{i}^{\dagger}a_{j}$. The first moment is the expectation value given by

[TABLE]

The second moment is obtained as

[TABLE]

where $\rho_{12}=\rho_{1}\rho_{2}(1-P_{12})$. Hence, the second central moment is given by

[TABLE]

The third moment of $A$ reads

[TABLE]

where $\rho_{123}=\rho_{1}\rho_{2}\rho_{3}(1-P_{12})(1-P_{13}-P_{23})$. The third central moment of $A$ reads

[TABLE]

where we used the fact that $n_{i}^{2}=n_{i}$ for a Slater determinant and $\Lambda_{ijk}$ is the symmetrized version of the term $n_{i}(1-3n_{j})(1-n_{j}n_{k})$ given by

[TABLE]

The fourth moment of $A$ reads

[TABLE]

where $\rho_{1234}=\rho_{1}\rho_{2}\rho_{3}\rho_{4}(1-P_{12})(1-P_{13}-P_{23})(1-P_{14}-P_{24}-P_{34})$. The fourth central moment of $A$ reads

[TABLE]

where the symmetrized terms are

[TABLE]

Note that the symmetrized terms $\Lambda^{(4a)}_{ijkl}$ and $\Lambda^{(4b)}_{ijkl}$ are obtained by realizing, from Eq. (60), the following equalities

[TABLE]

and

[TABLE]