Supersolid formation in a dipolar condensate by roton instability

PHYSICAL REVIEW A 108, 033316 (2023) Supersolid formation in a dipolar condensate by roton instability Aitor Alaña ,1,2Iñigo L. Egusquiza ,1,2and Michele Modugno 1,2,3 1Department of Physics, University of the Basque Country UPV/EHU, 48080 Bilbao, Spain 2EHU Quantum Center, University of the Basque Country UPV/EHU, 48940 Leioa, Biscay, Spain 3IKERBASQUE, Basque Foundation for Science, 48013 Bilbao, Spain (Received 31 May 2023; accepted 8 September 2023; published 26 September 2023) We characterize the role of 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, which provides a reasonably accurate prediction for the growth rate of the most unstable mode observed in the simulations. To analyze the behavior of the system, we employ the inverse participation ratio, which conveniently characterizes the different degree of localization in the superfluid and supersolid phases. By means of a suitable effective ansatz for the density, we derive a simple yet effective expression for the formation time of the supersolid. This expression provides valuable insights regarding its scaling behavior with respect to the s-wave scattering length. DOI: 10.1103/PhysRevA.108.033316 I. INTRODUCTION Ultracold dipolar gases, characterized by significant magnetic or electric dipole moments, have emerged as a fascinating area of research in the field of cold atoms and quantum gases. Unlike traditional Bose-Einstein condensates (BECs), with only short-range contact interactions, dipolar condensates exhibit also long-range, anisotropic interactions (see, e.g., Refs. [1–3]). Notably, their excitation spectrum is characterized by a rotonic mode [4–9], which may become unstable [10,11] and lead to exotic and complex behavior, including the formation of droplets [12] and supersolids [13–16], among others [17–19]. In particular, in recent years an intense research activity has been directed towards unraveling the properties of the supersolid phase of dipolar gases, both from the experimental [20–33] and theoretical point of view [34–47]. Proposed in the past century [48,49], supersolids are an exciting phase of matter combining a superfluid nature [50,51] with the translational symmetry-breaking characteristic of solid structures [52–55]. Notably, the supersolid (SS) phase of matter can be achieved through both a classical transition from a gas to a supersolid [30] and a quantum transition from an unmodulated superfluid (SF) to a supersolid. Experimentally and theoretically both discontinuous [13,14,22,56–58] and continuous [27,59] features have been observed, reminiscent of the first-order and second-order transitions predicted in the thermodynamic limit in two dimensions (2D) and one dimension, Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. respectively. Interestingly, the effective dimensionality of the system can be controlled by tuning the transverse confinement and the atom number, as experimentally observed by Biagioni et al. [31] and theoretically discussed in Ref. [44]. 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 [14,26,31,44]. This corresponds to the time required for intrinsic fluctuations— associated to the roton instability—to grow up and break the translational invariance of the system. This distinctive rotonic spectrum of a dipolar superfluid was recently observed and characterized in the groundbreaking experiment by Chomaz et al. in Ref. [7]. Here, we provide a theoretical characterization of the role of this instability in the dynamics of formation of a supersolid in a quasi-one-dimensional (quasi-1D) dipolar condensate, by considering a typical configuration of the experiment in Ref. [31]. 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 [44] (unlike what occurs in 2D [32]). 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 one-dimensional theory proposed by Blakie et al. [36], which provides a simplified yet accurate characterization of the simulation outcomes. Specifically, we investigate the evolution of the system by analyzing the inverse participation ratio (IPR), a valuable measure for characterizing the varying degrees of localization in the superfluid and supersolid phases. By employing an effective ansatz 2469-9926/2023/108(3)/033316(7) 033316-1 Published by the American Physical Society ALAÑA, EGUSQUIZA, AND MODUGNO PHYSICAL REVIEW A 108, 033316 (2023) for the density, we derive a simple expression containing a minimal set of fitting parameters, that provides a reliable description of the IPR behavior in the regime dominated by the growth of the roton mode. 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 Sec. II, we provide an overview of the system under consideration and briefly summarize the relevant formulas defining the extended GrossPitaevskii 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. [36] 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. [31] and theoretically analyzed in Ref. [44]. It consists of a dipolar condensate, at zero temperature, composed by N=3×104 magnetic atoms of 162Dy—with tunable s-wave scattering length asand dipolar length add =130a0(a0being 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 [60] and the Lee-Huang-Yang (LHY) correction accounting for quantum fluctuations, within the local density approximation [61–63]. The energy functional can be written as E=EGP +Edd +ELHY with EGP =¯h2 2m|∇ψ(r)|2+V(r)n(r)+g 2n2(r)dr, Edd =Cdd 2 n(r)Vdd(r−r)n(r)drdr,(1) ELHY =2 5γLHY n5/2(r)dr, where EGP =Ek+Eho +Eint is the standard GP energy functional including the kinetic, potential, and contact interaction terms, V(r)=(m/2)α=x,y,zω2 αr2 αis the harmonic trapping potential, n(r)=|ψ(r)|2represents the condensate density (normalized to the total number of atoms N), g=4π¯h2as/mis the contact interaction strength, Vdd (r)=(1 −3 cos2θ)/(4πr3) the interparticle dipole-dipole potential, Cdd ≡μ0μ2its strength, μthe modulus of the dipole moment μ,rthe distance between the dipoles, and θ the angle between the vector rand the dipole axis, cos θ=μ· r/(μr). As in Refs. [31,44], we consider the magnetic dipoles to be aligned along the zdirection by a magnetic field B.The LHY coefficient is γLHY =128√π¯h2a5/2 s/(3m)(1 +32 dd/2), with dd =μ0μ2N/(3g). As discussed in Refs. [31,44], the equilibrium configuration of the system corresponds to either a conventional SF state or a SS state. The transition from one phase to the other can be induced by tuning the s-wave scattering length as. For the present values of the number of atoms and trapping frequencies, the critical point is located at ac s≃94.4a0and the transition has a continuous character [31,44]. 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= ac s+1.5a0. For the sake of simplicity, this value is kept fixed throughout this work. Then, the condensate is quenched in the SS phase by a sudden change of the contact scattering length, to a final value a(fi) s=ac s−δas. In the present study, δasis varied in the range [1.0,6.5]a0. The dynamics of the system following the quench is obtained by solving the GP equation [51] i¯h∂tψ=δE[ψ,ψ∗]/δψ∗,(2) where the energy functional E[ψ,ψ∗] is the one in Eq. (1) [64]. 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 IPR, which 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 IPR =|ψ|4dr. In the present case, it turns out to be proportional to the contact interaction energy Eint [see Eq. (1)], IPR ∝Eint(t)=g 2n2(r,t)dr.(3) In order to facilitate the ongoing discussion, it is convenient to normalize the above expression as ¯ Eint(t)≡Eint(t)/Eint(0),(4) 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 telapsed after the quench is shown in Fig. 1(a), for different values of δas. It is noteworthy that all the cases presented exhibit the same qualitative behavior, which is exemplified in Fig. 1(b) for the case δas=2.0a0. 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 δas, as indicated by the blue arrows in Fig. 1(a). 033316-2 SUPERSOLID FORMATION IN A DIPOLAR CONDENSATE … PHYSICAL REVIEW A 108, 033316 (2023) (a) (b) FIG. 1. (a) Behavior of ¯ Eint(t) during the evolution after the quench, for different values of δasranging from 1 a0to 3 a0.The red dots correspond to the values ¯ E∗ 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 ¯ ESS int of the SS ground state at the same value of δas(same line type). (b) The case with δas=2a0in 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. After some time, the supersolid structure emerges, following a sudden increase of ¯ Eint. 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 as, due to the highly out-of-equilibrium nature of both the quench and the formation process (see also the discussion in Refs. [7,16,17]). 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 δas[up to δas=2a0, for the cases shown in Fig. 1(a)], the value ¯ E∗ int at which the instability saturates almost coincides with the equilibrium value ¯ ESS int for the ground state, indicated by horizontal lines in the figure. This is no longer the case at higher values of δas,forwhich IPR saturation takes place at lower values than IPR for the SS ground state, ¯ E∗ int <¯ ESS int . This effect becomes more prominent by increasing δas, 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. 1(b), the time frame during which the system exhibits instability is illustrated as a (blue) shaded area 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 well-formed 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 1(a) shows that this formation time gets reduced by increasing δas. 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. [36], 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 of factorizing the condensate wave function as ψ(r)=ϕ(x)χ(y,z), where the transverse wave function χ(y,z) is conveniently taken as a Gaussian, χ(y,z)=(√πl)−1e−(ηy2+z2)/2ηl2,(5) with l=lylz,ly(lz) being the 1/ehalf width of the transverse density along the yaxis (zaxis), and η=ly/lzthe transverse anisotropy of the density distribution. From a Gaussian fit of the transverse profile of the initial SF density distribution we obtain ly≃0.82 μm and lz≃2.19 μm[65]. Integrating out the transverse directions, one obtains an effective 1D GP model where the contact interaction strength and the LHY coefficient get renormalized as g1D =g/(2πlylz), γ(1D) LHY =γ⊥γLHY, with γ⊥=2/5π3/2l3. The dipole-dipole effective interaction potential can be conveniently approximated in momentum space by ˜ V(1D) dd (q)=1 1+η μ0μ2 Dy 2πlylzq2eq2Ei(−q2)+2−η 3,(6) with q≡η1/4lk/√2 and krepresenting the momentum component relative to the xaxis. The description can be further simplified by considering small deviations from uniformity, with linear density n0. The longitudinal wave function can be expressed as ϕ(x,t)=√n0+δϕ(x,t), with the latter term representing a small perturbation over the initial state ϕ0=√n0. The collective excitations of the condensate are then described by the associated Bogoliubov–de Gennes equations, see, e.g., Refs. [7,11,36,60]. By looking for solutions of the form δϕ(x,t)∝u(x)e−iωt+v∗(x)eiωtand considering that for a uniform system the excitations are plane waves of momentum ¯hk, namely, uk(x)=Uke−ikx and vk(x)=Vke−ikx, one finally obtains the following expression for the excitation energies [36]: ¯hωk=±     ¯h2k2 2m¯h2k2 2m+2n0˜ V(k)+3γ(1D) LHY n3/2 0,(7) 033316-3 ALAÑA, EGUSQUIZA, AND MODUGNO PHYSICAL REVIEW A 108, 033316 (2023) FIG. 2. Imaginary part of ωk, as a function of the s-wave scattering length asand of the excitation momenta k. The vertical (red) dashed line indicates the critical value of the scattering length, ac 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 ac s, where it is replaced by the most unstable mode (black dot-dashed line). The whole line is indicated as kR(as). where ˜ V(k) denotes the Fourier transform of the interaction potential, see Eq. (6), ˜ V(k)=g1D +˜ V(1D) dd (η1/4lk/√2).(8) 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 ac s[7,36]. 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 nonuniformity 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. [39,47]. Specifically, by defining n(x)≡ n(r)dydz,wesetn0=cn(0), where cis chosen to reproduce the critical value of the s-wave scattering length ac s=94.4a0 obtained from numerical simulations (see Sec. III). We find c≃0.5. In Fig. 2we show the imaginary part of the positive branch of the frequency ωkas a function of the s-wave scattering length asand 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, ac s=94.4a0. Above ac s the spectrum is purely real, corresponding to a stable SF state. Below ac sthe frequency of some modes becomes imaginary, as indicated by the colored area. This unstable region broadens by decreasing the value of as. The dot-dashed line corresponds to the position of the most unstable mode, for as<ac s, which connects to that of the roton mode, for as>ac s. We indicate its position as kR(as). It is also worth noting that the LHY correction in Eq. (7) has a negligible contribution to the spectrum in Fig. 2for the parameter values at hand. 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, in order to capture the essential features of the exponential growth of the roton mode, we make use of the following ansatz for the density, n(x,t)≃n0+nR(t) cos(kRx), with nR(t)=ηRe2t/τR, where ηRrepresents the initial population of the roton mode, and τ−1 R≡Im[ωkR]. Such ansatz is phenomenological, providing a first approximation towards the nonhomogeneous density distribution. Then, the normalized inverse participation ratio in Eqs. (3) and (4) can be approximated as ¯ Eint(t)≃1 Ln2 0n(x,t)2dx =1+β1 ηR n0 e2t/τR+β2 η2 R n2 0 e4t/τR,(9) where Lis the total length of the system and the parameters β1and β2are numerical factors that depend on the geometry of the system. Typically, in a uniform system, their determination is straightforward. However, in the present case, we treat them as free parameters to accommodate the inherent nonuniformity of the system. Based on Eq. (9), we introduce the following expression for fitting the numerical data, f(t)≡1+2γAe2t/τ +A2e4t/τ ,(10) where we have conveniently reabsorbed the term ηR/n0into the definition of A, namely, A≡√β2ηR/n0,γ=β1/(2√β2). The above expression will be used for discussing the scaling behavior of the instability. It depends on the three independent parameters: A,γ, and τ. In this respect, it is important to mention that while τis an intrinsic property of the system, determined by Eq. (7), Adepends through ηR/n0on the initial population of the unstable modes, namely, on the preparation of the initial state, whereas γis just a numerical factor. To illustrate the instability behavior, in Fig. 3we consider the case δas=5a0, as an example. In Fig. 3(a), the evolution of ¯ Eint(t) (see also Fig. 1) is compared with the fitting function f(t)inEq.(10). The fit is restricted to the shaded area, where the instability becomes manifest [66]. Obviously, the simplified model in Eq. (10) cannot describe the initial fluctuations in the SF phase, that stem from the nonuniform nature of the actual system, nor the collective oscillations of the SS state, which 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(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 033316-4 SUPERSOLID FORMATION IN A DIPOLAR CONDENSATE … PHYSICAL REVIEW A 108, 033316 (2023) FIG. 3. Time evolution of the system after a quench with δas= 5a0. (a) Plot of ¯ Eint(t) (thick line) along with the fitting function in Eq. (10) (thin red line). The gray area indicates the range in which the fit has been applied (see text). (b) Density plot of the momentum distribution along the xdirection as a function of time. The horizontal (green) dashed line corresponds to the nominal most unstable modes at ±kR,seeFig.2. Side peaks at integer multiples of the former are also visible (notice that the plot is in logarithmic scale). excitation regime, and Eq. (10) 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. 3(b), where we show a density plot representing the time evolution of the momentum distribution along the xdirection (see also Ref. [7]). The horizontal (green) dashed lines correspond to the nominally most unstable mode in Fig. 2, namely, the roton excitations at ±kR. Notably, it is in good agreement with the result of the numerical simulation (see also Ref. [7,13]). The same analysis can be repeated for the other values of δas considered in Sec. III.InFig.4we compare the inverse growth rate extracted from the fit, τfit, with the value τR=Im[ωkR]−1 corresponding to the most unstable mode in Fig. 2.Wealso show, as a light shaded area, how τRchanges by varying the linear density n0. 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. 4we also show the total formation time τSS,as red squares. Interestingly, τSS displays a scaling behavior as a function of δasthat closely resembles that of τR, namely, τSS ≃ατR(δas) (notice the log scale of the vertical axis), with αbeing a proportionality factor. In the present case we find α≃6.5, see the black dotted line in the figure. Actually, the FIG. 4. Plot of the different characteristic times τentering the formation of the supersolid. τfit is the inverse growth rate of the instability extracted from a fit of ¯ Eint(t) as in Fig. 3, to be compared with τRobtained from the stability analysis in Fig. 2. The light shaded area represents the variation of the latter upon a variation of ±5% 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). figure shows that there must also be a slight correction to the scaling, namely, α=αδas. This can be explained as follows. From the simplified model in Eq. (10) it is straightforward to express the supersolid formation time as τSS =τR 2ln ¯ E∗ int −1+γ2−γ+ln n0 √βηR.(11) This result confirms that the scaling of τSS as a function of δasis indeed mostly determined by the behavior of τR(δas), with the other parameters providing logarithmic corrections. Regarding the initial population ηRof the unstable modes, it is worth noting that—given a certain realization of the initial state (and of its momentum distribution)—it may slightly depend on asdue to the varying position of the roton momentum kR(as) (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 [14,31,44]. 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 the term ln(n0/ηR)inthe range of values of δasconsidered 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 superfluid-supersolid transition, as experimentally reported by Biagioni et al. [31] and theoretically investigated in Ref. [44]. 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 δasof the interaction quench. We have 033316-5 ALAÑA, EGUSQUIZA, AND MODUGNO PHYSICAL REVIEW A 108, 033316 (2023) also computed the corresponding excitation spectrum by using the effective one-dimensional approach of Blakie et al. [36], 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 δasis mostly determined by the scaling of the inverse growth rate of the roton mode, namely, τSS(δas)≃αδasτR(δas), with αδas being a proportionality factor fixed by the initial 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 than by the energy difference between the initial superfluid state and the final supersolid state, as it was speculated in Refs. [14,44]. The present analysis—which admits a straightforward extension also to other cases besides a quench (e.g., by means of a time rescaling as discussed in Ref. [7])—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. ACKNOWLEDGMENTS This work was supported by Grant No. PID2021126273NB-I00 funded by MCIN/AEI/10.13039/5011000 11033 by “ERDF A way of making Europe,” by the Basque Government through Grant No. IT1470-22, and by the European Research Council through the Advanced Grant “Supersolids” (Grant No. 101055319). [1] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Rep. Prog. Phys. 72, 126401 (2009). [2] M. A. Baranov, M. Dalmonte, G. Pupillo, and P. Zoller, Chem. Rev. 112, 5012 (2012). [3] V. Yukalov, Laser Phys. 28, 053001 (2018). [4] L. Santos, G. V. Shlyapnikov, and M. Lewenstein, Phys. Rev. Lett. 90, 250403 (2003). [5] D. H. J. O’Dell, S. Giovanazzi, and G. Kurizki, Phys.Rev.Lett. 90, 110402 (2003). [6] S. Ronen, D. C. E. Bortolotti, and J. L. Bohn, Phys.Rev.Lett. 98, 030406 (2007). [7] L. Chomaz, R. M. W. van Bijnen, D. Petter, G. Faraoni, S. Baier, J. H. Becher, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Nat. Phys. 14, 442 (2018). [8] D. Petter, G. Natale, R. M. W. van Bijnen, A. Patscheider, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev. Lett. 122, 183401 (2019). [9] J.-N. Schmidt, J. Hertkorn, M. Guo, F. Böttcher, M. Schmidt, K. S. H. Ng, S. D. Graham, T. Langen, M. Zwierlein, and T. Pfau, Phys.Rev.Lett.126, 193002 (2021). [10] S. Giovanazzi, D. O’Dell, and G. Kurizki, Phys.Rev.Lett.88, 130402 (2002). [11] S. M. Roccuzzo and F. Ancilotto, Phys. Rev. A 99, 041601(R) (2019). [12] H. Kadau, M. Schmitt, M. Wenzel, C. Wink, T. Maier, I. FerrierBarbut, and T. Pfau, Nature (London) 530, 194 (2016). [13] L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Phys. Rev. Lett. 122, 130405 (2019). [14] F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Phys.Rev.X9, 011051 (2019). [15] Y. Kora and M. Boninsegni, J. Low Temp. Phys. 197, 337 (2019). [16] L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen et al.,Phys. Rev. X 9, 021012 (2019). [17] M. Wenzel, F. Böttcher, T. Langen, I. Ferrier-Barbut, and T. Pfau, Phys.Rev.A96, 053630 (2017). [18] L. Klaus, T. Bland, E. Poli, C. Politi, G. Lamporesi, E. Casotti, R. N. Bisset, M. J. Mark, and F. Ferlaino, Nat. Phys. 18, 1453 (2022). [19] S. V. Andreev, Phys.Rev.B95, 184519 (2017). [20] L. Pollet, Nature (London) 569, 494 (2019). [21] T. Donner, Physics 12, 38 (2019). [22] L. Tanzi, S. M. Roccuzzo, E. Lucioni, F. Famà, A. Fioretti, C. Gabbanini, G. Modugno, A. Recati, and S. Stringari, Nature (London) 574, 382 (2019). [23] M. Guo, F. Böttcher, J. Hertkorn, J.-N. Schmidt, M. Wenzel, H. P. Büchler, T. Langen, and T. Pfau, Nature (London) 574, 386 (2019). [24] G. Natale, R. M. W. van Bijnen, A. Patscheider, D. Petter, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev. Lett. 123, 050402 (2019). [25] L. Tanzi, J. G. Maloberti, G. Biagioni, A. Fioretti, C. Gabbanini, and G. Modugno, Science 371, 1162 (2021). [26] J. Hertkorn, J.-N. Schmidt, F. Böttcher, M. Guo, M. Schmidt, K. S. H. Ng, S. D. Graham, H. P. Büchler, T. Langen, M. Zwierlein et al.,Phys.Rev.X11, 011037 (2021). [27] D. Petter, A. Patscheider, G. Natale, M. J. Mark, M. A. Baranov, R. van Bijnen, S. M. Roccuzzo, A. Recati, B. Blakie, D. Baillie et al.,Phys. Rev. A 104, L011302 (2021). [28] M. A. Norcia, C. Politi, L. Klaus, E. Poli, M. Sohmen, M. J. Mark, R. N. Bisset, L. Santos, and F. Ferlaino, Nature (London) 596, 357 (2021). [29] F. Böttcher, J.-N. Schmidt, J. Hertkorn, K. S. H. Ng, S. D. Graham, M. Guo, T. Langen, and T. Pfau, Rep. Prog. Phys. 84, 012403 (2021). [30] M. Sohmen, C. Politi, L. Klaus, L. Chomaz, M. J. Mark, M. A. Norcia, and F. Ferlaino, Phys.Rev.Lett.126, 233401 (2021). [31] G. Biagioni, N. Antolini, A. Alaña, M. Modugno, A. Fioretti, C. Gabbanini, L. Tanzi, and G. Modugno, Phys. Rev. X 12, 021019 (2022). [32] T. Bland, E. Poli, C. Politi, L. Klaus, M. A. Norcia, F. Ferlaino, L. Santos, and R. N. Bisset, Phys. Rev. Lett. 128, 195302 (2022). [33] J. Sánchez-Baena, C. Politi, F. Maucher, F. Ferlaino, and T. Pohl, Nat. Commun. 14, 1868 (2023). [34] Y.-C. Zhang, F. Maucher, and T. Pohl, Phys.Rev.Lett.123, 015301 (2019). [35] S. C. Schuster, P. Wolf, S. Ostermann, S. Slama, and C. Zimmermann, Phys. Rev. Lett. 124, 143602 (2020). 033316-6 SUPERSOLID FORMATION IN A DIPOLAR CONDENSATE … PHYSICAL REVIEW A 108, 033316 (2023) [36] P. B. Blakie, D. Baillie, and S. Pal, Commun. Theor. Phys. 72, 085501 (2020). [37] S. M. Roccuzzo, A. Gallemí, A. Recati, and S. Stringari, Phys. Rev. Lett. 124, 045702 (2020). [38] A. Gallemí, S. M. Roccuzzo, S. Stringari, and A. Recati, Phys. Rev. A 102, 023322 (2020). [39] P. B. Blakie, D. Baillie, L. Chomaz, and F. Ferlaino, Phys. Rev. Res. 2, 043318 (2020). [40] M. N. Tengstrand, D. Boholm, R. Sachdeva, J. Bengtsson, and S. M. Reimann, Phys. Rev. A 103, 013313 (2021). [41] Y.-C. Zhang, T. Pohl, and F. Maucher, Phys. Rev. A 104, 013310 (2021). [42] J. Hertkorn, J.-N. Schmidt, M. Guo, F. Böttcher, K. S. H. Ng, S. D. Graham, P. Uerlings, H. P. Büchler, T. Langen, M. Zwierlein et al.,Phys.Rev.Lett.127, 155301 (2021). [43] J. Hertkorn, J.-N. Schmidt, M. Guo, F. Böttcher, K. S. H. Ng, S. D. Graham, P. Uerlings, T. Langen, M. Zwierlein, and T. Pfau, Phys.Rev.Res.3, 033125 (2021). [44] A. Alaña, N. Antolini, G. Biagioni, I. L. Egusquiza, and M. Modugno, Phys. Rev. A 106, 043313 (2022). [45] S. M. Roccuzzo, S. Stringari, and A. Recati, Phys. Rev. Res. 4, 013086 (2022). [46] T. Ilg and H. P. Büchler, Phys.Rev.A107, 013314 (2023). [47] J. C. Smith, D. Baillie, and P. B. Blakie, Phys. Rev. A 107, 033301 (2023). [48] A. J. Leggett, Phys.Rev.Lett.25, 1543 (1970). [49] G. V. Chester, Phys.Rev.A2, 256 (1970). [50] E. P. Gross, J. Math. Phys. 4, 195 (1963). [51] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation and Superfluidity (Oxford University Press, New York, 2016), Vol. 164. [52] E. P. Gross, Phys. Rev. 106, 161 (1957). [53] D. Kirzhnitis and Y. A. Nepomnyashchii, Sov. Phys. JETP 32, 1191 (1971). [54] M. Boninsegni and N. V. Prokofev, Rev. Mod. Phys. 84, 759 (2012). [55] S. L. Sondhi and F. J. Burnell, Physics 2, 49 (2009). [56] Y. Pomeau and S. Rica, Phys. Rev. Lett. 72, 2426 (1994). [57] T. Macrì, F. Maucher, F. Cinti, and T. Pohl, Phys. Rev. A 87, 061602(R) (2013). [58] Z.-K. Lu, Y. Li, D. S. Petrov, and G. V. Shlyapnikov, Phys. Rev. Lett. 115, 075303 (2015). [59] N. Sepúlveda, C. Josserand, and S. Rica, Phys. Rev. B 77, 054513 (2008). [60] S. Ronen, D. C. E. Bortolotti, and J. L. Bohn, Phys. Rev. A 74, 013623 (2006). [61] R. Schützhold, M. Uhlmann, Y. Xu, and U. R. Fisher, Int. J. Mod. Phys. B 20, 3555 (2006). [62] F. Wächtler and L. Santos, Phys.Rev.A94, 043618 (2016). [63] M. Schmitt, M. Wenzel, F. Böttcher, I. Ferrier-Barbut, and T. Pfau, Nature (London) 539, 259 (2016). [64] For details of the numerical methods see Ref. [44]. [65] A different approach, based on the Thomas-Fermi approximation, is also possible (see Ref. [7]). In this work, we have chosen the Gaussian ansatz due to its generality, and we do not expect any qualitative differences to arise from the choice of the factorization model. [66] In this case the fit returns τ≃3.5ms,whichisverycloseto the expected value, see Fig. 4. We also have A≃2×10−5, indicating a very low initial population of the roton mode, and γ≃2.6. 033316-7