# First-principles quantum dynamics for fermions:

Application to
molecular dissociation

###### Abstract

We demonstrate that the quantum dynamics of a many-body Fermi-Bose system can be simulated using a Gaussian phase-space representation method. In particular, we consider the application of the mixed fermion-boson model to ultracold quantum gases and simulate the dynamics of dissociation of a Bose-Einstein condensate of bosonic dimers into pairs of fermionic atoms. We quantify deviations of atom-atom pair correlations from Wick’s factorization scheme, and show that atom-molecule and molecule-molecule correlations grow with time, in clear departures from pairing mean-field theories. As a first-principles approach, the method provides benchmarking of approximate approaches and can be used to validate dynamical probes for characterizing strongly correlated phases of fermionic systems.

## Introduction

The physics of interacting fermions is the basis of many of the most important phenomena in condensed matter physics, ultracold gases, and quantum chemistry. A fundamental issue is how the microscopic interactions at the quantum level give rise to collective and emergent effects in many-body systems. For many situations, particularly in condensed matter systems, static or equilibrium correlation functions are sufficient to connect theory and experiment, and sophisticated techniques have been developed to calculate and measure them [1].

Addressing similar questions in the domain of many-body dynamics, however, has limitations in condensed matter systems. Ultracold quantum gases, on the other hand, allow creation of highly controllable implementations of analogue many-body systems for which the dynamical evolution and correlations are directly accessible [2, 3, 4, 5, 6]. The purity and tunability of these ‘tailor-made’ analogue systems means that ultracold quantum gases are ideal for testing fundamental ideas in quantum many-body physics and are leading candidates for dynamical ‘quantum simulation’ [7, 10, 9, 8]. In order to make predictions from the underlying theory and to validate the potential simulators [11, 9], or to benchmark approximate approaches, a numerical simulation of the real-time dynamics is required. Similar requirements of exact simulation of many fermions arise in determining the quantum chemistry of complex molecular systems [12].

In this work we perform first-principles dynamical simulations of a fermion-boson model. We use a Gaussian stochastic method based on a generalized phase-space representation of the quantum density operator [13]. The fermion-boson model forms the underlying basis for a broad range of phenomena in condensed matter and ultracold atom physics. It was originally proposed in the context of high-temperature superconductivity [14], but in ultracold gases it corresponds to the theory of resonance superfluidity with Feshbach molecules [15]. The latter forms the basis of a two-channel model for describing the physics of the BCS-BEC crossover [16]. More recently, the fermion-boson model has been used for analyzing the decay of double occupancies (doublons) [17] in a driven Fermi-Hubbard system [18]. The particular situation that we simulate here corresponds to spontaneous dissociation of a Bose-Einstein condensate (BEC) of molecular dimers into fermionic atoms [19, 20], in which case the model provides the fermionic equivalent of parametric downconversion in quantum optics: the production of pairs of entangled particles.

The Gaussian phase-space method represents an extension to fermions of successful bosonic techniques [21, 22, 25, 26, 23, 24, 27, 28, 29, 30]. The essence of the method is the mapping of the density operator evolution onto a Fokker-Planck equation for a phase-space distribution, via a continuous Gaussian operator basis [13]. The evolving distribution is then sampled with stochastic differential equations (SDEs) for the phase-space variables [21, 22]. The mapping to the phase-space distribution is exact [21, 22, 23, 31] provided no boundary terms arise in deriving the the Fokker-Planck equation. In practice, such terms may develop after some simulation time, but in an easily verifiable way [23, 31], putting a well characterised upper limit to a useful simulation duration. Numerical signatures of systematic errors include: (i) the onset of spiking behaviour due to the presence of near-singular trajectories; (ii) sudden dramatic increase in the statistical uncertainties; and (iii) development of power-law tails in the probability distribution. All these signatures, well-known from the early studies of real-time dynamics of bosonic systems [23] and from equilibrium calculations for fermion systems using imaginary-time techniques [13, 31], carry over to the present simulations of fermion dynamics and are verified in the numerical examples that we present below.

As in any stochastic method, sampling error limits the precision of the results. However, unless the distribution develops power-law tails, indicated by the above-mentioned signatures, this uncertainty can be made arbitrarily small by increasing the number of trajectories.

From the physical point of view, the Gaussian phase-space method can be viewed as providing the quantum corrections, through additional stochastic terms, to different mean-field approaches. For example, with certain factorization assumptions [32], the method reduces to a time-dependent Hartree-Fock formalism. Furthermore, neglecting the stochastic terms recovers the approximate pairing mean-field theory (PMFT) [19, 33], to which we compare the phase-space results. While often accurate for determining particle number densities, the mean-field approach gives no direct information about higher-order correlations, and its accuracy is not known a priori. In contrast to this, the first-principles simulations presented here reveal significant development of higher-order correlations.

For the first application of the fermionic phase-space method to a multimode dynamical problem, we consider a uniform molecular BEC (MBEC) initially in a coherent state at zero temperature, with no atoms present. Assuming sufficiently low densities, we neglect -wave scattering interactions to simplify the treatment. The Hamiltonian of this fermion-boson model [14] is given by

(0) |

where labels the plane-wave modes and labels the effective spin state for the atoms. Even though we will present the numerical results for a one-dimensional (1D) system, we formulate the problem in the general case as the method is straightforward to use in higher dimensions. The fermionic number and pair operators are defined as and , with , while the bosonic molecular operator obeys . The atom-molecule coupling (invoked by a magnetic Feshbach resonance sweep or optical Raman transitions) is characterized by [33], where is the size of the quantization box, and mediates an effective interaction between the atoms. The first term, , contains the kinetic energy of the atoms (of mass ), while the detuning corresponds to the total dissociation energy imparted onto the system by the external fields.

Because of the symmetry between spins in the Hamiltonian, and the equal initial populations, we need only to consider the number operator for one of the spin states . An additional operator identity that follows from the Hamiltonian is

(0) |

which arises because the condensate to which the atom pairs are coupled is assumed to be homogeneous. One consequence of eq. (Introduction) is that the relative number of atoms with equal and opposite momenta is perfectly squeezed [20], i.e. with zero variance. It also means that the second-order atom-atom correlation function reduces to . Thus the atom-atom correlation function can be determined from the number density alone.

One effective approximate approach for treating the dynamics of dissociation is the PMFT [19, 33], which is obtained by assuming atom-molecule decorrelation and by replacing the molecular operator by a coherent mean-field amplitude, .

In this paper we solve the full Hamiltonian (Introduction) exactly, and in order to quantify deviations from the PMFT behavior we evaluate several correlation functions. The departures from Wick decorrelation are analyzed via the correlation coefficient

(0) |

which is unity within the PMFT.

To examine molecule-atom pair correlations and the second-order coherence of the molecular field, we define

(0) |

Again, within the PMFT, these will be unity. We may expect that, over time, correlations will develop between the molecular and atomic fields; the Gaussian phase-space simulations give exact quantitative accounts of these effects.

## Fermionic phase-space representation

As mentioned above, the essence of the Gaussian phase-space method is the mapping of the density operator evolution into a set of stochastic differential equations [13] that can be solved on a computer. For a large class of second-quantized system Hamiltonians, such as those containing no higher than quartic terms in field operators, such mapping introduces no additional approximations. For plane-wave modes, the quantum state governed by Hamiltonian (Introduction) can be mapped onto a complex phase-space of dimension : , with and . The phase-space equations are structurally similar to the Heisenberg equations for the corresponding operators, but the stochastic terms are not unique and can be modified by different choices of diffusion gauge [24, 13, 34]. This freedom allows the equations to be tailored to have different numerical properties. One specific set of Ito SDEs is [34]:

(0) |

where the derivative is with respect to a scaled time, , with . We have normalized the molecular field by its maximum (initial) value, , where is the initial number of molecules. The complex Gaussian noises obey . This form of eqs. (Fermionic phase-space representation) shows that with drift terms of order , the noise terms are , i.e. the noise and therefore non-mean-field corrections to correlations become more important for decreasing . In practice we convert the equations to Stratonovich form and integrate with a semi-implicit method. Stochastic averages of the variables give the first-order operator moments; normally ordered higher-order moments are obtained by averages of the corresponding Wick decomposition [13], e.g. . Note that the final, averaged moment will not satisfy Wick’s theorem for a general quantum state.

The stochastic sampling assumes a sufficiently bounded distributions, such that any boundary terms could have been neglected in obtaining the Fokker-Planck equation. Previous experience with bosons [23] and fermions [31] has shown that spikes and a sudden rapid growth of the statistical sampling errors in observables are seen when the tails of the probability distributions do not decay fast enough. Through comparisons with independent numerical solutions we confirm that the onset of such spiking behaviour and the rapidly growing sampling error signifies the limit of the useful simulation time. Although it can be controlled somewhat with gauges, [24, 13, 34], the finite simulation time is a well-known limitation of stochastic phase-space methods. In the numerical examples of figs. 1, 3 and 4 below, the simulation results are shown for time windows well below the spiking time. We simulate a sufficiently large number of stochastic trajectories to reduce the statistical errors to below the thickness of the relevant curves shown in the figures. In addition, we use identities such as (Introduction) and different gauges as a further demonstration that the simulations are exact up until the emergence of spiking behaviour.

Unlike quantum Monte Carlo approaches that are well-suited to calculation of exact ground-state properties or to simulation through imaginary time [35], the Gaussian stochastic method does not suffer from a ‘dynamical sign problem’ [36]. Other approaches for real-time simulations include the time-dependent density functional theory [37], although in practice this method is often restricted in accuracy by the need for exact functionals. Methods that use matrix-product-states based algorithms have been very successful for applications to one spatial dimension [38, 39, 40, 41, 42], however, as these methods require a truncated basis they do not fulfill the strict benchmarking criteria that a first-principles method can provide. An interesting direction in recent years has been the extension to fermionic systems of stochastic wavefunction approaches [43], which are similar in spirit to phase-space methods.

## Few-mode system

To confirm the validity of our numerical implementation of the phase-space method, we first independently solve a small system with molecules and atomic modes in a standard number-state basis. For this test system, with a bosonic number-basis truncation of , the Hilbert space has dimension . In fig. 1 we show the population in the momentum modes calculated using the phase-space method (with stochastic trajectories) and the number-state basis; we also illustrate the identity (Introduction) by calculating and plotting directly and comparing it with . The agreement between the two methods is excellent. The top two curves in fig. 1 illustrate the deviation of the PMFT prediction (dashed-dotted curve) from the exact calculation (solid curve) for the resonant mode .

In fig. 2, we show the population dynamics of the resonant mode, , together with the explicit behaviour of the statistical errors due to the stochastic sampling. The number of stochastic trajectories in this example (which otherwise is the same as in fig. 1) is small enough to render the sampling errors visible. At the same time, we have chosen the time window longer than the onset of spiking behaviour and the sudden dramatic growth of statistical errors to explicitly illustrate the signatures of systematic errors that limit the useful simulation time of the phase-space method.

To further evaluate the differences between treating the Hamiltonian (Introduction) exactly and using the approximate PMFT, we plot in fig. 3 (a) the correlation coefficient . Clear deviations are seen as time evolves; the deviations illustrate that in the exact treatment the following inequality holds, , whereas the PMFT prescribes an equality sign. Next, we consider the molecule-atom and molecule-molecule second-order correlations, and [see figs. 3 (b) and (c)]. Within the PMFT, both correlations are identically equal to 1. However, our exact results show that the molecule-atom correlation initially grows with time while the total atomic population grows. Then it changes to anti-correlation as the resonant-mode atoms start to re-associate. Meanwhile the effect on the molecular field from the atom interactions is revealed as it gradually loses its second-order coherence, albeit not by a significant amount for this few-mode system.

## Multi-mode system

We now use the phase-space method for simulating large 1D systems, with atomic modes and () molecules at densities m. In these cases, the number-state calculation is impossible as the dimension of the Hilbert space is enormous (). In fig. 4 we show the evolution of the number of molecules for three different cases. For the top curve, the initial number () is much larger than the number of available atomic modes, each of which hosts at most 1 atom due to the Pauli blocking. Accordingly, we see negligible depletion of the MBEC, which makes the relative size of the bosonic fluctuations very small. Hence, we do not observe significant deviations from the PMFT, including in the molecular second-order coherence, eq. (Introduction), which differed from 1 by less than in this case.

The situation changes for the bottom curve, for which is comparable with the number of atomic modes within the relevant width of the momentum distribution near ; we estimate this number [33] to be . In this case, we see strong molecular depletion and an increased role of bosonic quantum fluctuations so that the PMFT starts to show disagreement with the exact result. Admittedly the disagreement is still very small, implying that the predictions of the PMFT for total particle numbers can be rather accurate. The same is not true, however, for higher-order correlations, shown in the right column of fig. 3 for the same parameters as the lowest curve of fig. 4. Here, the large depletion of the MBEC and the increased role of quantum fluctuations are manifest – beyond the predictability of the PMFT – in strong higher-order correlations. The correlation coefficient clearly deviates from one, though to a lesser extend than in the few-mode system. The deviations of the molecule-atom and molecule-molecule correlations from and , on the other hand, are more dramatic.

The development of decoherence in the molecular field can largely be accounted for by the dephasing of atomic-molecular oscillations due to total number uncertainty. The frequency of the oscillations depends on the initial number of molecules; with a range of frequencies, the oscillations get out of phase and thus, for example, prevent complete disassociation of the molecular field from being seen in the average. As illustrated by the dashed grey curve in fig. 3 (f), this effect can be reproduced by an ensemble of mean-field trajectories with different initial numbers. Here we have used a Poissonian weighting to give the same number-distribution as the initial coherent molecular field in the exact calculation, with [44]. The large values of , which are possible with a state containing a superposition (or mixture) of a few low-occupation number states [45], occur at maximum depletion when the number-uncertainty is relatively large.

## Summary

We have demonstrated a successful application of a fermionic phase-space representation to first-principles quantum dynamics of a fermion-boson model. We simulated the coherent molecular dissociation to fermionic atoms and found significant higher-order correlations that cannot be accounted for by the approximate pairing mean-field theory. The knowledge of such correlations and the development of experimental probes to measure them provide the most accurate characterization of quantum many-body phases in strongly correlated systems.

The method is exact up until clear ‘spiking’ signatures emerge in the stochastic trajectories. For the present model, the useful simulation duration encompasses the time required for partial recombination of the atoms into molecules and for significant higher-order correlations to emerge. The accuracy during this useful simulation time was independently confirmed by comparison to number state calculations (for a small system) and by checking of conserved quantities.

Although we have here reported only on 1D simulations, we have also implemented 2D and 3D calculations and found that the method works reliably in higher dimensions. Extensions of the method to implement -wave scattering interactions will enable the study of non-equilibrium dynamics in a broader class of fermionic systems of current experimental interest, such as atomic Mott insulators in optical lattices and the BCS-BEC crossover problem.

The authors acknowledge support by the Australian Research Council.

## References

- [1] Vojta M., Rep. Prog. Phys., 66 (2003) 2069.
- [2] Greiner M. et al., Nature (London), 419 (2002) 51.
- [3] Dürr S., Volz T. Rempe G., Phys. Rev. A, 70 (2004) 031601(R).
- [4] Greiner M. et al., Phys. Rev. Lett., 94 (2005) 110401.
- [5] Öttl A. et al., Phys. Rev. Lett., 95 (2005) 090404.
- [6] Perrin A. et al., Phys. Rev. Lett., 99 (2007) 150405.
- [7] Bloch I. et al., Rev. Mod. Phys., 80 (2008) 885.
- [8] Buluta I. Nori F., Science, 326 (2009) 108.
- [9] Cazalilla M.A. Rigol M., New J. Phys., 12 (2010) 055006.
- [10] Rigol M., Dunjko V., Yurovsky V. Olshanii M., Phys. Rev. Lett., 98 (2007) 050405. Rigol M., Dunjko V., Olshanii M., Nature, 452 (2008) 854.
- [11] Trotzky S. et al., Nat. Phys., 6 (2010) 998.
- [12] Miller W. H., PNAS, 102 (2005) 6660.
- [13] Corney J. F. Drummond P. D., Phys. Rev. Lett., 93 (2004) 260401, Phys. Rev. B, 73 (2006) 125112, J. Phys. A: Math. Gen., 39 (2006) 269.
- [14] Friedberg R. Lee T. D., Phys. Rev. B, 40 (1989) 6745.
- [15] Kheruntsyan K. V. Drummond P. D., Phys. Rev. A, 61 (2000) 063816; Holland M. et al., Phys. Rev. Lett., 87 (2001) 120406; Timmermans E. et al., Phys. Lett. A, 285 (2001) 228.
- [16] Ohashi Y. Griffin A., Phys. Rev. Lett., 89 (2002) 130402,
- [17] Strohmaier N. et al., Phys. Rev. Lett., 104 (2010) 080401,
- [18] Jördens R. et al., Nature (London), 455 (2008) 204.
- [19] Jack M. W. Pu H., Phys. Rev. A, 72 (2005) 063625.
- [20] Kheruntsyan K. V., Phys. Rev. Lett., 96 (2006) 110401.
- [21] Drummond P. D. Gardiner C. W., J. Phys. A, 13 (1980) 2353.
- [22] Gardiner C. W., Handbook of Stochastic Methods (Springer, Berlin) 2008.
- [23] Gilchrist A., Gardiner C. W. Drummond P. D., Phys. Rev. A, 55 (1997) 3014.
- [24] Deuar P. Drummond P. D., Phys. Rev. A, 66 (2002) 033812, J. Opt. B, 5 (2003) S281.
- [25] Drummond P. D. Carter S. J., J. Opt. Soc. Am. B, 4 (1987) 1565.
- [26] Steel M. J. et al., Phys. Rev. A, 58 (1998) 4824.
- [27] Corney J. F. Drummond P. D., Phys. Rev. A, 68 (2003) 063822.
- [28] Savage C. M., Schwenn P. E. Kheruntsyan K. V., Phys. Rev. A, 74 (2006) 033620; Savage C. M. Kheruntsyan K. V. , Phys. Rev. Lett., 99 (2007) 220404.
- [29] Deaur P. Drummond P. D., Phys. Rev. Lett., 98 (2007) 120402.
- [30] Perrin A. et al., New J. Phys., 10 (2008) 045021.
- [31] Corboz P. et al., Phys. Rev. B, 77 (2008) 085108.
- [32] Rahav S. Mukamel S., Phys. Rev. B, 79 (2009) 165103.
- [33] Davis M. J. et al., Phys. Rev. A, 77 (2008) 023617,
- [34] Ögren M., Kheruntsyan K. V. Corney J. F., Comput. Phys. Commun. in press, arXiv:1008.0970v1., (2011)
- [35] von der Linden W., Physics Reports, 220 (1992) 53.
- [36] Egger R. Mak C. H., Phys. Rev. B, 50 (1994) 15210.
- [37] Runge E. Gross E. K. U., Phys. Rev. Lett., 52 (1984) 997; Verdozzi C., Phys. Rev. Lett., 101 (2008) 166401.
- [38] Cazalilla M. A. Marston J. B., Phys. Rev. Lett., 88 (2002) 256403.
- [39] Daley A. J. et al., J. Stat. Mech., (2004) P04005.
- [40] White S. R. Feiguin A. E., Phys. Rev. Lett., 93 (2004) 076401.
- [41] Kollath C. et al., Phys. Rev. A, 74 (2006) 041604.
- [42] Bañuls M. C. et al., Phys. Rev. Lett., 102 (2009) 240603.
- [43] Juillet O., Gulminelli F. Chomaz Ph., Phys. Rev. Lett., 92 (2004) 160401; Montina A. Castin Y., Phys. Rev. A, 73 (2006) 013618.
- [44] More specifically, we perform an ensemble of PMFT calculations, each one having a different initial average number of molecules , ranging over integer values between and . The results of each simulation are then averaged with a Poissonian weight factor . The resulting pair correlation is given by .
- [45] As a simple example, consider the state . It has a value of , corresponding to the peak in Fig. 3 (f) at , where .