Numerical calculation of dipolar quantum droplet stationary states

We describe and benchmark a method to accurately calculate the quantum droplet states that can be produced from a dipolar Bose-Einstein condensate. Our approach also allows us to consider vortex states, where the atoms circulate around the long-axis of the filament shaped droplet. We apply our approach to determine a phase diagram showing where self-bound droplets are stable against evaporation, and to quantify the energetics related to the fission of a vortex droplet into two non-vortex droplets.


I. INTRODUCTION
Dipolar condensates consist of highly magnetic atoms that interact with a long-ranged and anisotropic dipole-dipole interaction (DDI). Experiments with dipolar condensates of dysprosium [1-3] and erbium [4] atoms have prepared the system into one or several self-bound droplets that can preserve their form, even in the absence of any external confinement. These droplets occur in the dipole-dominated regime, where the short ranged s-wave interactions are tuned to be weaker than the DDIs [5,6]. In this regime the standard meanfield theory provided by the Gross-Pitaevskii equation is inadequate, as it predicts that the condensate is unstable to mechanical collapse. Here (beyond meanfield) quantum fluctuation corrections become important [7][8][9] and stabilize the droplets against collapse [10][11][12][13]. Droplets have been produced with 10 3 -10 4 atoms and peak densities roughly an order of magnitude higher than that of typical condensates. These droplets are still well within the dilute regime (i.e. na 3 1, where n is the density and a is the interaction length scale). In this regime the extended Gross-Pitaevskii equation (EGPE) is expected to provide a good description, and indeed the EGPE has been successfully used to model the equilibrium and dynamical properties of droplets. The excitation spectrum of the droplets has been the subject of theoretical and experimental studies [4,6,14], and more recently the possibility of preparing a droplet in an excited vortex state has been considered [15,16], although these have been found to be unstable.
So far little attention in the literature has been given to the details of EGPE calculations for stationary dipolar droplet states. The DDI needs to be treated with care since it is singular and long-ranged. Because the droplets are small, dense, and highly elongated (taking a filament-like shape), it becomes more difficult to treat the DDI accurately without taking large and dense numerical grids. Due to these technical challenges few groups outside of those already working on dipolar condensates have reported calculations for these droplets. In contrast we note that quantum droplets have also been realized in two component condensates [10,[17][18][19]. The absence of DDIs in these droplets makes the calculations more straightforward, and a rather large number of theory groups have reported work on this system.
In this paper we aim to outline a numerical technique to calculate self-bound quantum droplet states, including the case where the droplet has a vortex. We demonstrate the importance of using a truncated interaction potential to accurately evaluate the DDI energy, and introduce a simple gradient flow technique for obtaining the stationary solutions of the EGPE. We present benchmark results for the energy and chemical potential for vortex and non-vortex droplet stationary states. As applications we calculate a phase diagram for the stability of the droplet states, and quantify the tendency for vortex droplets to undergo fission into a pair of non-vortex droplets.
The outline of the paper is as follows. In Sec. II we introduce the EGPE for describing quantum droplet states. We also rewrite the formalism in cylindrical coordinates and introduce our dimensionless units. In Sec. III we introduce the Besselcosine based method that we use to discretize the EGPE and to enforce states to have a particular value of angular momentum. A variational approach is also presented and used to test the accuracy of the discretization when evaluating the various energy terms related to the EGPE. The gradient flow (or imaginary time) method is introduced in Sec. IV to obtain energy minimising solutions of the EGPE. The main results are presented in Sec. V including: benchmark results for droplet energies and chemical potentials, a phase diagram indicating where the droplets are stable to evaporation, and results quantifying the energetic instability of vortex droplets to fission. We then conclude in Sec. VI.

II. FORMALISM
A. Extended Gross-Pitaevskii Equation

General formulation
Several works [2- 6,[11][12][13][20][21][22][23][24] have established that the ground states and dynamics of a dipolar condensate in the droplet regime is well-described by the EGPE. In this formalism the time-independent stationary state wavefunction Ψ is a Here V tr describes any trapping potential, µ is the chemical potential and g s = 4πa s 2 /M is the s-wave coupling constant, arXiv:2012.11103v1 [cond-mat.quant-gas] 21 Dec 2020 with a s being the s-wave scattering length. The potential describes the effects of the long-ranged DDIs where the kernel is This is written for the case of dipoles polarized along the z axis by an external field, where θ is the angle between r and the z axis. The dipole coupling constant is g dd = 4πa dd 2 /M , where a dd = M µ 0 µ 2 m /12π 2 is the dipole length determined by the magnetic moment µ m of the particles. The leading-order quantum fluctuation correction to the chemical potential for a uniform system of density n is s π (1 + 3 2 2 dd ) and dd ≡ a dd /a s [2, 8,13]. The effects of quantum fluctuations are included in Eq. (1) using the local density approximation n → |Ψ(x)| 2 . Stationary EGPE states are also local minima of the energy functional where represent the kinetic, trap potential, interaction and quantum fluctuation energies, respectively.

Dimensionless cylindrical formulation
We now write the problem in a form utilizing the cylindrical symmetry and introducing natural units. While we do not present any results here that include a trapping potential (we focus on self-bound states in the absence of confinement), for generality we include a cylindrically symmetric trapping potential in our formulation. We take this to be of the form V tr = 1 2 M (ω 2 ρ ρ 2 + ω 2 z z 2 ), where ρ = x 2 + y 2 , and {ω ρ , ω z } are the trap frequencies. For this case the system is cylindrically symmetric (since we have also chosen the dipoles to be along z) and we can choose to look for stationary solutions of the form where ψ is real, φ = arctan(y/x) and s is an integer specifying the z-component angular momentum of the state. By separating variables, and introducing units of length x 0 = a dd and energy E 0 = 2 /M a 2 dd , we arrive at the effective cylindrical GPE Here is the single particle Hamiltonian, which features the Bessel differential operator We also note that for this choice of units the nonlinear coupling constants are g s → 4π −1 dd , g dd → 4π, and thus all specified in terms of dd . The convolution used to evaluate Φ is performed in threedimensions, but the resulting Φ is a cylindrically symmetric function: where the Fourier transformed density and DDI kernel arẽ n(k ρ , k z ) = 2π dρ dz ρJ 0 (k ρ ρ)e −ikzz |ψ(ρ, z)| 2 , (17) respectively. We also note that here we choose the condensate to be normalized to the number of atoms N , i.e.

III. BESSEL-COSINE NUMERICAL REPRESENTATION
To accurately treat the terms appearing in the EGPE operator we use a different discretization for the radial and axial dimensions, and we discuss these separately in the next sub-

sections.
A. Radial Bessel treatment

Bessel grid and quadrature
In the radial direction we consider N ρ points, that nonuniformly span the interval (0, R), given by which we refer to as the q-order Bessel grid. Here K q = α q Nρ+1 /R , and {α qi } are the ordered non-zero roots of the Bessel function J q (x) of integer order q. In this work we will need several such radial grids, each with the same number of points and range but of different orders. For this reason we need to adopt a notation that explicitly indicates the order. We also introduce the reciprocal space grid of the same order that spans the interval (0, K q ).
The radial grid points can be associated with a quadraturelike integration formula [25,26] where w qi are the real-space quadrature weights and we have used the notation g qi = g(ρ qi ) to denote the function g(ρ) sampled on the q-order grid. This integration requires that the functions of interest have limited spatial range, i.e. g(ρ > R) = 0.
Similarly, for the reciprocal k ρ -space we have a quadraturelike integration formula for a functiong(k ρ ) wherew qi are the k ρ -space quadrature weights andg qi =g(k qi ). Result (24) is only valid on our grid if the function is bandwidth limited, i.e.g(k ρ > K q ) = 0.

Hankel transformation
The Bessel grid is useful because it allows an accurate twodimensional (2D) Fourier transformation of functions of the form where we have used ρ to denote the planar position vector with polar coordinates (ρ, φ). The 2D Fourier transform of F isF wheref is the m-th order Hankel transform, arising because the angular integral of e imφ−ikρ·ρ yields the J m Bessel function. Here we have introduced k ρ to represent the 2D k-space vector, with polar coordinates (k ρ , φ k ).
For the case where the angular momentum of the function [i.e., m from Eq. (26)] and the order of the grid used to sample the function are the same (i.e. q = m) an accurate discrete Hankel transform can be implemented. For this case the transform yieldsf sampled on the k qi grid from f sampled on the ρ qi grid (see [27]). The explicit form of this discrete transform is obtained by evaluating (28) using the quadrature formula (22) where is the q-order Hankel transformation matrix and f qj = f (ρ qj ) etc.
Similarly, the inverse 2D transform is accomplished by the inverse Hankel transform with discrete form The Discrete Hankel transform matrices are not exact inverses of each other, but for typical grid sizes (N ρ ∼ 10 2 ) we have , which is adequate for our purposes.

Radial Laplacian operator
Hankel transforms are particularly useful for accurately evaluating the kinetic energy operator in radially symmetric cases (e.g. see [28]). To see this we note that by separating variable the 2D Laplacian ∇ 2 ρ acting on the function f (ρ)e isφ is equivalent to the Bessel differential operator D s (13) acting on f (ρ) [cf. Eq. (12)]. Using that J s are eigenfunctions of , we can utilize the Hankel transform to act on a radial function with D s : Allowing us to implement a spectrally accurate radial kinetic energy matrix [i.e. discretized version of T ρ = − 1 2 D s , c.f. Eq. (12)]: B. Axial cosine treatment

Trigonometric grids and quadrature
Our system has reflection symmetry along z so functions of interest will be of definite parity and here we will concern ourselves with the case of even parity. Utilizing this symmetry we use a half-grid on the interval (0, Z) spanned by N z equally spaced points where ∆z = Z/N z . The corresponding reciprocal grid spans the interval (0, K z ) where ∆k = π/Z and K z = π/∆z. The appropriate quadrature for this grid is the rectangular rule where g i = g(z i ).

Cosine transformations
The 1D Fourier transform for an even parity function is given by the cosine transform We can discretize this transform onto our chosen grid as where is the transformation matrix, This can be identified as the type-IV discrete cosine transformation (e.g. see [29]).

Axial Laplacian operator
Utilizing that the derivative operator is diagonal in k zspace, we can implement a discrete second derivative operator as This allows us to define a spectrally accurate axial kinetic

C. Treatment of EGPE operators
To find eigenstates of the effective cylindrical EGPE problem (10) we combine the radial and axial treatments described above to define a cylindrical (2D) mesh of points ρ ij ≡ (ρ si , z j ), choosing the radial order q to match the stationary state circulation s.

Local interaction terms
The contact interaction term in Eq. (11) is a nonlinear term that is local in position space and is evaluated as and similarly for the quantum fluctuation term Note: we have assumed that ψ ij is real, but modulus sign is needed on the quantum fluctuation term to properly deal with any case where ψ ij is negative.

DDI term
The DDI potential Φ can be evaluated using the convolution theorem in cylindrical coordinates, as given in Eq. (16). To evaluate this expression we first need to Fourier transform the density (17), which we obtain by applying the cosine transform along z and the quadrature rule to evaluate the radial transformñ Here we use the prime to indicate that the Fourier transformed density is evaluated on the cylindrical k-space mesh of grid points k ij = (k 0i , k j ), noting that the radial k-space points are defined for the 0-order Bessel grid 1 . We then evaluate the potential Φ from Eq. (16) as and thus the full DDI term The functionĨ(k) is of critical importance for accurate numerical calculations of the DDIs. For clarity we hereon refer to the analytic form of this function introduced in Eq. (18) as the bare k-space potentialĨ bare . This result is singular since the k → 0 limit does not exist, reflecting the long-ranged anisotropic character of the interaction. As such the quadrature in (48) will converge slowly as the density of k-space points increases, or equivalently as the spatial ranges R, Z are made larger (also see discussion in Refs. [30,31]). This slow convergence can be understood as the Φ having contributions from periodic copies (arising from using Fourier transforms on a finite interval) of the density distribution. These periodic copies, separated by twice the grid range, interact with each 1 The radial transform in Eq. (47) is identical to (29) other causing a shift in the interaction energy. The unphysical influence of these copies decays with the cube of the grid range, making it impractical to calculate accurate results using I bare .
One approach to deal with this problem is to use spherical coordinates in which the Jacobian introduced by the coordinate transform removes the singularity ofĨ bare , but this requires the use of a non-uniform Fourier transform [31,32] in the algorithm. An alternative approach, first proposed for dipolar BECs in Ref. [30] (also see [33,34]), is to introduce a truncation of the DDI, i.e. restrict the range of the DDI to the physical extent of the grids used, so that interactions between periodic copies are formally zero (also see [35], and related treatments for the Coulomb case [36,37]). We choose to follow this approach here since it allows us to work with the cylindrical grid and associated transforms that we have introduced. The issue now becomes how to obtain the appropriate truncated DDI potential in k-space.
First let us consider a spherical truncation derived and applied to dipolar BECs in Ref. [30] (also see [38]). Here the truncated interaction in position space is i.e. a cutoff radius r c such that the DDIs do not occur between any two points separated by a distance of more than 2r c . In practice choosing r c = min{R, Z} ensures that no interactions can occur between periodic copies of the density. Fortunately the Fourier transform of I sph can be calculated analyticallỹ and is seen to regularize the behaviour ofĨ bare (k) near k = 0. The spherically truncated potential is most useful for situations where the condensate is nearly spherical so it is natural to choose R ≈ Z. For highly anisotropic cases the natural grid choice is R Z or R Z for cigar or pancake shaped condensates, respectively. Here the spherical truncation is impractical as the range of the narrow dimension would have to be extended to match the range of the other dimension, introducing many redundant grid points. In these situations it is natural to introduce a cylindrical truncation The Fourier transform of this truncated kernel,Ĩ cyl , does not have an analytic result and needs to be calculated numerically.
We refer the interested reader to Ref. [39] for details about the numerical calculation.
In the testing we perform we evaluate the DDI term (49) using Φ ij evaluated according to (48) withĨ set to one of {Ĩ bare ,Ĩ sph ,Ĩ cyl }. If not specified otherwise, we useĨ cyl

D. EGPE operator and energy functional
We have now discussed all the operators needed to evaluate the EGPE operator i.e.
The expectation of the EGPE operator, normalized by the field norm gives the expected value of the chemical potential. This is evaluated numerically as where dV i = 4π∆zw si denotes the combined quadrature weights, and N [ψ] = ij ψ 2 ij dV i is the normalization. For a stationary state solution of the EGPE ψ with eigenvalue µ [see (10)] we have that µ EGP [ψ] = µ.
We can also evaluate the energy functional for the system [see Eq. (4)] which will be a local minimum for stationary solutions of the EGPE. Here we have introduced the labels E C , E D and E Q to refer to the standard GPE energy functional (including single particle and contact interactions), the dipole interaction energy and the quantum fluctuation energy, respectively.

E. Gaussian variational solution
It is useful to have a simple variational solution as an initial guess for the numerically calculated stationary state solutions and to validate the accuracy of our numerical methods. Here we consider a Gaussian state with angular circulation of s where the cylindrical amplitude is with width parameters {σ ρ , σ z }.
We can use this state to analytically evaluate the separate energy terms E C , E D , and E Q , as identified in Eq. (55). We obtain (also see [5,15]) where c s = (2s)!/4 s (s!) 2 , We have also used g s and g dd in place of their dimensionless values introduced earlier to make it easier to transform these results to other choices of units. The ansatz (56) can be used to provide a variational solution to the GPE. This involves minimizing the nonlinear function to determine {σ ρ , σ z }.

F. Gaussian test of numerical representation
We use the analytic results (58)-(60) to benchmark the accuracy of our numerical evaluation of the various terms appearing in the GPE. To do this we sample the variational Gaussian solution on a cylindrical grid and numerically evaluate the terms corresponding to the individual energy contributions as specified in Eq. (55), and compute the absolute value of the relative error with respect to the analytic results.
We show the results in Table I for various values of s, different grid choices and methods for evaluating the DDI term. For all our grid choices the single particle and contact interaction term (E C ) has an absolute relative error of less than 10 −13 , so we do not list it in the table. Instead we focus on the DDI and quantum fluctuation terms that tend to have larger errors.
For the DDI term we present results for the bare, spherically truncated, and cylindrically truncated interactions (see Sec. III C 3). Evaluating E D using the bare interaction is always inaccurate, and converges slowly to the exact result as the grid range increases [cf. cases (a) and (b)]. The spherically cutoff interaction works well when the grid has a similar radial and axial range [cf. cases (b) and (c)], and is thus most  useful for states where the density distribution has a similar radial and axial extent. The cylindrically cutoff interaction is seen to work well in all cases including highly anisotropic situations.
The accuracy of the quantum fluctuation term is reduced in the s = 1 case relative to the s = 0 and s = 2 cases for similar numbers of points. This finding is consistent with the analysis presented by Ogata [40] on the accuracy of numerical Bessel quadrature (22). Ogata shows that q-order Bessel quadrature converges exponentially for an integrand of the form |x| 2q+1 h(x)dx, if h(x) is analytic on the real axis (−∞, ∞). For our Gaussian ansatz (57) the quantum fluctuation term (46) including the Jacobian (and neglecting zcoordinates) is of the form ρ 5s+1 g(ρ), where g(ρ) represents the Gaussian part. Casting this integrand in the form Ogata analyses and taking s = q we have |ρ| 2s+1 [|ρ 3s |g(ρ)]. The term is square brackets is only analytic for s even. We note that by choosing s = q (e.g. taking q = 1/2) this integration can be performed more accurately, although we do not pursue this further here.

IV. GRADIENT FLOW SOLUTION OF THE GPE
Here we present a simple gradient flow solver based on our cylindrical discretization. This is an energy minimising scheme for finding ground states 2 . The gradient flow involves solving the time-dependent GPE in imaginary time, i.e. solving the flowψ = −L[ψ]. However, normalization of the field tends to decrease under this evolution, so it is necessary to renormalize during the evolution. We follow Ref. [41] (also see [42]) and discretize the evolution using a backwardsforwards Euler scheme. Here time is advanced in time steps ∆t n , as t n+1 = ∆t n+1 + t n , with n = 0, 1, 2, . . . and t 0 = 0. During such a step the updated wavefunction ψ + is obtained from the current wavefunction ψ(t n ) according to where (65) is the renormalization (projection), and we have introduced Notice that the term involving the kinetic energy operator is implemented with a backwards-Euler step, giving us good stability for dense grids (where the kinetic energy operator is large). By subtracting µ EGP in the V eff term, we ensure that to O(∆t 2 ) the field normalization is constant under the gradient flow (e.g. see [43]), which improves the performance of the algorithm. Hereon we suppress explicit notation of the position indices of the wavefunction for notational brevity, however terms appearing are to be evaluated as described in Secs. III C and III D.
The semi-implicit equation (64) is formally solved by inverting the spatial differential operator. This can be done using Fourier transformation F, yielding an explicit expression for ψ + : This can be efficiently implemented numerically using the operators and transforms we introduced earlier 3 . Using a forwards-Euler approach for the potential and nonlinear terms in Eq. (64) has the advantage that these terms are explicit, however care needs to be taken with the time-step to ensure that the algorithm is stable. In practice we choose the time step according to where || || ∞ denotes the maximum (over all spatial points) of the absolute value of the argument evaluated with ψ(t n ) and a ∆t 1 is a constant. For the results we present here we generally take a ∆t = 0.5 4 .
We are interested in obtaining the lowest energy stationary state subject to the imposed angular momentum s. We can quantify the backwards error through the residual r ≡ L[ψ] − µψ. In particular we use the measure where the N −1/2 factor is to make the measure independent of normalization choice 5 . We terminate the gradient flow once r ∞ decreases below a desired value (typically 10 −15 ). We mention that r ∞ has to be used with care. First, it depends on choice of units 6 , scaling as x −7/2 0 . Second, because µ can be negative (i.e. self-bound), or even zero, the sensitivity of this measure is also dependent on the case under consideration. In practice we can evaluate (69) for a particular solution ψ by applying the EGPE operator (53) and using (54) to obtain µ. Alternatively, we can approximately evaluate this as [see Eq. (64)]. As a second measure of solution quality we have developed a virial theorem for the EGPE. The virial relation is with the terms being the components of energy [see Eq. (4)], and where we have assumed that any trap is harmonic in form.
We have obtained this virial theorem by considering how the energy functional transforms under a scaling of coordinates (e.g. see [44]). The terms in Eq. (71) can be evaluated using the techniques described in Secs. III C and III D. In general we find that |Λ V | is a useful quantity to assess our numerical solutions, as it is sensitive to the residual r ∞ and the quality of the spatial discretisation.  Table II using a∆t = 0.5, showing (c) the error r∞, (d) time step ∆tn, (e) energy error, and (f) the absolute virial expectation error. In the energy error Eref is a reference energy calculated from a state using a larger more dense grid.

V. RESULTS
In Table II we present results for the energy, chemical potential and |Λ V | of stationary self-bound states with s = 0 and s = 1 obtained by the gradient flow method. The density profiles of the states are shown in Figs. 1(a) and (b). We also show the gradient flow evolution for case (b) of Table II in Figs. 1(c)-(f). We use the variational solution ψ G as the initial condition for the gradient flow and the time step ∆t n quickly settles to a value of ∆t n ≈ 6.4 (for a ∆t = 0.5) during the flow. The flow terminates at r ∞ = 5 × 10 −17 after 5591 steps, taking about 15 seconds to execute on a workstation computer. For this case the flow is stable for a ∆t 0.84, and takes a proportionately shorter amount of time to execute with a larger value of a ∆t , but becomes unstable and does not converge when a ∆t 0.84.
The results in Table II reveal how sensitive the energy and chemical potential are to the choice of grid [i.e. range (R, Z) and number of points (N ρ , N z )], showing that accurate results (greater than 9 significant figures) can be obtained when the grid ranges are approximately twice the spatial extent of this droplet and for sufficiently many points. This grid range dependence arises because the DDI term involves a convolution (e.g. see discussion in Ref. [45]). The s = 1 case typically requires more spatial points to have a similar level of accuracy as the s = 0 case. We expect that this arises from the poorer performance of the Bessel quadrature for the quantum fluctuation term when s = 1, as noted in Sec. III F. The results in the table also show that the absolute virial expectation |Λ V | is qualitative similar to the magnitude to the energy error (i.e. number of converged digits in E), suggesting that |Λ V | provides a useful additional method to characterize the solution accuracy.
We also consider the calculation of a phase diagram to predict where s = 0 and s = 1 self-bound droplets have a negative energy. This requirement ensures that the droplet is energetically stable against evaporating into the trivial E = 0 state where the atoms are dispersed over all space. This condition is adequate to ensure that s = 0 droplets are ground states that are dynamically stable. For the s = 1 case this requirement does not ensure dynamic stability as the droplet can decay by fission. We discuss this aspect later.
Our results for the regions where the droplet states have negative energy are shown in Fig. 2(a). Note by using coordinates of N and −1 dd this phase diagram has no remaining dependence on other parameters, and is in this sense universal. In general the s = 0 and s = 1 regions have similar shapes, although the s = 0 region is larger (extends to higher −1 dd ). This is because the s = 0 droplet has a considerably lower energy for the same parameters [e.g. see Fig. 2(b) and Table II]. This difference is from the large energy cost associated with hosting a vortex in the droplet. This arises from the kinetic energy of the vortex, but also because the s = 1 droplet is wider than the s = 0 droplet [cf. Figs. 1(a) and (b)], making the DDI energy less negative.
The results in Fig. 2 also compare the stationary EGPE solutions and variational Gaussian approach. For example, the markers in Fig. 2(a) show the boundary to the negative energy region determined by the EGPE solutions, while the boundary of the shaded region is determined variationally. The EGPE results are obtained using the gradient flow method to solve for a state at an initial (low) value of −1 dd stating from a variational Gaussian solution. The resulting EGPE solution for that case is then used as the starting point for the gradient flow at a slightly higher value of −1 dd , and so on. This process is stopped once the localized state disperses (unbinds and spreads out over the range of the grid). In Fig. 2(b) we show the energy of a sequence of states obtained this way for a particular case of N = 10 4 . We also indicate the point where E = 0, used to identify the boundary. The variational results are located as an energy minimum of the function given in Eq. (63). We follow this minimum as −1 dd increases, noting that eventually it becomes a local minimum (when E > 0) and then changes to a saddle point (causing the branch to terminate). We also mention that while the energy is positive near the end of the branches, the chemical potential remains negative [ Fig. 2(c)].
The phase diagram in Fig. 2(a) collates the results of many analyses of the type in Fig. 2(b) for different atom numbers N . In general the E = 0 boundary predicted by the variational Gaussian theory is close to that obtained from the EGPE solution. This is because the energies predicted by the two approaches are similar near the transition point. However, we note that the energy of the EGPE state tends to be significantly lower than the variational state for smaller values of −1 dd where the droplets are more deeply bound. In this regime [e.g. s = 0 state in Fig. 1(a)] the droplet density profile has a relative flattop axial density distribution, which is not well-captured by a Gaussian.
In Figs. 2(d) and (e) we show examples of the droplet states for cases with slightly positive energy. We can compare these states to the more deeply bound states (differing only in the values of −1 dd ) from Fig. 1. This comparison, particularly for the s = 0 case, emphasizes how much the droplet size can change with −1 dd . For example, over the range of −1 dd considered in Fig. 2(b) the axial length of the s = 0 state changes by approximately a factor of 5. This necessitates careful choice of numerical girds to ensure the states are calculated accurately as −1 dd changes. We can assess the stability of the s = 1 vortex droplet to undergoing fission, whereby it splits into two s = 0 droplets [see Fig. 3(a)]. This scenario has been observed in dynamical simulations of the s = 1 vortex droplet [15]. We can define a fission energy as where E s=1 N is the energy of a s = 1 vortex droplet with N atoms, and E s=0 N 2 is the energy of a s = 0 droplet with N/2 atoms. Thus ∆E N represents the excess energy of the vortex droplet compared to two independent (i.e. well-separated) s = 0 droplets. When this energy is positive, then the vortex is potentially unstable to fission. The angular momentum of the vortex droplet would need to be carried by the motion of the separating s = 0 droplets, meaning that kinetic energy will also be important in determining if fission can occur. Thus ∆E N > 0 is a necessary but not sufficient condition for fission.
In Fig. 3(b) we show some results for the fission energy computed using the EGPE and variational solutions. For the results with N = 10 4 both approaches predict that ∆E N > 0, and thus the vortex droplet is potentially unstable to fission. However, for the larger atom number of N = 10 5 the two predictions are different: the EGPE results predict that ∆E N > 0 at all values of −1 dd considered, while the variational results suggest stability (i.e. ∆E N < 0) for −1 dd 0.25. Here the variational approximation is failing because the droplets are not well described by the Gaussian ansatz. We note that this problem was originally considered in a preprint by Cidrim et al. [46], who presented EGPE results that appeared to support the variational prediction that ∆E N < 0 for large N and small −1 dd . Because ∆E N is determined by a difference in two calculations, it can be quite sensitive to the accuracy of the EGPE calculations. This example emphasizes the need for high accuracy calculations of dipolar droplets.

VI. CONCLUSIONS
In this paper we have presented a method to accurately solve for dipolar quantum droplets in a cylindrical geometry allowing for the inclusion of angular momentum. This work builds on the discretization introduced by Ronen et al. [30], extending it to include angular momentum in the stationary state, a cylindrical cutoff (truncation) of the DDI kernel, and the quantum fluctuation term. Using a simple gradient flow technique we demonstrate that this approach is able to obtain accurate results for the dense, highly elongated (filament shaped) quantum droplets that form in dipolar BECs in the regime where the DDIs dominate the contact interactions. We also show that without a careful treatment of the DDI term in this regime it may be difficult to obtain the droplet energy to better than ∼10% accuracy. Such errors would make calculations such as the fission energy (which depends in the difference of energies between two states) difficult to compute.
We have also presented benchmark energy and chemical potential calculations for self-bound droplets. There are very few such benchmark results in the literature and we expect these will be important for comparisons with different approaches that may be developed in the future. We present a generalization of the virial theorem for dipolar EGPE and find that this can be used to test the solution accuracy. As an application of our method we have also presented a phase diagram for the energetic stability of self-bound s = 0 and s = 1 droplets, and considered the stability of s = 1 droplets against fission.
There are many avenues for future development of this work. We note that the EGPE solutions we have presented were obtained using a simple but efficient backwardsforwards Euler gradient flow method. It would be of interest to consider other more efficient solvers such as conjugate gradient solvers (e.g. see [30,35,43]) or a fully implicit backwards-Euler method. Such solvers will be useful in situations where efficiency becomes more important, such as fully three-dimensional (3D) cases that cannot be reduced to cylindrical symmetry. Typically, non-cylindrically symmetric ground states occur when the confinement is not symmetric about the dipole polarization axis (e.g. [47]), or when the system favours the formation of a droplet array (including supersolid) [48,49]. This latter case is of significant interest in the community due to the recent observation of supersolid states [50][51][52][53][54][55]. Another area requiring attention in the fully 3D case is an efficient and accurate truncation scheme for the DDI kernel, which is a critical tool to enable an efficient grid choice to represent the state.