Variational theory for the ground state and collective excitations of an elongated dipolar condensat
本站小编 Free考研考试/2022-01-02
P Blair Blakie,1, D Baillie, Sukla PalDodd-Walls Centre for Photonic and Quantum Technologies, New Zealand and Department of Physics, University of Otago, Dunedin 9016, New Zealand
First author contact:1Author to whom any correspondence should be addressed. Received:2020-02-21Revised:2020-04-21Accepted:2020-05-15Online:2020-07-20
Abstract We develop a variational theory for a dipolar condensate in an elongated (cigar shaped) confinement potential. Our formulation provides an effective one-dimensional extended meanfield theory for the ground state and its collective excitations. We apply our theory to investigate the properties of rotons in the system comparing the variational treatment to a full numerical solution. We consider the effect of quantum fluctuations on the scattering length at which the roton excitation softens to zero energy. Keywords：dipolar Bose–Einstein condensate;ground states;collective excitations
PDF (857KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite Cite this article P Blair Blakie, D Baillie, Sukla Pal. Variational theory for the ground state and collective excitations of an elongated dipolar condensate. Communications in Theoretical Physics, 2020, 72(8): 085501- doi:10.1088/1572-9494/ab95fa
1. Introduction
Bose–Einstein condensates of highly magnetic atoms, such as chromium, erbium, and dysprosium [1–3], realise a dilute quantum system with long-ranged and anisotropic dipole–dipole interactions (DDIs). Recently experiments with such dipolar condensates have observed a roton excitation [4, 5], i.e.a local minimum in the excitation dispersion relation of the condensate at a non-zero wave vector. The roton excitation, originally introduced in the study of super fluid Helium [6], has been of significant theoretical interest in dipolar condensates where it emerges from an interplay between the DDIs and confinement (e.g. see [7–15]). Much of the theoretical attention has focused on the case of a system confined to a planar (pancake) geometry with the dipoles aligned along the tightly confined direction. Experiments instead have realised roton excitations in an elongated (cigar) geometry with the dipoles oriented along one of the tightly confined directions (e.g. see figure 1). While the planar case can have cylindrical symmetry, the elongated system does not and in general requires a full three-dimensional (3D) calculation. Such calculations for the system ground states and its collective excitations demand significant computational resources, mostly arising from the large and dense numerical grids required to carefully resolve the singular DDI potential.
Figure 1.
New window|Download| PPT slide Figure 1.Schematic of the system geometry we consider in this paper: a condensate of atoms with magnetic dipoles aligned along the y-axis by an external magnetic field, and with tighter confinement in the xy-plane relative to the z-axis.
The theoretical description of a dipolar condensate is usually provided by meanfield theory. However, recent developments in the field have revealed that quantum fluctuations can play an important role in the regime of strong DDIs, for example leading to the formation of stable quantum droplets and supersolids (e.g. see [16–23]). Extended meanfield theory includes the leading order effects of quantum fluctuations within a local density approximation [17, 18, 24]. In this formalism the stationary states of the system are described by the extended Gross–Pitaevskii equation (eGPE), with the collective excitations described by the associated Bogoliubov–de Gennes (BdG) equations [25, 26].
In this paper we develop an approximation that allows us to reduce the 3D extended meanfield theory to a tractable one-dimensional (1D) form. Our approach is to use a Gaussian ansatz to describe the tightly confined transverse directions, and by integrating this out we obtain an effective 1D form of the theory, albeit with some variational parameters from the Gaussian. Such an approach is a rather obvious path to take and has been used for non-dipolar condensates (e.g. see [27]). However, for the dipolar case the Gaussian cannot be integrated out against the DDI potential to yield an analytic result for the interaction term, except for the special case where the Gaussian is isotropic. An important result of this paper is that we introduce a useful approximate analytic result for general Gaussian. Based on this we develop a simple 1D effective theory for the stationary states and the collective excitations of the system. We emphasise that this description is for a 3D dipolar condensate, and is not applicable to the true 1D regime where interactions remain weak compared to the confinement and the quantum fluctuations take a different form [28]. Indeed, in the regime of interest (e.g. where rotons occur) the interaction energy scale is typically larger than the transverse confinement energy and the transverse degrees of freedom must be treated variationally. Furthermore, in this regime magnetostriction effects can be large, causing the Gaussian to distort appreciably from the geometry imposed by the confining potential. This occurs because the DDIs are anisotropic and the energy of the system is reduced by having more particles in a relative head-to-tail orientation.
We compare our results from the variational 1D theory to full numerical solutions of the 3D problem for both ground states and the excitation spectrum. We use our theory to predict the value of the s-wave scattering length (readily adjusted in experiments using Feshbach resonances) where the roton mode goes soft (i.e. to zero energy), and the associated value of the wave vector of the soft mode. By comparing to results that exclude the quantum fluctuation term we can assess the effect of quantum fluctuations on the roton properties.
The outline of the paper is as follows. In section 2 we introduce extended meanfield theory and our approach for simplifying it to an effective variational 1D theory. The main results are presented in section 3, before we conclude in section 4.
2. Formalism
2.1. Extended meanfield theory for dipolar condensates
2.1.1. 3D eGPE
The system of interest is a dilute Bose gas of atoms with a magnetic moment μm polarised along the y-axis. The extended meanfield theory for this system identifies stationary states of the matterwave field &PSgr;0 as solutions of the eGPE $\mu {{\rm{\Psi }}}_{0}={{ \mathcal L }}_{3{\rm{D}}}{{\rm{\Psi }}}_{0}$, where$ \begin{eqnarray}{{ \mathcal L }}_{3{\rm{D}}}{{\rm{\Psi }}}_{0}=\left[-\displaystyle \frac{{{\hslash }}^{2}{{\rm{\nabla }}}^{2}}{2m}+V({\boldsymbol{x}})+{\rm{\Phi }}({\boldsymbol{x}})+{\gamma }_{\mathrm{QF}}| {{\rm{\Psi }}}_{0}{| }^{3}\right]{{\rm{\Psi }}}_{0}.\end{eqnarray}$Here μ is the chemical potential and$ \begin{eqnarray}{\rm{\Phi }}({\boldsymbol{x}})=\int {\rm{d}}{\boldsymbol{x}}^{\prime} \,U({\boldsymbol{x}}-{\boldsymbol{x}}^{\prime} )| {{\rm{\Psi }}}_{0}({\boldsymbol{x}}^{\prime} ){| }^{2},\end{eqnarray}$is the interaction term, with interaction potential$ \begin{eqnarray}U({\boldsymbol{r}})={g}_{s}\delta ({\boldsymbol{r}})+\displaystyle \frac{3{g}_{{dd}}}{4\pi {r}^{3}}\left(1-3\displaystyle \frac{{y}^{2}}{{r}^{2}}\right),\end{eqnarray}$where the contact interaction coupling constant is ${g}_{s}=4\pi {{\hslash }}^{2}{a}_{s}/m$, as is the s-wave scattering length, the DDI coupling constant is ${g}_{{dd}}=4\pi {{\hslash }}^{2}{a}_{{dd}}/m$, and ${a}_{{dd}}\,\equiv m{\mu }_{0}{\mu }_{m}^{2}/12\pi {{\hslash }}^{2}$ is the dipole length. This theory includes the leading order quantum fluctuations in the local density approximation, where the coefficient of this term is [18, 24, 29]$ \begin{eqnarray}{\gamma }_{\mathrm{QF}}=\displaystyle \frac{32}{3}{g}_{s}\sqrt{\displaystyle \frac{{a}_{s}^{3}}{\pi }}(1+\tfrac{3}{2}{\epsilon }_{{dd}}^{2}),\end{eqnarray}$with ${\epsilon }_{{dd}}={a}_{{dd}}/{a}_{s}$. The atoms are taken to be confined by an external potential $V({\boldsymbol{x}})=\tfrac{1}{2}m{\sum }_{\nu =x,y,z}{\omega }_{\nu }^{2}{\nu }^{2}$. Here we focus our attention on the regime used in experiments to observe roton excitations: a cigar shaped trap with ${\omega }_{x},{\omega }_{y}\gg {\omega }_{z}$ (including the pure tube case with ωz=0). The energy functional associated with the 3D eGPE is$ \begin{eqnarray}E=\int {\rm{d}}{\boldsymbol{x}}\,{{\rm{\Psi }}}_{0}^{* }\left[-\displaystyle \frac{{{\hslash }}^{2}{{\rm{\nabla }}}^{2}}{2m}+V({\boldsymbol{x}})+\tfrac{1}{2}{\rm{\Phi }}({\boldsymbol{x}})+\tfrac{2}{5}{\gamma }_{\mathrm{QF}}| {{\rm{\Psi }}}_{0}{| }^{3}\right]{{\rm{\Psi }}}_{0}.\end{eqnarray}$
2.1.2. BdG theory of excitations
The collective excitations of this system are Bogoliubov quasiparticles, which can be obtained by linearising the time-dependent GPE ${\rm{i}}{\hslash }\dot{{\rm{\Psi }}}={{ \mathcal L }}_{3{\rm{D}}}{\rm{\Psi }}$ about a stationary state as$ \begin{eqnarray}{\rm{\Psi }}={{\rm{e}}}^{-{\rm{i}}\mu t/\hslash }\left[{{\rm{\Psi }}}_{0}+\displaystyle \sum _{\nu }\left({\lambda }_{\nu }{U}_{\nu }{{\rm{e}}}^{-{\rm{i}}{\epsilon }_{\nu }t/\hslash }-{\lambda }_{\nu }^{* }{V}_{\nu }^{* }{{\rm{e}}}^{{\rm{i}}{\epsilon }_{\nu }^{* }t/\hslash }\right)\right],\end{eqnarray}$where λν is a small complex amplitude. The quasiparticle modes Uν, Vν and energies εν satisfy the BdG equations [25]$ \begin{eqnarray}\left(\begin{array}{cc}{{ \mathcal L }}_{3{\rm{D}}}-\mu +X & -X\\ X & -({{ \mathcal L }}_{3{\rm{D}}}-\mu +X)\end{array}\right)\left(\begin{array}{c}{U}_{\nu }\\ {V}_{\nu }\end{array}\right)={\epsilon }_{\nu }\left(\begin{array}{c}{U}_{\nu }\\ {V}_{\nu }\end{array}\right),\end{eqnarray}$where X is the exchange operator given by$ \begin{eqnarray}{Xf}\equiv {{\rm{\Psi }}}_{0}\int {\rm{d}}{\boldsymbol{x}}^{\prime} U({\boldsymbol{x}}-{\boldsymbol{x}}^{\prime} )f({\boldsymbol{x}}^{\prime} ){{\rm{\Psi }}}_{0}^{* }({\boldsymbol{x}}^{\prime} )+\tfrac{3}{2}{\gamma }_{\mathrm{QF}}| {{\rm{\Psi }}}_{0}{| }^{3}f.\end{eqnarray}$
2.2. Reduction to an effective 1D eGPE
2.2.1. General approach
We approximate the 3D solutions in the elongated trap to be of the separable form ${{\rm{\Psi }}}_{0}({\boldsymbol{x}})={\psi }_{0}(z)\chi ({\boldsymbol{\rho }})$, where χ is the transverse mode function with ${\boldsymbol{\rho }}=(x,y)$ being the radial coordinate vector and $\int {\rm{d}}{\boldsymbol{\rho }}\,| \chi {| }^{2}=1$. Integrating out the transverse directions we obtain the 1D eGPE operator for the axial wavefunction ψ0:$ \begin{eqnarray}{{ \mathcal L }}_{z}=\displaystyle \int {\rm{d}}{\boldsymbol{\rho }}\,{\chi }^{* }{{ \mathcal L }}_{3{\rm{D}}}\chi ,\end{eqnarray}$$ \begin{eqnarray}={{ \mathcal E }}_{\perp }-\displaystyle \frac{{{\hslash }}^{2}}{2m}\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{z}^{2}}+\tfrac{1}{2}m{\omega }_{z}^{2}{z}^{2}+{{\rm{\Phi }}}_{z}(z)+{\gamma }_{\mathrm{QF}}{\gamma }_{\perp }| {\psi }_{0}{| }^{3},\end{eqnarray}$where ${\gamma }_{\perp }\equiv \int {\rm{d}}{\boldsymbol{\rho }}| \chi {| }^{5}$, and$ \begin{eqnarray}{{ \mathcal E }}_{\perp }=\int {\rm{d}}{\boldsymbol{\rho }}\,{\chi }^{* }\left[-\displaystyle \frac{{{\hslash }}^{2}{{\rm{\nabla }}}_{{\boldsymbol{\rho }}}^{2}}{2m}+\tfrac{1}{2}m({\omega }_{x}^{2}{x}^{2}+{\omega }_{y}^{2}{y}^{2})\right]\chi .\end{eqnarray}$The effective interaction term is$ \begin{eqnarray}{{\rm{\Phi }}}_{z}(z)={{ \mathcal F }}_{z}^{-1}\left\{{\tilde{U}}_{z}({k}_{z}){{ \mathcal F }}_{z}\{| {\psi }_{0}{| }^{2}\right\},\end{eqnarray}$where we have introduced the effective 1D k-space interaction kernel$ \begin{eqnarray}{\tilde{U}}_{z}({k}_{z})=\int \displaystyle \frac{{\rm{d}}{{\boldsymbol{k}}}_{\rho }}{{\left(2\pi \right)}^{2}}\tilde{U}({\boldsymbol{k}}){\left|{{ \mathcal F }}_{{\boldsymbol{\rho }}}\{| \chi {| }^{2}\}\right|}^{2}.\end{eqnarray}$In the above results $\tilde{U}({\boldsymbol{k}})$ is the Fourier transform of equation (3), given by$ \begin{eqnarray}\tilde{U}({\boldsymbol{k}})={g}_{s}+{g}_{{dd}}\left(3\displaystyle \frac{{k}_{y}^{2}}{{k}^{2}}-1\right),\end{eqnarray}$${{ \mathcal F }}_{z}\{f\}$ denotes the 1D Fourier transform $f(z)\to \tilde{f}({k}_{z})$ and ${{ \mathcal F }}_{{\boldsymbol{\rho }}}\{f\}$ denotes the 2D Fourier transform $f({\boldsymbol{\rho }})\to \tilde{f}({{\boldsymbol{k}}}_{\rho })$. The associated energy functional takes the form$ \begin{eqnarray}\begin{array}{rcl}E & = & \displaystyle \int {\rm{d}}z\,{\psi }_{0}^{* }\left[{{ \mathcal E }}_{\perp }-\displaystyle \frac{{{\hslash }}^{2}}{2m}\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{z}^{2}}+\tfrac{1}{2}m{\omega }_{z}^{2}{z}^{2}\right]{\psi }_{0}\\ & & +\displaystyle \int {\rm{d}}z\,{\psi }_{0}^{* }\left[\tfrac{1}{2}{{\rm{\Phi }}}_{z}(z)+\tfrac{2}{5}{\gamma }_{\mathrm{QF}}{\gamma }_{\perp }| {\psi }_{0}{| }^{3}\right]{\psi }_{0}.\end{array}\end{eqnarray}$
2.2.2. Anisotropic Gaussian approximation
Here we introduce a convenient analytic form for χ. Our choice is the Gaussian$ \begin{eqnarray}{\chi }_{\sigma }({\boldsymbol{\rho }})=\displaystyle \frac{{{\rm{e}}}^{-(\eta {x}^{2}+{y}^{2}/\eta )/2{l}^{2}}}{\sqrt{\pi }l},\end{eqnarray}$of mean width $l=\sqrt{{l}_{x}{l}_{y}}$, and anisotropy $\eta ={l}_{y}/{l}_{x}$, where lx (ly) is the 1/e half width of $| {\chi }_{\sigma }{| }^{2}$ along the x-axis (y-axis). We use σ to collectively denote the variational parameters σ={l, η}, which are determined by minimising the system energy.
Using ${\chi }_{\sigma }$ we can analytically evaluate key terms in the 1D theory. First we denote ${{ \mathcal E }}_{\perp }$ evaluated with χσ as$ \begin{eqnarray}{{ \mathcal E }}_{\sigma }(l,\eta )=\displaystyle \frac{{{\hslash }}^{2}}{4{{ml}}^{2}}\left(\eta +\displaystyle \frac{1}{\eta }\right)+\displaystyle \frac{{{ml}}^{2}}{4}\left(\displaystyle \frac{{\omega }_{x}^{2}}{\eta }+{\omega }_{y}^{2}\eta \right).\end{eqnarray}$Similarly, ${\gamma }_{\perp }\to {\gamma }_{\sigma }=\tfrac{2}{5{\pi }^{3/2}{l}^{3}}$. We are unaware of a general analytic result for ${\tilde{U}}_{z}({k}_{z})$ evaluated using χσ, which we denote as ${\tilde{U}}_{\sigma }({k}_{z})$. However, we have obtained the useful approximate result$ \begin{eqnarray}{\tilde{U}}_{\sigma }({k}_{z})=\displaystyle \frac{{g}_{s}}{2\pi {l}^{2}}+\displaystyle \frac{{g}_{{dd}}}{2\pi {l}^{2}}\left\{\displaystyle \frac{3[{Q}_{\sigma }^{2}{{\rm{e}}}^{{Q}_{\sigma }^{2}}\mathrm{Ei}(-{Q}_{\sigma }^{2})+1]}{1+\eta }-1\right\},\end{eqnarray}$with Ei being the exponential integral2This can also be written in terms of the incomplete Gamma function Γ using $\mathrm{Ei}(-x)=-{\rm{\Gamma }}(0,x)$ for x>0 (see [30]). , and ${Q}_{\sigma }\equiv \tfrac{1}{\sqrt{2}}{k}_{z}{\eta }^{1/4}l$.
2.2.3. Justification for equation (18)
For the particular case of an isotropic Gaussian (i.e. η=1) equation (18) is exact (see [30–33]). While a general analytic result for $\eta \ne 1$ is unavailable, we can calculate the limiting behaviour$ \begin{eqnarray}{\tilde{U}}_{\sigma }({k}_{z})=\displaystyle \frac{{g}_{s}}{2\pi {l}^{2}}+\left\{\begin{array}{ll}\tfrac{{g}_{{dd}}}{2\pi {l}^{2}}\tfrac{2-\eta }{1+\eta }, & {k}_{z}\to 0\\ -\tfrac{{g}_{{dd}}}{2\pi {l}^{2}}, & {k}_{z}\to \infty \end{array}\right..\end{eqnarray}$We have arrived at result (18) by inspection and numerical experiment: it satisfies the required limiting behaviour and reduces to the exact isotropic result at η=1.
In figure 2 we compare the accuracy of equation (18) to a full calculation of the kernel obtained by numerically evaluating (13) with $\chi \to {\chi }_{\sigma }$. To ensure the numerical calculation is accurate it is performed using a large and dense two-dimensional transverse grid of points and using a cutoff k-space DDI potential to avoid finite size boundary effects (e.g. see [34]). The results in figure 2 show that while our approximate analytic result (18) is not identical to the numerical result, it is generally in very good agreement over a wide range of η values (i.e. $0.2\lesssim \eta \lesssim 20$). Note that ground states with η<1 can occur due to the confinement (i.e. when ωy > ωx), and are also favoured by the interactions when gdd<0, which can be arranged by rotationally tuning the dipoles [35, 36]. We expect that for the regimes of interest the error associated with making the Gaussian approximation is much more significant than any additional error introduced by using equation (18) to describe its interactions.
Figure 2.
New window|Download| PPT slide Figure 2.(a) Comparison of the (dotted lines) analytic result (18) to (solid lines) numerically calculated ${\tilde{U}}_{\mathrm{num}}$ (obtained by numerically evaluating equation (13) using $\chi \to {\chi }_{\sigma })$ for the effective k-space kernel. Results shown for several values of η and for gs=0. The exact result for η=1 is also shown (dashed line). (b) The maximum absolute error of the approximation ${\tilde{U}}_{\sigma }$ compared to the ${\tilde{U}}_{\mathrm{num}}$ over the kz-range shown in (a), (in units of ${g}_{{dd}}/4\pi {l}^{2}$).
2.2.4. Variational theory
Here we summarise the results developed in section 2.2.2 and succinctly present the variational 1D eGPE theory that forms the main formalism result of this paper.
The axial orbital ψ0 satisfies the 1D eGPE:$ \begin{eqnarray}\mu {\psi }_{0}={{ \mathcal L }}_{\sigma }{\psi }_{0},\end{eqnarray}$where$ \begin{eqnarray}{{ \mathcal L }}_{\sigma }={{ \mathcal E }}_{\sigma }-\displaystyle \frac{{{\hslash }}^{2}}{2m}\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{z}^{2}}+\tfrac{1}{2}m{\omega }_{z}^{2}{z}^{2}+{{\rm{\Phi }}}_{\sigma }(z)+{g}_{\mathrm{QF}}| {\psi }_{0}{| }^{3},\end{eqnarray}$with ${g}_{\mathrm{QF}}={\gamma }_{\mathrm{QF}}{\gamma }_{\sigma }$, and Φσ is evaluated according to equation (12) but using ${\tilde{U}}_{\sigma }$ in place of ${\tilde{U}}_{z}$.
Since ${{ \mathcal L }}_{\sigma }$ depends on χσ we also need a procedure to obtain the parameters {l, η}. To do this we consider the total system energy per particle:$ \begin{eqnarray}{ \mathcal E }[{\psi }_{0};l,\eta ]={{ \mathcal E }}_{\sigma }(l,\eta )+{{ \mathcal E }}_{z}[{\psi }_{0};l,\eta ],\end{eqnarray}$where$ \begin{eqnarray}\begin{array}{c}\begin{array}{l}{{ \mathcal E }}_{z}[{\psi }_{0};l,\eta ]=\displaystyle \frac{2}{5N}\displaystyle \int {\rm{d}}z\,{g}_{{\rm{QF}}}| {\psi }_{0}{| }^{5}\,\\ \,+\,\displaystyle \frac{1}{N}\displaystyle \int {\rm{d}}z\,{\psi }_{0}^{\ast }\left(-\displaystyle \frac{{\hslash }^{2}}{2m}\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{z}^{2}}+\tfrac{1}{2}m{\omega }_{z}^{2}{z}^{2}+\tfrac{1}{2}{{\rm{\Phi }}}_{\sigma }(z)\right){\psi }_{0}\end{array}\end{array}\end{eqnarray}$is the energy functional for the ψ0 orbital, and $N=\int {\rm{d}}z\,| {\psi }_{0}{| }^{2}$.
For ωz=0 the axial wavefunction can be uniform ${\psi }_{0}\to \sqrt{n}$, where n is the linear density. In this regime the energy per particle (22) reduces to$ \begin{eqnarray}{{ \mathcal E }}_{u}(l,\eta )={{ \mathcal E }}_{\sigma }(l,\eta )+\tfrac{1}{2}n{\tilde{U}}_{\sigma }(0)+\tfrac{2}{5}{g}_{\mathrm{QF}}{n}^{3/2},\end{eqnarray}$i.e.the ground state is determined by minimising a simple nonlinear function.
2.2.5. Quasi-1D theory
Predictions for ψ0 can be made within the quasi-1D approximation using the procedure outlined for the variational theory, but with l and η held fixed to the values for the harmonic oscillator ground state of the transverse confinement, i.e.for ${l}^{2}={\hslash }/m\sqrt{{\omega }_{x}{\omega }_{y}}$ and $\eta =\sqrt{{\omega }_{x}/{\omega }_{y}}$, with ${{ \mathcal E }}_{\sigma }={\hslash }({\omega }_{x}+{\omega }_{y})/2$. We denote the harmonic oscillator state as χho and the associated k-space kernel as ${\tilde{U}}_{\mathrm{ho}}$. This quasi-1D approximation will only be accurate when the interaction terms in ${{ \mathcal E }}_{z}$ remain small compared to ${{ \mathcal E }}_{\sigma }$.
2.2.6. Effective 1D form of the excitations
Making the same shape approximation for the transverse form of the excitations [14] we set ${U}_{\nu }({\boldsymbol{x}})={u}_{\nu }(z){\chi }_{\sigma }({\boldsymbol{\rho }})$ and ${V}_{\nu }({\boldsymbol{x}})={v}_{\nu }(z){\chi }_{\sigma }({\boldsymbol{\rho }})$ and integrating out ${\chi }_{\sigma }({\boldsymbol{\rho }})$ the BdG equations (7) reduce to$ \begin{eqnarray}\left(\begin{array}{cc}{{ \mathcal L }}_{\sigma }-\mu +{X}_{\sigma } & -{X}_{\sigma }\\ {X}_{\sigma } & -({{ \mathcal L }}_{\sigma }-\mu +{X}_{\sigma })\end{array}\right)\left(\begin{array}{c}{u}_{\nu }\\ {v}_{\nu }\end{array}\right)={\epsilon }_{\nu }\left(\begin{array}{c}{u}_{\nu }\\ {v}_{\nu }\end{array}\right),\end{eqnarray}$where ${X}_{\sigma }f={\psi }_{0}{{ \mathcal F }}_{z}^{-1}\left\{{\tilde{U}}_{\sigma }({k}_{z}){{ \mathcal F }}_{z}\{{\psi }_{0}f\right\}+\tfrac{3}{2}{g}_{\mathrm{QF}}{\psi }_{0}^{3}f$.
In general the BdG equations need to be discretised and solved numerically, however for the case of a uniform ground state an analytic solution can be obtained. Here the excitations are plane waves of momentum ${\hslash }{k}_{z}$, i.e.${u}_{\nu }(z)\to {{\rm{u}}}_{{k}_{z}}{{\rm{e}}}^{{\rm{i}}{k}_{z}z}$, ${v}_{\nu }(z)\to {{\rm{v}}}_{{k}_{z}}{{\rm{e}}}^{{\rm{i}}{k}_{z}z}$, ${\epsilon }_{\nu }\to \epsilon ({k}_{z})$, with excitation energy$ \begin{eqnarray}\epsilon ({k}_{z})=\sqrt{{\epsilon }_{0}({k}_{z})\left[{\epsilon }_{0}({k}_{z})+2n{\tilde{U}}_{\sigma }({k}_{z})+3{g}_{\mathrm{QF}}{n}^{3/2}\right]},\end{eqnarray}$where ${\epsilon }_{0}({k}_{z})={{\hslash }}^{2}{k}_{z}^{2}/2m$.
3. Results
3.1. Numerical methods
In this subsection we briefly outline the various numerical methods used to solve for the results we present later.
3.1.1. Uniform cases
For cases without axial trapping (ωz=0) we restrict our attention to the regime where the ground state is uniform and specified by the linear density n.
The variational 1D eGPE theory reduces to minimising the nonlinear function (24) for l and η. The BdG excitation energies are then directly given by evaluating equation (26).
The 3D eGPE reduces to the determining the transverse mode $\chi ({\boldsymbol{\rho }})$. We do this by discretising $\chi ({\boldsymbol{\rho }})$ on a two-dimensional numerical grid and apply discrete Fourier transformations to apply the kinetic energy operator with spectral accuracy, and to evaluate the interaction term Φ. For high accuracy the 3D k-space kernel is cutoff in the transverse direction to the range of the numerical grid (e.g. see [34, 37]). The eGPE is solved using a gradient flow technique [38]. The excitations for this case are of the form of plane waves along z, reducing the BdG equations to a 2D form that can be solved using large-scale eigensolvers (i.e. the implicitly restarted Arnoldi method).
3.1.2. Fully trapped cases
For ${\omega }_{z}\ne 0$ the variational theory (including the quasi-1D theory) involves solving for ψ0 on a 1D numerical grid. We use a set of equally spaced points allowing us to use discrete Fourier transformations to evaluate the kinetic energy operator and the interaction term Φσ. To improve accuracy we implement an axial cutoff of the k-space kernel ${\tilde{U}}_{\sigma }$: this is obtained by Fourier transforming the real-space interaction potential [30] restricted to the z-spatial range of the grid used for the numerical calculation. The orbital ψ0 is solved using a gradient flow technique for given values of l and η, thus determining a minimum energy solution of ${{ \mathcal E }}_{z}$ (23). An optimisation scheme is used to adjust l and η, then ψ0 is solved with the new parameters, and this procedure iterates until the minimum of the full energy functional (22) is found.
The 3D ground states are obtained using 3D numerical grids and discrete Fourier transforms A cylindrically cutoff k-space kernel is used to improve accuracy of the Φ evaluation. The ground states are found using a conjugate gradient technique to minimise the energy functional (also see [34, 39]).
3.2. Uniform ground states
In figures 3(a)–(d) we compare results obtained from the 3D and variational theories for the transverse density profile of a 164Dy condensate at a linear density of $n=2.5\times {10}^{3}/\mu {\rm{m}}$ for two values of as. The lower value of scattering length considered (${a}_{s}=95{a}_{0}$) is close to where the roton excitation softens to zero energy and becomes dynamically unstable (see section 3.4). For reference the harmonic oscillator ground state (i.e. quasi-1D result) is also shown. Here we observe that both the 3D and variational eGPE solutions have a much larger transverse width than the harmonic oscillator ground state since the system parameters are outside quasi-1D regime3We also note that the case in figures 3(b) and (d) has nas=12.6, well-satisfying the requirement ${{na}}_{s}\gg 1$ established in [28] for the quantum fluctuations to be described by the 3D formalism we use here. . We also note that while the confining potential is isotropic the condensates exhibits magnetostriction, i.e.significantly elongates in the y-direction compared to the x-direction.
Figure 3.
New window|Download| PPT slide Figure 3.Comparison of the variational (red lines) and 3D eGPE (black lines) solutions for a uniform infinite system. The 1/e density contours of the transverse modes of the 3D eGPE χ and the variational approach χσ for (a) ${a}_{s}=120{a}_{0}$ and (b) ${a}_{s}=95{a}_{0}$. The harmonic oscillator ground state ${\chi }_{\mathrm{ho}}$ is shown for reference (blue lines). In (c) and (d) we compare the transverse mode profiles along the x (dash–dot) and y (lines) axes for the cases given in (a) and (b), respectively. (e) The effective 1D k-space interaction kernel obtained from the various transverse functions for ${a}_{s}=120{a}_{0}$ (dashed lines) and ${a}_{s}=95{a}_{0}$ (solid lines). The 3D eGPE result ${\tilde{U}}_{z}$ is obtained by evaluating equation (13) using χ. The variational ${\tilde{U}}_{\sigma }$ and the harmonic oscillator ${\tilde{U}}_{\mathrm{ho}}$ results are obtained from equation (18). Results for 164Dy using ${a}_{{dd}}=130.8\,{a}_{0}$, with ${\omega }_{x,y}=2\pi \times 150\,\mathrm{Hz}$, ${\omega }_{z}=0$, and $n=2.5\times {10}^{3}\ \mu {{\rm{m}}}^{-1}$.
In figure 3(e) we compare the effective 1D k-space interaction kernel for the various theories. While the uniform ground state only depends on the value of the kernel at kz=0 (24), the excitations are sensitive to its non-zero kz behaviour (26). Our results show that our approximate kernel ${\tilde{U}}_{\sigma }$ closely matches that obtained from the full 3D eGPE solution. In comparison, the quasi-1D kernel based on the harmonic oscillator ground state (see section 2.2.5) is a poor approximation.
3.3. Fully trapped ground states
In figure 4 we present results for ground states with axial confinement. The line density profiles of the solutions reveal that the variational theory is in reasonable agreement with the 3D eGPE solution, although it generally tends to have a lower peak density. Except for small atom numbers N, for which the interaction effects are negligible, the quasi-1D case is in poor agreement with the other theories. A significant difference between the variational and the 3D result arises because the variational solution has the separable form ${\rm{\Psi }}({\boldsymbol{x}})\,={\psi }_{0}(z){\chi }_{\sigma }({\boldsymbol{\rho }})$, and thus has the same transverse profile for all z. At higher N we can see that this not a good approximation to the 3D solution: the transverse profile at z=0 (where the density is highest) is more strongly affected by interactions (larger average width and anisotropy) than it is for higher values of $| z| $ (see inset to figure 4(a)). For the case of condensates with contact interactions an effective 1D theory has been developed that allows the transverse profile χσ to vary slowly with z (see [27]). Such a theory could be developed for the dipolar case, although we do not pursue this here (also see [40]). In figures 4(b) and (c) we compare the ground state energy, observing that over the wide parameter regime considered the variational prediction for the energy is typically within a few percent of the full 3D solution.
Figure 4.
New window|Download| PPT slide Figure 4.Comparison of line density profiles $n(z)=\int {\rm{d}}{\boldsymbol{\rho }}| {\rm{\Psi }}({\boldsymbol{x}}){| }^{2}$ and energies for a 164Dy condensate in a trap with ${\omega }_{x,y}=2\pi \,\times 150\,\mathrm{Hz}$, ${\omega }_{z}=2\pi \times 20\,\mathrm{Hz}$, for ${a}_{s}=100{a}_{0}$ and various atom numbers N. (a) The line density along the z-axis calculated using the 3D eGPE (black), variational (red) and the quasi-1D (blue) theories. Inset shows the 1/e density contours (relative to the density at ${\boldsymbol{\rho }}={\bf{0}}$) in the ${\boldsymbol{\rho }}$-plane for the various theories for N=5×104. The contour of the 3D eGPE solution is evaluated at z=0 (black) and z=15 μm (grey). (b) Energy per particle of the three theories and (c) error in the energy of the variational and the quasi-1D theories relative to the 3D eGPE results.
3.4. Uniform system excitations: roton softening
In figures 5(a) and (b) we compare the predictions of the variational and 3D theories for the spectrum of a uniform case as as is varied. In these results we see that a roton (i.e. a local minimum in the excitation dispersion relation) appears for ${a}_{s}\approx 95{a}_{0}$ and lowers in energy as as is further decreased. Our calculations predict that the roton hits zero energy at the critical value of scattering length ${a}_{s}^{* }$ with ${a}_{s}^{* }=91.6{a}_{0}$ (${a}_{s}^{* }=93.3{a}_{0}$) according the variational (3D) theory for the density considered. In general we find that the variational theory predicts a lower value of ${a}_{s}^{* }$ than the 3D result (also see figures 5(c) and (e)). For ${a}_{s}\lt {a}_{s}^{* }$ the uniform state is dynamically unstable.
Figure 5.
New window|Download| PPT slide Figure 5.Excitation dispersion relations obtained from the (a) variational (dark blue) and (b) 3D (light blue) BdG theories for various values of as as labelled. In (a) the black circle indicates the roton coordinates (${k}_{\mathrm{rot}},{\epsilon }_{\mathrm{rot}}$) for ${a}_{s}=95{a}_{0}$. In (b) the higher excitations bands for ${a}_{s}=150{a}_{0}$ are shown (black dotted lines). The (c) roton energy and (d) roton wavevector as as changes for the variational (dark blue line) and 3D (light blue dots) theories. In (c) the critical value ${a}_{s}^{* }$ at which the roton energy goes to zero is indicated for each theory with an arrow and the fit function $\alpha \sqrt{{a}_{s}-{a}_{s}^{* }}$ (with α a fitting parameter) is also shown (dashed lines). In (d) the value of the critical roton wavevector (${k}_{\mathrm{rot}}^{* }$) when ${\epsilon }_{\mathrm{rot}}=0$ for ${a}_{s}={a}_{s}^{* }$ is indicated by an arrow for each theory. The roton critical values (e) ${a}_{s}^{* }$ and (f) ${k}_{\mathrm{rot}}^{* }$ as the system density changes. We compare the variational (lines) and 3D (symbols) theories both including (dark blue line for variational, light blue circles for 3D) and neglecting (red line for variational, magenta crosses for 3D) the quantum fluctuation term. Results for 164Dy with ${\omega }_{x,y}=2\pi \times 150\,\mathrm{Hz}$, and in (a)–(d) the density is n=2.5×103μm−1.
Identifying the local minimum in the dispersion relation with the roton energy εrot and wavevector krot (see figure 5(a)), we can monitor the behaviour of the roton as as varies in figures 5(c) and (d). The roton energy is seen to soften to zero as ${({a}_{s}-{a}_{s}^{* })}^{1/2}$ for as above but close to ${a}_{s}^{* }$ [4]. The roton wave vector tends to increase as as decreases, and obtains the value ${k}_{\mathrm{rot}}^{* }$ at the critical point ${a}_{s}^{* }$. We observe that the variational theory predicts ${k}_{\mathrm{rot}}^{* }$ to be larger than that obtained from the 3D theory, however the krot values of both theories are similar at the same as values (see figure 5(d)). We note that our results show that the roton wave vector occurs at a value slightly higher than the inverse harmonic oscillator length along the dipole direction, i.e.$\sqrt{m{\omega }_{y}/{\hslash }}=1.56/\mu $m, similar to the observations of experiments [4].
In figures 5(e) and (f) we examine the values of ${a}_{s}^{* }$ and ${k}_{\mathrm{rot}}^{* }$ for a range of system densities. We also include results without the quantum fluctuation term. For small n, the quantum fluctuations have a small effect and the theories make similar predictions4The results for $n\lesssim 500/\mu {\rm{m}}$ on figures 5(e) and (f) have ${{na}}_{s}\lesssim 1$, and the 3D treatment of quantum fluctuations is inappropriate. A quantitative treatment of this regime is outside the scope of this work. . However, in general we find that ${a}_{s}^{* }$ is larger when the quantum fluctuation term is neglected. We understand this arises because the quantum fluctuations effectively act as a repulsive interaction and tend to stabilise (i.e. lift the energy) of the roton. Thus with quantum fluctuations a lower as value is needed to destabilise the condensate. We also find that ${a}_{s}^{* }$ as a function of n has a maximum when quantum fluctuations are included. In contrast without quantum fluctuations ${a}_{s}^{* }$ monotonically increases with n, slowly approaching the value add in the large density limit.
4. Conclusions and outlook
In this paper we have reported the development of a simple theory for a dipolar condensate in an elongated confining potential. Our main result is the effective 1D variational eGPE for stationary states, and the associated BdG theory of its collective excitations. This theory is practical to solve with modest computational resources, yet provides a good quantitative description of the full 3D solution.
In the application of our theory we have focused on the typical density, interaction and trap parameter regimes used in current experiments. For example, the rotons observed in the experiments of Chomaz et al [4] occurred (prior to structure formation occurring) in an elongated system with linear densities in the range $2\times {10}^{3}-4\times {10}^{3}\,\mu {\rm{m}}$−1. This regime is well-beyond where the quasi-1D approximation is valid. Using our theory we have predicted the scattering length ${a}_{s}^{* }$ and roton wave vector ${k}_{\mathrm{rot}}^{* }$ at the point where roton softens to zero energy as a function of system density. Our results show that quantum fluctuations lower the value of ${a}_{s}^{* }$ and cause it to have a non-monotonic dependence on density. These predictions could be investigated in future experiments and may relate to the non-monotonic dependence of the value of as where the supersolid transition was observed (see figure 1(g) of [41]).
Here we have restricted our focus to the regime where the system does not develop density modulations, which tends to occur at lower values of as (e.g. after the roton softens and causes a dynamic instability). Such modulations can indicate the onset of a supersolid state, as has been observed in three recent experiments working with dipolar condensates in elongated trapping potentials. So far theoretical studies of the ground states and their excitations of these modulated states required large scale numerical methods for cases with (see [22, 23, 41–44]) and without (see [15]) axial confinement. Our theory can treat such modulated states, but a full and systematic treatment of this is beyond the current scope and will be examined in future work.
Acknowledgments
We acknowledge the contribution of NZ eScience Infrastructure (NeSI) high-performance computing facilities, support from the Marsden Fund of the Royal Society of New Zealand, and valuable discussions with F Ferlaino and L Chomaz.
PasquiouBBismutGMaréchalEPedriPVernacLGorceixOLaburthe-TolraB2011 Spin relaxation and band excitation of a dipolar Bose–Einstein condensate in 2D optical lattices Phys. Rev. Lett. 106 015301 DOI:10.1103/PhysRevLett.94.160401 [Cited within: 1]
LuMBurdickN QYounS HLevB L2011 Strongly dipolar Bose–Einstein condensate of dysprosium Phys. Rev. Lett. 107 190401 DOI:10.1103/PhysRevLett.107.190401
LuMBurdickN QLevB L2012 Quantum degenerate dipolar Fermi gas Phys. Rev. Lett. 108 215301 DOI:10.1103/PhysRevLett.107.190401
ChomazLvan BijnenR M WPetterDFaraoniGBaierSBecherJ HMarkM JWächtlerFSantosLFerlainoF2018 Observation of roton mode population in a dipolar quantum gas Nat. Phys. 14 442 DOI:10.1038/s41567-018-0054-7 [Cited within: 4]
PetterDNataleGvan BijnenR M WPatscheiderAMarkM JChomazLFerlainoF2019 Probing the roton excitation spectrum of a stable dipolar Bose gas Phys. Rev. Lett. 122 183401 DOI:10.1103/PhysRevLett.122.183401 [Cited within: 1]
SantosLShlyapnikovG VLewensteinM2003 Roton-maxon spectrum and stability of trapped dipolar Bose–Einstein condensates Phys. Rev. Lett. 90 250403 DOI:10.1103/PhysRevLett.90.250403 [Cited within: 1]
RonenSBortolottiD C EBohnJ L2007 Radial and angular rotons in trapped dipolar gases Phys. Rev. Lett. 98 030406 DOI:10.1103/PhysRevLett.98.030406
BlakieP BBaillieDBissetR N2012 Roton spectroscopy in a harmonically trapped dipolar Bose–Einstein condensate Phys. Rev. A 86 021604 DOI:10.1103/PhysRevA.86.021604
CorsonJ PWilsonR MBohnJ L2013 Stability spectroscopy of rotons in a dipolar Bose gas Phys. Rev. A 87 051605 DOI:10.1103/PhysRevA.87.051605
CorsonJ PWilsonR MBohnJ L2013 Geometric stability spectra of dipolar Bose gases in tunable optical lattices Phys. Rev. A 88 013614 DOI:10.1103/PhysRevA.88.013614
Jona-LasinioMŁakomyKSantosL2013 Time-of-flight roton spectroscopy in dipolar Bose–Einstein condensates Phys. Rev. A 88 025603 DOI:10.1103/PhysRevA.88.025603
BissetR NBlakieP B2013 Fingerprinting rotons in a dipolar condensate: super-poissonian peak in the atom-number fluctuations Phys. Rev. Lett. 110 265302 DOI:10.1103/PhysRevLett.110.265302
BissetR NWilsonR MBaillieDBlakieP B2016 Ground-state phase diagram of a dipolar condensate with quantum fluctuations Phys. Rev. A 94 033619 DOI:10.1103/PhysRevA.94.033619 [Cited within: 2]
ChomazLBaierSPetterDMarkM JWächtlerFSantosLFerlainoF2016 Quantum-fluctuation-driven crossover from a dilute Bose–Einstein condensate to a macrodroplet in a dipolar quantum fluid Phys. Rev. X 6 041039 DOI:10.1103/PhysRevX.6.041039
SchmittMWenzelMBöttcherFFerrier-BarbutIPfauT2016 Self-bound droplets of a dilute magnetic quantum liquid Nature 539 259 DOI:10.1038/nature20126
BöttcherFSchmidtJ-NWenzelMHertkornJGuoMLangenTPfauT2019 Transient supersolid properties in an array of dipolar quantum droplets Phys. Rev. X 9 011051 DOI:10.1103/PhysRevX.9.011051 [Cited within: 1]
TanziLLucioniEFamàFCataniJFiorettiAGabbaniniCBissetR NSantosLModugnoG2019 Observation of a dipolar quantum gas with metastable supersolid properties Phys. Rev. Lett. 122 130405 DOI:10.1103/PhysRevLett.122.130405 [Cited within: 2]
SalasnichLParolaAReattoL2002 Effective wave equations for the dynamics of cigar-shaped and disk-shaped Bose condensates Phys. Rev. A 65 043614 DOI:10.1103/PhysRevA.65.043614 [Cited within: 2]
EdlerDMishraCWächtlerFNathRSinhaSSantosL2017 Quantum fluctuations in quasi-one-dimensional dipolar Bose–Einstein condensates Phys. Rev. Lett. 119 050403 DOI:10.1103/PhysRevLett.119.050403 [Cited within: 2]
DeuretzbacherFCremonJ CReimannS M2010 Ground-state properties of few dipolar bosons in a quasi-one-dimensional harmonic trap Phys. Rev. A 81 063616 DOI:10.1103/PhysRevA.81.063616
DeuretzbacherFCremonJ CReimannS M2013 Erratum: ground-state properties of few dipolar bosons in a quasi-one-dimensional harmonic trap Phys. Rev. A 87 039903 DOI:10.1103/PhysRevA.87.039903
GiovanazziSO'DellD H J2004 Instabilities and the roton spectrum of a quasi-1D Bose–Einstein condensed gas with dipole–dipole interactions Eur. Phys. J. D 31 439445 DOI:10.1140/epjd/e2004-00146-7 [Cited within: 1]
BaoWCaiYWangH2010 Efficient numerical methods for computing ground states and dynamics of dipolar Bose–Einstein condensates J. Comput. Phys. 229 7874 DOI:10.1016/j.jcp.2010.07.001 [Cited within: 1]
AntoineXLevittATangQ2017 Efficient spectral computation of the stationary states of rotating Bose–Einstein condensates by preconditioned nonlinear conjugate gradient methods J. Comput. Phys. 343 92 DOI:10.1016/j.jcp.2017.04.040 [Cited within: 1]
KnightM JBlandTParkerN GMartinA M2019 Improved low-dimensional wave equations for cigar-shaped and disk-shaped dipolar Bose–Einstein condensates arXiv:1908.02395 [cond-mat.quant-gas] [Cited within: 1]
TanziLRoccuzzoS MLucioniEFamàFFiorettiAGabbaniniCModugnoGRecatiAStringariS2019 Supersolid symmetry breaking from compressional oscillations in a dipolar quantum gas Nature 574 382 DOI:10.1038/s41586-019-1568-6
GuoMBöttcherFHertkornJSchmidtJ-NWenzelMBüchlerH PLangenTPfauT2019 The low-energy goldstone mode in a trapped dipolar supersolid Nature 574 386 DOI:10.1038/s41586-019-1569-5
NataleGvan BijnenR M WPatscheiderAPetterDMarkM JChomazLFerlainoF2019 Excitation spectrum of a trapped dipolar supersolid and its experimental evidence Phys. Rev. Lett. 123 050402 DOI:10.1103/PhysRevLett.123.050402 [Cited within: 1]
,1, D Baillie, Sukla PalDodd-Walls Centre for Photonic and Quantum Technologies, New Zealand and Department of Physics, 
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
