Numerical model of the Gross–Pitaevskii equation for rotating Bose–Einstein condensates using smoothed-particle hydrodynamics

Vortex lattice phenomena observed in rotating Bose–Einstein condensates (BECs) have garnered attention in many fields related to condensed matter physics. Clarifying the essential mechanism of lattice formation in BECs can significantly contribute to these disciplines. The numerical analyses of the lattice formation process in rotating superfluid helium-4 or ultracold atomic gases have been major topics in this field. Several computational schemes for the direct simulation of the Gross–Pitaevskii (GP) equation, a nonlinear Schrödinger equation for interacting bosons, have been developed over the past decades. Specifically, mesh-based approaches, such as the finite difference (FD)^1–3 or finite element (FE)^4–6 methods, have frequently been adopted for the spatial discretization of the GP equations. Employing a high-order time-integrating scheme is essential to ensure numerical stability for the iterative computation over a long period; hence, the implicit time-integration schemes using the Crank–Nicolson^2,7,8 or alternating direction implicit (ADI)^9,10 have been employed. Because the implicit approaches result in matrix computations for solving eigenvalue problems, they introduce the unit of the imaginary time and algebraically update the state of the system. Alternative methods for solving the eigenvalue problem for the time-dependent Schrödinger equation include the imaginary-time propagation methods^11–15 and the Fourier spectra method.¹⁶ Till date, these numerical approaches have achieved a certain success in capturing the complex dynamics of lattice formation in BECs.

Because the GP equation is a nonlinear Schrödinger equation for interacting bosons, its Hamiltonian has the form of a many-particle interacting system.^17–19 By describing the GP equation in numerical simulations as a many-particle interaction system similar to the picture of BEC physics, one can expect a more accurate reproduction of microscopic interactions among atomic particles using analytical particles in the computational domain. However, this is challenging with Eulerian models, such as FE or FD, because discretized points have no physical meaning as particles. Therefore, Lagrangian models have an advantage over these Eulerian models. Unfortunately, to the best of our knowledge, no study has reported on the Lagrangian particle approximation to discretize the GP equation to reproduce quantum lattice. The vortex filament method (VFM) is employed for local interaction calculations among several quantum vortices on the basis of the Biot–Savart law;^20–24 however, currently, it is not suitable for the collective dynamics simulations of the condensates because of the difference in physical scales. In contrast, smoothed particle hydrodynamics (SPH),²⁵ which originated in astrophysics, is used for fluid dynamics calculations and for general purposes, such as a finite particle approximation for the partial differential equation.^26–28 SPH can be expected to work as a Lagrangian particle numerical scheme for the direct simulations of the GP equation. In addition, a recent study²⁹ suggested that a two-fluid model for superfluid helium-4 can be solved using SPH to produce vortex lattices under specific conditions. Another study¹⁹ reported a theoretical proposal that the motion equation for inviscid fluid in the two-fluid model becomes equal to the quantum fluid equation derived from the GP equation under specific conditions, in the case of SPH formalism. Therefore, simulating the behavior of the GP equation in SPH form may contribute to finding certain connection to these studies.

This study proposed a new numerical scheme for vortex lattice formation in a rotating BEC using SPH with an explicit real-time integration scheme. Specifically, the Gross–Pitaevskii (GP) equation was proposed in complex representation to obtain a pair of time-dependent equations, which were obtained in the real and imaginary parts, respectively. These equations were then discretized and solved based on the SPH particle approximation. We adopted the fourth-order Runge–Kutta method for time evolution; our scheme has the accuracy of the second order in space and the fourth order in time evolution. The results of SPH simulations of a rotating Bose gas trapped in a harmonic potential were qualitatively consistent with previously reported experiments and simulations. Notably, we succeeded in reproducing the geometric patterns of formed lattices for several cases, for example, the hexagonal lattice observed in the experiments of rotating BECs. We also report that the simulation began with the periodic oscillation of the condensate, which attenuated and maintained a stable rotation with slanted elliptical shapes; however, the surface was excited to be unstable and generated ripples, which grew into vortices and then penetrated inside the condensate, forming a lattice. We confirmed that each branch point of the phase of wavefunctions corresponded to each vortex. Meanwhile, the Lagrangian approaches, such as SPH, incur huge computational costs. Thus, we developed an explicit real-time integration scheme for SPH to be compatible with computational accelerators, such as graphics processing units (GPU). Consequently, this study demonstrated that the Lagrangian numerical models with explicit time integration can reproduce the dynamics of vortex-lattice formation by selecting the appropriate scheme.

The remainder of this paper is structured as follows: In Sec. II, we explain the targeted governing equation and its discretized form in SPH. We present the several numerical tests for vortex lattice formation in rotating BECs in Sec. III and discuss the simulation results in Sec. IV. Finally, Sec. V concludes this paper.

In the first equation of Eq. (2), we name the first and second terms as the rotational and oscillatory terms, respectively. We refer to $( C | ψ | 2 + V − μ )$ in the parentheses of the oscillatory term as the potential term. The same terminology is applied to the second equation of Eq. (2). In the repeated process of the first equation of Eq. (2), we first update the potential term using the positions and the current values of u and w. Next, we calculate the oscillation term by summing the product of the potential term with w and the Laplacian of w. Here, the Laplacian of w has been obtained by substituting w for $ϕ$ in Eq. (7). We also calculate the rotational term by using the position and gradient of u, which was obtained by substituting u for $ϕ$ in Eq. (6), where $∂ a ϕ$ represents the a-component of the gradient $∇ ϕ$⁠. We update the variable X after obtaining the right-hand side, the difference between the rotational and oscillatory terms. We update the variable Y through a similar process for the second equation of Eq. (2). We then obtain the updated u and w from X and Y using the third and fourth equations of Eq. (2). We iteratively repeat this process in an explicit time-integrating scheme using the fourth-order Runge–Kutta method. Then, the values of u and w are renormalized to satisfy Eq. (3).

To ensure the numerical accuracy of the Laplacian calculation, Eq. (7) can be rewritten to the form of the Laplacian operator of the moving particle semi-implicit (MPS). This is a particle method where the gradient at a point is estimated by the weighted arithmetic mean of the local gradients defined between the point and other neighboring points in space.^34,35 In brief, MPS is one form of the finite difference method (FDM) generalized in the Lagrangian form.²⁷ For the detailed explanation of the reformulation of the Laplacian operator, refer to the appendix in Ref. 36. In simulations, we applied Eqs. (6) and (7) to the first term and the first part in the square bracket of the second term in Eq. (2) after the reformulation. In contrast to the SPH scheme for fluid dynamics, other stabilization techniques or improved SPH operators are not required to compute Eq. (2).

Numerical simulations of the vortex lattice formations were performed using our SPH scheme. The parameters $( γ , C , N , λ , a , a h , ω ⊥ )$ introduced in Sec. II A were set as $( 0.03 , 500 , 3 × 10 5 , 9.2 , 5.77 , 0.728 , 2 π × 108.56 )$⁠, respectively, to replicate the previous numerical tests based on the FDM-based scheme reported in Ref. 30. Here, the units of the last three parameters a, a_h, and $ω ⊥$ are $nm , μ m$⁠, and Hz, respectively. We estimated the chemical potential μ as $μ ≈ ( 15 / 8 N λ a / a h ) 2 / 5$ based on the Thomas–Fermi (TF) model. For the anisotropy parameters ϵ_x and ϵ_y, we set ϵ_y to be zero during simulations while rapidly increasing ϵ_x from zero to 0.051 282 051 28 in 20 ms at the beginning of the simulation; this condition corresponds to the case that the entire anisotropy $ϵ : = ( 1 + ϵ x ) − ( 1 + ϵ y ) / ( 1 + ϵ x ) + ( 1 + ϵ y )$ changes from zero to 0.025 as per the numerical tests in Ref. 30 and experiments by Madison et al.³⁷ Simulations were performed for different values of Ω: 0.57, 0.7, 0.86, and 0.9 to compare the resulting lattices.

For computational conditions, we set the simulation domain (L_x, L_y) to be (40, 40). An initial profile was assigned to the real component of the wavefunction according to the profile given by TF approximation around the z-axis with a radius of 8.855 817 453 25, which was determined from chemical potential μ. The resolutions in the x and y directions (N_x, N_y) were set as (380, 380). The same SPH parameters optimized in the numerical tests for a wave propagation problem reported in Sec. 5.2 in Ref. 36 were used; the kernel radius h was $0.4 d$⁠, where d is the distance between the two nearest discrete points arranged along the orthogonal grids of the two-dimensional space in the initial state. We adopted the Gaussian kernel function. Here, the Gaussian kernel does not exhibit the compact support property, where the kernel exhibits a value of zero at a finite distance from its center point. Thus, a cutoff distance with a radius of 3h was introduced. In addition, we imposed the Dirichlet boundary condition $u = w = 0$ on the edge of the simulation domain. As a boundary treatment for the rotational and oscillatory forces, we introduced another cutoff distance d_c at around the outer area of the simulation domain; the forces computed at the distance of d_c or more significant from the center point were attenuated using a Gaussian filter to be smoothly connected to the edge of the domain where the value of wavefunction was always zero. Thus, we set the value of d_c to 15 to ensure the sufficient space for the condensates.

Equation (1) assumes the existence of one particle in the ground state in the z-direction, and, therefore, the z-component of the wavefunction can be described by Gaussian distribution.³⁰ Thus, Eq. (1) postulates a quasi-2D situation. We can replicate this situation in our SPH model using a 3D kernel function. As the cutoff length of the kernel function was set to 3h and $h = 0.4 d$⁠, the kernel function contained three layers of discrete points in one direction. Thus, we placed the positive and negative layers parallel to the xy-plane at z = 0 such that the kernel function could hold three layers in the z-direction. We then projected their contributions onto the xy-plane at z = 0. These layers can be virtually set because the 2D and 3D Gaussian kernels only differ in the scale of $C sph$⁠, the normalization coefficient of the SPH kernel (Sec. II B). The theoretical value of $C sph$ is obtained as $( π h 2 ) − 1$ for 2D and $( π 3 / 2 h 3 ) − 1$ for 3D calculations.²⁵ However, because actual calculations include discretization and summation errors, the direct use of these analytical values does not guarantee sufficient accuracy of normalization. We calibrated $C sph$ by referring to the correspondence table between the number of vortices and Ω reported in Ref. 30. Specifically, we assumed the linear relationship between the value of Ω and the number of vortices between $Ω = 0.57$ and $Ω = 0.9$⁠. Subsequently, $C sph$ was estimated by referring to the results for $Ω = 0.9$⁠. We then used this value in the remaining three cases of $Ω = 0.57 , Ω = 0.7$⁠, and $Ω = 0.86$⁠. We set the discrete time dt to be $1.0 × 10 − 4 ω ⊥ − 1$ and computed $1000 ω ⊥ − 1$ in total; this is approximately equal to 1466 ms on the original timescale. We implemented our scheme on a single GPU, NVIDIA GeForce RTX2080 Ti, using C/C++ and CUDA.³⁸ We utilized several computational techniques to speed up the simulations; neighbor-particle lists are adopted to reduce the computational cost of finding particles within the radius of interactions from $O ( N 2 )$ to $O ( N )$⁠. We combined the linked list technique^39,40 with the neighbor-particle list to reduce memory usage. Although we developed our original calculation code, several open-source frameworks should serve as references for implementation.^40,41

Figure 1 shows the snapshots of an SPH simulation of the vortex lattice formation obtained by solving the GP equation in Eq. (2) for $Ω = 0.7$⁠. The top row in Fig. 1 displays the density profiles $| ψ | 2$ approximately at (a) 0.5, (b) 39, (c) 244, (d) 391, and (e) 977 ms. The lower rows (f)–(j) show the phase profiles obtained as $θ = atan ( u / w )$ for the corresponding density profiles in (a)–(e), respectively. A video of the snapshots in Fig. 1 is provided as integral multimedia in Fig. 2 (Multimedia view) (ancillary files for the preprint version). Figure 3 shows the snapshots of vortex lattice formations under the same conditions as in Fig. 1, for different cases of (a) $Ω = 0.57$⁠, (b) $Ω = 0.7$⁠, (c) $Ω = 0.86$⁠, and (d) $Ω = 0.9$⁠, after a sufficient time longer than 1400 ms. The lower panels (e)–(h) show the phase profiles for the corresponding density profiles in (a)–(d), respectively. Figure 4 shows the simulation results for $Ω = 0.706 25$ under the same conditions as in Fig. 3, evidently demonstrating that we could reproduce a hexagonal lattice, which is a typical shape for vortex lattices in rotating BECs.

The following observations are made in Fig. 1. The simulation began with the periodic oscillation of condensate, which gradually attenuated, and the condensate maintained a stable rotation with slanted elliptical shapes. However, concurrently, the surface of condensates was gradually excited to be unstable and generated ripples, which grew into vortices and then penetrated into the inside condensates, eventually forming a lattice. Notably, we observed a regular single pentagon after a sufficient time longer than 977 ms. It was confirmed that the geometric patterns of lattices for all the comparable cases of $Ω = 0.57 , Ω = 0.7$⁠, and $Ω = 0.86$ in Fig. 3 were consistent with the numerical tests in Ref. 30, respectively. In addition, the increase in the number of vortices with an increase in Ω was observed, as shown in Fig. 3. Figure 5(a) shows Ω-dependence of the number of vortices for the simulations displayed in Fig. 3; the dependence was quantitatively consistent with the numerical results for all the reported cases of $Ω = 0.57 , Ω = 0.7$⁠, and $Ω = 0.86$ in Ref. 30. Furthermore, the phase patterns were confirmed to be well reproduced in all the cases. Here, u and w were the real and imaginary parts of the wavefunction, and hence, the phase discontinuities in the lower rows of Figs. 1 and 3 represent the branch cuts in the complex plane. Through comparisons of the geometrical locations of the quantum vortices in the density profiles presented in the upper rows, it was confirmed that these branch points corresponded to the locations of vortices. The shape of the hexagonal lattice in Fig. 4 is confirmed to agree with that of lattices reported in the related work.^37,42 In summary, the key observations specific to the formation of vortex lattices were consistent with numerical^30,31,42 and experimental results,³⁷ demonstrating that our approach exhibited a certain degree of accuracy.

Nevertheless, it is still significant that we could reproduce the dynamics of vortex lattice formation using the explicit time-integrating scheme and SPH, which is a Lagrangian calculation scheme. As mentioned in Sec. I, the GP equation is a nonlinear Schrödinger equation for interacting bosons, and its Hamiltonian has the form of a many-particle interacting system. Now that we have succeeded in describing the GP equation in SPH, which has a form of many-particle interaction systems similar to the quantum mechanical picture, we expect to accurately model the microscopic interactions among atomic particles in future work. In addition, the explicit time-integration scheme is frequently preconceived as less accurate than the implicit time-integration one. Similarly, the Lagrangian calculation schemes are believed to be less accurate than the Euler approaches, for example, the finite-difference or finite-element methods. However, this study demonstrated that the Lagrangian numerical models with explicit time integration could reproduce the dynamics of vortex-lattice formation by selecting the appropriate scheme.

Furthermore, our time-integrating scheme, which describes the GP equation in complex representation and simultaneously solves a pair of two time-dependent equations obtained for the respective parts of wavefunctions, can be intuitively understood in contrast to the case of conventional methods owing to the real-time integration. As it does not require a special mathematical operation system to perform algebraic operations involving imaginary units, its implementation on a computational accelerator, such as a GPU, is straightforward. In fact, all calculations in this study were performed on an NVIDIA GeForce RTX2080 Ti. Thus, as our SPH scheme is a computational scheme that is compatible with high-performance computing, we can expect to explore the nonlinear dynamics of vortex lattice formation using HPC resources like GPU-rich supercomputers in the near future.

This study proposed a static SPH scheme wherein all computation points were fixed in space. This offers many advantages in view of computational physics, such as the practical ability to place calculation points anywhere in space. However, it is significant from the perspective of condensed matter physics that we reproduced vortex lattices in SPH form. A recent study suggested that a two-fluid model for superfluid helium-4 can be solved using SPH to produce vortex lattices under specific conditions; their simulations showed that a gathering of low-density components became a spinning vortex, and they form a vortex lattice eventually.²⁹ In contrast, this study confirmed that the branch points of the phase of wavefunctions, where density cannot be defined owing to singularities, corresponded to the holes in the vortex lattice using SPH simulations. Therefore, the results in this study and Ref. 29 are similar in that both confirmed in SPH formalism that the vortices correspond to the points where densities are vacant. Accordingly, SPH-based approaches can be helpful in determining the connectivity between these phenomena, which was observed on micro- and macroscopic scales. In reality, another study theoretically presents that the motion equation for inviscid fluid in the two-fluid model becomes equal to the quantum fluid equations derived from the GP equations, in the SPH form, given that fluid exhibits moderate density change and has sufficiently slow fluid speed compared to the critical speed of superfluid helium-4.¹⁹ Such connectivity between classical and quantum mechanics can be explored by developing a dynamic SPH model for the GP equation. By allowing the computational points to flow in our SPH scheme proposed in this study, we may expect to simultaneously solve the quantum fluid equations, that is, the macroscopic motion equations derived from the GP equations and the original GP equations in the SPH form. The higher-order scheme would be necessary to this end; nevertheless, our SPH scheme has great potential to advance both physical and computational studies on physics revolving around the dynamics of vortex lattice formation.

This study proposed a new numerical scheme for vortex lattice formation in a rotating BEC using smoothed particle hydrodynamics (SPH) with an explicit real-time integration scheme. Specifically, we described the Gross–Pitaevskii (GP) equation in complex representation to obtain a pair of time-dependent equations, which were solved simultaneously after discretization based on SPH particle approximation. We adopted the fourth-order Runge–Kutta method for the update; our numerical scheme has the accuracy of the second order in space and the fourth order in time. We then performed the simulations of rotating Bose gas trapped in a harmonic potential. Consequently, our simulations are qualitatively consistent with both experiments and simulations reported in the previous studies regarding the phenomena characteristic of quantum lattices: the dynamic processes to form lattices, and the relationship between the angular frequency and the properties of the number of vortices and the geometric pattern of formed lattices. Notably, we succeeded in reproducing the geometric patterns of formed lattices for several cases; for example, the hexagonal lattice observed in the experiments of rotating BECs. We confirmed that the simulation began with the periodic oscillation of the condensate, which attenuated and maintained a stable rotation with slanted elliptical shapes. However, the surface was excited to be unstable and generated ripples, which grew into vortices and then penetrated the inside the condensate, eventually forming a lattice. We confirmed that each branch point of the phase of wavefunctions corresponds to each vortex. The same observations were reported in experiments and the simulations using the higher-order finite difference scheme, demonstrating that our approach exhibited a certain degree of accuracy. In conclusion, we reproduced the dynamics of vortex lattice formation using SPH with the explicit time-integrating scheme and real-time update; we successfully developed a new Lagrangian method for the simulations of vortex lattice formation in rotating BECs.

This study was supported by JSPS KAKENHI (Grant No. 22K14177). The author would like to thank Editage (www.editage.jp) for English language editing. The author, in particular, acknowledges Professor Katsuhiro Nishinari and the administrative staff at the Nishinari Laboratory and the RCAST, University of Tokyo. The author is also grateful to his family for their warm encouragement.

The authors have no conflicts to disclose.

Satori Tsuzuki: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.