Welcome to visit Communications in Theoretical Physics,
Mathematical Physics

High precision solutions to quantized vortices within Gross–Pitaevskii equation

  • Hao-Hao Peng , 1 ,
  • Jian Deng , 2, ,
  • Sen-Yue Lou , 3 ,
  • Qun Wang , 1
Expand
  • 1Department of Modern Physics, University of Science and Technology of China, Hefei, Anhui 230026, China
  • 2Institute of Frontier and Interdisciplinary Science, Key Laboratory of Particle Physics and Particle Irradiation (MOE), Shandong University, Qingdao, Shandong 266237, China
  • 3School of Physical Science and Technology, Ningbo University, Ningbo 315211, China

Author to whom any correspondence should be addressed.

Received date: 2022-04-27

  Revised date: 2022-08-04

  Accepted date: 2022-08-04

  Online published: 2022-08-23

Copyright

© 2022 Institute of Theoretical Physics CAS, Chinese Physical Society and IOP Publishing

Abstract

The dynamics of vortices in Bose–Einstein condensates of dilute cold atoms can be well formulated by Gross–Pitaevskii equation. To better understand the properties of vortices, a systematic method to solve the nonlinear differential equation for the vortex to very high precision is proposed. Through two-point Padé approximants, these solutions are presented in terms of simple rational functions, which can be used in the simulation of vortex dynamics. The precision of the solutions is sensitive to the connecting parameter and the truncation orders. It can be improved significantly with a reasonable extension in the order of rational functions. The errors of the solutions and the limitation of two-point Padé approximants are discussed. This investigation may shed light on the exact solution to the nonlinear vortex equation.

Cite this article

Hao-Hao Peng , Jian Deng , Sen-Yue Lou , Qun Wang . High precision solutions to quantized vortices within Gross–Pitaevskii equation[J]. Communications in Theoretical Physics, 2022 , 74(9) : 095002 . DOI: 10.1088/1572-9494/ac86bd

1. Introduction

Vortices are strongly nonlinear excitations in superfluids, which are quantized as topologic defects from the long range quantum phase coherence [13]. The formation, stability and dynamical properties of vortices have been intensively studied experimentally [48] and theoretically [9, 10]. Vortices play important roles not only in many-body quantum systems, but also in dark matter [1113] and the phase transition of the early universe [1417]. The structure of nonlinear dynamics for vortices is in connection with the gravitational field equation for the metric of black holes, which makes it possible to study the Penrose process of rotating black holes in laboratories [18, 19]. Vortices are also necessary degrees of freedom in turbulence and the intermediate state of an over-populated off-equilibrium system. The generation and clustering of vortices and the annihilation of vortex–antivortex pairs are closely related to the scaling law of a non-thermal fixed point [2026].
The static and dynamical properties of vortices and vortex lattices at zero temperature limit can be well described by the Gross–Pitaevskii equation (GPE), which is a nonlinear Schrödinger equation with the nonlinearity being determined by interactions through a mean-field approximation [27, 28]. This equation emerges in various nonlinear phenomena, and has been extensively studied in connection with nonlinear optics, plasma physics and fluid dynamics. The GPE helps the prediction and description on experimental observations of nonlinear effects, such as vortices and interaction among them [29, 30].
As an individual object, the vortex has an internal structure and can also be involved in an external evolution. How the wave function of a vortex responds to a perturbation or to the nonuniformity of the condensate is an interesting and complicated problem. The structure of the vortex core is determined by the properties of the vortex itself, but it should also match the outer wave function determined by the collective motion and interaction with the environment. Through a unique coordinate transform and matching the asymptotic expansion inside and outside the vortex core, a set of differential equations have been derived to describe the overall structure and collective motion of the vortex in a perturbative way [31, 32].
To construct a perturbation theory, the ideal starting point is a system with an exact solution. If the disturbance is not large, a series of perturbative solutions can be obtained systematically, so that a reliable understanding of the complex system can be achieved. The various physical quantities associated with the perturbed system can be expressed as corrections to the original one. For example, the contribution from the trace anomaly of the QED energy momentum tensor to the mass of the hydrogen atom, the simplest QED bound state, is calculated in [33], where the Lamb shift from the vacuum polarization and radiative corrections to the electron mass are perturbative corrections to the ground state with an exact wave function.
To study the vortex dynamics in a perturbative way, the stationary wave function of a vortex in a uniform condensate is the starting point. But due to the nonlinear nature of this problem, an exact solution to the nonlinear ordinary differential equations (ODE) has not been found yet. The best alternative we can think of is to look for an approximate function that is simple, compact and as precise as possible. This can be done with the method of two-point Padé approximants. The basic idea of Padé approximants is to construct a rational function which can reproduce the Taylor expansion series of the target function within a given order. The rational function can mimic the singularity of the target function since it contains pole structure, so the convergence radius of the original Taylor expansion can be extended (see an example in [34]). If the target function is constrained by two boundary conditions, the rational function should reproduce two Taylor expansion series simultaneously, which is called two-point Padé approximants [35, 36]. We found that this scheme is particularly suitable for finding approximate solutions to the nonlinear equation for the vortex. The main purpose of this paper is to explore the precision limit of this scheme, hoping to provide a clue for the exact solution of the vortex.
The outline of this paper is as follows. In section 2, the static and rotationally symmetric equation for the voetex is derived. In section 3, a very high precision numerical solution to the vortex equation is obtained with a relaxation method, which can be used to test the accuracy of the approximate solutions. In section 4, the results from two-point Padé approximants are presented. In section 5, the errors and challenge of this approach are discussed. The conclusions and discussions are given in section 6.

2. Vortices within GPE

The Bose–Einstein condensate (BEC) in cold atom systems at zero temperature can be well described by the Gross–Pitaevskii (GP) theory as the mean-filed approximation of quantum field theories [1]. In GP theory, the ground state and weakly excited states of the condensate are described by the complex wave function ψ(r, t), which satisfies GPE,
$\begin{eqnarray}{\rm{i}}{\hslash }\displaystyle \frac{\partial }{\partial t}\psi ({\text{}}r,t)=-\displaystyle \frac{{{\hslash }}^{2}}{2m}{{\rm{\nabla }}}^{2}\psi ({\text{}}r,t)+\lambda | \psi ({\text{}}r,t){| }^{2}\psi ({\text{}}r,t),\end{eqnarray}$
where ℏ is the reduced Planck constant, m is the particle mass, and the coupling constant λ can reproduce the s-wave scattering length, λ = 4πas2/m in the Born approximation. The investigation on the exact and numerical solutions to GPE improves our understanding on the nonlinear dynamics of matter waves in BEC [30]. In this work, we will focus on the static and rotationally symmetric solution to GPE, the quantum vortex, a kind of topological excitation in a superfluid in presence of the local orbital angular momentum. The wave function can be written in the following form
$\begin{eqnarray}\psi ({\text{}}r,t)=\sqrt{{n}_{0}}f\left(\eta =\displaystyle \frac{| {\text{}}r| }{\xi }\right){{\rm{e}}}^{{\rm{i}}s\varphi }{{\rm{e}}}^{-{\rm{i}}\mu t},\end{eqnarray}$
where n0 is the bulk number density, $\xi =1/\sqrt{2m\lambda {n}_{0}}$ is the healing length, φ is the azimuthal angle, μ is the chemical potential, s is the winding number which must be an integer to keep the wave function single-valued. The normalized amplitude function f(η) satisfies the nonlinear ordinary differential equation (ODE)
$\begin{eqnarray}\displaystyle \frac{1}{\eta }\displaystyle \frac{{\rm{d}}}{{\rm{d}}\eta }\left[\eta \displaystyle \frac{{\rm{d}}f(\eta )}{{\rm{d}}\eta }\right]+\left(\displaystyle \frac{\mu }{\lambda {n}_{0}}-\displaystyle \frac{{s}^{2}}{{\eta }^{2}}\right)f(\eta )-f{\left(\eta \right)}^{3}=0.\end{eqnarray}$
Physically, the area modified by a vortex is rather limited, the density profile should return to the bulk value in a region far away from the vortex core, which requires μ/λn0 = 1. So the ODE can be put into the form
$\begin{eqnarray}f^{\prime\prime} (\eta )+\displaystyle \frac{1}{\eta }f^{\prime} (\eta )+(1-\displaystyle \frac{{s}^{2}}{{\eta }^{2}})f(\eta )-f{\left(\eta \right)}^{3}=0,\end{eqnarray}$
with the boundary conditions f(η → 0) = ksηs and f(η →∞ ) = 1. Here ks is the connecting parameter, which plays an important role in the vortex solution.
The purpose of this work is to solve equation (4). Without the cubic term, the ODE becomes Bessel's differential equation with Bessel functions as the formal solution, but the boundary condition at infinity cannot be satisfied. With the cubic term, however, the equation is a nonlinear ODE. It is very hard, if not impossible, to obtain the exact analytical solution to it. So far, we can only find a similar example with an exact solution. It can be easily checked that $f=\tanh (\eta /\sqrt{2})$ is one solution to the nonlinear ODE, f″ + f(1 − f2) = 0 with the same boundary conditions. This solution is a direct extension of the Riccati equation. In appendix A, we will demonstrate another way to obtain this exact solution which may provide some hints for solving equation (4). The example gives us a little hope that a sufficiently precise numerical solution can guide the guesswork about the possible form of the analytical solution.
For the numerical solution, the differential equation with two-point boundary conditions can be solved with shooting method, which requires a fine tuning of ks to achieve the right limit f(ηηC) → 1 and $f^{\prime} (\eta \to {\eta }_{C})\to 0$ in the range η < ηC. The shooting results of f(η) for the winding numbers s = 1, 2, 3 are shown in figure 1. It seems that shooting works well, but the results are not stable for a large ηC, because a rigorous precision in ks and an inaccessible small step size in the finite difference are required. Fortunately, we find a relaxation method which can achieve a very high precision in the whole region of 0 < η < ∞ . The technical details will be discussed in the next section.
Figure 1. The profile function f(η) are obtained with shooting method in the range η < ηC = 10. Black solid line, red dashed line and blue dash-dotted line stand for the results with s = 1, 2, 3 respectively.

3. Relaxation iteration for a high precision numerical solution

A high precision numerical solution to equation (4) can be achieved by a relaxation iteration. The method is built up with three steps. Firstly, we take a replacement for the target function, f(η) → ηs/(ηs + 1) + h(η). The auxiliary function ηs/(ηs + 1) is to make the unknown function h(η) satisfy the boundary condition h(0) = h( ∞ ) = 0, which is convenient for solving the ODE. The function h(η) satisfies the nonlinear ODE as follows
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{{{\rm{d}}}^{2}h}{{\rm{d}}{\eta }^{2}}+\displaystyle \frac{1}{\eta }\displaystyle \frac{{\rm{d}}h}{{\rm{d}}\eta }+\left[1-\displaystyle \frac{{s}^{2}}{{\eta }^{2}}-\displaystyle \frac{3{\eta }^{2s}}{{\left(1+{\eta }^{s}\right)}^{2}}\right]h\\ =\,{h}^{3}+\displaystyle \frac{3{\eta }^{s}}{1+{\eta }^{s}}{h}^{2}+\displaystyle \frac{{\eta }^{s}\left({s}^{2}{\eta }^{s}({\eta }^{s}+3)-{\eta }^{2}(2{\eta }^{s}+1)\right)}{{\eta }^{2}{\left(1+{\eta }^{s}\right)}^{3}},\end{array}\end{eqnarray}$
where we have put all linear terms of h(x) to the left-hand-side. Secondly, we replace the variable to x = η/(R + η), so that the limit η → ∞ is shifted to x → 1. The positive number R sets a finite scale, we choose R = 1 without loss of generality. The nontrivial change is
$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{d}}}^{2}h}{{\rm{d}}{\eta }^{2}}+\displaystyle \frac{1}{\eta }\displaystyle \frac{{\rm{d}}h}{{\rm{d}}\eta } & \to & \displaystyle \frac{{\left(1-x\right)}^{4}}{{R}^{2}}\displaystyle \frac{{{\rm{d}}}^{2}h}{{\rm{d}}{x}^{2}}-2\displaystyle \frac{{\left(1-x\right)}^{3}}{{R}^{2}}\displaystyle \frac{{\rm{d}}h}{{\rm{d}}x}\\ & & +\displaystyle \frac{{\left(1-x\right)}^{3}}{{R}^{2}x}\displaystyle \frac{{\rm{d}}h}{{\rm{d}}x}.\end{array}\end{eqnarray}$
In this way, the ODE is ready for the discretization and iteration.
Finally we discretize x as xi = iΔx, where i = 0, ⋯ ,N and the step size is Δx = 1/N. So the derivatives become
$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{d}}}^{2}h}{{\rm{d}}{x}^{2}} & = & \displaystyle \frac{{h}_{i+1}+{h}_{i-1}-2{h}_{i}}{{{\rm{\Delta }}}_{x}^{2}},\\ \displaystyle \frac{{\rm{d}}h}{{\rm{d}}x} & = & \displaystyle \frac{{h}_{i+1}-{h}_{i-1}}{2{{\rm{\Delta }}}_{x}}.\end{array}\end{eqnarray}$
Then the nonlinear ODE is changed to a set of linear algebraic equation
$\begin{eqnarray}{a}_{i}{h}_{i-1}+{b}_{i}{h}_{i}+{c}_{i}{h}_{i+1}={w}_{i}({h}_{i}),\end{eqnarray}$
for i = 1, ⋯ ,N − 1, where ai, bi, ci, wi(hi) are well defined at each point. Because of the well-designed boundary condition, h0 = hN = 0 does not interfere with the iteration, then equation (8) can be rewritten as
$\begin{eqnarray}\begin{array}{l}\left(\begin{array}{ccccccc}{b}_{1} & {c}_{1} & 0 & \cdots & 0 & 0 & 0\\ {a}_{2} & {b}_{2} & {c}_{2} & \cdots & 0 & 0 & 0\\ 0 & {a}_{2} & {b}_{3} & \cdots & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & {b}_{N-3} & {c}_{N-3} & 0\\ 0 & 0 & 0 & \cdots & {a}_{N-2} & {b}_{N-2} & {c}_{N-2}\\ 0 & 0 & 0 & \cdots & 0 & {a}_{N-1} & {b}_{N-1}\end{array}\right)\left(\begin{array}{c}{h}_{1}\\ {h}_{2}\\ {h}_{3}\\ \vdots \\ {h}_{N-3}\\ {h}_{N-2}\\ {h}_{N-1}\end{array}\right)\\ =\left(\begin{array}{c}{w}_{1}\\ {w}_{2}\\ {w}_{3}\\ \vdots \\ {w}_{N-3}\\ {w}_{N-2}\\ {w}_{N-1}\end{array}\right).\end{array}\end{eqnarray}$
We may start by choosing a set of reasonable values for hi and obtain wi(hi), then a new set of hi can be determined by solving the tridiagonal algebraic equation (9). In this way, a high precision solution can be achieved within dozens of iterations. In the last iteration, the maximum error to hi is smaller than 10−25, which means that the algorithm has a fast speed of convergence. The main source to the error comes from the discretization, which can be improved by increasing the number of points N. We find the maximum difference between the results with N = 225 and N = 226 is of the order 10−16, which sets the precision of our numerical solution. The main challenge to go beyond larger N is to handle a larger dimension vector with high precision on a personal computer, that is the reason we stopped at N = 226. For the vector with such a dimension, the iteration to achieve the required precision can be done in a few hours.
For applications, the numerical solution in the form of data sample is not convenient. It is better to look for an approximate analytical solution. The present numerical solution can be used to verify the accuracy of the approximate analytical solution.
The most important input for the approximate analytical solution is the value of connecting parameter ks. We can get it from the numerical solution in the following way. To determine ${k}_{s}={\mathrm{lim}}_{\eta \to 0}f(\eta )/{\eta }^{s}$ with a higher precision, the behavior of f(η) near η = 0 can be parameterized by
$\begin{eqnarray}f(\eta )\to \left(\displaystyle \frac{\eta }{\sqrt{1+{\eta }^{2}}}+\tilde{h}(\eta )\right){\tanh }^{s-1}(\eta ),\end{eqnarray}$
for example, for winding number s = 1, we find that f(η)/η can be best fitted by RN(1 − η2/8) for small η, where RN is the fitting parameter corresponding to the solution with N grid points. The left panel in figure 2 shows the numerical result for f(η)/ηRN with N = 226. The dependence of RN on N is in table 1 and visualized in the right panel of figure 2.
Figure 2. Left: numerical results for f(η)/ηRN with N = 226. Right: numerical results for ${R}_{N}-{k}_{1}^{\mathrm{ext}}=\beta {\gamma }^{{\mathrm{log}}_{2}N}$ as a function of ${\mathrm{log}}_{2}N$.
Table 1. The numerical values of RN and ${k}_{1}^{\mathrm{ext}}$ extrapolated from RN.
N = 2n RN
n = 18 0.583 189 495 872
n = 19 0.583 189 495 8634
n = 20 0.583 189 495 8611
n = 21 0.583 189 495 860 52
n = 22 0.583 189 495 860 37
n = 23 0.583 189 495 860 341
n = 24 0.583 189 495 860 332
n = 25 0.583 189 495 860 3300
n = 26 0.583 189 495 860 329 47
${k}_{1}^{\mathrm{ext}}$ 0.583 189 495 860 329 28
We may extract the value ${k}_{1}={\mathrm{lim}}_{N\to \infty }{R}_{N}$ by an extrapolation because these points can be best fitted in the form ${R}_{N}-{k}_{1}^{\mathrm{ext}}=\beta {\gamma }^{{\mathrm{log}}_{2}N}$ with fitting parameters ${k}_{1}^{\mathrm{ext}}$,β and γ. Then ${k}_{1}^{\mathrm{ext}}$ can be regarded as a high-precision approximation to the exact value of k1. Its reliability will be discussed in the following sections.
In the same way, the extrapolated values of ks for s = 2, 3 can be obtained
$\begin{eqnarray}{k}_{2}^{\mathrm{ext}}=0.153\,099\,102\,859\,539,\end{eqnarray}$
$\begin{eqnarray}{k}_{3}^{\mathrm{ext}}=0.026\,183\,420\,721\,62.\end{eqnarray}$
The precision of the connecting parameters is reduced for larger s, because it is harder to extract the small η behavior of f(η) ∝ ηs. But in any way, the present results are precise enough for our purpose. We note that these connecting parameters can be calculated in a semi-analytical scheme [37]. A similar method is discussed in section 5, in which the connecting parameters are traced by choosing the proper root of a polynomial with increasing order. The biggest problem of such semi-analytical schemes is to find the right root of a high-order polynomial, which is computationally expensive, meanwhile the improvement on the precision of ks is very limited in comparison with the results from the extrapolation method.

4. Two-point Padé approximants for the vortex solution

Analytically, we may try to mimic the solution to equation (4) in terms of Taylor expansion series from both points η = 0 and η = ∞ simultaneously. This is the basic idea of two-point Padé approximants. Without loss of generality, we may define g(η) by taking a shift 1/2 in f(η), g(η) = f(η) − 1/2, and then expand g(η) at η = 0 and η = ∞ in the form
$\begin{eqnarray}g(\eta )=\displaystyle \frac{{c}_{0}}{2}+\sum _{l=1}^{\infty }{c}_{l}{\eta }^{l},\,\,\,\,\,\,\,\,\eta \to 0,\end{eqnarray}$
$\begin{eqnarray}g(\eta )=-\displaystyle \frac{{c}_{0}}{2}-\sum _{l=1}^{\infty }{c}_{-l}{\eta }^{-l},\,\,\,\eta \to \infty ,\end{eqnarray}$
where c0 = −1 and other coefficients ci, ⋯ ,ci can be derived from equation (4). For η = 0, by inserting equation (13) into equation (4), we find cl = 0 for 0 < l < s and a recursive relation for cl for ls,
$\begin{eqnarray}\begin{array}{l}{c}_{l+2}\left[{\left(l+2\right)}^{2}-{s}^{2}\right]+{c}_{l}-\sum _{{k}_{1},{k}_{2},{k}_{3}\geqslant s}^{{k}_{1}+{k}_{2}+{k}_{3}=l}{c}_{{k}_{1}}{c}_{{k}_{2}}{c}_{{k}_{3}}\\ =\,0,\,\,\,\mathrm{for}\,\,\,l\geqslant s.\end{array}\end{eqnarray}$
The last term on the left-hand-side of equation (15) comes from the nonlinearity of equation (4). Without this term, the coefficients are just those of Bessel's functions of the first kind. For the complete form, we can read out that only the coefficients with even/odd indices can be nonzero for the even/odd s, and their values are solely determined by the first nonzero one, which is cs = ks. Take s = 1 as an example, we have
$\begin{eqnarray}\begin{array}{rcl}{c}_{3} & = & -\displaystyle \frac{1}{8}{c}_{1},\\ {c}_{5} & = & \displaystyle \frac{1}{192}({c}_{1}+8{c}_{1}^{3}),\\ {c}_{7} & = & -\displaystyle \frac{1}{9216}({c}_{1}+80{c}_{1}^{3}),\\ {c}_{9} & = & \displaystyle \frac{1}{737280}({c}_{1}+656{c}_{1}^{3}+1152{c}_{1}^{5}).\end{array}\end{eqnarray}$
For η = ∞ , we insert equation (14) into (4) and obtain c−1 = 0 and a recursive relation for cl with l ≥ 2,
$\begin{eqnarray}\begin{array}{l}{c}_{-l+2}\left[{\left(l-2\right)}^{2}-{s}^{2}\right]+{c}_{-l}-\sum _{{k}_{1},{k}_{2},{k}_{3}\geqslant 0}^{{k}_{1}+{k}_{2}+{k}_{3}=l}{c}_{-{k}_{1}}{c}_{-{k}_{2}}{c}_{-{k}_{3}}\\ =0,\,\,\,\mathrm{for}\,\,\,l\geqslant 2.\end{array}\end{eqnarray}$
It can be read out that all coefficients with odd indices are vanishing, the values of coefficients with even indices are solely determined by s. For the case s = 1, we have
$\begin{eqnarray}{c}_{-2}=\displaystyle \frac{1}{2},\,{c}_{-4}=\displaystyle \frac{9}{8},\,{c}_{-6}=\displaystyle \frac{161}{16},\,{c}_{-8}=\displaystyle \frac{24661}{128}.\end{eqnarray}$
In figure 3, we show the values of coefficients cl for s = 1, where c1 takes the value of ${k}_{1}^{\mathrm{ext}}$. The absolute values of nonzero coefficients are shown in a logarithmic scale, the sign of each coefficient is indicated by the color of the point: red/blue stands for the positive/negative value. The magnitude of nonzero cl with l ≤ 0 is fitted by 1.92 × 2l/2Γ[ −l − 1/2] as the green solid line. The absolute values of nonzero cl with l > 0 is fitted by 1.125 × 2.512 1034l as the black solid line. So that we can read out the convergence radius for the series in equation (13) is about 2.5, while the convergence radius for equation (14) is zero if we take 1/η as the expansion variable. It seems that polynomials of Taylor expansion cannot approximate the solution to equation (4). We have to seek a new type of function to incorporate the boundary conditions from both sides and describe the numerical solution in the whole range 0 < η <+ ∞ .
Figure 3. The coefficients cl for s = 1. The absolute values of cl are shown in a logarithmic scale, the sign of cl is indicated by the color of the point: red/blue stands for the positive/negative value. The green solid line stands for 1.92 × 2l/2Γ[ −l − 1/2]. The black solid line stands for 1.125 × 2.512 1034l.
The two-point Padé approximants in terms of rational functions have the form
$\begin{eqnarray}{\hat{g}}_{i,j}(\eta )=\displaystyle \frac{{P}_{i,j}(\eta )}{{Q}_{i,j}(\eta )}=\displaystyle \frac{\sum _{l=0}^{m}{\alpha }_{l}{\eta }^{l}}{\sum _{l=0}^{m}{\beta }_{l}{\eta }^{l}},\,\,\,\mathrm{with}\,\,m=\displaystyle \frac{1}{2}(i+j).\end{eqnarray}$
In this form, Pi,j(η) and Qi,j(η) are polynomials with the truncation order m, and αl and βl are the coefficients determined by requiring that the expansion of ${\hat{g}}_{i,j}(\eta )$ at η → 0+ and 1/η → 0+ agrees with equations (13) and (14) up to and including ci−1ηi−1 and cjηj respectively. So the residues of the approximants are of the order i and j + 1 at the boundaries, i.e.
$\begin{eqnarray}{\hat{g}}_{i,j}(\eta )-g(\eta )={ \mathcal O }({\eta }^{i},{\eta }^{-j-1}).\end{eqnarray}$
With cl in equations (13, 14), the polynomials in the numerator and denominator in equation (19) can be constructed in a determinant representation [36],
$\begin{eqnarray}{Q}_{i,j}(\eta )=\left|\begin{array}{cccc}1 & \eta & \cdots & {\eta }^{m}\\ {c}_{i-1} & {c}_{i-2} & \cdots & {c}_{i-m-1}\\ {c}_{i-2} & {c}_{i-3} & \cdots & {c}_{i-m-2}\\ \vdots & \vdots & \ddots & \vdots \\ {c}_{i-m} & {c}_{i-m-1} & \cdots & {c}_{-j}\end{array}\right|,\end{eqnarray}$
and Pi,j(η) is obtained from Qi,j(η) by replacing the first row with the vector
$\begin{eqnarray}\begin{array}{l}\left\{{S}_{m-1}(\eta ),\,\eta {S}_{m-2}(\eta ),\,\cdots ,\,{\eta }^{m-1}{S}_{0}(\eta ),{\eta }^{m}{T}_{0}(\eta )\right\}\\ \quad \mathrm{for}\,\,\,i\geqslant m,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\left\{{T}_{0}(\eta ),\,\eta {T}_{1}(\eta ),\,\cdots ,\,{\eta }^{m-1}{T}_{m-1}(\eta ),{\eta }^{m}{T}_{m}(\eta )\right\}\\ \quad \mathrm{for}\,\,\,i\leqslant m,\end{array}\end{eqnarray}$
where Sk(η) and Tk(η) are defined as
$\begin{eqnarray}{S}_{k}(\eta )=\displaystyle \frac{{c}_{0}}{2}+\sum _{l=1}^{k}{c}_{l}{\eta }^{l},\end{eqnarray}$
$\begin{eqnarray}{T}_{k}(\eta )=-\displaystyle \frac{{c}_{0}}{2}-\sum _{l=1}^{k}{c}_{-l}{\eta }^{-l}.\end{eqnarray}$
Note that equations (24) and (25) are just the series of equations (13) and (14) up to k respectively.
In such a way, the two-point Padé approximants for the solution to the ODE (4) can be easily constructed. The truncation index (i,j) sets the size and structure of Pi,j and Qi,j, while the numerical value of each matrix element can be determined by the recursive relations in equations (15) and (17) with the connecting parameter cs = ks. But in practice, the expected precision in comparison with the numerical solution is not guaranteed by the constructed function (19), although ks can be determined with extraordinary precision. The quality of the result depends on the choice of (i,j), so we have to check whether the required accuracy can be achieved for particular values of i and j. We regard the high precision numerical solution discussed in the previous section as an ‘exact' solution, from which the accuracy of the approximants can be estimated. For the primary vortex with s = 1, an economic approximate function is given by setting i = 9 and j = 3, which is the best choice for m = (i + j)/2 ≤ 6. The constructed function is parametrized with 12 coefficients in the form
$\begin{eqnarray}\tilde{f}(\eta )=\displaystyle \frac{0.000536403\,{\eta }^{6}+0.00499313\,{\eta }^{5}+0.0279403\,{\eta }^{4}+0.109179\,{\eta }^{3}+0.297503\,{\eta }^{2}+0.583189\,\eta }{0.000536403\,{\eta }^{6}+0.00499313\,{\eta }^{5}+0.0282085\,{\eta }^{4}+0.111676\,{\eta }^{3}+0.312211\,{\eta }^{2}+0.510131\,\eta +1}.\end{eqnarray}$
The Taylor expansion of this fractional polynomial can reproduce c1, ⋯ ,c8 for η → 0 and c−1, ⋯ ,c−3 for $\tfrac{1}{\eta }\to 0$. The accuracy of the function with i = 9 and j = 3, defined by $\max [| \tilde{f}(\eta )-f(\eta )| ]$, is about 1.2 × 10−3, as shown in the left panel of figure 4. It is much better than one general-purpose interpolation routine with limited input data points. The accuracy can be significantly improved for larger m, as shown in the middle panel of figure 4, the accuracy is about 1.5 × 10−8 for the function with i = 26 and j = 10. The coefficients in this function are listed in appendix B. This approximate function with 36 parameters is good enough for the precise simulation of vortices within GPE.
Figure 4. The difference between the approximate functions and the ‘exact' solution for three typical truncation indices.
In the same way, the two-point Padé approximants for the profile functions of vortices with larger winding numbers can be achieved. The coefficients of the fractional polynomials for s = 2, 3 are given in appendix B. In both cases, the accuracy is better than 10−6 in the whole range of η. It is interesting to see that the approximants with even s contain only even terms in Pi,j and Qi,j because the coefficients cl with odd l are all vanishing in equations (13) and (14).
Before we work on even higher order functions, we must keep in mind that the coefficients will all change if we choose different i and j. The reason is that the polynomials in the numerator or denominator result from matrix determinants, any change to the matrix will give very different polynomials. So we cannot expect the coefficients of polynomials will be determined order by order as in a perturbation theory.
To explore the limit of this method, we try to calculate Pi,j(η) and Qi,j(η) for large i and j. The problem is that the matrices of Pi,j(η) and Qi,j(η) involve coefficients cl that differ dramatically in magnitudes as shown in figure 3, it is a numerical challenge to keep the precision in calculating the determinant of such matrices of large dimension. We found Mathematica can handle these matrices with i ≤ 100 and m ≤ 60. The best approximate function that can be accessible corresponds to the matrix for i = 89 and j = 23, the accuracy is of the order 10−14, as shown in the right panel of figure 4. The structure and coefficients of the rational functions may indicate some hints for the exact solution to the nonlinear ODE. In appendix A, we give an example of how the Padé approximants help obtain an exact solution to a nonlinear ODE. But for equation (4), we fail to find any hint for an exact solution except approximate functions at very high precision.

5. Error estimation and challenge

The precision of Padé approximants cannot be improved by simply increasing the truncation order. The asymptotic errors for both boundaries can be estimated as
$\begin{eqnarray}\begin{array}{rcl}\left|g(\eta )-{\hat{g}}_{i,j}(\eta )\right| & = & \left|\displaystyle \frac{{D}_{i-m,m}}{{D}_{i-m-1,m-1}}\right|{\eta }^{i}(1+{ \mathcal O }(\eta ))\\ & & \mathrm{for}\,\,\,\eta \to 0,\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}\left|g(\eta )-{\hat{g}}_{i,j}(\eta )\right| & = & \left|\displaystyle \frac{{D}_{i-m-1,m}}{{D}_{i-m,m-1}}\right|{\eta }^{-j-1}(1+{ \mathcal O }({\eta }^{-1}))\\ & & \mathrm{for}\,\,\,\eta \to \infty ,\end{array}\end{eqnarray}$
with Dr,t being the determinant of the coefficient matrix
$\begin{eqnarray}{D}_{r,t}=\left|\begin{array}{cccc}{c}_{r} & {c}_{r+1} & \cdots & {c}_{r+t}\\ {c}_{r-1} & {c}_{r} & \cdots & {c}_{r+t-1}\\ \vdots & \vdots & \ddots & \vdots \\ {c}_{r-t} & {c}_{r-t+1} & \cdots & {c}_{r}\end{array}\right|.\end{eqnarray}$
The error for small η is controlled by the ratio ∣Dim,m/Dim−1,m−1∣. To have an intuitive idea about the ratio, we estimate it as a function of m for the case i = j = m. As shown in the left panel of figure 5, the ratio can be approximated as 102 × 5.5m. This feature reveals that the convergence radius of ${\hat{g}}_{i,j}$ has been extended from 2.5 to 5.5 by Padé approximants for small η. For large η, the dependence of ∣D−1,m/D0,m−1∣ on m is shown in the right panel of figure 5 and can be described by 10−10 × 18m. This means that the error is under control for η > 18, or the convergence radius for 1/η is extended from 0 to 1/18. Although the change of the convergence radius is small, the coverage range of the approximants is expanded significantly. Actually, this is a remarkable progress of Padé approximants which aim to reproduce a complicated function from both boundaries simultaneously. Inside the gap between two convergence regions, the rational function connects both regions smoothly, which should reproduce the target function with a similar or slightly worse precision if the target function does not contain a singularity in the gap. The error can be further reduced by fine-tuning of i and j if the target function is known with higher precision. As we have already shown in figure 4, the vortex profile function can be reproduced with high accuracy, and the behaviors of errors can be well understood as discussed above.
Figure 5. Left: the ratio ∣D0,m/D−1,m−1∣ (points) varies with m, the red line shows the function 102 × 5.5m. Right: the ratio ∣D−1,m/D0,m−1∣ (points) varies with m, the blue line shows the function 10−10 × 18m.
In the description of Padé approximants for vortex profiles, we find that Pi,j, Qi,j and Dr,t are sensitive to the connecting parameter ks. Since the corresponding matrices have a similar structure, we take Dr,t as an example to discuss its dependence on ks. For simplicity, we look at D0,m and k1. From equation (29), D0,m can be expressed as a polynomial of c1 and the order increases fast with the dimension of the matrix. We find that D0,m has some real roots, denoted as c1i, which are very close to k1. The distance between each c1i and k1 is shown in figure 6, where the magnitude is plotted in logarithm scale and the sign is marked by color. A bunch of roots near k1 indicate that D0,m varies dramatically when c1 approaches k1. So it can be understood that a small deviation from k1 will lead to an uncontrollable variation of Dr,t when the dimension of the matrix becomes large. Similar behavior is also observed in Pi,j and Qi,j as functions of c1. So the accuracy of the Padé approximants for large m strongly depends on the precision of k1. As shown in figure 7, for the best result that we find (the right panel of figure 4), a very small change in the value of k1 leads to a significant increase in the error of the approximants. So the precision of k1 is crucial for that of Padé approximants for large m.
Figure 6. The relative distance between k1 and the roots of D0,m for different m. The magnitude is shown in logarithm scale, while the sign is indicated by color: red points stand for ${c}_{1}^{i}\gt {k}_{1}$, blue points stand for ${c}_{1}^{i}\lt {k}_{1}$. The bottom line indicates the chain of the closest roots, which is technically traceable without knowing k1 at first.
Figure 7. The maximum deviation of the approximants from the ‘exact' numerical solution as a function of the variation of k1.
The minimal distance in figure 6 can be extremely small for large m, it seems that this feature can be employed to find the asymptotic value of k1 from c1i with a carefully chosen initial value and iteration of increasing m. Technically, we can set the only possible root of D0,4 as the initial value, and replace it with one root of D0,5 which is real and closest to the original one. Repeat this procedure, a chain of the closest root can be found, approaching the true value of k1. As shown in figure 6, the bottom line is traceable without knowing k1 at first. This idea is workable for m < 80, with the precision at the level of 10−12. But for even larger m, it is difficult to find the right root of a high-order polynomial, because the density of roots near k1 increases while the minimal distance between them decreases exponentially, which raises the risk of finding wrong roots. This is the reason why we do not seek the value of k1 from roots of Dr,t but import the extrapolated value from a series of numerical solutions with a reliable precision. At the present stage, with Mathematica running on a laptop, we can confirm the value of k1 with 17 significant digits, which helps us to promote the precision of the two-point Padé approximants to 10−14. The method may help us to search for the exact solution of the vortex profile.

6. Summary and conclusion

Understanding the vortex dynamics is important for the study of superfluids and many other physical systems. Due to the nonlinear nature of the vortex dynamics, the evolution of the vortex can only be studied in a perturbative way, which requires a high precision solution to the stationary wave function of the vortex. The nonlinear ODE for the vortex profile function can be solved numerically with high precision.
For practical applications, we have constructed semi-analytical solutions to this problem with two-point Padé approximants. These solutions are presented with simple rational functions, which have a high accuracy and are ready for use in the precise simulation of the vortex dynamics. The coefficients of rational functions strongly depend on the value of the connecting parameter which can be obtained in a high precision from an extrapolation of the numerical solution or a systematic root-finding in a series of polynomials.
The accuracy of the approximate solutions can be improved significantly with a reasonable extension in the order of rational functions. With a systematic scan of the truncation orders, the best accuracy is found to be at the order of 10−14. The errors of the approximate functions and the limitation of two-point Padé approximants are discussed, which can extend our understanding on the nonlinearity of the vortex dynamics. This investigation may provide the clue for an exact solution to the vortex profile function.
The methods and algorithms developed in this work for high precision solutions to the vortex equation can be applied to other nonlinear systems, such as the Schrödinger–Poisson problem [38, 39], the shape and properties of Bose stars formed by dark matter through the universal gravitational interaction and so on.

J D is supported by the Natural Science Foundation of Shandong Province under Grant No. ZR2020MA099. Q W is supported in part by the National Natural Science Foundation of China (NSFC) under Grants No. 12135011, 11890713 (a subgrant of 11890710), and by the Strategic Priority Research Program of the Chinese Academy of Sciences (CAS) under Grant No. XDB34030102.

We consider the nonlinear ODE

$\begin{eqnarray}\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{x}^{2}}f(x)+f(x)-{f}^{3}(x)=0,\end{eqnarray}$
with the boundary conditions f(0) = 0 and f( ∞ ) = 1. At the zero point, we may insert the formal solution
$\begin{eqnarray}f(x)=\sum _{0}^{\infty }{c}_{i}{x}^{i},\end{eqnarray}$
into the ODE. It is easy to check that all even terms are vanishing, i.e. c2k = 0 for k = 0, 1, ⋯ . The odd terms satisfy the recurrence relation
$\begin{eqnarray}{c}_{n+2}(n+2)(n+1)+{c}_{n}-\sum _{{k}_{1}+{k}_{2}+{k}_{3}=n}{c}_{{k}_{1}}{c}_{{k}_{2}}{c}_{{k}_{3}}=0.\end{eqnarray}$
The only unsettled parameter is c1, by which other coefficients can be expressed, for example, we have the following expressions for c3, ⋯ ,c9 in terms of c1
$\begin{eqnarray}\begin{array}{rcl}{c}_{3} & = & -\displaystyle \frac{1}{3!}{c}_{1},\\ {c}_{5} & = & \displaystyle \frac{1}{5!}({c}_{1}+6{c}_{1}^{3}),\\ {c}_{7} & = & -\displaystyle \frac{1}{7!}({c}_{1}+66{c}_{1}^{3}),\\ {c}_{9} & = & \displaystyle \frac{1}{9!}({c}_{1}+612{c}_{1}^{3}+756{c}_{1}^{5}).\end{array}\end{eqnarray}$
We may represent the solution in the form of Padé approximants
$\begin{eqnarray}{\tilde{f}}^{[N]}(x)=\displaystyle \frac{\sum _{i=0}^{N}{a}_{i}{x}^{i}}{\sum _{j=0}^{N}{b}_{j}{x}^{j}}.\end{eqnarray}$
Because the solution is determined by the ratio, it is natural to demand b0 = 1. Because c0 = 0, we can set a0 = 0. Other coefficients should satisfy
$\begin{eqnarray}{a}_{i}=\sum _{j=0}^{N}{b}_{j}{c}_{i-j},\,\,\,\,\,\mathrm{for}\,\,1\leqslant i\leqslant N.\end{eqnarray}$

For the boundary condition at x → ∞ , we can change the variable to z = 1/x, so that equation (A1) can be rewritten as

$\begin{eqnarray}{z}^{4}\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{z}^{2}}f(z)+2{z}^{3}\displaystyle \frac{{\rm{d}}}{{\rm{d}}z}f(z)+f(z)-{f}^{3}(z)=0.\end{eqnarray}$
Then we can take the formal solution in the form
$\begin{eqnarray}f(z)=\sum _{0}^{\infty }{d}_{i}{z}^{i}.\end{eqnarray}$
It is easy to obtain d0 = 1 and di = 0 for i > 0. So there is another constraint for coefficients
$\begin{eqnarray}{a}_{i}=\sum _{j=0}^{N}{b}_{j}{d}_{j-i}={b}_{i},\,\,\,\,\,\mathrm{for}\,\,1\leqslant i\leqslant N.\end{eqnarray}$
Combining equations (A6) and (A9), we obtain a set of linear equations for (a1, ⋯ ,aN) or (b1, ⋯ ,bN). We may write the equations into a matrix form. As an example, the equations for N = 5 read
$\begin{eqnarray}\left(\begin{array}{ccccc}1 & 0 & 0 & 0 & 0\\ -{c}_{1} & 1 & 0 & 0 & 0\\ 0 & -{c}_{1} & 1 & 0 & 0\\ \displaystyle \frac{{c}_{1}}{6} & 0 & -{c}_{1} & 1 & 0\\ 0 & \displaystyle \frac{{c}_{1}}{6} & 0 & -{c}_{1} & 1\end{array}\right)\left(\begin{array}{c}{b}_{1}\\ \\ {b}_{2}\\ \\ {b}_{3}\\ \\ {b}_{4}\\ \\ {b}_{5}\end{array}\right)=\left(\begin{array}{c}{c}_{1}\\ 0\\ -\displaystyle \frac{{c}_{1}}{6}\\ 0\\ \displaystyle \frac{{c}_{1}^{3}}{20}+\displaystyle \frac{{c}_{1}}{120}\end{array}\right).\end{eqnarray}$
Because the matrix is already in a triangle form, the solution can be easily obtained
$\begin{eqnarray}\left(\begin{array}{}\\ \\ {a}_{1}\\ \\ \\ {a}_{2}\\ \\ \\ {a}_{3}\\ \\ \\ {a}_{4}\\ \\ \\ {a}_{5}\\ \end{array}\right)=\left(\begin{array}{}\\ {b}_{1}\\ \\ \\ {b}_{2}\\ \\ \\ {b}_{3}\\ \\ \\ {b}_{4}\\ \\ \\ {b}_{5}\end{array}\right)=\left(\begin{array}{c}{c}_{1}\\ {c}_{1}^{2}\\ {c}_{1}^{3}-\displaystyle \frac{1}{6}{c}_{1}\\ {c}_{1}^{4}-\displaystyle \frac{1}{3}{c}_{1}^{2}\\ {c}_{1}^{5}-\displaystyle \frac{9}{20}{c}_{1}^{3}+\displaystyle \frac{1}{120}{c}_{1}\end{array}\right).\end{eqnarray}$
To determine the value of c1, we may demand one more constraint on the solution by taking aN+1 = bN+1 = 0, which can be written as
$\begin{eqnarray}\sum _{j=1}^{N}{b}_{j}{c}_{N+1-j}+{c}_{N+1}=0.\end{eqnarray}$
It gives an algebra equation for c1. As an example, for N = 5, the equation reads
$\begin{eqnarray}\displaystyle \frac{{c}_{1}^{2}}{90}\left(90{c}_{1}^{4}-51{c}_{1}^{2}+4\right)=0.\end{eqnarray}$
The meaningful root to the above equation is $\sqrt{(51+3\sqrt{123})/180}\approx 0.687\,481$. In this way, we can get a series of c1 values 0.69976, 0.704 416, 0.706143, 0.706768, 0.706989 for N = 6, 7, 8, 9, 10 respectively. We may check that the approximate value of c1 becomes closer to $\sqrt{2}/2\approx 0.707\,107$ for large N. If we set ${c}_{1}=\sqrt{2}/2$, we find the solution in equation (A11) can be simplified to ${a}_{i}={b}_{i}={\left(\sqrt{2}\right)}^{i}/(2\cdot i!)$. So the solution can be put into the form
$\begin{eqnarray}\begin{array}{rcl}{\tilde{f}}^{[N]}(x) & = & \displaystyle \frac{\sum _{i=1}^{N}{\left(\sqrt{2}x\right)}^{i}/(2\cdot i!)}{1+\sum _{i=1}^{N}{\left(\sqrt{2}x\right)}^{i}/(2\cdot i!)}\\ & = & \displaystyle \frac{\sum _{i=0}^{N}{\left(\sqrt{2}x\right)}^{i}/i!-1}{\sum _{i=0}^{N}{\left(\sqrt{2}x\right)}^{i}/i!+1}.\end{array}\end{eqnarray}$
It happens to give $f(x)=\tanh (x/\sqrt{2})$ if we send N to ∞ .

In this appendix we list the coefficients in approximate functions in (19) with which one can reproduce the vortex profiles in high accuracy.

s = 1 s = 3
l αl βl αl βl
0 0 1.000 000 000 0 1.000 000 000
1 5.831 894 959 × 10−1 5.873 264 261 × 10−1 0 4.062 542 828 × 10−1
2 3.425 226 024 × 10−1 5.455 940 613 × 10−1 0 2.100 817 734 × 10−1
3 2.452 860 386 × 10−1 2.468 318 072 × 10−1 2.618 342 072 × 10−2 6.454 678 403 × 10−2
4 1.011 343 919 × 10−1 1.145 178 304 × 10−1 1.063 712 681 × 10−2 2.069 992 436 × 10−2
5 3.831 447 420 × 10−2 3.898 275 190 × 10−2 3.864 195 663 × 10−3 5.236 196 331 × 10−3
6 1.137 856 193 × 10−2 1.179 179 362 × 10−2 1.025 235 177 × 10−3 1.313 346 436 × 10−3
7 2.909 882 040 × 10−3 2.971 398 890 × 10−3 2.391 152 072 × 10−4 2.871 142 311 × 10−4
8 6.323 721 506 × 10−4 6.423 893 504 × 10−4 4.809 356 698 × 10−5 5.537 192 270 × 10−5
9 1.180 693 171 × 10−4 1.194 567 590 × 10−4 8.539 789 317 × 10−6 9.539 350 138 × 10−6
10 1.905 344 120 × 10−5 1.922 168 115 × 10−5 1.358 658 441 × 10−6 1.482 637 465 × 10−6
11 2.672 381 372 × 10−6 2.690 240 172 × 10−6 1.934 240 689 × 10−7 2.074 012 815 × 10−7
12 3.271 123 471 × 10−7 3.287 644 644 × 10−7 2.463 814 537 × 10−8 2.607 377 411 × 10−8
13 3.497 660 717 × 10−8 3.510 863 602 × 10−8 2.842 462 008 × 10−9 2.973 581 152 × 10−9
14 3.254 331 773 × 10−9 3.263 312 251 × 10−9 2.975 340 172 × 10−10 3.082 927 325 × 10−10
15 2.613 493 838 × 10−10 2.618 418 071 × 10−10 2.742 356 229 × 10−11 2.830 506 106 × 10−11
16 1.786 780 526 × 10−11 1.788 474 199 × 10−11 2.276 667 859 × 10−12 2.335 377 575 × 10−12
17 9.848 466 014 × 10−13 9.848 466 014 × 10−13 1.958 886 160 × 10−13 1.958 886 160 × 10−13
18 3.387 345 177 × 10−14 3.387 345 177 × 10−14 1.304 660 342 × 10−14 1.304 660 342 × 10−14

s = 2
l αl βl
0 0 1.000 000 000
2 1.530 991 029 × 10−1 4.398 930 447 × 10−1
4 5.458 897 193 × 10−2 6.936 791 316 × 10−2
6 5.406 591 644 × 10−3 6.042 107 443 × 10−3
8 2.685 749 628 × 10−4 2.859 576 645 × 10−4
10 7.837 336 242 × 10−6 8.153 179 530 × 10−6
12 1.472 555 945 × 10−7 1.513 096 732 × 10−7
14 1.934 727 229 × 10−9 1.970 129 182 × 10−9
16 1.693 952 735 × 10−11 1.724 410 709 × 10−11
18 1.522 898 667 × 10−13 1.522 898 667 × 10−13

1
Pitaevskii L P Stringari S 2003 Bose–Einstein Condensation Oxford Oxford University Press

2
Fetter A L 2009 Rotating trapped Bose–Einstein condensates Rev. Mod. Phys. 81 647 691

DOI

3
Sonin E 2016 Dynamics of Quantised Vortices in Superfluids Cambridge Cambridge University Press

4
Matthews M R 1999 Vortices in a Bose–Einstein Condensate Phys. Rev. Lett. 83 2498 2501

DOI

5
Leanhardt A E 2002 Imprinting vortices in a Bose–Einstein condensate using topological phases Phys. Rev. Lett. 89 190403

DOI

6
Inouye S 2001 Observation of vortex phase singularities in Bose–Einstein condensates Phys. Rev. Lett. 87 080402

DOI

7
Abo-Shaeer J R Raman C Vogels J M Ketterle W 2001 Observation of vortex lattices in Bose–Einstein condensates Science 292 476 479

DOI

8
Scherer D R Weiler C N Neely T W Anderson B P 2007 Vortex formation by merging of multiple trapped Bose–Einstein condensates Phys. Rev. Lett. 98 110402

DOI

9
Carretero-Gonzalez R 2008 Dynamics of vortex formation in merging Bose–Einstein condensate fragments Phys. Rev. A 77 033625

DOI

10
Deng J Wang Q Zhang H Vortex formation in collisions of Bose–Einstein condensates arXiv:2202.13300 [nucl-th]

11
Berezhiani L Khoury J 2016 Dark matter superfluidity and galactic dynamics Phys. Lett. B 753 639 643

DOI

12
Berezhiani L Khoury J 2015 Theory of dark matter superfluidity Phys. Rev. D 92 103510

DOI

13
Hui L Joyce A Landry M J Li X 2021 Vortices and waves in light dark matter JCAP 01 011

DOI

14
Kibble T W B 1976 Topology of cosmic domains and strings J. Phys. A 9 1387 1398

DOI

15
Zurek W H 1985 Cosmological experiments in superfluid helium? Nature 317 505 508

DOI

16
Hindmarsh M B Kibble T W B 1995 Cosmic strings Rep. Prog. Phys. 58 477 562

DOI

17
Zurek W H 1996 Cosmological experiments in condensed matter systems Phys. Rep. 276 177 221

DOI

18
Garay L J Anglin J R Cirac J I Zoller P 2000 Black holes in Bose–Einstein condensates Phys. Rev. Lett. 85 4643 4647

DOI

19
Solnyshkov D D Leblanc C Koniakhin S V Bleu O Malpuech G 2019 Quantum analogue of a Kerr black hole and the Penrose effect in a Bose–Einstein condensate Phys. Rev. B 99 214511

DOI

20
Berges J Rothkopf A Schmidt J 2008 Non-thermal fixed points: effective weak-coupling for strongly correlated systems far from equilibrium Phys. Rev. Lett. 101 041603

DOI

21
Berges J Hoffmeister G 2009 Nonthermal fixed points and the functional renormalization group Nucl. Phys. B 813 383 407

DOI

22
Berges J Boguslavski K Schlichting S Venugopalan R 2015 Universality far from equilibrium: from superfluid Bose gases to heavy-ion collisions Phys. Rev. Lett. 114 061601

DOI

23
Nowak B Schole J Sexty D Gasenzer T 2012 Nonthermal fixed points, vortex statistics, and superfluid turbulence in an ultracold Bose gas Phys. Rev. A 85 043627

DOI

24
Schole J Nowak B Gasenzer T 2012 Critical dynamics of a two-dimensional superfluid near a non-thermal fixed point Phys. Rev. A 86 013624

DOI

25
Karl M Gasenzer T 2017 Strongly anomalous non-thermal fixed point in a quenched two-dimensional Bose gas New J. Phys. 19 093014

DOI

26
Deng J Schlichting S Venugopalan R Wang Q 2018 Off-equilibrium infrared structure of self-interacting scalar fields: universal scaling, vortex-antivortex superfluid dynamics and Bose–Einstein condensation Phys. Rev. A 97 053606

DOI

27
Gross E P 1963 Hydrodynamics of a superfluid condensate J. Math. Phys 4 195 207

DOI

28
Pitaevskii P P 1961 Vortex lines in an imperfect bose gas J. Exp. Theor. Phys. 13 451 454

29
Frantzeskakis D J Carretero-González R 2008 Emergent Nonlinear Phenomena in Bose–Einstein Condensates: Theory and Experiment[M] Berlin Springer

30
Liu W M Kengne E 2019 Schrödinger Equations in Nonlinear Systems 1st edn Singapore Springer

31
Svidzinsky A A Fetter A L 2000 Dynamics of a vortex in a trapped Bose–Einstein condensate Phys. Rev. A 62 063617

DOI

32
Koens L Martin A M 2012 Perturbative behavior of a vortex in a trapped Bose–Einstein condensate Phys. Rev. A 86 013605

DOI

33
Sun B D Sun Z H Zhou J 2021 Trace anomaly contribution to hydrogen atom mass Phys. Rev. D 104 056008

DOI

34
Tian Y H Li W D Shen Y Dai W S 2021 Padé approximant approach to singular properties of quantum gases: the ideal cases Commun. Theor. Phys. 73 065602

DOI

35
Mccabe J H Murphy J A 1976 Continued fractions which correspond to power series expansions at two points IMA J. Appl. Math. 17 233

DOI

36
Sidi A 1980 Some aspects of two-point Padé approximants J. Comput. Appl. Math. 6 9

DOI

37
Boisseau B Forgacs P Giacomini H 2007 An analytical approximation scheme to two point boundary value problems of ordinary differential equations J. Phys. A 40 215

DOI

38
Braaten E Mohapatra A Zhang H 2016 Dense axion stars Phys. Rev. Lett. 117 121801

DOI

39
Levkov D G Panin A G Tkachev I I 2018 Gravitational Bose–Einstein condensation in the kinetic regime Phys. Rev. Lett. 121 151301

DOI

Outlines

/