Quantum and Classical Behavior in Interacting Bosonic Systems

It is understood that in free bosonic theories, the classical field theory accurately describes the full quantum theory when the occupancy numbers of systems are very large. However, the situation is less understood in interacting theories, especially on time scales longer than the dynamical relaxation time. Recently there have been claims that the quantum theory deviates spectacularly from the classical theory on this time scale, even if the occupancy numbers are extremely large. Furthermore, it is claimed that the quantum theory quickly thermalizes while the classical theory does not. The evidence for these claims comes from noticing a spectacular difference in the time evolution of expectation values of quantum operators compared to the classical micro-state evolution. If true, this would have dramatic consequences for many important phenomena, including laboratory studies of interacting BECs, dark matter axions, preheating after inflation, etc. In this work we critically examine these claims. We show that in fact the classical theory can describe the quantum behavior in the high occupancy regime, even when interactions are large. The connection is that the expectation values of quantum operators in a single quantum micro-state are approximated by a corresponding classical ensemble average over many classical micro-states. Furthermore, by the ergodic theorem, a classical ensemble average of local fields with statistical translation invariance is the spatial average of a single micro-state. So the correlation functions of the quantum and classical field theories of a single micro-state approximately agree at high occupancy, even in interacting systems. Furthermore, both quantum and classical field theories can thermalize, when appropriate coarse graining is introduced, with the classical case requiring a cutoff on low occupancy UV modes. We discuss applications of our results.


I. INTRODUCTION
Interacting bosons play an important role in many regimes, such as condensed matter systems of 4 He, superfluidity, to particle physics of Higgs, W-bosons, and to various cosmological applications, such as inflation, preheating, etc. One especially interesting application is to dark matter axions [1][2][3]. Since axions are very light (a QCD axion might have mass m a ∼ 10 −5 eV or so) their density and hence occupancy numbers are huge, so their bosonic character is important.
In general the behavior of interacting quantum particles is extraordinarily difficult to calculate, except in special circumstances. One important regime is the ultraquantum regime of bosons at very high occupancy. In this regime, one would certainly not use classical particle physics. Instead one appears to require the full treatment of the many particle Schrödinger equation. However, it is well understood, at least for free or nearly free systems, that in such a regime, a different type of classical approximation emerges: namely the classical field theory. Such a treatment is often employed for interacting systems, for example to understand the basic dynamics of Bose-Einstein condensates [4], to run lattice simulations of preheating [5,6], and to compute the evolution of the whole universe governed by scalar dark matter [7].
One might be concerned that this treatment is a little Electronic address: mark.hertzberg@tufts.edu too naive. It is possible that the classical field theory is no longer applicable on some time scale associated with interactions. Indeed there is inevitably at least one very noticeable difference between interacting quantum and classical systems: In the classical system, the state evolves uniquely. While in the quantum system, even states with rather well defined initial values for, say, the field and its conjugate momentum, will inevitably have these quantities become less and less well defined as wavefunctions spread. This spread is most dramatic in interacting theories, especially among systems that exhibit chaos. So even for bosons at high occupancy, several modes can interact with one another causing a huge spread in the wavefunction and calling into question any possible role of the classical approximation. Although the above argument was not quite the motivation, related concerns appeared in the recent work of Ref. [8] (earlier discussion appears in Refs. [9,10]), where it was claimed that quantum and classical interacting theories deviate on a dynamical time scale τ , even at high occupancy. At first sight this claim seems implausible, since the time scale τ is a property that can be defined purely within the classical theory. So it is very strange that classical physics should fail on a time scale independent of Planck's constant . In fact there exists an interesting literature on this subject, including the work of Refs. [11,12] where it is shown that agreement between quantum and classical thermal systems is in fact correct to order ∼ 2 in anharmonic systems. Other notable work includes Ref. [13] where it is shown that the Boltzmann equation and classical field theory are related at high occupancy. Thermalization in classical field theory has been studied in Ref. [14] and work on simulating quantum fields using classical physics includes [15].
Nevertheless numerical studies in Ref. [8] appear to indeed justify the claim that classical physics generically fails in interacting systems. The evidence presented was to consider a system that begins in a state of definite particle number |{N i } for a set of modes labelled i. This state was evolved by the Heisenberg equation of motion for a specific toy model for a choice of N i . On the other hand, a related classical problem was also studied, where the annihilation operatorsâ i were replaced by C numbers ψ i (as is usual for classical field theory), whose initial magnitudes were set to √ N i and phases set to θ i = 0. The classical values of N i (t) were shown to deviate spectacularly from the quantum expectation values of N i (t) on time scales longer than ∼ τ . Furthermore it was shown that the classical values of N i (t) kept oscillating significantly throughout the simulation, without settling down (see Figure 1 top left panel), while the quantum expectation values N i (t) settled down to near constant values at late times. It was thus suggested that the quantum system has relaxed to thermal equilibrium, while the classical system has not.
In this paper, we take a critical view of these conclusions. We point out that the appropriate comparison between quantum and classical is not to compare a quantum expectation value in a quantum state to the evolution of one very special classical micro-state (one with θ i = 0 initially), but to an ensemble of classical states. Indeed a classical micro-state, which of course oscillates forever in a closed frictionless system, should not be compared to a quantum coarse grained quantity, namely an expectation value. In fact both classical and quantum micro-states oscillate wildly forever, while averaged values in both theories can settle down.
We provide numerical and analytical evidence that the quantum expectation values of the occupancy numbers are approximately given by the classical ensemble average of classical micro-states with initial phases drawn randomly from a uniform distribution (for a sample, see Figure 1 top right panel). This is appropriate as the initial quantum states are chosen to be states of definite occupancy number and hence they have completely unspecified phases. This is in contrast to the work of Ref. [8] which focusses only on the special θ i = 0 case, which is not connected to the quantum state in any meaningful way.
We show that expectation values of classical states also settle down and approach the equilibrium values (see Figure 1 bottom right panel) in the same way the quantum expectation values do. Hence both classical and quantum treatments can exhibit thermalization if the thermodynamic limit is taken and appropriate coarse graining is introduced.
Finally we comment on an application to continuum field theory. In this case the classical ensemble average can be replaced by a spatial average of a single classi-cal micro-state by the ergodic theorem (assuming an underlying translationally invariant distribution). Hence in this way, even a single classical micro-state can approximate correlation functions of the quantum theory, despite the quantum spreading of the wavefunction. This "spreading" is captured by the ensemble or ergodicity. As an application, this means that the correlation length of dark matter axions is captured, at least approximately, by the classical theory. (Of course in the continuum field theory, one should be concerned about the UV behavior of the classical theory involving low occupancy modes, which require artificial regulation. But this is not relevant to the work of Ref. [8] which focusses only on a handful of finite frequency modes, all at high occupancy.) Our paper is organized as follows: In Section II we introduce a class of interacting models of bosons. In Section III we provide numerical results from evolving and ensemble averaging the classical evolution. In Section IV we compute the thermal averages and compare to numerics. In Section V we show analytically that the expectation values match. In Section VI we discuss the implications for correlation functions and correlation lengths in local field theories. Finally, in Section VII we present a discussion.

II. BOSONIC MODELS
Our primary motivation comes from systems of N tot bosons with a conserved particle number. This usually emerges in the non-relativistic limit (important examples include laboratory studies of 4 He due to conservation of baryon number and dark matter axions due to small annihilation cross sections). Furthermore, we will focus on 2 → 2 scattering processes as these tend to dominate in the non-relativistic theory. A set of relevant interactions are of the standard form Since we are interested in bosons at high occupancy, it is useful to pass to the second quantized language using creation and annihilation operatorsâ † i ,â i (where index i labels each mode with wave-vector k i ). For a discrete set of momenta the above Hamiltonian can be re-written aŝ where ω i = |k i | 2 /(2m) is the frequency of each mode. The potential is re-organized into some collection of coefficients encoded in the couplings Λ kl ij (related to the Fourier transform of V ). The conservation of momentum requires Λ kl ij to only be non-zero when the momenta associated with each index satisfy k i + k j = k k + k l . The hermiticity of the Hamiltonian requires Λ ji * lk = Λ kl ij . Note that the convention of eq. (2) is to allow overcounting of indices, so we can take without loss of generality Λ kl ij = Λ kl ji = Λ lk ij = Λ lk ji . In principle we would like to study a huge set of wavevectors to properly address the continuum theory. However, the numerical evolution is very difficult. It suffices for the present purposes to simply study a toy problem built out of only a handful of oscillators. We will follow the interesting and useful model of Ref. [8] (earlier introduced in Ref. [9]).
The details of the model are as follows: Only 5 modes are included, labelled by i = 1, . . . , 5, with k i ∝ i just a scalar. The non-relativistic dispersion relation is replaced by a linear dispersion relation ω = k. So the frequencies of the 5 oscillators are taken to be integer multiples of a fundamental frequency ω 0 as The only non-zero independent couplings are taken to all equal some overall strength of interaction Λ 0 as Λ 23 14 = Λ 24 15 = Λ 34 25 = Λ 13 22 /2 = Λ 24 33 /2 = Λ 15 33 /2 = Λ 35 44 /2 = Λ 0 . As a concrete example of initial conditions, Ref. [8] took the initial quantum state to have definite occupancy numbers of |{N i } = |12, 25, 4, 12, 1 . The coupling is chosen to be Λ 0 = ω 0 /10, while the inverse frequency ω −1 0 can be used a unit of time.
In Ref. [8] the system was evolved under the Heisenberg equations of motion following from the Hamiltonian in (2). They then outputted the expectation value of the occupancy operators N i (t) = â † i (t)â i (t) over time, finding that on a fairly short time scale the occupancy number expectation values settled towards their equilibrium values.

III. CLASSICAL BEHAVIOR
Let us now study this problem carefully within the framework of classical physics. The classical theory arises from replacing the annihilation operatorsâ i by C numbers ψ i . The corresponding classical equation of motion comes from Hamilton's equations using the fact that i ψ * is the momentum conjugate to ψ. This gives the set of ODEs The complex variables ψ i are specified by a magnitude and a phase that we can write as where N i (t) has an interpretation as an "occupancy number" in the quantum theory.
According to Ref. [8] the classical analogue of the initial quantum state of definite occupancy numbers |Ψ(t = 0) = |{N i } is to choose ψ i (t = 0) = √ N i with all phases vanishing θ i (t = 0) = 0. However, this is not the right analogue. Since the quantum state has an unspecified phase, then we should not try to connect it to a classical state of a specific value of θ i = 0. Instead the classical analogue is an ensemble of states with starting values of θ i drawn independently and randomly from the uniform distribution on the domain In Section V we will explain why this is the appropriate ensemble of initial states. These initial states can then each be evolved under the classical equations of motion (4). Finally we can output the expectation value of the modulus square of the field after ensemble averaging over s initial sets of phases, which we denote N i (s) ens . In the classical problem, we can scale out the overall number of particles N tot , which is only important in the quantum problem. We can define the fractional occupancy numbers n i (t) ≡ N i (t)/N tot , which satisfy the conservation law i n i (t) = 1, and we can define a dimensionless coupling parameter where h = 7 is the number of interaction terms in the Hamiltonian.
We have solved this system of classical equations numerically with results presented in Figure 1. In the top left panel we output the special case in which all the phases are set θ i (t = 0) = 0. This state is highly nongeneric, but was used as representing the classical theory in Ref. [8]. In the top right panel we output a much more generic case in which the θ i (t = 0) are chosen randomly. This evolution exhibits considerably more chaos than the special case. In the bottom panels we then pick s = 30 and s = 30000 random sets of initial θ i and average the solutions. Even for s = 30 in lower left panel we see somewhat less variation compared to the top right panel. For s = 30000 we are essentially in the limit in which we have achieved the true ensemble average  The ensemble averaged evolution is seen to be similar to the evolution of the quantum expectation value N i of the quantum state that was computed numerically in Ref. [8] (we refer the reader to those figures for comparison). Although the quantum case can only be computed efficiently for finite N i , we will show in Section V that in the high N i regime, the two answers will approximately agree Furthermore, at a fixed time t, we expect these to converge in a limit in which we take N tot → ∞ and Λ 0 → 0, while keepingΛ finite.

IV. THERMAL AVERAGES
The ensemble averaged classical occupancy numbers roughly approach some equilibrium value (presumably one needs to include a large number of oscillators to truly be in the thermodynamic limit and reach true equilibrium). One can enquire whether they approximate the thermal equilibrium valuesN i . Computing the exact thermal equilibrium is ordinarily difficult as we would need to perform statistical mechanics of a nonlinear system. However, in this special toy model, we can focus on the free theory HamiltonianĤ 0 which is the conserved momentum P tot in a theory with dispersion ω i = k i . (In other cases, one sometimes just approximates equilibrium by using the free theory. This can fail, such as for attractive interactions in the continuum field theory [16]).
In Ref. [8] it was claimed that thermal equilibrium in the classical theory would mean equipartition of (free theory) energy into each of the 5 oscillators:Ē 0i =N 1 ω 1 = N 2 ω 2 =N 3 ω 3 =N 4 ω 4 =N 5 ω 5 . Using eq. (3) this immediately gives a set of values forN i which disagree considerably with the late time values seen in Figure 1 lower right panel. However it is incorrect to use equipartition of energy in this case because it ignores the fact that the number of particles is fixed.
Instead the correct treatment of classical thermal equilibrium in this framework is to use the micro-canonical ensemble with momentum and number of particles fixed. So the macro-state is specified by M = {P tot , N tot }. The entire set of allowed micro-states are any µ = {N i , θ i } that satisfy the 2 constraints Thermal averages for the occupancy numbers are then given by the following integrals and similarly for other moments such asN 2 i .
For the example studied in Figure 1, the input momentum to particle number ratio is P tot /N tot = 127ω 0 /54. By carrying out these integrals we obtain the following set of thermal equilibrium values in the classical theory: {0.325, 0.293, 0.175, 0.118, 0.089}, (13) {σ ni } ≈ {0.145, 0.211, 0.132, 0.088, 0.066}, (14) wheren i ≡N i /N is the mean fractional occupancy in each mode and σ 2 ni ≡n 2 i −n 2 i is the variance. The thermal averages {n i } are indicated by arrows in the lower right panel of Figure 1. We see that they match the late time ensemble averages of the simulation quite well. We have also checked that they match a long time temporal average of a single micro-state quite well too. Also note that the size of the σ ni is of the same order asn i , so the fluctuations are large.
In the quantum theory, the thermal averages come from an almost identical calculation to eq. (12), but with integrals replaced by discrete sums. Hence it is obvious that these two approaches agree at high occupancy where the discrete sum may be approximated by an integral. So the quantum fluctuations are equally large. This is due to the spreading of the (occupancy basis) wavefunction and is captured by the spread in members of the classical ensemble.

V. QUANTUM TO CLASSICAL CONNECTION
Here we would like to explain why the classical ensemble average reproduces the quantum expectation values. To begin, consider the classical equations of motion (4). We denote the initial values as ψ i (t = 0) = φ i and we can solve this system of equations as a Taylor series in time. If we form the modulus square of ψ i , this provides the Taylor series of the classical occupancy N i (t), which takes the form where the coefficient of t p in the Taylor expansion is a polynomial in φ i , φ * j of order 2p + 2. On the other hand we can also compute the time evolution in the quantum theory. The Heisenberg equation of motion is similar to eq. (4), but with the replacement ψ i →â i and ψ * i →â † i . Let us denote the initial values for these operators asâ i (t = 0) =b i . In principle we can solve this system of equations, but it is much more difficult due to the fact that the creation and annihilation operators do not commutê However there is tremendous simplification in the high occupancy regime. In this case the typical values of these operators are large, in the sense that expectation values â † iâ i are large. So in this regime we do not need to be concerned about the δ ij correction of (16). This is a relative error of O(1/N i ). Hence we can freely commute these operators, which means that the structure of the solution reduces to precisely eq. (15), at each order, with where, as above, the coefficient of t p in the Taylor expansion is a polynomial inb i ,b † j of order 2p + 2, with coefficients matching the classical case (15). Now we would like to compute the expectation value ofN i (t) in an initial state of definite occupancy |Ψ(t = 0) = |{N i } . We can compute this expectation value term by term in the series (17) by using the standard ways in which creation and annihilation operators act on stateŝ Lets illustrate by taking the expectation value of a representative term in the expansion (17), namely a collection of 4 operators. The orthogonality of the states |{N i } leads to (20) In the classical case, the analogue is to perform an ensemble average over phases of the initial ψ i (t = 0) = φ i = √ N i e iθi(t=0) , so the corresponding classical term has expectation value Carrying out the integral gives exactly eq. (20). (It turns out that this particular contribution vanishes when inserted into eq. (15) or (17) for the toy model, but this illustrates the basic idea). This correspondence carries over to all the various higher order terms in the series (15) and (17) (many of which do not vanish). Hence the final results are indeed expected to agree in the large N i limit where we take N tot → ∞ and Λ → 0 holding their product finite. For finite N i , it is not obvious what is the time scale for departure of these theories, as the errors may grow in the higher order terms in the Taylor expansion. But the classical numerics presented earlier, with the comparison to the quantum numerics in Ref. [8], indicates that they quantitatively agree for times much longer than τ , and we believe they roughly agree for extremely long times as they exhibit very similar equilibria.

VI. CORRELATION FUNCTIONS
If we consider the continuum field theory, this connection allows us to express position space correlation functions of local operators in terms of the corresponding classical averages. For example, we may be interested in the two-point correlation function {N i }|ψ † (x, t)ψ(y, t)|{N i } , which ordinarily encodes the correlation length as the scale over which the correlations fall off. At high occupancy we can express this as Now in some sense the quantum and classical theories still differ in so far as the left hand side involves the expectation value of a single state, while the right hand side involves ensemble averaging many states. However, this is not necessarily so. The classical ensemble average may be replaced by a spatial average of a single classical micro-state ψ µ by the ergodic theorem (assuming an underlying translationally invariant distribution, as provided by inflation, for example) where the initial value of the classical field in k-space is taken to be ψ µ (k i , t = 0) = √ N i e iθi , where θ i is randomly chosen. Together eqs. (22,23) provide an important result: For a quantum micro-state, initially specified by |{N i } , and a classical micro-state, initially specified by µ = {N i , θ i }, the quantum and classical correlation functions (and correlation lengths) approximately agree at high occupancy for long times. Also, under certain circumstances, a temporal form of the ergodic theorem may be applicable too.

VII. DISCUSSION
In the case of axion dark matter, it was shown in Ref. [16] that the correlation length according to the classical theory is small because the interactions (gravity and self-interactions) are attractive rather than repulsive. The current analysis shows that the classical result for attractive interactions carries over directly to the quantum theory too. Hence axion dark matter does not lead to long range correlations. Instead it can (at least partially) thermalize and lead to the formation of Bose-Einstein condensate clumps, such as Bose stars. For details see Refs. [16][17][18][19][20][21][22] and for related discussions see Refs. [23][24][25][26][27].
Our new results indicate that the classical description of bosonic fields can be entirely adequate, even though quantum wavefunctions do spread appreciably in interacting systems. This has application to not only axions, but to preheating simulations, etc. Also, this behavior can be seen in several other familiar contexts. For example, if one considers interacting billiard balls on a frictionless table, the wavefunctions spread, especially after each collision, so the expectation values of each ball's position x i (t) settle down at late times, while the classical micro-state x i (t) oscillates wildly. It is understood that this doesn't prevent classical physics from remaining a useful description of billiard balls.
Another familiar example is that of fluids governed by the Navier-Stokes equation. At the level of the effective field theory, one can, in principle, quantize the fluid's density and momentum density to formulate a quantum theory. Here there is interesting nonlinear behavior, such as turbulence, that is captured accurately by the classical field theory. Again one does not need to be concerned that the quantum wavefunction of the fluid has spread out on time scales longer than some dynamical time.
In all these cases, one can, in principle, perform an ensemble averaging over some space of initial conditions and use the classical evolution to mimic the quantum expectation values. Moreover, when an ergodic theorem applies, some spatial or temporal average can simply be performed to capture this.
In practice, even this step is often unnecessary, however, since decoherence provides an effective collapsing of the wavefunction. So one can essentially utilize the classical theory with a single history, bearing in mind that one should not attempt to predict the future trajectory with detailed precision in chaotic systems, but only to represent the basic character of what an individual observer might see. Furthermore, for certain special states, such as a BEC, a single classical field configuration is usually accurate in describing its behavior, as the fluctuations around the condensate δψ are often small [4].
Finally, let us remark on a special class of initial quantum states, namely coherent states, which are often thought of as the most classical. In this case, we expect a single classical micro-state to match the quantum expectation value on a time scale that is parametrically ∼ τ lnN , as one expects ∼ lnN collisions for the small initial quantum uncertainty ∼ 1/

√N
to grow to be O(1) [28] as the system of interacting oscillators is chaotic. That there should be improved agreement between classical and quantum in the high occupancy limit is essentially guaranteed by the Ehrenfest theorem. Contradicting this well established theorem, Ref. [8] claimed that the time scale for agreement is still only ∼ τ . We believe this is an artifact of running simulations with occupancy numbers that are too small to see the ∼ lnN enhancement. Indeed in order to study the coherent state, Ref. [8] used mean initial occupancy numbers of |0, 12, 16, 0, 0 , whose average value isN = 5.6. This is not a particularly high occupancy number and so the parametric enhancement of ∼ lnN is only an O(1) change to the ∼ τ estimate. Instead one would need to study much higher occupancy numbers to clearly establish the logarithmic enhancement beyond the dynamical time scale τ . In any case, as is the main point of this paper, we believe ensemble averaging (or ergodicity) is still essential to mimic the quantum thermalization for times t τ lnN . We know that at late times these simple systems thermalize and exhibit the same thermal distribution in the high occupancy regime. So we are then assured to have agreement for both early and intermediate times, due to Ehrenfest theorem, and at late times, due to similar thermalization after averaging. For coherent states, we can imagine some procedure of drawing the starting values of N i and θ i from a distribution of relative widths O(1/ √ N i ) around their starting mean values and then ensemble averaging.