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. 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:3 ), given by
$ \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 ( $ \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 integral2 , 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 behaviour18 ) by inspection and numerical experiment: it satisfies the required limiting behaviour and reduces to the exact isotropic result at η = 1.
$ \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 (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. (a) Comparison of the (dotted lines) analytic result ( |
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:12 ) but using ${\tilde{U}}_{\sigma }$ in place of ${\tilde{U}}_{z}$.
$ \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 (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.
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 regime3 . 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. 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 ( |
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. 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. 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 predictions4 . 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.