Thermalization in a Quantum Harmonic Oscillator with Random Disorder

We propose a possible scheme to study the thermalization in a quantum harmonic oscillator with random disorder. Our numerical simulation shows that through the effect of random disorder, the system can undergo a transition from an initial nonequilibrium state to a equilibrium state. Unlike the classical damped harmonic oscillator where total energy is dissipated, total energy of the disordered quantum harmonic oscillator is conserved. In particular, at equilibrium the initial mechanical energy is transformed to the thermodynamic energy in which kinetic and potential energies are evenly distributed. Shannon entropy in different bases are shown to yield consistent results during the thermalization.


Introduction
Microscopic description of the thermalization has been a longstanding question in isolated quantum systems. The biggest enigma is due to the discordance between reversible microscopic laws and irreversible macroscopic phenomena. Among many fundamental issues, one common question is whether an isolated quantum system can reach thermal equilibrium, i.e., the state with maximum entropy [1][2][3][4][5][6][7][8]. Fortunately, ultracold quantum gases, which are both pure and controllable, provide an excellent platform to study the nonequilibrium dynamics for isolated quantum systems. In atomic Bose-Einstein condensate (BEC) experiments, Kinoshita et al. [9] showed no evidence of thermalization by pairwise collision, from the Tonks-Girardeau limit to the intermediate coupling regime. However, the dissipative motion of oscillating BEC in a disordered trap, done by Dries et al. [10], did manifest the thermalization [11].
The quantum harmonic oscillator is one of the simplest systems in physics and plays a central role in a wide variety of fields [36]. For example, quantum harmonic oscillator appears everywhere in quantum optics and its properties have been seen regularly in experiments. In ultracold atoms, one can use the Feshbach resonance technique to tune the s-wave scattering length to zero [37]. Inspired by the experiment of an oscillating BEC in disordered trap [10], we aim to study the thermalization in a disordered quantum harmonic oscillator. It will be shown in this system that the reversible microscopic quantum mechanics actually conceals the irreversible macroscopic phenomena of thermodynamics.
As mentioned earlier, works from over a decade ago have already discussed the related phenomena in disordered BEC in the mean-field regime [19,20,[30][31][32][33]35]. We nevertheless try to present an exact result for a "noninteracting" quantum oscillator in a disorder trap. In principle, the exact noninteracting harmonic oscillator could still be very different from the disordered dipole-oscillating BEC in the mean-field regime, even in the noninteracting limit. Our results show that equilibrium properties are very different between the current disordered noninteracting quantum oscillator and the disordered dipole-oscillating BEC in the mean-field regime [38]. More interestingly, it will be shown in the current quantum system that the equilibrium properties are very similar to those for a classical harmonic oscillator in the microcanonical ensemble. The kind of quantum-to-classical transition indicates strongly the "intrinsic decoherence" of the isolated system [39,40].
A possible scheme to this problem will first be proposed. In the theoretical simulation, we alternatively prepare an initial coherent state with a centroid velocity v 0 , which is an out-of-equilibrium state. Owing to the multiple scattering with the disorder, more and more (initial) mechanical energy will transform to the thermodynamic one. Once the mechanical energy is fully transformed to the thermodynamic one, the system reaches the equilibrium. During the whole process from nonequilibrium to equilibrium, total energy is conserved. Moreover, at equilibrium, total energy is evenly distributed between the kinetic energy and the potential energy associated with harmonic trap.
The paper is organized as follows. In Section 2, we propose the possible scheme to study the thermalization in a quantum harmonic oscillator with random disorder. Theoretical approach for the simulation is also introduced. In Section 3, we show the results of real-space density and momentum distributions when the system reaches the equilibrium. In Section 4, due to the effect of disorder, we show that mechanical energy is transformed to the thermodynamic one. In addition, at equilibrium thermodynamic energy is evenly distributed between kinetic energy and harmonic potential energy. In Section 5, Shannon entropy in different bases are studied for the thermalization process. Section 6 is a conclusion.

The Approach
To study the thermalization of a quantum harmonic oscillator with random disorder, our approach is to solve the time-evolution wavefunction ψ of the particle for the whole process from nonequilibrium to equilibrium. The corresponding properties can be accessed based on the solved ψ.
For simplicity, we consider a one-dimensional quantum harmonic oscillator in a disordered trap. The time-dependent Schrödinger equation (TDSE) describing the system is where m is the particle mass, ω is the trapping frequency, and V dis (x) is the Gaussian correlated disorder potential which satisfies the autocorrelation function where V D and σ D correspond to the strength and correlation length of the random disorder, respectively. Figure 1 illustrates the set-up of a possible experiment. Before the disorder potential is turned on, the system is in the ground state with the Gaussian wavefunction. At t = 0, the trap is abruptly displaced to the left, so the system gains energy and starts to move to the left. Afterwards, it oscillates. Note that at the same (t = 0), the disorder potential is turned on and due to the effect of it, the system will eventually come to equilibrium. Alternatively, for convenience in simulation, at t = 0 an initial coherent state with a velocity v 0 at the center is prepared for the system, where ψ g (x) = (mω/πh) 1/4 exp(−mωx 2 /2h) is the Gaussian ground state. After release, owing to the multiple scattering from the disorder, the wave packet starts to redistribute. When evolution time is long enough (t t th , t th is the thermalization time defined in Section 5), and the system approaches the equilibrium. Together with (2) and (3), we have numerically solved (1) for a long enough time to see the phenomenon of thermalization. In our simulation, ε 0 ≡hω and the trapping length l 0 ≡ √h /mω are taken as the units of energy and length.
Here, we detail the simulation method to solve the TDSE (1). We performed the fast Fourier transform (FFT) to calculate the integration and differentiation in space. The calculation domain is fixed at [−64π, 64π] with a grid of N = 12, 288 points. Moreover, the time integration is done by adaptive Runge-Kutta method of orders 4 and 5, built-in in Matlab software. The random Gaussian correlated disorder potential is taken as with x i the grid points. The generated random numbers A i are subject to A i = 0 and A 2 i = 1, and the normalization ∞ −∞ f (x) 2 dx = 1 is finally made to satisfy Equation (2). We have been assured that the wave function has not been contaminated by numerical noise (round-off error) even after long-time computation. One can also consider the non-Gaussian correlated disorders such as the sinc(x) = sin x/x correlated function considered in [30,34]. We have also done the calculation based on the sinc correlated disorder function and the results are found to be the same as the presented ones. It is also worth noting that all results obtained are from a single (fixed) random potential. The result does not average over noise at any points.

Equilibrium Distribution
It is important to note that "equilibrium" has a precise meaning in thermodynamics, namely, the systems is described with a Boltzmann distribution when it is coupled with a reservoir at temperature T. There is no "reservoir" in the current isolated system. Here, "equilibrium" means a state of balance or a stable situation in which different parts of energies will no longer exchange internally. As will be seen clearly, disorder plays the role to assist the energy exchange and lead to the "equilibrium" of the system. In this regard, the thermalization discussed in the present paper can be viewed as a kind of "intrinsic thermalization".
It is also important to note that a system connected to a random disorder cannot in general be regarded as an isolated one. The isolated system under consideration refers to the combination of the harmonic oscillator and the random disorder. The real and time-independent random Gaussian correlated disorder considered does fulfill the time reversibility of Equation (1) and, at the same time, results in the maximization of Shannon entropy (see Section 5).
With the solved dynamic wavefunction ψ(x, t), one is ready to see the real-space density distribution At t = 0 before release, the real-space density distribution is given by Once the system starts to oscillate, due to multiple scattering with the random disorder, the system will eventually reach equilibrium from nonequilibrium. Figure 2 shows the equilibrium real-space density distributions at t t th . Three cases of disorder strength V D = 50ε 0 , 100ε 0 , and 200ε 0 are shown. The disorder correlation length σ D = 0.01l 0 and the initial velocity v 0 = 50l 0 ω are fixed for all three cases. As shown, the three equilibrium distributions have almost overlapped with each other. It means that different finite disorder potential strengths lead to the same equilibrium distribution for the isolated quantum oscillator. Of equal importance, the distribution is very close to the microcanonical distribution for a classical harmonic oscillator. In a microcanonical ensemble of a classical harmonic oscillator with a given energy E 0 , the probability distribution in phase space is and thus the corresponding real-space density distribution is The singularities at the classical turning points correspond to E 0 = mv 2 0 /2 = 1250ε 0 . The transition from a quantum distribution to a classical one provides a clear evidence of quantum (intrinsic) decoherence for the current system.
It is also of interest to study the momentum distribution ρ(k, t). Momentum distribution is given by where ψ(k, t) is the Fourier transform of ψ(x, t). As at t = 0 before release, the real-space wavefunction is ψ( (3)), and the corresponding initial momentum distribution is then ρ(k, , the Fourier transformation of ψ g (x). In fact, ψ g (k) is also a Gaussian function. The initial momentum distribution is shown by the red curve in Figuer 3. Figure 3 shows only the case with V D = 50ε 0 .
When the system reaches equilibrium (t t th ), the numerically solved momentum distribution is shown by the blue curve in Figure 3, featuring two peaks at the high-k turning points. When comparing the t = 0 and t t th momentum distributions in Figure 3, it is important to note the following. (i) The case of t = 0 is asymmetric which implies a nonequilibrium state; the case of t t th is symmetric which implies an equilibrium state. (ii) Through the effect of random disorder, initial mechanical energy E 0 = mv 2 0 /2 has been transformed to thermodynamic energy at equilibrium, which is distributed among different k states in a wide range up to the cut-off turning points.  . Equilibrium momentum distributions ρ(k, t) for a quantum harmonic oscillator with random disorder. Red and blue curves correspond to the distribution at t = 0 (nonequilibrium) and t t th (equilibrium), respectively. Parameters used are V D = 50ε 0 , σ D = 0.01l 0 , and v 0 = 50l 0 ω.

Energy Distribution
In the present system, total energy is the combination of three parts: where K(t) = |h∂ x ψ| 2 dx/(2m) is the kinetic energy, U = (mω 2 x 2 /2)|ψ| 2 dx is the potential energy associated with harmonic trapping, and V dis = V dis (x)|ψ| 2 dx is the potential energy associated with random potential. For random potential, the associated V dis is typically small compared to K and U. Figure 4 shows the evolutions of K and U, as well as the total energy E for the case of V D = 50ε 0 , σ D = 0.01l 0 , and v 0 = 50l 0 ω.
Unlike the classical damped harmonic oscillator where total energy is dissipated, the total energy E of the current disordered quantum harmonic oscillator is conserved. The disorder does not store energy during the thermalization process. It plays the role only to redistribute the energy and result in the final energy distribution. Both K and U are oscillating at the beginning. When t just passes t th , the oscillation becomes relatively smaller, and when t t th , they approach the static limit. As a matter of fact, energy distribution is another indication for the system going from a nonequilibrium to an equilibrium state.
One important feature shown in Figure 4 is that when the system reaches equilibrium, thermodynamic energy E is evenly distributed in K and U. Thus, at equilibrium the virial theorem is satisfied: K = U = E/2. This is analogous to the case for a classical harmonic oscillator, the time-averaged K = U = E/2. Moreover, for the present system, the equilibrium temperature T eq can be obtained as k B T eq /2 = K(t t th ) = U(t t th ). In view of the static value shown in Figure 4, K(t t th ) U(t t th ) 625ε 0 , it is identified that T eq 1250ε 0 /k B . Figure 4. Energy evolution of a quantum harmonic oscillator with random disorder. Total energy is conserved in the whole process and kinetic energy and potential energy associated with harmonic trapping are evenly distributed at equilibrium. The parameters used are V D = 50ε 0 , σ D = 0.01l 0 , and v 0 = 50l 0 ω.

Entropy Evolution
Studies of thermalization often follow the von Neumann entropy whereρ = |ψ(t) ψ(t)| is the density operator. Thus, the production rate of S v is whereρ ≡ ∂ρ/∂t and the second term in (11) vanishes due to the conservation of total states. It is well known that S v (t) as well as dS v (t)/dt are basis-independent, i.e., invariant under unitary transformations. In terms of arbitrary bases, However, if both basis |i and |j are chosen to be the eigenstates ofρ, then or S v (t) is constant, which implies von Neumann entropy is not eligible to describe the time direction of an isolated quantum system [41,42]. Alternatively we consider the Shannon entropy [43][44][45] where the basis |i is arbitrary but should rule out the eigenstates of eitherρ orĤ. Unlike the basis-independent von Neumann entropy, Shannon entropy is basis-dependent. To acknowledge that the choice is arbitrary, we consider three choices. When the basis |i is chosen to be the position eigenstates, the corresponding Shannon entropy where ρ(x, t) = |ψ(x, t)| 2 is the real-space density distribution (see also (5)). Alternatively, if the basis |i is chosen to be the momentum eigenstates, the corresponding Shannon entropy is where ρ(k, t) = | ψ(k, t)| 2 is the momentum distribution (see also (9)). Figure 5 shows the evolutions of both S x (t) and S p (t). Taking into account all possible phase spaces, the result of sum of the two, S x (t) + S p (t), is also shown. While the parameters used, V D = 50ε 0 and v 0 = 50l 0 ω, are just a numerical trial, they are close to those of a real experimental set-up with V D = 50.9ε 0 and v 0 = 37.5l 0 ω [10]. In addition, we consider a shorter disorder correlation length σ D = 0.01l 0 . All three cases consistently show the increase of entropy for the thermalization process, and they maximize (saturate) when the system approaches equilibrium. Therefore, the use of Shannon entropy is suitable to describe the thermalization process in the current isolated quantum system. In Figure 5, the black dotted line corresponds to the Shannon entropy for a classical harmonic oscillator in a microcanonical ensemble. It is obtained by substituting the real-space density distribution ρ c (x) in Equation (8) into S x in (15). In view of Figure 5, the thermalization time in this case after which the system reaches equilibrium can be unambiguously identified to be t th = 50(1/ω). To have a better understanding on how the thermalization time depends on the experimental set-up, in Figure 6 we study the effects of V D and v 0 on S(t) = S p (t) + S x (t). Figure 6a shows the evolution of S(t) for three disorder strengths: V D = 50ε 0 , 100ε 0 , and 200ε 0 with fixed v 0 = 50l 0 ω and σ D = 0.01l 0 . The thermalization times for the three cases are identified to be t th = 50(1/ω), 25(1/ω), and 12.5(1/ω), respectively. Therefore, the larger the V D is, the faster the system thermalizes. More interestingly, the effect of V D on the thermalization rate (1/t th ) is linear, 1/t th ∝ V D . The linear behavior seems to be quite general in all regimes (noninteracting and interacting) in an isolated quantum system.  Figure 6b shows the evolution of S(t) for three initial velocities: v 0 = 50l 0 ω, 25l 0 ω, and 12.5l 0 ω with fixed V D = 50ε 0 and σ D = 0.01l 0 . The thermalization times for the three cases are found to be t th = 50(1/ω), 25(1/ω), and 17(1/ω), respectively. Thus, t th goes monotonically with v 0 . In an interacting Bose condensate with disorder, it has been found that when v 0 > c (c being the sound velocity), thermalization will develop, while when v 0 ≤ c, thermalization is hardly developed or will not develop [46]. It thus suggested in an interacting system that the thermalization is intimately related to the generation of elementary excitations [10]. In the present noninteracting system, in contrast, the situation is very different. There are only free real particles without the quasiparticles. Due to the effect of disorder, the system tends to thermalize faster with a smaller v 0 as in this case lesser initial mechanical energy needs to be transformed to the thermodynamic one.
To check the applicability of the Shannon entropy, we also consider the basis of the generalized coherent states, |i = |α, n , where [47] x|α, n = 1 Here, H n (x) are Hermite polynomials,x = α, n|x|α, n = √ 2l 0 |α| cos(ωt − θ), p = α, n|p|α, n = −( √ 2h|α|/l 0 ) sin(ωt − θ), and θ is some arbitrary phase. The corresponding Shannon entropy is given by For simplicity, for this basis we choose the case with t = 0,x = 0, andp = √ 2h|α|/l 0 ≡ mv 0 . In this case, Equation (17) results in which is identical to our initial wavefunction shown in (3). Thus, the initial entropy is expected to be zero in the case of a single basis function. Figure 7 shows the simulation result of S α (t) with v 0 = 12.5l 0 ω. We only show the case of relatively smaller v 0 , as it takes much longer computing time to calculate the cases of higher v 0 . For comparison, we also show S(t), which was already shown in Figure 6b. While the Fermi-Pasta-Ulam-Tsingou (FPUT) recurrence effect [48,49] for the current oscillating system is significant, the maximization process of S α (t) is still evident. The key is to follow local minima of S α (t). Of more importance, the thermalization time is predicted to be t th 17(1/ω), consistent with that predicted in S(t). This concludes that Shannon entropy is an eligible one to study the thermalization in an isolated quantum system. Figure 7. Evolutions of Shannon entropy S α (t) for a quantum harmonic oscillator with random disorder (red curve). The parameters are v 0 = 12.5l 0 ω, V D = 50ε 0 , and σ D = 0.01l 0 . As the FPUT recurrence effect is significant, the key is to follow local minima of S α (t). For comparison, S(t) (blue curve) is also shown for the same set of parameters (see also in Figure 6b). The thermalization time is identified to be t th 17(1/ω), consistent in both cases.

Conclusions
In conclusion, we propose a simple, possible scheme to study the thermodynamics in a quantum harmonic oscillator with random disorder. We have numerically shown that due to the effect of random disorder, the system can undergo from nonequilibrium to equilibrium state, i.e., thermalization. During the thermalization process, total energy of the system is conserved and at equilibrium the kinetic and potential energies are evenly distributed. Shannon entropy in different bases are studied and shown to maximize during the thermalization process.

Conflicts of Interest:
The authors declare no conflict of interest.