Improved low-dimensional wave equations for cigar-shaped and disk-shaped dipolar Bose-Einstein condensates

TL;DR

This paper derives improved low-dimensional equations for dipolar Bose-Einstein condensates with arbitrary polarization, accurately predicting ground states and outperforming standard reductions.

Contribution

It introduces a variational ansatz-based method for deriving effective 1D and 2D equations that better approximate 3D results for dipolar BECs.

Findings

Accurately predicts ground state densities.

Strong agreement with 3D results even outside strict regimes.

Significant improvement over standard reductions.

Abstract

Within the formalism of the Gross-Pitaevskii equation, we derive effective one- and two-dimensional equations for cigar- and pancake-shaped dipolar Bose-Einstein condensates with arbitrary polarization angle. These are based on an ansatz for the condensate wavefunction whose width in the tightly-confined direction/s is treated variationally. The equations constitute a coupled partial differential for the low-dimensional wavefunction and algebraic equations for the width parameters. This approach accurately predicts the ground state densities of cigar-shaped and pancake-shaped dipolar Bose-Einstein condensates, and gives strong agreement with the three-dimensional results, even as the trapping is relaxed away from the strict quasi-one- and quasi-two-dimensional regimes. This approach offers a significant improvement over the standard one- and two-dimensional reduction.

Equations

$$i\hbar\partial_{t}\psi(\bm{x},t)=\Bigg{\{}-\frac{\hbar^{2}}{2m}\nabla^{2}+\mathcal{V}_{\text{trap}}(\bm{x})+gN|\psi(\bm{x},t)|^{2}\\ +\Phi_{\text{dd}}(\bm{x})\Bigg{\}}\psi(\bm{x},t),$$

$${\Phi }_{\text{dd}}(\mathbf{x})=N\int \text{d}{\mathbf{x}}^{\prime }{\mathcal{V}}_{\text{dd}}(\mathbf{x}-{\mathbf{x}}^{\prime })|\psi ({\mathbf{x}}^{\prime },t){|}^2,$$

$${\mathcal{V}}_{\text{dd}}(\mathbf{x})=\frac{C_{\text{dd}}}{4\pi }\frac{1-3{\cos }^2\theta }{|\mathbf{x}{|}^3}.$$

$$\mathcal{S}=\int\text{d}t\int\text{d}\bm{x}\,\psi^{*}(\bm{x},t)\Bigg{\{}i\hbar\partial_{t}+\frac{\hbar^{2}}{2m}\nabla^{2}-\mathcal{V}_{\text{trap}}(\bm{x})\\ -\frac{1}{2}gN|\psi(\bm{x},t)|^{2}-\frac{1}{2}\Phi_{\text{dd}}(\bm{x})\Bigg{\}}\psi(\bm{x},t),$$

$$\Phi_{\text{dd}}(\bm{x})=\mathcal{F}^{-1}\Big{\{}\tilde{\mathcal{V}}_{\text{dd}}(\bm{k})\,\mathcal{F}[|\psi(\bm{x},t)|^{2}]\Big{\}},$$

$$\frac{\tilde{\mathcal{V}}(\mathbf{k})}{C_{\text{dd}}}=\frac{k_x^2{\sin }^2\alpha +k_x{k}_z\sin (2\alpha )+k_z^2{\cos }^2\alpha }{k_x^2+k_y^2+k_z^2}-\frac{1}{3},$$

$${\mathcal{V}}_{\text{trap}}(\mathbf{x})=\frac{1}{2}m{\omega }_{\perp }^2(x^2+y^2)+\frac{1}{2}m{\omega }_z^2{z}^2,$$

$$U_{\text{dd}}^{{\rho }_c,Z_c}(\mathbf{r})=\{\begin{array}{ll}\frac{C_{\text{dd}}}{4\pi }\frac{1-3{\cos }^2\theta }{r^3},& \rho <{\rho }_c\text{ and }|z|<Z_c,\\ 0,& \text{otherwise}.\end{array}$$

$${\tilde{U}}_{\text{dd}}^{{\rho }_c,Z_c}(\mathbf{k})=\frac{C_{\text{dd}}}{3}\left(3{\cos }^2{\theta }_k-1\right)+C_{\text{dd}}{e}^{-Z_c{k}_{\rho }}F({\theta }_k,z_k)$$

$$-C_{\text{dd}}{\int }_{{\rho }_c}^{\infty }\text{d}\rho \rho {\int }_0^{Z_c}\text{d}z\cos (k_{z}z)\frac{{\rho }^2-2z^2}{{\left({\rho }^2+z^2\right)}^{5/2}}{J}_0(k_{\rho }\rho ),$$

$$\psi (\mathbf{x},t)$$

$$\phi (x,y,t;\sigma (z,t))=\frac{\exp \left[\frac{-(x^2+y^2)}{2\sigma (z,t{)}^2}\right]}{\sqrt{\pi }\sigma (z,t)}$$

$$\mathcal{S}_{1\text{D}}=\int\text{d}z\int\text{d}t\,f^{*}\Bigg{\{}i\hbar\partial_{t}+\frac{\hbar^{2}}{2m}\partial_{z}^{2}-\frac{1}{2}m\omega_{z}^{2}z^{2}-\frac{\hbar^{2}}{2m\sigma^{2}}\\ -\frac{1}{2}m\omega_{\perp}^{2}\sigma^{2}-\frac{gN}{4\pi}\frac{|f|^{2}}{\sigma^{2}}-\frac{1}{2}\Phi_{\text{dd}}^{1\text{D}}\Bigg{\}}f.$$

$$\Phi_{\text{dd}}^{1\text{D}}\Bigg{|}_{\alpha=0}=\frac{C_{\text{dd}}N}{6\pi\sigma^{2}}\mathcal{F}_{z}^{-1}\Big{\{}(3qe^{q}\Gamma[0,q]-1)\mathcal{F}_{z}\{|f|^{2}\}\Big{\}},$$

$$i\hslash {\partial }_{t}f$$

$$\displaystyle\hskip 80.29999pt+\frac{gN|f|^{2}}{2\pi\sigma^{2}}+\Phi_{\text{dd}}^{1\text{D}}\Big{|}_{\alpha=0}\Bigg{\}}f,$$

$$0$$

$$\tilde{t}=t{\omega }_z,\tilde{z}=\frac{z}{{\ell }_z},{\tilde{a}}_s=\frac{a}{{\ell }_z},\tilde{\sigma }=\frac{\sigma }{{\ell }_z},{\tilde{k}}_z=k_z{\ell }_z.$$

$$i{\partial }_{\tilde{t}}\tilde{f}$$

$$0$$

$$\psi (\mathbf{x},t)=\frac{f(z,t)}{\sqrt{\pi {\sigma }_x{\sigma }_y}}\exp \left\{-\frac{x^2}{2{\sigma }_x^2}-\frac{y^2}{2{\sigma }_y^2}\right\},$$

$$i{\partial }_{\tilde{t}}\tilde{f}$$

$$+\frac{3\sqrt{2\pi }{\tilde{\beta }}_{1\text{D}}{\widetilde{ϵ}}_{\text{dd}}}{{\tilde{\sigma }}_x}{\mathcal{F}}_{\tilde{z}}^{-1}\left[\int \text{d}{\tilde{k}}_y{e}^{-2{\pi }^2{\tilde{k}}_y^2{\tilde{\sigma }}_y^2}\left(\tilde{A}(q){\sin }^2\alpha +\tilde{B}(q){\cos }^2\alpha -\frac{1}{3}\right){\mathcal{F}}_{\tilde{z}}\{|\tilde{f}{|}^2\}\right]\}\tilde{f},$$

$$0$$

$$-3\sqrt{2\pi }{\tilde{\beta }}_{1\text{D}}{ϵ}_{\text{dd}}{\tilde{\sigma }}_x{\mathcal{F}}_{\tilde{z}}^{-1}\left[\int \text{d}{\tilde{k}}_y{e}^{-2{\pi }^2{\tilde{k}}_y^2{\tilde{\sigma }}_y^2}\left(\tilde{C}(q){\sin }^2\alpha +\tilde{D}(q){\cos }^2\alpha +\frac{1}{3}\right){\mathcal{F}}_{\tilde{z}}\{|\tilde{f}{|}^2\}\right],$$

$$0$$

$$+\frac{12\sqrt{2}{\pi }^{5/2}{\tilde{\sigma }}_y^4{\tilde{\beta }}_{1\text{D}}{ϵ}_{\text{dd}}}{{\tilde{\sigma }}_x}{\mathcal{F}}_{\tilde{z}}^{-1}\left[\int \text{d}{\tilde{k}}_y{e}^{-2{\pi }^2{\tilde{k}}_y^2{\tilde{\sigma }}_y^2}\left(\tilde{A}(q){\sin }^2\alpha +\tilde{B}(q){\cos }^2\alpha -\frac{1}{3}\right){\mathcal{F}}_{\tilde{z}}\{|\tilde{f}{|}^2\}\right],$$

$$\tilde{A}(q)$$

$$\tilde{B}(q)$$

$$\tilde{C}(q)$$

$$\tilde{D}(q)$$

$$\Gamma =\frac{\max \{{\tilde{\sigma }}_y\}}{\max \{{\tilde{\sigma }}_x\}}.$$

$$\psi (\mathbf{x},t)$$

$$=f(x,y,t)\frac{\exp \left\{-\frac{z^2}{2\eta (x,y,t{)}^2}\right\}}{{\pi }^{1/4}\sqrt{\eta (x,y,t)}},$$

$$\begin{gathered} {\mathcal{S}}_{2\text{D}}=\int \text{d}x\text{d}y\text{d}t{f}^{\ast }\{i\hslash {\partial }_t+\frac{{\hslash }^2}{2m}{\nabla }_{\perp }^2-\frac{{\hslash }^2}{4m{\eta }^2}\frac{1}{4}m{\omega }_z^2{\eta }^2 \\ -\frac{1}{2}m{\omega }_{\perp }^2(x^2+y^2)-\frac{gN|f{|}^2}{2\sqrt{2\pi }\eta }-\frac{1}{2}{\Phi }_{\text{dd}}^{2\text{D}}\}f, \end{gathered}$$

$$\Phi_{\text{dd}}^{2\text{D}}=\frac{C_{\text{dd}}N}{\sqrt{2\pi}\eta}\mathcal{F}_{x,y}^{-1}\Bigg{\{}\Big{[}\sin^{2}\alpha\,G(\xi)+\cos^{2}\alpha\,F(\xi)\\ -\frac{1}{3}\Big{]}\mathcal{F}_{x,y}\left\{|f|^{2}\right\}\Bigg{\}}f,$$

$$G(\xi )$$

Full text

I Introduction

Gaseous Bose-Einstein condensates (BECs) are today routinely created by experimental groups over the world. In the decade following the creation of the first BEC in 1995 Anderson et al. (1995), typical BECs were made from alkali metals with small magnetic dipole moments, such as lithium-7 Bradley et al. (1997), sodium-23 Davis et al. (1995), and potassium-41 Modugno et al. (2001). In such BECs, the atom-atom interactions are short-range and isotropic, arising from van der Waals forces, and are well accounted for theoretically as a contact pseudo-potential Pitaevskii and S. (2003). In more recent times, BECs have also been produced using atoms with significant magnetic dipole moments. This began with the condensation of chromium-52 in 2005 Griesmaier et al. (2005), followed by dysprosium-164 Lu et al. (2011), erbium-168 Aikawa et al. (2012), and erbium-166 Chomaz et al. (2016). In these so-called dipolar BECs, the magnetic dipoles are polarized in a common direction by an externally-applied magnetic field. Alongside this experimental effort there has been a renewed interest in the theoretical understanding of dipolar BECs, as well as dipolar quantum gases more generally Lahaye et al. (2009); Baranov et al. (2012).

In contrast to the van der Waals/contact interactions, the dipole-dipole interaction (DDI) between the atoms is anisotropic (part attractive and part repulsive depending on the relative orientation of the dipoles) and long-range. These characteristics underpin a range of striking physical phenomena, such as magnetostriction Yi and You (2001); Santos et al. (2000); Góral et al. (2000); Fischer (2006); Lahaye et al. (2007), the existence of a roton dip in the excitation spectrum O’Dell et al. (2003); Santos et al. (2003); Chomaz et al. (2018), and roton-driven density corrugations Dutta and Meystre (2007); Wilson et al. (2008); Lu et al. (2010a) and recently-reported supersolid phases Tanzi et al. (2019); Böttcher et al. (2019); Chomaz et al. (2019).

There is strong interest in dipolar BECs in elongated and flattened geometries, for a variety of reasons. Firstly, such geometries can be exploited to enhance the role of the DDIs, as shown in Fig. 1: in an elongated dipolar BEC with dipoles polarised along the primary axis, the attractive head-to-tail alignment is favoured, whereas in a flattened geometry, the repulsive side-by-side configuration is favoured. Secondly, strong trapping in one or two dimensions is required to support the presence of the roton and its rich manifestations mentioned above. Thirdly, such geometries allow access to dipolar physics in lower dimensions, including one-dimensional bright solitons Cuevas et al. (2009); Baizakov et al. (2015); Edmonds et al. (2017), two-dimensional bright solitons Pedri and Santos (2005); Jiang et al. (2016); Huang et al. (2017); Tikhonenkov et al. (2008); Raghunandan et al. (2015); Chen et al. (2016), one-dimensional dark solitons Pawłowski and Rzążewski (2015); Bland et al. (2015, 2017); Edmonds et al. (2016); Umarov et al. (2016), vortices and vortex lattices Cooper et al. (2005); Zhang and Zhai (2005); Yi and Pu (2006); Wilson et al. (2009); Abad et al. (2009); Malet et al. (2011); Mulkerin et al. (2013); Gautam (2014); Martin et al. (2017); Kumar et al. (2017); Jia et al. (2018); Cai et al. (2018); Prasad et al. (2019), and near-integrability Tang et al. (2018).

The workhorse model for BECs in general is the Gross-Pitaevskii equation (GPE), a zero-temperature mean-field model for the macroscopic wavefunction of the BEC Pitaevskii (1961); Gross (1961); Pitaevskii and S. (2003); Pethick and Smith (2008); Barenghi and Parker (2016). In the absence of DDIs, it has the form of a cubic nonlinear Schrödinger equation (NLSE), while it can be extended to dipolar BECs through the inclusion of an additional, non-local term Lahaye et al. (2009). Theoretical studies of low dimensional BECs typically employ the reduced 1D or 2D forms of the GPE, in which the condensate wavefunction is written as the product of the low dimensional wavefunction and a static profile (the non-interacting harmonic oscillator state) in the tightly confined direction(s). By integrating out the latter, one arrives at the effective 1D/2D GPE. These equations are strictly valid in the so-called quasi-1D and quasi-2D regimes. Assuming the condensate to by cylindrically symmetric, with radius $R_{\perp}$ and axial half-length $R_{z}$, then the quasi-1D regime requires $R_{z}>\xi>R_{\perp}$, where $\xi$ is the healing length of the condensate. Similarly, the quasi-2D regime requires $R_{\perp}>\xi>R_{z}$. The form and validity of these equations have been established for conventional BECs Olshanii (1998); Petrov et al. (2000a, b); Menotti and Stringari (2002); Muñoz Mateo and Delgado (2006); Pitaevskii and S. (2003); Barenghi and Parker (2016), and later extended to dipolar BECs Parker and O’Dell (2008); Cai et al. (2010).

The widespread use of the reduced 1D and 2D GPEs is motivated by their simplicity and ease of solution. However, their predictive power in capturing realistic 3D behaviour has its limits: the strict quasi-1D and quasi-2D regimes in which they are valid can be challenging to reach experimentally and they are not quantitatively accurate in more relaxed geometries. Moreover, 1D and 2D models can fail to capture important physical effects, with the notable example being the collapse instability of attractively-interacting BECs, which is not captured by the 1D GPE Pérez-García et al. (1997); Salasnich et al. (2002a); Parker et al. (2007).

For conventional (non-dipolar) BECs, low dimensional equations have been pursued which capture some of the residual three-dimensional physics and overcome the above challenges Salasnich et al. (2002b). These so-called non-polynomial Schrödinger equations take the profile of the wavefunction in the tightly-confined direction(s) as a Gaussian ansatz with variable width, rather than fixed width. They provide excellent agreement with the 3D solutions in cigar-shaped and pancake-shaped condensates, even if the tight confinement is relaxed away from the strict quasi-1D and quasi-2D regime. Moreover, the 1D non-polynomial GPE successfully captured the 3D physics of the collapse instability Salasnich et al. (2002a).

In the dipolar case, the strict 1D/2D GPEs have been derived and shown to capture the 3D behaviour when the system is deep in the 1D/2D regime, but depart away from this regime Cai et al. (2010); Baillie and Blakie (2015). To date, the above approach has not been extended to the dipolar BECs. One can expect that such equations would allow for accurate, convenient prediction of elongated and flattened dipolar BECs over a greater parameter regime, and may be particularly important in capturing the three-dimensional physics of the collapse and rotons.

Motivated by this, in this work we use the variable-width approach to derive effective lower-dimensional equations for cigar-shaped and pancake-shaped dipolar BECs, in the manner of Salasnich et al. Cai et al. (2010). Because both the lower-dimensional wavefunction and the width are variational, we minimise the energy of the system with respect to both, and arrive at a wave equation for the wavefunction, and a polynomial equation for the width/s. Solving these coupled equations in tandem allows for an accurate prediction of the system. We will show that these equations accurately predict the ground state solutions, and provide a significant improvement over the standard 1D/2D effective equations. In Section II we describe the over-arching theoretical model of the 3D dipolar GPE, while in Sections III and IV we proceed to derive and test the variable-width 1D and 2D equations, comparing to solutions of the full 3D dipolar GPE. In Section V we conclude our work.

II The 3D Model

Consider a dipolar BEC of $N$ atoms at zero temperature, confined in an external trapping potential $\mathcal{V}_{\text{trap}}(\bm{x})$. The time-evolution of the macroscopic wavefunction, $\psi(\bm{x},t)$, is given by the 3D dipolar GPE (dGPE) Lahaye et al. (2009)

[TABLE]

where $\nabla^{2}$ is the Laplacian, $g=4\pi\hbar^{2}a_{s}/m$ is the $s$-wave coupling constant with $a_{s}$ being the $s$-wave scattering length and $m$ being the mass of the constituent bosons. $|\psi(\bm{x},t)|^{2}$ is the probability density of the atoms, satisfying the normalisation condition $\int\text{d}\bm{x}\,|\psi(\bm{x},t)|^{2}=1$. $\Phi_{\text{dd}}$ is the mean-field potential generated by the DDIs, which is given by the convolution

[TABLE]

where $\mathcal{V}_{\text{dd}}(\bm{x})$ is the DDI interaction term

[TABLE]

Here $\theta$ is the angle between the polarisation direction of the dipoles and the inter-particle vector between two interacting dipoles (see Fig. 2) and $C_{\text{dd}}=\mu_{0}\mu_{\text{d}}^{2}$, where $\mu_{0}$ is the magnetic vacuum permeability and $\mu_{\text{d}}$ is the magnetic dipole moment of the bosons. We introduce the dimensionless parameter $\epsilon_{\text{dd}}=C_{\text{dd}}/(3g)$ which parametrises the strength of the DDIs relative to the contact interactions.

The 3D dGPE can be recast into a more useful form as the Euler-Lagrange equation of the action functional of the system. Such an approach results in a clear method of dimension reduction, as shown by Salasnich et al. Salasnich et al. (2002b). The action functional of the system is

[TABLE]

where $\psi^{*}(\bm{x},t)$ is the complex conjugate of $\psi(\bm{x},t)$. The 3D dGPE in Eq. (1) is the Euler-Lagrange equation of Eq. (4), taken with respect to $\psi^{*}(\bm{x},t)$. We recast the mean-field dipolar potential through the convolution theorem

[TABLE]

where $\mathcal{F}$ ($\mathcal{F}^{-1}$) denotes the Fourier (inverse Fourier) transform, and $\tilde{\mathcal{V}}_{\text{dd}}(\bm{k})$ is the Fourier transform of the DDI, given by Lahaye et al. (2009)

[TABLE]

where $\alpha$ is the angle of polarisation of the dipoles measured from the positive $z$-axis in the $xz$-plane (see Fig. 2).

The dipolar BEC is confined in a cylindrically symmetric harmonic oscillator trap with trap frequency $\omega_{z}$ in the $z$ direction, and $\omega_{\perp}$ in the $x$- and $y$-directions,

[TABLE]

For convenience we define the trap ratio $\gamma=\omega_{\perp}/\omega_{z}$: a cigar-shaped trap corresponds to $\gamma\gg 1$ and a pancake-shaped trap to $0<\gamma\ll 1$. We also introduce the harmonic oscillator lengths, $\ell_{z}=\sqrt{\hbar/m\omega_{z}}$ and $\ell_{\perp}=\sqrt{\hbar/m\omega_{\perp}}$.

The accuracy of our low-dimensional equations is found by comparing their ground state solution to the 3D result. The 3D ground state solutions are found via imaginary-time propagation of the 3D dGPE using the split-step Fourier method Fleck et al. (1976, 1978). We employ a 2563 grid with discretization $(\Delta x,\Delta y,\Delta z)=(0.04,0.04,0.0008)\ell_{\perp}$ for pancake geometries, and $(\Delta x,\Delta y,\Delta z)=(0.008,0.008,0.04)\ell_{\perp}$ in cigar geometries. Since fast Fourier transform algorithms are naturally periodic, we introduce a cylindrical cut-off to the dipolar potential, which restricts the range of the DDI to radius $\rho=\sqrt{x^{2}+y^{2}}<\rho_{c}$ and length $|z|<z_{c}$ in real-space in order to reduce the effects of alias copies, giving a real-space DDI

[TABLE]

The Fourier transform of this is semi-analytic Lu et al. (2010b), given by

[TABLE]

where $J_{0}$ is the zeroth-order Bessel function of the first kind and $F(\theta_{k},z_{k})=\sin^{2}\theta_{k}\cos\left(Z_{c}k_{z}\right)$ $-\sin\theta_{k}\cos\theta_{k}\sin\left(Z_{c}k_{z}\right)$. The final line requires numerical integration, which is computationally costly but only required once as part of the simulation setup.

III Variable-Width One-Dimensional Equations

We assume a prolate harmonic trap $\gamma\gg 1$ giving rise to a cigar-shaped dipolar BEC. Following Salasnich et al. Salasnich et al. (2002b), we write the 3D wavefunction in the separable form,

[TABLE]

where $f$ is the 1D wavefunction along $z$ and $\varphi$ is taken to have a Gaussian profile in $x$ and $y$, whose width $\sigma$ varies in time and with $z$,

[TABLE]

Note that both $f$ and $\varphi$ are normalised to unity. The strict 1D limit derived previously Cai et al. (2010) corresponds to the fixed width $\sigma=\ell_{\perp}$. However, here we consider variable width, treating $\sigma$ as a variational parameter.

To derive the variable-width 1D equations, we substitute the above form of the wavefunction into the action functional of Eq. (4), along with the DDI potential (5) and trapping potential (7), and integrate out the $x$ and $y$ dimensions. We also assume that $\varphi$ is slowly varying along the $z$-direction, allowing us to write $\nabla^{2}\varphi\approx\nabla_{\perp}^{2}\varphi$ where $\nabla_{\perp}^{2}=\partial_{x}^{2}+\partial_{y}^{2}$. We then arrive at the 1D action functional,

[TABLE]

Here, $\Phi_{\text{dd}}^{1\text{D}}$ is the 1D DDI potential which is derived from Eq. (6) by evaluating the Fourier (inverse Fourier) transform in the $x$ and $y$ ($k_{x}$ and $k_{y}$) directions, the result of which depends upon the polarisation direction $\alpha$. $\Phi_{\text{dd}}^{1\text{D}}$ is only analytic for $\alpha=0$; for $\alpha\neq 0$, the evaluation of the Fourier and inverse Fourier transforms have to be completed numerically.

III.1 $z$-Polarisation Case

If $\alpha=0$, the dipoles are aligned along the positive $z$-axis and share an axis with the primary axis of the cigar-shaped BEC. The 1D DDI potential then becomes,

[TABLE]

where $q=2\pi^{2}k_{z}^{2}\sigma^{2}$, $\mathcal{F}_{z}$ ($\mathcal{F}_{z}^{-1}$) denotes the 1D Fourier (inverse Fourier) transform and $\Gamma[0,\cdot]$ is the incomplete gamma function. With this DDI potential in the 1D action functional in Eq. (12), we apply the Euler-Lagrange equations to the conjugate of the wavefunction, $f^{*}$, and the variational width, $\sigma$. This leads to the following system of equations, an integro-differential equation for $f(z)$ coupled to an equation for $\sigma$,

[TABLE]

These are the variable-width GPE equations for the cigar-shaped dipolar BEC with $\alpha=0$.

For convenience, we scale our variables into dimensionless form using harmonic oscillator units based on the trapping along $z$: the timescale is $\omega_{z}^{-1}$ and the lengthscale is $\ell_{z}$. The dimensionless variables (denoted by tilde) are then,

[TABLE]

It follows that the dimensionless wavefunction is $\tilde{f}=f\sqrt{\ell_{z}}$ and the harmonic oscillator length in the $xy$ plane is $\tilde{\ell}_{\perp}=1/\sqrt{\gamma}$. With these transformations, the variable-width 1D equations for $\alpha=0$ become,

[TABLE]

where $\tilde{\beta}_{1\text{D}}=2\tilde{a}_{s}N$ and $\mathcal{F}_{\tilde{z}}$ ($\mathcal{F}_{\tilde{z}}^{-1}$) is the scaled Fourier (inverse Fourier) transform. These equations have a similar form to the 3D dGPE given in Eq. (1). The first and second terms in Eq. (17) are the kinetic energy, and the third and fourth terms are the potential energy. The final two terms are due to the interactions: the first of these is a contact interaction term combining the van der Waals interactions and the short-range component of the DDIs, while the second is the long-range interactions of the DDIs.

For $\epsilon_{\text{dd}}=0$, Eq. (18) can be solved analytically for $\tilde{\sigma}$, and substitution into Eq. (17) results in the NPSE of Salasnich et al. Salasnich et al. (2002b). Such a simplification cannot be implemented for $\epsilon_{\text{dd}}\neq 0$, and therefore Eq. (17) and Eq. (18) must be solved numerically in a self-consistent manner.

In the top panels of Fig. 3 we plot the ground state density $|f(z)|^{2}$ as predicted by the variable-width 1D model [Eqs. (17) and (18)], the fixed-width model [Eq. (17) with $\tilde{\sigma}=\ell_{\perp}$], and the 3D dGPE, for a highly-elongated system ($\gamma=80$) and a moderately-elongated system ($\gamma=10$). In both cases, there is close agreement between the 3D dGPE and the variable-width approach, while the fixed-width approach has a noticeable deviation. The lower panels of Fig. 3 show the variational width $\sigma(z)$ according to Eq. (18). For both cases, the departure of $\sigma$ from the harmonic oscillator length $\ell_{\perp}=\ell_{z}/\sqrt{\gamma}$ is significant, varying by up to 10-20%.

III.2 Arbitrary Polarisation

If the dipole moments are polarised away from the $z$-axis ($\alpha\neq 0)$, the 1D DDI potential is not be analytic and must be computed numerically. Furthermore, the anisotropic nature of the DDIs will cause the cylindrical symmetry of the dipolar BEC to be broken. To accommodate this, we extend the ansatz for the wavefunction to include an anisotropic Gaussian profile in the $xy$ plane

[TABLE]

where $\sigma_{x}(z,t)$ and $\sigma_{y}(z,t)$ are the variational widths. Using this wavefunction, we evaluate the 3D DDI potential in Eq. (5), the Fourier transforms in $x$ and $y$, and the inverse Fourier transforms in $k_{x}$ and $k_{y}$ direction. Only one of these integrals can be completed analytically. As such, the final equations will involve an inverse Fourier transform in either the $k_{x}$ or $k_{y}$ coordinate in Fourier space, and this will need to evaluated numerically.

Following the above approach, we integrate out the $x$ and $y$ directions in the action functional, and employ the Euler-Lagrange equations. This leads to three coupled equations, an intrego-differential equation for $\tilde{f}$ and two integral equations for $\sigma_{x}$ and $\sigma_{y}$. In dimensionless form, these are,

[TABLE]

where

[TABLE]

and $\tilde{k}=\sqrt{\tilde{k}_{y}^{2}+\tilde{k}_{z}^{2}}$ and $q=\sqrt{2}\pi\sqrt{\tilde{k}_{y}^{2}+\tilde{k}_{z}^{2}}\,\tilde{\sigma}_{x}$, and $\text{erfc}(\cdot)$ is the complimentary error function. These are the variable-width GPE equations for the cigar-shaped dipolar BEC with arbitrary polarization angle $\alpha$. For $\alpha=0$, for which the cylindrical symmetry of the BEC is preserved, the integrals can be completed analytically and the equations reduce to our previous $\alpha=0$ result, Eqs. (17) and (18).

We proceed with $\alpha=\pi/2$, that is, dipoles polarized perpendicular to the primary axis of the trap, so as to break the symmetry about this axis. In the top panels of Fig. 4 we compare the ground state density according to the variable-width approach, fixed-width approach and the 3D dGPE, for $\gamma=10$ and $80$. Again, the variable-width approach more closely agrees with the 3D dGPE than the fixed-width approach. For both values of $\gamma$, the widths again deviate significantly from $\ell_{\perp}$ (bottom panels), emphasizing that the fixed-width assumption is invalid. Furthermore, for $\gamma=10$ (bottom-left panel) we see strong anisotropy in the $xy$ plane, due to the dipole-induced magnetostriction of the condensate in the $x$ direction. This emphasizes the importance of allowing the anisotropy in the transverse profile of the BEC. In the case for $\gamma=80$ (bottom-right panel), for which the system is more deeply in the 1D regime and the transverse trapping dominates over the interactions, the BEC approaches cylindrical symmetry.

These results are particularly useful, since the variable-width approach matches very closely with the 3D dGPE in the moderate trapping regime, $\gamma\sim 10$, where the transverse trap isn’t extreme enough to diminish the effects of the DDIs. This excellent agreement only occurs because we take into account the breaking of the cylindrical symmetry. Without the variational widths, the approach would only work in the case where the effect of the DDIs are effectively nullified by the extreme trapping ratio $\gamma$, as seen in Cai et al. (2010).

To characterise the anisotropy of the BEC in the $xy$ plane, we evaluate the transverse aspect ratio at $z=0$ (where $\sigma_{x}$ and $\sigma_{y}$ are maximal),

[TABLE]

In Fig. 5 we plot $\Gamma$ as a function of $\epsilon_{\rm dd}$ for two trap ratios, $\gamma=5$ and $10$. We see the $\Gamma$ decreases monotonically in both cases, corresponding to the increasing magnetostriction along $x$ as the DDI strength $\epsilon_{\rm dd}$ is increased (relative to the contact interactions). The transverse aspect ratio is closer to unity for the higher trap ratio $\gamma=10$. This is because for large $\gamma$, the transverse energy effectively diminishes the interaction energy, such that transverse wavefunction becomes closer to the isotropic harmonic oscillator state. Remarkably, we see excellent agreement between the variable-width predictions (markers) and the 3D dGPE (dashed lines) throughout. This demonstrates a marked improvement over the approach with isotropic fixed-width, for which $\Gamma$ would be fixed to unity throughout.

IV Variable-Width Two-Dimensional Equations

We now consider a flattened dipolar BEC confined in an oblate harmonic potential with $\gamma\ll 1$. Following Salasnich et al. Salasnich et al. (2002b) we write the wavefunction in the following separable form,

[TABLE]

where $f(x,y,t)$ is the two-dimensional wavefunction, and $\eta(x,y,t)$ is the variational width, which describes the extent of the BEC in the confined $z$-direction. In the standard 2D reduction Cai et al. (2010), $\eta$ is fixed to the oscillator length in the $z$-direction, $\ell_{z}=\sqrt{\hbar/(m\omega_{z})}$.

To derive the variable-width 2D equations, we follow the same process as for the 1D case, except we now integrate out the $z$-component of the 3D action functional, and the $z$ and $k_{z}$ components of the Fourier and inverse Fourier transforms, respectively, in the DDI potential given by Eq. (6). Fortunately, the integrals can be completed analytically, so there is no particular need to separate out the $\alpha=0$ equations as a special case. Also, we assume that the axial wavefunction, $\varphi(z,t;\eta(x,y,t))$ is slowly varying in the $xy$-plane, allowing us to write $\nabla^{2}\varphi(z,t)\approx\partial_{z}^{2}\varphi(z,t)$.

After the integration, the 2D action functional is,

[TABLE]

where $\nabla_{\perp}^{2}=\partial_{x}^{2}+\partial_{y}^{2}$. $\Phi_{\text{dd}}^{2\text{D}}$ is the 2D DDI term, derived from Eq. (6) by integrating out the $z$ and $k_{z}$-components of the Fourier transforms, and has the form,

[TABLE]

where $\mathcal{F}_{x,y}$ ($\mathcal{F}^{-1}_{x,y}$) is the 2D Fourier (inverse Fourier) transform, and

[TABLE]

where $\xi=\sqrt{2}\pi\eta k_{\perp}$, $k_{\perp}=\sqrt{k_{x}^{2}+k_{y}^{2}}$, and $\text{erfc}(\cdot)$ is the complimentary error function. We can see that this 2D DDI potential has a clear dependence on the polarisation direction $\alpha$.

An important aspect of this DDI potential is the asymmetric dependence on the Fourier space coordinates $k_{x}$ and $k_{y}$, which is in contrast to the 1D DDI potential, which had a dependence only on $k_{z}$ (after integration) and hence maintained symmetry. This asymmetry takes into effect magnetostriction already, and as such, a fixed-width approach to calculating the ground state for dipoles polarised away from the axis of symmetry will prove to be quite accurate as compared to the 3D dGPE. Furthermore, a calculation of the transverse aspect ratio can be reached without the inclusion of a variational width. The inclusion of $\eta(x,y,t)$, which we recall is a measure of the width in the $z$-direction, simply increases the accuracy of these predictions. This was not the case for the 1D BEC, as the fixed-width approach set the aspect ratio to unity, and no magnetrostriction was observed unless we let the orthogonal widths $\tilde{\sigma}_{x}$ and $\tilde{\sigma}_{y}$ vary.

For this 2D case, we employ units defined by the trapping in the transverse plane, with a time-scale $1/\omega_{\perp}$ and a length scale $\ell_{\perp}=\sqrt{\hbar/(m\omega_{\perp})}$, and transform our dimensions to the dimensionless forms according to,

[TABLE]

The dimensionless 2D wavefunction becomes $\tilde{f}=\ell_{\perp}f$. The Euler-Lagrange equations with respect to the 2D wavefunction and variable width $\tilde{\eta}$ giving the following equations,

[TABLE]

where $\tilde{\beta}_{2\text{D}}=2\sqrt{2\pi}\tilde{a}_{s}N$, and

[TABLE]

Here $\mathcal{F}_{\tilde{x},\tilde{y}}$ ($\mathcal{F}_{\tilde{x},\tilde{y}}^{-1}$) denotes the Fourier (inverse Fourier) in the $\tilde{x}$ an $\tilde{y}$ ($\tilde{k}_{x}$ and $\tilde{k}_{y}$) directions. These equations have a similar form to the 1D equations for the arbitrary polarization direction (Eqs. (20),  (21), and (22)), maintaining a direct dependence upon $\alpha$. There is only one algebraic equation since $\tilde{\eta}$ is a measure of the physical extent of the BEC in the $z$-direction only.

First we consider the dipoles to be polarized in $z$, $\alpha=0$. In Fig. 6 (top panels) we plot the ground state density profile along the $x$-axis, $|\tilde{f}(\tilde{x},\tilde{y}=0)|^{2}$, for trap ratios $\gamma=1/10$ and $1/80$. The variable-width results (black dotted line) are in excellent agreement with those of the 3D dGPE (blue solid line) for both trap ratios. For $\gamma=1/10$, the fixed-width approach noticeably underestimates the correct density. This discrepany is reduced for $\gamma=1/80$, as the axial trapping energy dominates over the interaction energy of the BEC.

We next consider the dipoles to be polarized along $x$, $\alpha=\pi/2$, with the results shown in Fig. 7. Again, the variable width approach accurately reproduces the results of the 3D dGPE, while there is a noticeable deviation with the fixed-width approach. Finally, we consider the transverse aspect ratio of the BEC, $\Gamma$, measured from the physical extent of the wavefunction, $\tilde{f}(\tilde{x},\tilde{y})$ in the $xy$ plane. The variation of $\Gamma$ with $\epsilon_{\rm dd}$ is shown in Fig. 8 for two trap ratios. The aspect ratio decreases as the dipoles are increased in strength, due to the increasing magnetostriction in the $x$ direction. The variable-width predictions closely match 3D dGPE throughout.

V Conclusions

In this paper, we have derived effective 1D and 2D mean-field equations describing cigar-shaped and pancake-shaped dipolar Bose-Einstein condensates with arbitrary polarization angle. Previous works in this direction have obtained effective 1D and 2D equations by assuming that the condensate profile is fixed to the ground Gaussian harmonic oscillator state in the tightly-confined direction; this “fixed-width" approach is valid for strictly quasi-1D and quasi-2D BECs. Here we extend on this by allowing the Gaussian profile to have a variational width which is solved self-consistently; this “variable-width" approach leads to effective 1D and 2D equations for the low-dimensional wavefunction, coupled to equations for the widths of the tightly-confined directions.

We compute the ground state of cigar-shaped and pancake-shaped dipolar BECs, and show that our effective equations give excellent agreement with the full 3D dGPE solutions over a range of trap ratios, polarization angles and interaction strengths. In particular, they provide a significant improvement over the fixed-width approach in that they remain accurate and valid as the system geometry is relaxed away from the strict quasi-1D and quasi-2D regimes. In addition to this, they are more efficient to numerically solve than the 3D dGPE.

These approaches may provide a convenient means to model the high-dimensional physics of the roton and dipolar supersolid phases Dutta and Meystre (2007); Wilson et al. (2008); Lu et al. (2010a); Tanzi et al. (2019); Böttcher et al. (2019); Chomaz et al. (2019) within a low-dimensional model; in order to study the supersolid phases, our approach would need to be extended to incorporate the role of quantum fluctuations. A further interesting case would be to use this approach to study the dynamics of dark solitons in cigar-shaped BECs. A significant discrepancy has been noted between the prediction of the 3D GPE and standard 1D reduction Bland et al. (2017), and the success of the NPSE in capturing observed dark soliton oscillations in non-dipolar BECs Weller et al. (2008) suggests that this approach could well reconcile this discrepancy.

Acknowledgments

M. J. K. is supported by an Australian Government Research Training Program (RTP) Scholarship and A. M. M. would like to thank the Australian Research Council (Grant No. LE180100142) for support. T. B. and N. G. P. thank the Engineering and Physical Sciences Research Council (Grant No. EP/M005127/1) for support.