Supersolid formation in a dipolar condensate by roton instability

We characterize the role of the roton instability in the formation of a supersolid state of an elongated dipolar condensate, following a quench of the contact interactions across the superfluid-supersolid transition, as observed in recent experiments. We perform dynamical simulations by means of the extended Gross-Pitaevskii equation including quantum corrections, for different final values of the $s-$wave scattering length. The corresponding excitation spectrum is computed using an effective one-dimensional description, revealing that the calculated growth rates of the unstable roton mode accurately reproduce the observed behavior. Our results provide valuable insights regarding the formation time of the supersolid and its scaling behavior with respect to the $s$-wave scattering length.

A remarkable feature of this transition is that whereas one can relax an initial SS state onto a SF state by crossing the transition almost adiabatically, in the opposite direction the supersolid requires a finite time to form [13,23,28,41].This corresponds to the time required for intrinsic fluctuations -associated to the roton insta-bility -to grow up and break the translational invariance of the system.
Here, we provide a theoretical characterization of the role of this instability in the dynamics of formation of a supersolid in a quasi-1D dipolar condensate, by considering a typical configuration of the experiment in Ref. [28].In particular, we focus on the scenario in which the system crosses the SF-SS phase transition by undergoing a quench of the contact interactions which, for quasi-1D systems, can still provide a smooth connection between an unmodulated BEC and a supersolid array [41] (unlike what occurs in 2D [29]).We perform a systematic analysis by means of numerical simulations of the extended Gross-Pitaevskii equation, considering various amplitudes of the interaction quench.The corresponding roton spectrum is computed by means of the 1D theory proposed by Blakie et al. [33], which provides a simplified yet accurate characterization of the simulation outcomes.These results are then used to discuss the general features that determine the formation time of the supersolid and its scaling behavior with respect to the s-wave scattering length.
The paper is organized as follows.In Section II, we provide an overview of the system under consideration and briefly summarize the relevant formulas defining the extended Gross-Pitaevskii theory for dipolar condensates.In Sec.III, we discuss the general protocol employed to induce the interaction quench across the SF-SS phase transition, presenting the general phenomenology observed through numerical simulations.In Sec.IV we briefly review the 1D effective theory of Ref. [33] used to study the excitation spectrum of an elongated dipolar condensate.The latter is computed for our specific case in Sec.V, where it is used to characterize systematically the role of the roton instability in the formation of the supersolid.There, we also explore the implications of this instability on the scaling behavior of the supersolid formation time with respect to the s-wave scattering length.Finally, we present a summary of our findings and concluding remarks in Sec.VI.

II. SYSTEM
In the following analysis, we consider the quasi 1D configuration investigated in the experiment of Ref. [28] and theoretically analyzed in Ref. [41].It consists of a dipolar condensate, at zero temperature, composed by N = 3 × 10 4 magnetic atoms of 162 Dy -with tunable swave scattering length a s and dipolar length a dd = 130a 0 (a 0 being the Bohr radius) -trapped by a harmonic potential with frequencies (ω x , ω y , ω z ) = 2π × (15, 101, 94) Hz.While the choice of these parameters is motivated by their experimental feasibility in a specific case, it is worth noting that the following analysis is conceptually general, which enables its extension to other scenarios.
This system can be described in terms of an extended Gross-Pitaevskii (GP) theory including dipolar interactions [57] and the Lee-Huang-Yang (LHY) correction accounting for quantum fluctuations, within the local density approximation [58][59][60].The energy functional can be written as where E GP = E k + E ho + E int is the standard GP energy functional including the kinetic, potential, and contact interaction terms, V (r) = (m/2) α=x,y,z ω 2 α r 2 α is the harmonic trapping potential, n(r) = |ψ(r)| 2 represents the condensate density (normalized to the total number of atoms N ), g = 4π 2 a s /m is the contact interaction strength, V dd (r) = (1 − 3 cos 2 θ)/(4πr 3 ) the inter-particle dipole-dipole potential, C dd ≡ µ 0 µ 2 its strength, µ the modulus of the dipole moment µ, r the distance between the dipoles, and θ the angle between the vector r and the dipole axis, cos θ = µ • r/(µr).As in Refs.[28,41], we consider the magnetic dipoles to be aligned along the z direction by a magnetic field B.
As discussed in Refs.[28,41], the equilibrium configuration of the system corresponds to either a conventional superfluid (SF) state or a supersolid (SS) state.The transition from one phase to the other can be induced by tuning the s-wave scattering length a s .For the present values of the number of atoms and trapping frequencies, the critical point is located at a c s ≃ 94.4a 0 and the transition has a continuous character [28,41].

III. PROTOCOL
In order to study the formation of a supersolid, we adopt the following approach.The system is initially prepared in an equilibrium configuration within the SF phase, at a (in) s = a c s + 1.5 a 0 .For the sake of simplicity, this value is kept fixed throughout this work.Then, the condensate is quenched into the SS phase by a sudden change of the contact scattering length, to a final value a In the present study, δa s varied in the range [1.0, 6.5]a 0 .
The dynamics of the system following the quench is obtained by solving the GP equation [48] where the energy functional E[ψ, ψ * ] is the one in Eq. ( 1) [61].To uncover the physics involved in the formation process, we simplify the analysis by excluding dissipation mechanisms and particle losses, which are not essential for the present discussion.
A useful quantity for characterizing the formation process of the supersolid is represented by the inverse participation ratio (IPR), that measures the degree of localization (high IPR) or spread (low IPR) of a certain quantum state.For a continuous system, it can be defined as IP R = |ψ| 4 dV .In the present case, it turns out to be proportional to the contact interaction energy E int [see Eq. ( 1)], In order to facilitate the ongoing discussion, it is convenient to normalize the above expression as which is also equivalent to normalizing the inverse participation ratio to its initial value at t = 0.The behavior of this quantity as a function of the time t elapsed after the quench is shown in Fig. 1a, for different values of δa s .
It is noteworthy that all the cases presented exhibit the same qualitative behavior, which is exemplified in Fig. 1b for the case δa s = 2.0 a 0 .Just after the quench, the system retains its initial character of a SF condensate for a certain period of time, during which it undergoes mostly breathing-like oscillations.These oscillations exhibit a consistent pattern among the different values of δa s , as indicated by the blue arrows in Fig. 1a.After some time, the supersolid structure emerges, following a sudden increase of Ēint .As we shall see later on, this is a typical signature of an underlying formation process driven by the roton instability, namely the exponential growth of the corresponding mode, which becomes unstable.In Sec.IV, we will provide a detailed quantitative analysis and characterization of this behavior.
It should be emphasized that the SS state generated by this mechanism can display significant deviations from the supersolid ground state at the same value of a s , due to the highly out-of-equilibrium nature of both the quench and the formation process (see also the discussion in Refs.[6, 14,15]).In addition, since dissipation processes are not included in the present analysis, the system remains in such a highly excited state even long after the formation of the supersolid.However, it is worth noticing that for sufficiently low values of δa s (up to δa s = 2 a 0 , for the cases shown in Fig. 1a), the value Ē * int at which the instability saturates almost coincides with the equilibrium value ĒSS int for the ground state, indicated by horizontal lines in the figure.This is no longer the case at higher values of δa s , for which IPR saturation takes place at lower values than IPR for the SS ground state, Ē * int < ĒSS int .This effect becomes more prominent by increasing δa s , owing to the increase of nonlinear effects by moving deeper into the SS phase (having started in all cases in the same initial configuration).
In Fig. 1b, the time frame during which the system exhibits the instability is illustrated as a (blue) shaded area in between two consecutive oscillation maxima represented by blue and red dots.The density distributions of the corresponding states are shown in the insets.The first one corresponds to a state that can be still associated with the SF phase, despite displaying a weak density modulation.The other clearly corresponds to a wellformed SS state.Considering the above scenario, the time it takes for the supersolid to form from the instant of the quench, denoted as τ SS and named formation time in what follows, can be therefore conveniently defined by referring to the position of the red dot(s), as illustrated in the same figure.Figure 1a shows that this formation time gets reduced by increasing δa s .

IV. EFFECTIVE 1D DESCRIPTION
To gain a quantitative understanding of the superfluid formation process, we will employ the effective 1D model for an elongated dipolar condensate described in Ref. [33], simplifying the subsequent analysis.Owing to the strong transverse confinement of the present configuration, ω x /ω y,z ≈ 0.15, and the small asymmetry between the transverse trapping frequencies, ω y /ω z ≈ 1 (see Sec. II), this approximation is expected to be reasonably accurate.
In summary, this approach consists in factorizing the condensate wave function as ψ(r) = ϕ(x)χ(y, z), where the transverse wave function χ(y, z) is conveniently taken as a Gaussian, with l = l y l z , l y (l z ) being the 1/e half width of the transverse density along the y-axis (z-axis), and η = l y /l z the transverse anisotropy of the density distribution.
From a Gaussian fit of the transverse profile of the initial SF density distribution we obtain l y ≃ 0.82µm and l z ≃ 2.19µm.
Integrating out the transverse directions, one obtains an effective 1D GP model where the contact interaction strength and the LHY coefficient get renormalized as g 1D = g/(2πl y l z ), γ (1D) LHY = γ ⊥ γ LHY , with γ ⊥ = 2/5π 3/2 l 3 .The dipole-dipole effective interaction potential can be conveniently approximated in momentum space by with q ≡ η 1/4 lk/ √ 2 and k representing the momentum component relative to the x−axis.
The description can be further simplified by considering small deviations from uniformity, with linear density n 0 .The longitudinal wave function can be expressed as ϕ(x, t) = √ n 0 + δϕ(x, t), with the latter term representing a small perturbation over the initial state ϕ 0 = √ n 0 .The collective excitations of the condensate are then described by the associated Bogoliubov-de Gennes equations, see, e.g., [6,10,33,57].By looking for solutions of the form δϕ(x, t) ∝ u(x)e −iωt + v * (x)e iωt and considering that for a uniform system the excitations are plane waves of momentum k, namely u k (x) = U k e −ikx and v k (x) = V k e −ikx , one finally obtains the following expression for the excitation energies [33] where Ṽ (k) denotes the Fourier transform of the interaction potential, see Eq. ( 6), The above spectrum is characterized by a roton excitation (that is, a local minimum in the excitation dispersion relation) that softens to zero energy and becomes dynamically unstable when the s−wave scattering length is tuned below a certain critical value a c s [6,33].

V. STABILITY ANALYSIS
The quasi-1D effective formulation presented above allows a straightforward stability analysis of our system after the quench.In order to do so, we need to establish a criterion to account for the non-uniformity of the system.In particular, we use the fact that for a continuous transition, the critical point for the SF-SS transition is expected to coincide with the roton instability, see Refs.[36,44].Specifically, by defining n(x) ≡ n(r)dydz, we set n 0 = cn(0), where c is chosen to reproduce the critical value of the s−wave scattering length a c s = 94.4a0 obtained from numerical simulations (see Sec. III).We find c ≃ 0.5.
In Fig. 2 we show the imaginary part of ω k , as a function of the s−wave scattering length a s and of the excitation momenta k.We recall that the presence of imaginary frequencies in the spectrum is associated with exponentially growing modes that make the system modulationally unstable, if they are initially populated.In the figure, the vertical (red) dashed line corresponds to the critical value of the scattering length, a c s = 94.4a0 .Above a c s the spectrum is purely real, corresponding to a stable SF state.Below a c s the frequency of some modes becomes imaginary, as indicated by the colored area.This unstable region broadens by decreasing the value of a s .The dot-dashed line corresponds to the position of the most unstable mode, for a s < a c s , which connects to that of the roton mode, for a s > a c s .We indicate its position as k R (a s ).It is also worth noting that the LHY correction in Eq. ( 7) has a negligible contribution to the spectrum in Fig. 2 for the parameter values at hand.

A. Supersolid formation
Now that we have determined the properties of the excitation spectrum, we can revisit the dynamical behavior of the condensate after the quench, discussed in Sec.III, and examine the formation of the supersolid in terms of the emergence of exponentially growing modes.In particular, we use the fact that in the present case the exponentially growing modes take the form u k (x, t) = U k e −ikx e t/τ k , with τ −1 k ≡ Im[ω k ] indicating the growth rate.We also make the assumption that instability is dominated by the most unstable mode, k → R.Then, the normalized inverse participation ratio in Eqs. ( 3), ( 4) can be approximated by where we have defined . The above expression, which depends on the two independent parameters A and τ , fixes the scaling behavior of the instability.In this respect, it is important to mention that while τ is an intrinsic property of the system, determined by Eq. ( 7), A depends through U k on the initial population of the unstable modes, namely on the preparation of the initial state.
To illustrate the instability behavior, in Fig. 3 we consider the case δa s = 5 a 0 , as an example.In the top panel (a), the evolution of Ēint (t) (see also Fig. 1) is compared with the fitting function 1 + f A,τ (t) in Eq. ( 9), using A and τ as fitting parameters.The fit is restricted to the shaded area, where the instability becomes manifest.Obviously, the simplified model in Eq. ( 9) cannot describe the initial fluctuations in the SF phase, that stem from the non uniform nature of the actual system, nor the collective oscillations of the SS state, that clearly fall outside the regime of linear excitation of the initial state.However, it provides a clear explanation of the two different regimes observed after the quench, before the formation of the supersolid.The initial regime, in which the condensate retains its SF character, corresponds to the one in which the population of the unstable modes remains negligible, namely f A,τ (t) ≪ 1.Then, an evident exponential behavior emerges, leading to the formation of the supersolid, at t = τ SS .At this point, the system has already left the linear excitation regime, and Eq. ( 9) no longer applies.From τ SS on, the dynamics of the system once again become predominantly driven by nonlinear effects, with quantum fluctuations being instrumental for stabilizing the system.The growth of unstable modes is clearly visible in Fig. 3b, where we show a density plot representing the time evolution of the momentum distribution along the x−direction (see also Ref. [6]).The horizontal (green) dashed lines correspond to the nominally most unstable mode in Fig. 2, namely the roton excitations at ±k R .Notably, it is in good agreement with the result of the numerical simulation (see also Ref. [6,12]).
The same analysis can be repeated for the other values of δa s considered in Sec.III.In Fig. 4 we compare the inverse growth rate extracted from the fit, τ f it , with the value τ R ≡ Im[ω R ] −1 corresponding to the most unstable mode in Fig. 2. We also show, as a light shaded area, how τ R changes by varying the linear density n 0 .As one may expect, the instability mechanism is more effective for larger densities.This observation is consistent with the results of the GP simulations, which show that the formation of the SS pattern first occurs at the center of the condensate, where the density is highest, and subsequently spreads outward toward the lower density region in the tail.Overall, the agreement between the model predictions and the numerical data is remarkably good, indicating that the model captures the essential physics of the system despite the inherent simplifications leading to Eq. ( 7), such as the assumption of uniformity.
In Fig. 4 we also show the total formation time τ SS , as red squares.Interestingly, the dependence of τ SS on δa s follows a behavior similar to that of τ R , namely τ SS ≃ ατ R (δa s ) (notice the log scale of the vertical axis), ) δa s (units of a 0 ) Plot of the different characteristic times τ entering the formation of the supersolid.τ f it is the inverse growth rate of the instability extracted from a fit of Ēint(t) as in Fig. 3, to be compared with τR obtained from the stability analysis in Fig. 2. The light shaded area represents the variation of the latter upon a variation of +10% (downwards) and 5% (upwards) of the linear density n0.τSS is the formation time of the SS state defined in Fig. 1.The black dotted line corresponds to ατR(δas), with α ≃ 6.5 (see text).
with α being a proportionality factor.In the present case we find α ≃ 6.5, see the black dotted line in the figure .Actually, the figure shows that there must be also a slight correction to the scaling, namely α = α δas .This can be explained as follows.From the simplified model in Eq. ( 9) it is straightforward to express the supersolid formation time as This result confirms that the scaling of τ SS as a function of δa s is indeed mostly determined by the behavior of τ R (δa s ), with Ē * int (δa s ) and U R (δa s ) providing a logarithmic correction.Regarding the initial population of the unstable modes (U R ), it is worth noting that -given a certain realization of the initial state (and of its momentum distribution) -it may slightly depend on a s due to the varying position of the roton momentum k R (a s ) (see Fig. 3).Although we do not include temperature in our considerations, we signal that if there is an increase in the initial population at the roton mode because of thermal excitation then the formation time of the supersolid would be shortened with respect to the T = 0 prediction, as has been observed in experiments [13,28,41].It is also important to remark that the logarithmic terms, though they only provide a correction to the scaling, are essential for fixing the value of the proportionality constant α δas .In particular, we find that the major contribution comes from ln √ n 0 /U R in the range of values of δa s considered here.

VI. CONCLUSIONS
We have characterized the role of the roton instability in the formation of a supersolid state of an elongated dipolar condensate across a continuous superfluidsupersolid transition, as experimentally reported by Biagioni et al. [28] and theoretically investigated in Ref. [41].By means of numerical simulations based on the extended Gross-Pitaevskii equation, we have systematically analyzed the dynamic behavior of the system after a quench of the contact interactions considering various amplitudes δa s of the interaction quench.We have also computed the corresponding excitation spectrum by using the effective one-dimensional approach of Blakie et al. [33], finding that the calculated growth rates of the unstable roton modes provide an accurate characterization of the behavior observed in the simulations.
Our main finding is that the scaling behavior of the formation time τ SS of a supersolid state as a function of δa s is mostly determined by the scaling of the inverse growth rate of the roton mode, namely τ SS (δa s ) ≃ α δas τ R (δa s ), with α δas being a proportionality factor fixed by the ini-tial population of the (most) unstable mode.Remarkably, this implies that the supersolid formation time τ SS is determined by the properties of the initial state rather then by the energy difference between the initial superfluid state and the final supersolid state, as it was speculated in Refs.[13,41].The present analysis -that admits a straightforward extension also to other cases besides a quench (e.g., by means of a time rescaling as discussed in Ref. [6]) -offers a deeper insight into the formation dynamics of the supersolid structures in dipolar condensates, thereby providing a valuable framework for future experiments and theoretical studies.
FIG. 1.(a) Behavior of Ēint(t) during the evolution after the quench, for different values of δas ranging from 1 a0 to 3 a0.The red dots correspond to the values Ē * int at which the exponential growth saturates; the blue arrows the oscillations maxima before the formation of the SS pattern.The horizontal lines indicate the equilibrium value of ĒSS int of the SS ground state at the same value of δas (same line type).(b) The case with δas = 2 a0 in detail.Density plots (color scale weighted by each density distribution) show the SF configuration and the SS configuration corresponding to the blue and red dots, respectively.

k (µm − 1 )FIG. 2 .
FIG. 2. Imaginary part of ω k , as a function of the s−wave scattering length as and of the excitation momenta k.The vertical (red) dashed line indicates the critical value of the scattering length, a c s = 94.4a0.The black area corresponds to Im[ω k (as)] = 0.The (yellow) dot-dashed line represents the position of the roton mode, which disappears at a c s , where it is replaced by the most unstable mode (black dot-dashed line).The whole line is indicated as kR(as).

2 FIG. 3 .
FIG. 3. Time evolution of the system after a quench with δas = 5 a0.(a) Plot of Ēint(t) (thick line) along with the fitting function in Eq. (9) (thin red line).The grey area indicates the range in which the fit has been applied (see text).(b) Density plot of the momentum distribution along the x−direction as a function of time.The horizontal (green) dashed line corresponds to the nominal most unstable modes at ±kR, see Fig. 2. Side peaks at integer multiples of the former are also visible (notice that the plot is in logarithmic scale).