Oscillations of a quasi-one-dimensional dipolar supersolid

The spatial period of emerging density modulations in a dipolar condensate, which can transform in suitable conditions into SS, can be defined from the dispersion curve of elementary excitations. To derive the dispersion relation we consider the ground state wave function of the form

where $n_{0}$ is the constant background density, $\theta$ is the uniform phase

Next, we introduce weak perturbation to the ground state $\psi(x,t)=\psi_{0}+\delta\psi(x,t)$, assuming $\delta\psi(x,t)=\eta(x,t)e^{i\theta t}$ and $|\eta(x,t)|^{2}\ll n_{0}$. Substituting $\psi(x,t)$ into Eq. (2) and linearizing around $\psi_{0}(x,t)$ one obtains

Splitting the perturbation field into real and imaginary parts $\eta(x,t)=u(x,t)+iv(x,t)$, and taking the Fourier transform of the resulting equation, yields the following coupled system for the transformed components

where the overdot denotes the derivative with respect to time, while the tilde stands for the Fourier transform $\tilde{f}(k)=\int_{-\infty}^{\infty}f(x)e^{ikx}dx$, and $\tilde{R}(k)$ is given by

with ${\rm Ei}(z)$ being the exponential integral function [37]. Notice that the above equations are equivalent to a single harmonic oscillator equation

with the frequency

For small values of the wave vector $k$ the above frequency defines the energy spectrum $E(k)=v_{s}k$ of the elementary excitations obtained from the Bogoliubov-de Gennes analysis with a sound velocity

For larger values of $k$ and suitable parameters the dispersion relation (8) features a local minimum at $k=k_{r}$, corresponding to quasi-particles known as rotons, whose existence in dipolar condensates was theoretically predicted [38] and experimentally confirmed [39, 40]. Although initially the roton was linked to local vorticity in superfluid ⁴He [41], nowadays it is viewed as the precursor of a crystallization instability [42].

Figure 1 illustrates the dispersion relation (8) for different sets of parameter values. The real and positive values of the frequency $\Omega(k)>0$ (blue dashed and brown dot-dashed lines in Fig. 1) signify that the density modulations in the condensate are of standing wave type. A small energy gap near $k\sim k_{r}$ implies that the roton quasi-particles are more effectively generated in the condensate compared to other types of excitations [34]. The supersolid emerges as the frequency (8) becomes imaginary (red and black solid lines in Fig. 1) and roton instability sets in [43]. The steady growth of the amplitude of density waves, caused by the roton instability, is counterbalanced by the repulsive LHY term and/or three-body atomic scattering effects. It should be stressed that the balance solely between the attractive two-body ($g>0$) and repulsive dipolar ($g_{dd}<0$) interactions is unstable (corresponding to the curve (a) in Fig. 1). Therefore, the presence of repulsive LHY and/or three-body interaction terms is essential for the existence of stable SS in dipolar condensates. Either of these two higher-order nonlinear terms or a combination of them can stabilize the supersolid, producing very similar stationary density profiles. The spatial period of the supersolid represents its most important characteristic which is defined by the wave vector of rotons $\lambda=2\pi/k_{r}$. Another important parameter is the amplitude of the SS, expressed through its field contrast as $a_{f}={\rm max}(|\psi(x)|)-{\rm min}(|\psi(x)|)$, where ${\rm max}$ (${\rm min}$) denotes the maximum (minimum) over the space.

Below we study the existence of SS in dipolar quantum gases and its dynamics employing the variational approach and numerical integrations of the eGPE (2).

It was pointed out earlier that supersolids in dipolar condensates can exist due to the stabilizing effect of repulsive quantum fluctuations and/or three-body atomic interactions. The existence conditions and some characteristics of the supersolid can be evaluated using the variational approach. A variational theory based on minimization of the GPE energy functional was developed for the former [28, 29] and latter [44] cases, using the cosine modulated trial function $\psi(x)=\sqrt{n}\left(\cos\theta+\sqrt{2}\sin\theta\cos(kx)\right)$, with $n$ being the average linear density of the condensate. The variational parameters $\theta$ and $k$ characterize the amplitude and period of the density modulations, respectively. This ansatz is valid for weakly modulated density waves in the absence of trapping potential.

Below we develop a variational approach for supersolids in dipolar quantum gases using an alternative trial function

where the variational parameters ($A,a,k$) are introduced, meaning the amplitude of the background, amplitude and wave vector of density modulations, respectively. The normalization to the average linear density of the condensate

is adopted. Similar (sine-modulated) trial function was previously employed in Ref. [32]. These trial functions are suitable when the axial confinement is absent ($\alpha=0$). By substitution of the ansatz (10) into the energy density, corresponding to Eq. (2) without the harmonic trap

and integration over the space variable $E=\lim_{L\to\infty}\left(L^{-1}\int_{-L/2}^{L/2}{\cal E}dx\right)$, one obtains the following GP energy functional

The response function in the Fourier space $\tilde{R}(k)$, defined by Eq. (6), has the following properties

where the prime stands for the derivative with respect to $k$. Using these relations the expression for the amplitude of density modulations $a$ can be derived from minimization of the energy ($dE/dk=0$)

and the vanishing of the roton gap ($\Omega(k)=0$) at some $k>0$. The supersolid amplitude Eq. (15) implicitly depends on the system parameters through $k=k_{r}=f(g,p,\gamma,g_{dd},n)$. Calculations according to Eq. (8) and Eq. (15) for parameter values specified for the black line (d) in Fig. 1 gives for the roton minimum $k_{r}=4.1$, and SS amplitude $a=0.4$ at average linear density $n=12.12$, which is in good agreement with the results of eGPE ($a=0.46$).

The validity of the variational Eq. (15) is restricted to weak (sinusoidally modulated) supersolids. When it is composed of strongly localized spikes of density, the trial function (10) is no longer appropriate.

The majority of theoretical research concerns the infinite supersolid, which can be emulated in numerical simulations using periodic boundary conditions. Figure 2 illustrates the ground state of the supersolid and its oscillations induced by slight variation of the strength of dipolar interactions.

It should be pointed out, that in case of periodic boundary conditions, variation of the strength of atomic interactions leads to excitation of the amplitude mode only, leaving the width Eq.(16) unaffected. Owing to a superfluid nature of the system and equal heights of density peaks, the SS performs undamped oscillations of its amplitude with a single frequency $\omega\simeq 1.9$ (see Fig. 2b).

Oscillations of the effective width can be induced using the external potential providing a small shift of the density peaks from their equilibrium positions, meantime not violating the periodic boundary conditions. In numerical simulations we employ the following potential

with parameter values $V_{0}=0.5$, $\beta=2.24$, $w=3$. In addition to the displacement of the peak positions, their heights also become unequal. In Fig. 3 we show the ground state of the supersolid in presence of the potential (17), and its dynamics caused by subsequent elimination of the potential. The excitation of both the superfluid and crystalline modes is evident. The Fourier analysis of the time dependent width Eq. (16) reveals two distinct frequencies, associated with the superfluid ($\omega_{s}\sim 1$) and crystalline ($\omega_{c}\sim 2.4$) components of the system.

If the density peaks of the supersolid become unequal due to the external potential, a notable mixing of oscillations occurs. As a consequence, additional small peaks arise in the Fourier spectrum (Fig. 3d).

Supersolids in confined geometries exhibit more complex dynamics. Below we consider the SS confined to a soft-wall box potential

where the parameters $V_{0},h,w$ characterize its strength, spatial extent, and width of the transition region.

In experiments the potential (18) can be created by laser beams [47, 48]. To study oscillations of the supersolid first we create its ground state in the box potential, as shown in Fig. 4. The peculiarity of this case is the existence of density spikes near the trap edges [49]. The ground state with a density peak (drop) at $x=0$ emerges if the trap can accommodate an odd (even) number $m$ of supersolid periods $2h=m\lambda$.

Sudden reduction of the repulsive dipolar interactions gives rise to flow of superfluid matter into the central region of the trap, which shows up as raising (and subsequent oscillations) of the central density peaks (see Fig. 4b). This is a clear manifestation of the key property of supersolids, namely free flow of the superfluid component through its crystalline phase.

To identify the frequencies of particular modes we perform Fourier analysis of the time-dependent width $W(t)$ calculated according to Eq.(16), with $\psi(x,t)$ being the numerical solution of the eGPE (2). The Fourier transform reveals three frequencies (see Fig. 4d). The lower frequency is associated with the superfluid component, while other two higher frequencies are linked to the crystalline component of the SS. The additional frequency originates from two pronounced density peaks (edge states) of dipolar BEC confined to a box-potential (see Fig. 4a). These density peaks can be viewed as a separate crystal component, oscillating at a different frequency than the main crystalline phase, occupying the central part of the box-trap.

The parabolic potential is a frequently used setting in experiments with ultra-cold quantum gases. The ground state of the supersolid confined to a parabolic trap, numerically produced by the Pitaevskii phenomenological damping procedure, is depicted in Fig. 5. The figure gives evidence that the stronger is the parabolic trap, the more amount of superfluid matter is expelled to the trap periphery. A similar effect was noticed for the box trap where the edge states with increased matter density emerge even for parameter regimes, where the SS phase is absent. As a consequence, the lateral density peaks in the parabolic trap become less pronounced.

To study the collective dynamics of the SS in a parabolic trap we introduce its ground state (red solid line in Fig. 5a) into the eGPE as initial condition and reduce the trap strength from $\alpha=0.5$ to $\alpha=0.4$. The resulting collective dynamics of the supersolid is shown through the density plot in Fig. 5b. Notable periodic variation of the inter-peak separation can be seen. The average width of the supersolid $W(t)$ also shows complex dynamics, giving evidence of the superposition of few oscillations (Fig. 5c). The Fourier transform of the effective width ${\cal F}(W)$ reveals two main frequencies (Fig. 5d), corresponding to oscillations of the superfluid component ($\omega_{s}\sim 1$) and crystalline phase ($\omega_{c}\sim 2.3$), respectively.

Figure 5 represents the example of SS with dominating crystalline component, which shows a stronger peak in the Fourier spectrum than that of the superfluid component. Now we consider the opposite situation with dominating superfluid phase, which is shown in Fig. 6. The ground state is numerically obtained by the Pitaevskii phenomenological damping procedure starting from the uniform condensate with linear density $n_{0}=A^{2}=8.128$. For selected parameter values the density modulations are less pronounced in the central region of the trap, while notable density peaks show up near the borders, similar to edge states in the box-trap (see Fig. 4a).

To induce oscillations of the condensate initially residing in the ground state of the parabolic trap, we slightly change its strength $\alpha$. Figure 6 c illustrates the variation of the effective width of the condensate caused by a 5 % reduction of the parameter $\alpha$. Here a superposition of two frequencies can be noticed. The Fourier transform shown in Fig. 6 d reveals these two frequencies, associated with the superfluid and crystalline components. The crystalline phase, although quite weakly, shows up through the density peaks near trap edges. By comparing the results shown in figures 5 and 6 we conclude that when the SS is perturbed, its crystalline (superfluid) component oscillates with higher (lower) frequency, in accordance with experimental findings [8].

Now we estimate the parameters of our model using the quantities relevant to experiments [8, 33]. Consider ¹⁶²Dy atoms, whose magnetic dipole moment, background $s$-wave scattering length and atomic mass are $\mu_{Dy}=10\,\mu_{B}$, $a_{bg}=-130\,a_{0}$, $m_{Dy}=2.7\times 10^{-25}$ kg, respectively, with $\mu_{B}$, $a_{0}$ being the Bohr magneton and Bohr radius. The roton instability sets in, thus supersolid phase emerges, when the $s$-wave scattering length is tuned to $a_{s}=-98a_{0}\simeq 0.005\mu$m. The condensate of $N=4\times 10^{4}$ atoms is held in a tight quasi-1D trap with radial confinement frequency $\omega_{\bot}=2\pi\times 110$ Hz. The corresponding radial harmonic oscillator length, which is the adopted length scale in this work, is $a_{\bot}=0.75\,\mu$m. The unit of time is $\omega_{\bot}^{-1}=1.4\times 10^{-3}$ s. The ratio between the strengths of dipolar and contact interactions, computed with above defined parameters, is found to be $\epsilon_{d}=\mu_{0}\mu_{Dy}^{2}m_{Dy}/12\pi\hbar^{2}a_{s}\simeq 1.34$, therefore we have the dipolar interaction dominated regime. The coefficients of two-body contact interactions and long-range dipolar interactions in physical units are $g_{0}=4\pi\hbar^{2}a_{s}/m_{Dy}=-2.68\times 10^{-51}\,{\rm kg}\cdot{\rm m}^{5}/{\rm s}^{2}$ and $C_{d}=\mu_{0}\mu_{Dy}^{2}=1.08\times 10^{-50}\,{\rm kg}\cdot{\rm m}^{5}/{\rm s}^{2}$, yielding the dimensionless parameter $g_{dd}=C_{d}/g_{0}\simeq-4$. The strength of three-body atomic interactions $K_{3}=K_{r}+i\,K_{i}$ in dysprosium condensate was reported to be in the range $K_{r}=(3.91-19.6)\times 10^{-39}\,\hbar\cdot{\rm m}^{6}/{\rm s}$ for real conservative part and $K_{i}=7.8\times 10^{-42}\,\hbar\cdot{\rm m}^{6}/{\rm s}$ for imaginary part, characterizing the three-body recombination rate [15]. Since the imaginary part of this parameter is smaller than its real part by three orders of magnitude, we have omitted it in the eGPE (2), leaving only the conservative part $g_{2}=K_{r}$. The coefficient of quintic nonlinearity, computed using this value and above defined parameters is found to be in the range $p=-(0.007-0.036)$. The coefficient of quantum fluctuations corresponding to above mentioned parameter values is estimated as $\gamma\sim-0.05$.

Numerical simulations are performed using the integration domain of length ${\cal L}=(9-12)\lambda\,a_{\bot}\simeq(10-14)\,\mu$m. The amplitude of the background wave in dimensionless units is defined as $A=\sqrt{2\,|a_{s}|\,N/{\cal L}}\simeq 6$, where $N$ stands for the number of atoms in the condensate. The linear density of the condensate is $n_{l}=N/{\cal L}\sim 4\times 10^{3}\mu m^{-1}$.

The above estimates for dimensionless parameters obtained from experimentally possible quantities are in the range of values used in our numerical simulations. Since the radial confinement length is much greater than the atomic scattering length ($a_{\bot}\gg a_{s}$), usage of parameter values relevant to 3D geometry for numerical simulations in quasi-1D setting is justified. Some deviations can be adjusted by changing the tunable parameters of the system, such as the $s$-wave scattering length, frequency of radial confinement, and the strength of dipolar interactions. As for the roles of the repulsive three-body atomic collisions and quantum fluctuations, these are two possible factors responsible for the existence of stable supersolids [44, 8], and for a complete description of the problem, both terms should be taken into account.