Parametric Excitation of a Bose-Einstein Condensate: From Faraday Waves to Granulation

We explore, both experimentally and theoretically, the response of an elongated Bose-Einstein condensate to modulated interactions. We identify two distinct regimes differing in modulation frequency and modulation strength. Longitudinal surface waves are generated either resonantly or parametrically for modulation frequencies near the radial trap frequency or twice the trap frequency, respectively. The dispersion of these waves, the latter being a Faraday wave, is well-reproduced by a mean-field theory that accounts for the 3D nature of the elongated condensate. In contrast, in the regime of lower modulation frequencies we find that no clear resonances occur, but with increased modulation strength, the condensate forms an irregular granulated distribution that is outside the scope of a mean-field approach. We find that the granulated condensate is characterized by large quantum fluctuations and correlations, which are well-described with single-shot simulations obtained from wavefunctions computed by a beyond mean-field theory at zero temperature, the multiconfigurational time-dependent Hartree for bosons method.

I.

INTRODUCTION
Spatial patterns frequently emerge in driven fluids in a variety of contexts, including chemistry, biology, and nonlinear optics [1]. Instabilities in these systems can generally be categorized as Rayleigh-Bénard convection, Taylor-Couette flow, or parametric surface waves. One of the earliest and best known examples of the latter type are the surface waves found by Faraday when a vessel containing a fluid was shaken vertically [2]. The resulting standing wave patterns that appear on the fluid surface arise from parametric excitation of collective modes of the fluid. The Faraday experiment has been repeated in various geometries where complex patterns were observed for small driving amplitudes [3]. Chaotic behavior, such as sub-harmonic bifurcation, is seen when the drive amplitude is strong [3][4][5][6] and this behavior has been connected to the onset of turbulence [7].
A model of the Faraday instability has been developed for an inviscid fluid in which the underlying hydrodynamic equations have been linearized [8]. The linearized dynamics are described by a Mathieu equation,ẍ + p(t)x = 0, where x is the displacement, p(t) = Ω 2 (1 + cos(ωt)) is the drive, representing a parametrically driven (undamped) harmonic oscillator with a natural frequency Ω, drive frequency ω, and drive amplitude . Solving the equations using a Floquet analysis results in a series of resonances at ω = 2Ω/n, where n is an integer [9].
Superfluids are particularly interesting in the context of Faraday waves because the damping of collective modes can be much smaller than in normal fluids, and because patterns may dissipate by the formation of quantized vortices in two or three dimensions. Several theoretical works have investigated Faraday waves in Bose-Einstein condensates (BECs) of atomic gases [10][11][12][13][14][15][16]. To our knowledge, only three experiments on Faraday waves in superfluids have been performed, one in which a vessel containing liquid 4 He is vertically shaken in a way similar to the original Faraday experiment [17], a pioneering experiment in which Faraday waves were excited by modulation of the transverse trap frequency, ω r , of an elongated BEC of Rb atoms [18], and another in which a non-destructive imaging technique was used to observe Faraday waves in a BEC of Na atoms [19]. In the BEC experiments, the transverse breathing mode, excited at a frequency of 2ω r , strongly couples to the density, and hence, to the nonlinear interactions of the condensate. This coupling produces the longitudinal sound waves responsible for creating Faraday waves [18,19]. The spatial period of the Faraday waves was measured as a function of ω, and the response to the strength of the drive was investigated [18]. In a related BEC experiment, modulation of the scattering length in a regime of large modulation amplitude and frequency resulted in the stimulated emission of matter-wave jets from a 2D BEC of Cs atoms [20].
In this paper, we report measurements characterizing the response of an elongated BEC to direct modulation of the interaction parameter using a Feshbach resonance [21][22][23]. For drive frequencies near the first parametric resonance (ω near 2ω r ), we observe robust linear spatial patterns characterized by a spatial period λ F (ω) consistent with Faraday waves. We also observe the response of the gas to the next lowest "resonant" mode (ω near ω r ) [15].
We have also investigated how λ F depends on the interaction strength. These measurements are compared with a theory that fully incorporates radial, as well as axial dynamics using a variational method [15], and, as we will show, the agreement is excellent.
We also explore a different modulation regime, both experimentally and theoretically, where ω is far from any trap frequency. The behavior in this regime is distinctly different; no clear resonances are observed, and much larger and modulation times are needed to obtain a significant response. The response is not regular in this regime, and no clear patterns emerge; rather, modulation produces a series of irregular grains.
Granulation is found in a variety of systems extending over many length and energy scales [24,25]. In quantum gases, granular states have been discussed previously in the context of perturbed atomic BECs and explored theoretically using a mean-field approach [26,27]. Granular states have been defined to have the following properties [26]: i) they are dynamical quantum states where particles cluster in higher density grains interleaved by regions of very low density, ii) the spatial distribution of grains is random, and iii) the grain size is variable and of a multiscale nature.
Our theoretical description uses the multiconfigurational time-dependent Hartree method for bosons (MCTDHB) [28,29]. MCTDHB captures many of the salient experimental observations and goes systematically beyond a mean-field description obtained from the Gross-Pitaevskii equation. The discrepancies between the Gross-Pitaevskii mean-field description and both the experimental observations, and our MCTDHB results hint that granulation emerges concurrently with many-body correlations.

II. FARADAY WAVES
In our experiment, we confine a gas of up to 8 × 10 5 7 Li atoms in a single-beam optical dipole trap and cool them to well-below T c , the transition temperature for Bose-Einstein condensation [22]. This configuration results in a highly elongated cylindrical trapping geometry whose corresponding axial and radial harmonic frequencies are ω z = (2π)7 Hz and ω r = (2π)475 Hz, respectively. The atoms are optically pumped into the lowest ground state hyperfine level, |F = 1, m F = 1 , where their s-wave scattering length may be controlled using a broad Feshbach resonance located at 737.7 G [30][31][32][33]. The magnetic field is sinusoidally modulated according to B(t) =B + ∆B sin(ωt), resulting in a modulated scattering length, a(t). The modulation amplitude ∆B, modulation time t m and hold time t h following t m are varied for each value of the modulation frequency ω, as necessary to produce a Faraday pattern with similar contrast. After t h , we take a polarization phase contrast image [34] with a probe laser propagating along the x-axis, perpendicular to the cylindrical z-axis of the trap. These images provide column density distributions that we integrate along the y-axis to obtain line density profiles. We apply a fast-Fourier transform (FFT) to these profiles in order to determine the spectrum of spatial frequencies exhibited by the BEC following modulation.
A typical image of a single experimental run is shown in Fig. 1(a). In this example, ω = (2π)950 Hz is resonant with the Faraday mode at ω = 2ω r . A surface wave is generated after t m = 5 ms of modulation followed by t h = 20 ms. The FFT, shown in Fig. 1(b), features a single dominant peak corresponding to a spatial period of λ = 10 µm. Figure 2 shows the spatial period of the observed structure as a function of ω. Typically, t h = 0 ms and 20 < t m < 40 ms, with the exception of ω = ω r and ω = 2ω r . Near these resonances, the modulation time was kept short, t m = 20 ms and t m = 5 ms, respectively, followed by t h = 20 ms. The blue data points correspond to the spatial period of the primary peak in the FFT spectrum. Except for the point at ω = (2π)475 Hz, the period monotonically increases with decreasing ω. The blue line in Fig. 2 is the result of a 3D variational calculation of λ F [15], which fits the data well. We have verified that the standing wave surface wave amplitude oscillates at ω/2 for ω near 2ω r , consistent with its identification as a Faraday wave, which is excited parametrically.
The excitation at ω = (2π)475 Hz = ω r is not a sub-harmonic of the Faraday mode at 2ω r , dominates when ω is tuned to resonance at ω r , producing the observed primary peak.
but rather the next lowest mode in the infinite series of modes, identified as the "resonant" mode in Ref. 15. In addition to having a different dispersion relation, this mode is also weaker, and therefore more difficult to excite, except exactly on resonance, ω = ω r , where the growth rate of the resonant mode exceeds that of the Faraday mode [15]. A similar excitation at ω r was previously reported [18]. The theoretical calculation of the period of this mode is indicated in Fig. 2 by the red line, λ R [15].
We find that as ω is tuned away from 2ω r = (2π)950 Hz, a larger modulation amplitude ∆B and modulation time t m are required to obtain a pattern with similar contrast. For example, Fig. 3 displays the spectrum for ω = (2π)200 Hz, for which ∆B = 35 G, t m = 20 ms, and t h = 20 ms. Two peaks dominate the spectrum: the primary peak at lower spatial frequency, and a secondary peak at roughly twice this spatial frequency. These secondary peaks only appear for ω (2π)400 Hz, and are identified by the red data points in Fig. 2.
The appearance of the next lowest mode depends on being sufficiently near its resonance frequency at ω = ω r , and far enough off-resonant with the Faraday mode at ω = 2ω r that it does not dominate the FFT spectrum. We have looked for additional modes in the data, but the FFT spectrum is dominated by the off-resonant response to the 2ω r and ω r resonances, and we are unable to observe any resonances below ω r . A comparison of the period of these secondary peaks with the theoretically calculated solid red line indicates that they correspond to the resonant mode λ R .
We also explored a more impulsive regime, with short t m , and where ω is kept within 10% of the Faraday resonance at ω = 2ω r . In this case, with short t m , we find that the wavelength of the resulting Faraday pattern is constant, independent of ω.
The Faraday period also depends on the strength of the nonlinearity, as shown in Fig. 4, where both the measured and calculated [15] values of λ F are plotted vs. the interaction parameterāρ, whereρ is the line density obtained by integrating the column density along the transverse direction. The measured period is consistent with the 3D theory from Ref [15].
We have also explored the dynamics for the emergence of the Faraday pattern and its persistence following a short modulation time interval of t m = 5 ms near 2ω r . Figure 5 measurements of the condensate length vs. t h . Figure 5(b) shows the axial Thomas-Fermi radius during the same t h interval. It shows that a low frequency collective mode is excited by the coupling to the modulated nonlinearity. The parameters of this condensate place it between the 1D mean-field and the 3D cigar regimes [35]. In the 3D Thomas-Fermi limit, the lowest m = 0 quadrupolar mode for an elongated condensate has a frequency of 5/2 ω z while in the 1D limit the collective mode oscillates at √ 3ω z [35][36][37]. For ω z = (2π)7 Hz, the corresponding period for this mode is, therefore, ∼90 ms, which is close to the observed oscillation period of 95 ms. We find that the Faraday pattern is suppressed during axial compression, but subsequently revives as the condensate returns to its original size. The phase of the two oscillations, the FFT amplitude and the Thomas-Fermi radius, do not exactly coincide. We attribute this observation to the delay in the initial growth of the Faraday pattern. We have determined experimentally that the frequency of the collapse and revival of the Faraday pattern scales with the axial trap frequency. A similar collapse and revival of the Faraday wave was previously observed [19].

III. GRANULATION
A Faraday pattern is not observed for low frequency modulation, for which ω ω r . We find that as ω is reduced both modulation time t m and modulation amplitude ∆B must be increased in order to observe any change. As these parameters are increased, more spatial frequencies contribute (see Fig. 3), and as t m and ∆B are increased further, we observe random patterns spanning a broad spatial frequency range, resembling grains [26,27]. We do not observe a significant thermal fraction before, nor after modulation, and therefore we attribute the observed granular patterns to quantum fluctuations and use a theory applicable to pure states.
In Fig. 6(a) we show experimental images and compare them to Gross-Pitaevskii (GP) simulations. Note also that the axial and radial trap frequencies in this section are ω z = (2π)8 Hz and ω r = (2π)254 Hz, respectively. We observe that granulation is remarkably persistent in time after the modulation is turned off, and that its structure is random between different experimental runs. GP simulations for similar parameters are shown in  The GP ansatz is a product of one single-particle state φ GP : Ψ GP ∼ N k=1 φ GP (r k ). This is a "mean-field state" because all particles in the many-body system occupy the single-particle state φ GP (r). A GP product state cannot describe correlations, where the properties of one or several particles in the many-body system depend on the properties of other particles in it. We go beyond the mean-field GP theory by employing the multiconfigurational timedependent Hartree for bosons method (MCTDHB or MB), which can account for manybody correlations. The MCTDHB ansatz incorporates all possible configurations (n 1 , ..., n M ) of N particles in M single-particle states, |Ψ = n 1 ,n 2 ,...,n M C n 1 ,n 2 ,...,n M |n 1 , ..., n M . The MCTDHB ansatz can therefore self-consistently describe correlations in the many-body We simulate the in-situ single-shot images [39,40] from the wavefunctions obtained with MCTDHB for the various experimental parameters and for M = 2 modes (see Supplemental    In both the experiment (Fig. 8(a)) and MB theory ( Fig. 8(b)), we find that when ω < ω c the condensate is practically uncorrelated, as evidenced by C (2) (z, z ) ≈ g (2) (z, z ) ≈ 1.
However, when ω > ω c we find that the relatively constant correlation plane evolves into smaller correlated and anti-correlated regions, as shown Fig. 8(c) for the experiment and respectively (see Supplemental Materials [38]), which are plotted in Fig. 9(b). These may be used as a measure of the departure of our MB model from mean-field states. Many-body systems, where multiple eigenvalues of the 1 st order RDM are macroscopic (ie. of order N ), are termed fragmented [47,50]. At zero excitation only n  Both observations, the emergence of fragmentation and the loss of 2 nd order coherence, underscore that the granulation of Bose-Einstein condensates is a many-body effect. The system thus cannot be described by a mean-field product state any longer and has left the realm of GP theory. Although the transition to fragmentation is not sharp -the natural occupations n (1),(2) i take on continuous values -it is well established at sufficiently large ω. Granulation features randomly-distributed variably-sized grains of atoms which can be observed in single shot images. Fragmentation, or depletion, on the other hand, is characterized by the reduced density matrix and its (macroscopic) eigenvalues and is not necessarily accompanied by granulation of the density [47,50]. In our close-to-one-dimensional setup, we observe granulation to emerge side-by-side with fragmentation.
The dynamical evolution, as calculated from the MB theory, of the density is shown in Fig. 10(a) and 10(d) for ω < ω c and ω > ω c , respectively. In both cases, the modulation of the Thomas-Fermi radius follows the external perturbation. Once the modulation is turned off, the radius oscillates at its natural quadrupolar frequency. The 1 st -order spatial coherence is shown in Fig. 10(b,e) for the same parameters. The patterns that emerge and persist in g (1) (z, z ) demonstrate that spatial correlations between particles at distinct and distant locations in the granular state are present [ Fig. 10(e)]. The length-scale of the patterns in g (1) (z, z ) is similar to what is seen in Fig. 8 for g (2) (z, z ). We infer that the process of granulation in a BEC is accompanied by the emergence of non-local correlations in the many-body state. Fig. 10(f) shows the emergence of two macroscopic eigenvalues of the reduced one-body density matrix for ω > ω c . While these so-called natural occupations are unaffected by modulation for ω < ω c , as seen in Fig. 10(c), ω > ω c results in the second   The onset of granulation observed experimentally is shown in Fig. 11. The condensate was Granular states feature random patterns and lack periodicity in their distributions, distinguishing them from Faraday and shock waves [51]. The multi-characteristic nature of quantum grains is supported by our observation of additional anomalous features in real and momentum space. Indeed, we find signatures of different co-existing phases of perturbed quantum systems such as quantum turbulence and localization in granulated states.
We verified that the density in momentum space (as calculated from the MB theory) of the granulated state shows clear signs of a k −2 power-law scaling (see Supplemental Materials [38] and Fig. S2 therein) which indicates a connection to turbulent BECs [52][53][54].

IV. CONCLUSIONS
We have explored the response of a BEC to modulated interactions. In the regime where the drive frequency ω ω r , the drive couples to parametric and resonant modes that result in 1D spatial pattern formation. For ω near resonant with 2ω r or ω r , very little modulation time and amplitude are required to produce a significant response. Near these resonances the condensate undergoes breathing oscillations that persist for a long time, resulting in the formation of Faraday and resonant mode patterns for t h > 0. A pattern is also observed offresonance, but only with increased modulation amplitude and modulation time. Due to the long modulation time, the resulting pattern can be seen at t h = 0, and is a direct consequence of the applied modulation. The dispersion relation of both Faraday and resonant modes is well-represented by a mean-field theory that accounts for the 3D nature of the elongated condensate.
For lower drive frequencies, the modulated interactions only weakly couple to the condensate. Significant response is achieved only by increasing the modulation amplitude and time, and then, only above a critical modulation frequency ω c . Fluctuating and irregular spatial patterns, that we define as grains, may then emerge and persist for long periods of time. A theoretical description of granulation requires approaches that go beyond mean-field theory, indicating that quantum granulation is characterized by non-local many-body correlations and quantum fluctuations.

EXPERIMENTAL DETAILS
A pair of coils in Helmholtz configuration is used to produce a homogeneous magnetic field, B, which allows us to vary the interatomic interactions. For a given value of B the corresponding scattering length is determined from: where a bg = −24.5 a 0 , B ∞ = 736.8 G, ∆ = 192.3 G [30] and a 0 the Bohr radius. An oscillation of the bias field, B(t) =B + ∆Bsin(ωt), whereB is the mean and ∆B is the modulation amplitude. This produces an asymmetric a(t) since a is a non-linear function.
Thus, the mean isā, the maximum is a + , and minimum is a − .
To solve the time-dependent Schrödinger equation for many interacting particles, we apply the Multiconfigurational Time-Dependent Hartree theory for Bosons (MCTDHB) [28,29] and use the MCTDH-X numerical solver [41][42][43] for 1D and 3D simulations. The MCTDHB theory assumes a general ansatz Ψ = Ψ(R, t) for the N −particle problem and expands it on a many-body basis Ψ( At ω/2π = 30 Hz we have seen resonant behavior: the energy increases up to ≈ 10 times after 500 ms and the density is found to occupy all available space. Convergence checks are beyond the computational capacities and the point at ω/2π = 30 Hz has not been included in the plots. We attribute the resonant behavior at ω/2π = 30 Hz to its proximity to 2ω Q , where ω Q = √ 3ω z is the 1D quadrupolar frequency. We have also performed calculations for ω Q = 13.9 Hz and 2ω Q = 27.9 Hz, and have observed similar behavior.

2.
A three-dimensional system with ω z = 1, ω y = ω x = 32, N = 1000, M = 1 and g (3D) = 4πN exp a/l z = 222, using also a delta-type interaction pseudopotential. The computational grid was 512 × 64 × 64 wide. All other parameters are set as in paragraph 1. The modulation frequency was set to ω = 8.75ω z , that corresponds to the experimental value ω/2π = 70Hz. The amplitude of modulation of the interaction is, as before, always positive (results plotted in Fig. 6). commonly employed and gives the p-particle probability densities. The eigenbasis of the RDMs gives information on the p th -order coherence of the system. In particular, if there is more than one macroscopic eigenvalues of the first (second) order RDM then the system is fragmented and first (second) order coherence is lost.
Specifically, the p th -order reduced density matrix (RDM) is defined as [45]: where n   m is macroscopically occupied, or, n m /N ∼ 1 for some m, while n j /N ∼ 0 for j = m. If more than one natural orbital is macroscopically occupied then the system is called fragmented [47]. The diagonal we simply call density. The eigenfunctions φ (2) k (z 1 , z 2 ) of the 2 nd order RDM are known as natural geminals (NG). Their occupations satisfy j=1 n (2) j = N (N − 1) and are plotted in Fig. 9(c) (normalized to 1).
The pth order correlation function is: The skew diagonal (antidiagonal) gives the degree of correlation of the density at a point z with its antipodal at point [48] z ≡ −z (see Fig. 9). Similarly, the normalized pth order correlation function in momentum space can be defined, via the Fourier transformρ (p) (k 1 , . . . , k p |k 1 , . . . , k p ; t) of ρ (p) (z 1 , . . . , z p |z 1 , . . . , z p ; t). Note that |g (1) |, the spatial correlation function, is bounded like 0 ≤ |g (1) | ≤ 1 for any two points (z, z ). For Bose condensed and hence non-fragmented states, |g (1) | takes its maximal value everywhere in space and the state is first-order coherent. Moreover, if |g (2) | < 1 we term the state anticorrelated while for |g (2) | > 1 we term it correlated.
The 2 nd order correlation function of Fig. 8(a) and Fig. 8(c) for some observed distributions n(z) is given by: C (2) (z, z ) = n(z)n(z ) n(z) n(z ) .
We emphasize that in the expression for C (2) , the notation · corresponds to an average across experimental realizations.
The single-shot simulations plotted in Fig. 7(b) and Fig. S1 have been obtained with the method to obtain random deviates of the N -particle probability density |Ψ| 2 that is prescribed in Refs. [39,40]. In brief, the procedure relies on sampling the many-body probability density |Ψ| 2 as follows: one calculates the density ρ 0 (z), from the obtained solution |Ψ (0) ≡ Ψ of the MCTDHB equations. A random position z 1 is drawn from ρ 0 (z).
In continuation, one particle is annihilated at z 1 , the reduced density ρ 1 of the reduced system |Ψ (1) is calculated and a new random position z 2 is drawn. The procedure continues for N − 1 steps and the resulting distribution of positions (z 1 , z 2 , . . . , z N ) is a simulation of an experimental single-shot image.
The contrast parameter D quantifies the deviation of some spatial distribution n(z) = n(z; t 0 ) of a single shot at a given time t 0 from the parabolic (Thomas-Fermi-like) best fit n bf (z) = n bf (z; t 0 ) at the same time and is defined as: where i runs over all n gp pixels/grid points. The cutoff requirement C cutoff = 0.20 n bf (0) is set so that small (zero-excitation) fluctuations are wiped out and only values with large deviations are considered (see Fig. S1). Therefore, the resulting contrast parameter reflects only the large deviations of a given density from its parabolic best fit. To determine the best fits we used the gnuplot software to fit the polynomial p(z) = −a(z − b) 2 + c, where a, b, c ∈ R, to the obtained experimental or numerical distributions n(z) along z. The twodimensional experimental column densities have been integrated along y. The experimental data were also interpolated to a number of points along z so as to equal the grid used for the numerical simulations. An example of a processed image is shown in Fig. S1. t ≤ 250ms and after for t > 250ms). In the granulated case the momentum distribution scales like k −2 (straight line to guide the eye) for almost two decades, behavior that is characteristic of quantum turbulence. Contrary to the regular gas, this scaling remains even 250ms after the modulation.
FIG. S3. Experimental column densities exponentially fitted. Close to the threshold frequency ω/2π = 40Hz, where the system transitions from regular to granulated states, anomalous spatial distributions are seen (here, two experimental shots for the same initial conditions). These might bear resemblance to localized states, that have been shown to exist in BECs in optical lattices [49]. We fit our observed density distributions to C + A exp − |x−x0| α d and obtain α = 1.25 and 1.75 for the two shots. The transition from a regular to a localized states happens as α → 1. For comparison, we plot the parabolic Thomas-Fermi (TF) fit (blue).