Non-equilibrium evolution of Bose-Einstein condensate deformation in temporally controlled weak disorder

We consider a time-dependent extension of a perturbative mean-field approach to the dirty boson problem by considering how switching on and off a weak disorder potential affects the stationary state of an initially homogeneous Bose-Einstein condensate by the emergence of a disorder-induced condensate deformation. We find that in the switch on scenario the stationary condensate deformation turns out to be a sum of an equilibrium part and a dynamically-induced part, where the latter depends on the particular driving protocol. If the disorder is switched off afterwards, the resulting condensate deformation acquires an additional dynamically-induced part in the long-time limit, while the equilibrium part vanishes. Our results demonstrate that the condensate deformation represents an indicator of the generically non-equilibrium nature of steady states of a Bose gas in a temporally controlled weak disorder.


Introduction
All realistic physical systems inevitably involve a certain level of disorder due to the generic presence of an environment or a random distribution of imperfections. This is commonly understood as a nuisance which has to be coped with, especially in solid-state materials. However, a complementary and a more nuanced viewpoint is recently gaining the ground [1]. Namely, the disorder can lead to exciting, qualitatively novel phenomena that do not have any clean counterparts, such as the well investigated phenomenon of Anderson localization [2,3]. Thus, one has a potential of utilizing the disorder in order to engineer new states of matter, where the many-body localized phase is a recent prominent example [4,5], which is currently being debated [6]. This is in line with an ongoing quest to exploit the system environment as a tunable knob for quantum control [7][8][9][10][11][12][13].
Seminal studies of superfluid helium in porous vycor glass [14,15] instigated a theoretical interest in the so-called dirty boson problem that deals with the emergence of different phases of ultracold bosonic atoms in the presence of a static random potential. Starting with the experimental realization of Bose-Einstein condensates (BECs) in 1995 [16,17], there was a renewed theoretical effort in understanding the influence of disorder on samples of cold atoms [18]. Huang and Meng made significant progress in describing the thermodynamic equilibrium by developing a Bogoliubov treatment of a weakly interacting homogeneous Bose gas in a static delta-correlated random potential, where quantum, thermal and disorder fluctuations were assumed to be small [19]. After accounting for the disorder perturbatively they found that superfluidity persists in spite of it, although the superfluid depletion turned out to be larger than the condensate depletion. Their work was subsequently expanded by many others, using the original approach of second quantization [20][21][22][23] or employing the functional integral framework and the replica method [24][25][26]. The consideration of scenarios when the disorder correlation function has a non-zero correlation length σ showed that both condensate and superfluid depletions decrease with increasing σ, in particular for the case of a Gaussian [21,[27][28][29], Lorentzian [30], or laser speckle disorder [31,32]. Quite recently, the predictions of the Huang-Meng theory were experimentally confirmed by analyzing the cloud shape of a molecular 6 Li BEC in a disordered trap [33]. Furthermore, building on an equilibrium inhomogeneous Bogoliubov theory [34], the results presented in [35] clarified the relation between the condensate depletion and condensate deformation in the presence of a weak, possibly random, external potential. We adopt such a notion of the condensate deformation here, use it as a measure of the average effect of the disorder and generalize it to the nonequilibrium context.
Quench dynamics is an active field of research encompassing condensed matter physics and quantum information, with many intriguing applications to ultracold atomic gases [36]. In a recent experiment [37] quantum quenches of disorder in an ultracold bosonic gas in a lattice were used to dynamically probe the superfluid-Bose glass quantum phase transition at non-zero temperature. It is well known that quenches may lead to the appearance of nontrivial steady states. For instance, in our context, the authors of [38] studied an interaction quench of a 3D BEC in a static disordered potential and concluded that in the long-time limit the superfluid density tends to zero while the condensate density remains finite, i.e., that the quench dynamics enhances the ability of disorder to deplete the superfluid more than to deform the condensate. However, a quench is just a limiting case of a more general switch on and off protocol, which is a ubiquitous basic element of any experiment. In this work we consider a more generic scenario of switching the disorder on (off) during a certain time interval, while keeping the interaction strength constant, and to investigate the resulting dynamical behavior of the condensate deformation. Here we are, in particular, interested in the long-time limit of the emerging condensate deformation and in its tunability upon changing the respective system parameters.
To this end, the paper is organized as follows. In Sec. 2 we present a general mean-field theory for calculating the condensate deformation of a homogeneous Bose gas in a temporally controlled weak disorder. The relevant case studies involving the disorder switch on (off) scenarios are worked out and discussed in Sec. 3, followed by concluding remarks in Sec. 4.

General theory
We consider N identical weakly interacting ultracold bosons in a 3D box of large volume V , under the appropriate periodic boundary conditions. The weakly interacting theory of dirty bosons [19][20][21][22] demonstrates that Bogoliubov quasiparticles and disorder-induced fluctuations decouple in the lowest order. Thus, the leading correction due to the presence of a disorder potential is derivable from a mean-field theory [28,30]. Therefore, we suppose that at time t = 0 the system is in its equilibrium ground state, which is described by the macroscopic wave function Ψ 0 (x) that solves the stationary Gross-Pitaevskii equation The particle density n = N/V = |Ψ 0 (x)| 2 determines the equilibrium chemical potential µ 0 = gn. In the following we take, without loss of generality, a real-valued clean-case homogeneous wave function Ψ 0 (x) = √ n. The thermodynamic limit N → ∞, V → ∞ with n = const will be implicitly assumed at the final stage of all calculations.
At times t ≥ 0 an external disorder potential is switched on where u(x) is a random potential with the corresponding ensemble average . . . , and f (t) is a deterministic driving function such that f (0) = 0 and 0 ≤ f (t) ≤ 1. We assume that the disorder potential has zero ensemble average at every point, u(x) = 0, in order to eliminate the effects of a simple shift of the chemical potential. The two-point correlation function is supposed to be of the form so that homogeneity is restored after performing the spatial ensemble average. The disorder strength R(0) is related to the average height/depth of the hills/valleys of the random potential landscape. The disorder correlation length σ measures the average width of these hills/valleys, and R(x − x ) typically decays down to zero for distances |x − x | larger than several σ. In the k-space we have, correspondingly, The system evolution at t ≥ 0 is described by the time-dependent Gross-Pitaevskii equation, which reads in the frame rotating with the frequency We work in the regime where u(x) is small in comparison with all other energy scales of the system, and therefore we can treat its influence in a perturbative manner. The disorder slightly deforms the condensate and we make a perturbative ansatz for the wave function where |Ψ α (x, t)| = O(|u(x)| α ) denote perturbative corrections due to the disorder and obey the initial condition Ψ α (x, 0) = 0 for all α ≥ 1. Following Ref. [39], the difference between the disorder-averaged particle density |Ψ(x, t)| 2 and the density of the disorder-averaged stands out as a possible dynamical extension of the Bose-glass order parameter [38], in close analogy to the well-known Edwards-Anderson order parameter for spin glasses [40]. The last expression in (7) enables us to interpret it as the average particle density associated with disorder-induced condensate fluctuations. In agreement with the in-depth discussion of [35], it represents the ensemble average of the condensate deformation due to the disorder. Henceforth, we simply call it the condensate deformation. As we are going to show, its value even represents an indicator for the non-equilibrium feature of the system's stationary states.
Using the perturbative expansion (6), we obtain for the average particle density where the terms are grouped according to their respective perturbative order. On the other hand, the density of the disorder-averaged condensate is so that, up to the second order, the condensate deformation becomes Therefore, for calculations to that order, we only need the first perturbative correction Ψ 1 (x, t). The first-order coupled equations follow from (5) and read With the help of the Fourier transform and its inverse we get where ω k = 2 k 2 /(2m) is the free particle dispersion. In order to automatically incorporate the initial conditions Ψ ( * ) 1 (k, t = 0) = 0, we apply the Laplace transform and make the identification L[f ](s) ≡ f (s), for brevity of notation. With this, the differential equations (13) reduce to a system of algebraic relations which are solved by with the Bogoliubov dispersion Ω k = ω k ( ω k + 2µ 0 ). Using the inverse Laplace transform we finally find We note that, since u(k) = 0, it follows Ψ where we used (4) and K ( * ) (−k, t−t ) = K ( * ) (k, t−t ). In the next section we demonstrate that the last expression represents a time-dependent generalization of the seminal Huang and Meng equilibrium result [19].

Switch on -switch off case studies
In this section we apply the developed theory to two particular relevant scenarios. First, we focus our attention on an exponential switching on of the disorder potential, which leads to two physically distinct contributions to the stationary condensate deformation, which are generically of a non-equilibrium nature. Afterwards, we analyze the fate of both contributions upon switching the disorder off.

Switch on scenario
Let us first consider the case of an exponential introduction of the disorder potential via where τ > 0 determines the characteristic time scale. In this way one has f (0) = 0 and stationarity occurs for t τ since f (t) → 1 at large times. Hence, from (17) and (19) we obtain In the long-time limit we can safely neglect the exponentially decaying terms. In addition, the k-space integral of the remaining terms, that involve trigonometric functions, effectively represents a sum of infinitely many rapidly oscillating functions having incommensurate periods, which turns out to vanish as t → ∞. This leads to the stationary condensate deformation which notably consists of two physically different terms. The first one is the same for any switch on time and, thus, is not of dynamical origin. In fact, in the case of an adiabatic switching on of the disorder potential, which corresponds to the limit τ → ∞, only that term remains, so that we find This is exactly the part of the condensate deformation that corresponds to the equilibrium value [23,28,30,41]. The presence of an excess deformation in (22) signals the non-equilibrium feature of the achieved stationary state in general. Therefore, we arrive at a straightforward interpretation of the total stationary condensate deformation q τ in (22): the first term is the equilibrium part (23), while the second one is dynamically induced and depends on the switch on time scale τ and, thus, on the respective driving protocol. The dynamically induced contribution clearly ranges from zero in the adiabatic case up to the maximal value attained for a sudden quench of the disorder potential, which corresponds to the limit τ → 0. In the latter case, we have The above findings, which are applicable to any disorder correlation R(x), will now be illustrated for a relevant special choice of the spatially delta-correlated disorder Any effect of such a disorder represents the upper bound for every other correlated disorder scenario, because the influence of the disorder tends to decrease with increasing its correlation length. Such a disorder can even be realized experimentally via a random distribution of many neutral atomic impurities trapped in a deep optical lattice [42][43][44]. From the Huang-Meng theory [19] we know that the equilibrium condensate deformation (23) for the delta-correlated disorder (25) reads The general time-dependent result for the condensate deformation (21) is presented in Fig. 1 for three switch on times τ . The full blue curve corresponds to a rapid quench, i.e. for τ → 0. Following an initial increase the condensate deformation displays overshooting and eventually settles in a new stationary state. The green dot-dashed curve depicts the approach to the steady state in the opposite adiabatic regime, i.e. for τ → ∞. Here, the overshooting of the condensate deformation is absent for a slow enough switch on driving. Note that the stationary long-time limit can even be determined analytically from (22) using (25) for arbitrary time scale τ , and is given by where τ MF ≡ /µ 0 represents the characteristic mean-field time scale. The corresponding function is plotted in Fig. 2. From (27) we straight-forwardly deduce the special cases We again see that, after the adiabatic introduction of the disorder potential, the new deformed stationary state is reached, which corresponds precisely to the equilibrium result (26) [19]. Moreover, after a sudden quench of the disorder the system does not equilibrate at all, since the stationary condensate deformation is 5/2 times larger than the equilibrium one. The additional contribution of 3q HM /2 comes from the sudden quench itself and it represents the maximally possible dynamically generated condensate deformation for weak delta-correlated disorder.

Switch on -switch off scenario
Having identified the equilibrium and the purely dynamical contribution to the condensate deformation, we now aim at examining their fate upon switching the disorder off. Thus, we Figure 3: Exemplary switch on -switch off driving function (29) for τ = 100τ 1 = 10τ 2 .  (29), either adiabatically or through a rapid quench.
next focus on the following switch on − switch off scenario, which is given by the driving function and is depicted in Fig. 3. It corresponds to an exponential switching on with the characteristic time τ 1 , followed by a switch off with the decrease time τ 2 . The disorder disappears at t = 2τ . The time τ has to be much longer than τ 1 , τ 2 so that about t τ a transient stationary state, like in the previous scenario, is effectively reached. Hence, τ approximately determines the disorder duration time. For t 2τ a final stationary state is achieved. Using similar methodology as before, we find the final state condensate deformation First of all, we observe that the equilibrium contribution to the condensate deformation vanishes, as expected, once the disorder has finally disappeared. Furthermore, the resulting condensate deformation is dynamically accumulated during both the switching on and the switching off parts of the driving protocol. In the considered case, the dependence on the respective characteristic time scales τ 1 , τ 2 is identical for both parts. Hence, the time-reversed protocol would have led to the same end result. With this we conclude that the equilibrium condensate deformation is the only part that can be remedied in a reversible fashion in the switch on -switch off scenario. The previous finding is illustrated in Tab. 1 for the example of the delta-correlated disorder. Thus, the adiabatic switch on followed by the sudden quench down and the sudden quench up followed by the adiabatic switch off scenarios lead both to the same final condensate deformation of 3q HM /2. Furthermore, the largest dynamically generated deformation is three times larger than the equilibrium one (26) and it occurs by suddenly switching the disorder on and off.

Finite correlation length case
Finally, we analyze the effect of a finite correlation length σ on the condensate deformation. For this purpose we assume a Gaussian correlation function of the disorder,  (31) and on its correlation length σ. The dot-dashed, full, and dashed curves correspond to q τ /q HM being equal to 2, 1, and 1/2, respectively.
Experimental realization of such a disorder can be achieved by using laser speckles [33,45,46]. Since the dynamical effects are revealed already during switching on, we simply consider the driving function (20). The resulting condensate deformation (22) is plotted in Fig. 4 for a range of switch on times τ and correlation lengths σ. As expected, for σ ξ, where ξ = 2 /(2mµ 0 ) denotes the healing length, the behavior is approaching the delta-correlated case discussed in subsection 3.1. The increase of the correlation length σ leads to a drop of the stationary condensate deformation q τ below q HM and further towards zero for σ ξ, even for very rapid quenches of the disorder. This is in accordance with the general wisdom that the disorder influence decreases when σ becomes larger.

Conclusion and outlook
Our results clearly demonstrate that the condensate deformation is an indicator of the nonequilibrium feature of steady states of a Bose gas in a temporally controlled weak disorder. However, the exact nature of the observed steady states after switching the disorder on, and potentially off afterwards, still remains an elusive point, which warrants further investigation. Additional understanding could be gained, for instance, by proceeding along the route of Ref. [38]. Namely, it is a natural follow-up question to ask how much of the superfluid density would survive in the switching on and off scenarios considered here. Indeed, previous studies showed that superfluidity is more affected by disorder than the condensate itself both in equilibrium [19-22, 24, 47, 48] and in non-equilibrium [38]. Additional characterization of the emerging steady states would be gained by examining the time scales necessary to reach the respective stationary values of the condensate deformation and the superfluid fraction, in particular how they depend on the used protocols.
Our findings potentially have broader implications. The equilibrium and the dynamical contributions to the induced condensate deformation suggest that the same could generically hold for other physical quantities of interest, at least in the perturbative weak disorder regime. Even though we were dealing with the homogeneous Bose gas, we expect that in a harmonic trap the same qualitative pattern would emerge. As it was demonstrated that the Huang-Meng theory can be verified quantitatively by utilizing the homogeneous results within a local density approximation [33], the same might be achievable for the corresponding dynamical counterpart. Moreover, it is both an interesting and intriguing question how the dynamical condensate deformation would behave in the non-perturbative strong disorder regime, where the two aforementioned contributions might not be additive anymore. As this problem is analytically quite challenging, we hope that numerical simulations and experiments could shed the first light upon that situation.
Furthermore, the present work can be extended in multiple directions as the disorder driving could be of different types. One can, for instance, envision varying temporally some other feature of the disorder, such as its correlation length, instead of its strength. Alternatively, time-periodic driving has a great potential to expand the scope of the achievable non-equilibrium asymptotic states [49][50][51][52][53]. Also, there is an emerging and quite promising research direction of utilizing stochastic driving with a tunable time scale [54]. Moreover, one can ask, in the spirit of quantum control theory [55][56][57], whether a particular driving protocol could be designed in order to maximize or minimize the final condensate deformation.