This study proposed a new numerical scheme for vortex lattice formation in a rotating Bose–Einstein condensate (BEC) using smoothed particle hydrodynamics (SPH) with an explicit real-time integration scheme. Specifically, the Gross–Pitaevskii equation was described as a complex representation to obtain a pair of time-dependent equations, which were then solved simultaneously following discretization based on SPH particle approximation. We adopt the fourth-order Runge–Kutta method for time evolution. We performed simulations of a rotating Bose gas trapped in a harmonic potential, showing results that qualitatively agreed with previously reported experiments and simulations. The geometric patterns of formed lattices were successfully reproduced for several cases, for example, the hexagonal lattice observed in the experiments of rotating BECs. Consequently, it was 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 inside the condensate, forming a lattice. We confirmed that each branch point of the phase of wavefunctions corresponds to each vortex. These results demonstrate our approach at a certain degree of accuracy. In conclusion, we successfully developed a new SPH scheme for the simulations of vortex lattice formation in rotating BECs.

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–Nicolson2,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 methods11–15 and the Fourier spectra method.16 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),25 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 study29 suggested that a two-fluid model for superfluid helium-4 can be solved using SPH to produce vortex lattices under specific conditions. Another study19 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.

The dynamics of a rapidly rotating Bose gas trapped in a harmonic potential can be described by the following two-dimensional (2D) GP equation on condition that the z (vertical) component of the condensed wavefunction can be separated from the other components as follows:30,
( i γ ) ψ t = [ 2 + 1 4 { ( 1 + ϵ x ) x 2 + ( 1 + ϵ y ) y 2 } + C | ψ | 2 μ + i Ω ( x y y x ) ] ψ ,
(1)
where ϵx and ϵy indicate the anisotropy parameters of the harmonic potential, μ is the chemical potential, and γ is the phenomenological dissipation.31 Let us denote the angular frequency in the x or y (horizontal) direction as ω and that in the z (vertical) direction as ωz. When the ratio of ω to ωz, λ, is sufficiently larger than one ( λ = ω / ω z 1), the parameter C is expressed as C = 4 π λ N a / a h, where N is the number of bosonic particles projected to the two-dimensional plate, a is the scattering length, and ah represents the unit of length. Equation (1) is the dimensionless version of the original GP equation; the length, time, and wavefunction are scaled as x ¯ = a h x , t ¯ = ω 1 t, and ψ ¯ = N a h 1 ψ, respectively. Here, x ¯ , t ¯, and ψ ¯ are the variables of the original GP equation. Moreover, Ω in Eq. (1) indicates the scalar multiple of ω .
Let us describe Eq. (1) in the complex representation as a preparation for the SPH calculation with explicit time-integrating scheme. By substituting the complex representation of wavefunction ψ, ψ = u + i w, in Eq. (1) and comparing each part of the complex expression on the left and right sides, the following simultaneous equations are obtained after simple calculations:
X t = Ω ( x y u y x u ) [ 2 w ( C | ψ | 2 + V μ ) w ] , Y t = Ω ( x y w y x w ) + [ 2 u ( C | ψ | 2 + V μ ) u ] , X : = u γ w , Y : = γ u + w , V : = 1 4 { ( 1 + ϵ x ) x 2 + ( 1 + ϵ y ) y 2 } .
(2)
Here, | ψ | 2 = u 2 + w 2, where ψ must satisfy the following normalization conditions all the times:
| ψ | 2 d xdy = 1.
(3)
Equation (3) represents the law of quantum mechanics, which states that the integral of the square of wavefunction over the entire spatial domain is always one.32 To perform the simulations, we adopted the fourth-order Runge–Kutta method for the time integration of variables X and Y in Eq. (2). In addition, we normalized the wavefunctions according to Eq. (3) at every time step of the simulation. The SPH representation of the operators on the right-hand side of Eq. (2) is explained in Sec. II B.
The SPH was first described in astrophysics25,33 and is now widely used as a finite particle approximation method for continuous distributions. It exploits the fact that a discrete physical quantity φ can be expressed as a continuous quantity using the Dirac delta function δ [ φ ( r ) = φ ( r ) δ ( r r ´ ) d r ´]. In SPH, the δ function is approximated using a distribution function W, referred to as the kernel function as
φ ( r ) φ ( r ) W ( r r ´ , h ) d r ´ .
(4)
Here, W exhibits the properties of smoothness, symmetry around an axis [ W ( r ) = W ( r )], and normalization ( W d r = 1). Moreover, it converges to the δ function when the kernel radius h, or distribution width, approaches zero [ lim h 0 W ( r r ´ , h ) = δ ( r r ´ )]. A straightforward example of W satisfying the aforementioned conditions is the Gaussian kernel25 represented as W ( r r ´ ) = C sph / h d exp [ | r r ´ | 2 / h 2 ], where C sph denotes a normalization constant, h denotes the kernel radius, and d represents the dimension.
Equation (4) can be written in the discretized form considering summation approximation as
φ ( r i ) = j = 1 N p φ ( r j ) Δ V j W i j ,
(5)
where W i j = W ( | r i r j | , h ), and Δ V j represents the discretized small volume. After operating to Eq. (5) from the left, a simple calculation using vector analysis and Gauss' divergence theorem yields the gradient and Laplacian of φ as follows:33,
ϕ ( r i ) = j N p Δ V j [ ϕ ( r i ) + ϕ ( r j ) ] W i j ,
(6)
2 ϕ ( r i ) = j N p Δ V j ϕ ( r i ) ϕ ( r j ) | r i r j | 2 ( r i r j ) · W i j ,
(7)
where W i j denote the gradient of Wij. Equations (6) and (7) were derived assuming that each fragment has the same uniform volume from the standard forms in Ref. 33, wherein the volumes of the ith and jth fragments are expressed separately. This simplification is valid for this study because all the fragments, which are the discrete points of volume, are fixed at the same locations as in the initial state during simulations owing to that the governing equation in Eq. (1) has no advection term. In other words, this study presents a static model that describes the system in SPH form and simultaneously allows us to ensure numerical accuracy.

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.27 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, ah, 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.37 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 (Lx, Ly) 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 (Nx, Ny) 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 dc at around the outer area of the simulation domain; the forces computed at the distance of dc 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 dc 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.30 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.25 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.38 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 technique39,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.

FIG. 1.

The snapshots of the vortex lattice formation obtained by solving the GP equation in Eq. (2) for Ω = 0.7. The top row shows 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 given by θ = atan ( u / w ) for the corresponding density profiles in (a)–(e), respectively.

FIG. 1.

The snapshots of the vortex lattice formation obtained by solving the GP equation in Eq. (2) for Ω = 0.7. The top row shows 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 given by θ = atan ( u / w ) for the corresponding density profiles in (a)–(e), respectively.

Close modal
FIG. 2.

A video of (A) the phase and (B) density profiles of the snapshot in Fig. 1. Multimedia view: https://doi.org/10.1063/5.0143556.1

FIG. 2.

A video of (A) the phase and (B) density profiles of the snapshot in Fig. 1. Multimedia view: https://doi.org/10.1063/5.0143556.1

Close modal
FIG. 3.

The snapshots of vortex lattice formation for the same conditions as in Fig. 1, for (a) Ω = 0.57, (b) Ω = 0.7, (c) Ω = 0.86, and (d) Ω = 0.9, for 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.

FIG. 3.

The snapshots of vortex lattice formation for the same conditions as in Fig. 1, for (a) Ω = 0.57, (b) Ω = 0.7, (c) Ω = 0.86, and (d) Ω = 0.9, for 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.

Close modal
FIG. 4.

The snapshots of (A) the phase and (B) density profiles of a hexagonal lattice for Ω = 0.706 25.

FIG. 4.

The snapshots of (A) the phase and (B) density profiles of a hexagonal lattice for Ω = 0.706 25.

Close modal

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 numerical30,31,42 and experimental results,37 demonstrating that our approach exhibited a certain degree of accuracy.

FIG. 5.

Detailed breakdown of the simulations presented in Fig. 3: (a) Ω-dependence of the number of vortices and (b) time dependence of the deformation parameter α for different cases Ω = 0.57 , Ω = 0.7 , Ω = 0.86, and Ω = 0.9.

FIG. 5.

Detailed breakdown of the simulations presented in Fig. 3: (a) Ω-dependence of the number of vortices and (b) time dependence of the deformation parameter α for different cases Ω = 0.57 , Ω = 0.7 , Ω = 0.86, and Ω = 0.9.

Close modal
Our simulations were qualitatively consistent with both experiments and simulations in the previous studies30,31,37,42 in terms of the phenomena characteristic of quantum lattices: the dynamic processes to form lattices, and the relationship between the angular frequency and the properties of number of vortices and the geometric pattern of formed lattices. In addition, similar to these previous studies, we measured the deformation parameter α, defined below, to more precisely discuss the morphological aspects of vortex lattice formation,
α : = Ω x 2 y 2 x 2 + y 2 ,
(8)
where the bracket symbol ξ represents the expected value of the physical quantity ξ, which can be calculated as ξ = ξ | ψ | 2 d xdy if ξ is a scalar quantity. Figure 5(b) shows the time dependence of the parameter α for Ω = 0.57 , Ω = 0.7 , Ω = 0.86, and Ω = 0.9, respectively. In all cases of Ω, we observed that parameter α gets stabilized after experiencing oscillations for a certain time. Similar behaviors were observed in the experiments37 and the simulations that utilize the higher-order finite difference method.30 The parameter α was confirmed to drop abruptly to a value below 0.05 as the vortices enter the condensates when Ω = 0.7, similar to that reported in Ref. 30. However, an in-depth examination revealed several differences compared to the simulation reported in Ref. 30. In particular, the physical time elapsed for completing the vortex lattice formation process was approximately 735 ms for Ω = 0.7 in Ref. 30, whereas it was approximately 489 ms in our simulations. Thus, the vortex lattice formation process proceeded approximately 1.5 times faster in our SPH simulation. In addition, the physical time for the condensates to form an ellipsoid shape when Ω = 0.7 was approximately 39 ms in our case, whereas it was approximately 67 ms in Ref. 30. Further, regarding the initial oscillations, it was reported that the transition from sinusoidal to aperiodic oscillations occurs at approximately Ω = 0.75 in Ref. 30; however, in our simulation, such sinusoidal oscillations were observed only when Ω = 0.57. In fact, they were already non-sinusoidal in nature when Ω = 0.7. A possible reason for these discrepancies is that the estimation of C sph was inaccurate because it was roughly calibrated using the measured relationship between the number of vortices and Ω, as previously described. However, in our SPH scheme, C sph is a critical factor because the magnitude of the rotational force generated by the first term in Eq. (2) is proportional to C sph. In summary, the time dependence of α in our scheme is still under discussion and requires further research.

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.29 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.19 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.

1.
M. M.
Cerimele
,
M. L.
Chiofalo
,
F.
Pistella
,
S.
Succi
, and
M. P.
Tosi
, “
Numerical solution of the Gross-Pitaevskii equation using an explicit finite-difference scheme: An application to trapped Bose-Einstein condensates
,”
Phys. Rev. E
62
,
1382
(
2000
).
2.
S. K.
Adhikari
, “
Numerical study of the spherically symmetric Gross-Pitaevskii equation in two space dimensions
,”
Phys. Rev. E
62
,
2937
(
2000
).
3.
M. L.
Chiofalo
,
S.
Succi
, and
M. P.
Tosi
, “
Ground state of trapped interacting Bose-Einstein condensates by an explicit imaginary-time algorithm
,”
Phys. Rev. E
62
,
7438
(
2000
).
4.
G.
Vergez
,
I.
Danaila
,
S.
Auliac
, and
F.
Hecht
, “
A finite-element toolbox for the stationary Gross–Pitaevskii equation with rotation
,”
Comput. Phys. Commun.
209
,
144
(
2016
).
5.
P.
Henning
and
A.
Målqvist
, “
The finite element method for the time-dependent Gross–Pitaevskii equation with angular momentum rotation
,”
SIAM J. Numer. Anal.
55
,
923
(
2017
).
6.
P.
Heid
,
B.
Stamm
, and
T. P.
Wihler
, “
Gradient flow finite element discretizations with energy-based adaptivity for the Gross-Pitaevskii equation
,”
J. Comput. Phys.
436
,
110165
(
2021
).
7.
Y. Y.
Choy
,
W. N.
Tan
,
K. G.
Tay
, and
C. T.
Ong
, “
Crank-Nicolson implicit method for the nonlinear Schrodinger equation with variable coefficient
,”
AIP Conf. Proc.
1605
,
76
(
2014
).
8.
P.
Muruganandam
and
S.
Adhikari
, “
Fortran programs for the time-dependent Gross–Pitaevskii equation in a fully anisotropic trap
,”
Comput. Phys. Commun.
180
,
1888
(
2009
).
9.
P.
Wang
and
C.
Huang
, “
Split-step alternating direction implicit difference scheme for the fractional Schrödinger equation in two dimensions
,”
Comput. Math. Appl.
71
,
1114
(
2016
).
10.
L. Z.
Li
,
H.-W.
Sun
, and
S.-C.
Tam
, “
A spatial sixth-order alternating direction implicit method for two-dimensional cubic nonlinear Schrödinger equations
,”
Comput. Phys. Commun.
187
,
38
(
2015
).
11.
L.
Lehtovaara
,
J.
Toivanen
, and
J.
Eloranta
, “
Solution of time-independent Schrödinger equation by the imaginary time propagation method
,”
J. Comput. Phys.
221
,
148
(
2007
).
12.
P.
Bader
,
S.
Blanes
, and
F.
Casas
, “
Solving the Schrödinger eigenvalue problem by the imaginary time propagation technique using splitting methods with complex coefficients
,”
J. Chem. Phys.
139
,
124117
(
2013
).
13.
A.
Goldberg
and
J. L.
Schwartz
, “
Integration of the Schrödinger equation in imaginary time
,”
J. Comput. Phys.
1
,
433
(
1967
).
14.
D. L.
Feder
,
C. W.
Clark
, and
B. I.
Schneider
, “
Vortex stability of interacting Bose-Einstein condensates confined in anisotropic harmonic traps
,”
Phys. Rev. Lett.
82
,
4956
(
1999
).
15.
D. L.
Feder
,
C. W.
Clark
, and
B. I.
Schneider
, “
Nucleation of vortex arrays in rotating anisotropic Bose-Einstein condensates
,”
Phys. Rev. A
61
,
011601
(
1999
).
16.
P.
Muruganandam
and
S. K.
Adhikari
, “
Bose–Einstein condensation dynamics in three dimensions by the pseudospectral and finite-difference methods
,”
J. Phys. B
36
,
2501
(
2003
).
17.
J.
Rogel-Salazar
, “
The Gross–Pitaevskii equation and Bose–Einstein condensates
,”
Eur. J. Phys.
34
,
247
(
2013
).
18.
L.
Salasnich
, “
Bright solitons in ultracold atoms
,”
Opt. Quantum Electron.
49
,
409
(
2017
).
19.
S.
Tsuzuki
, “
Theoretical framework bridging classical and quantum mechanics for the dynamics of cryogenic liquid helium-4 using smoothed-particle hydrodynamics
,”
Phys. Fluids
34
,
127116
(
2022
).
20.
O. C.
Idowu
,
D.
Kivotides
,
C. F.
Barenghi
, and
D. C.
Samuels
,
Numerical Methods for Coupled Normal-Fluid and Superfluid Flows in Helium II
(
Springer
,
Berlin/Heidelberg
,
2001
), pp.
162
176
.
21.
A. W.
Baggaley
and
S.
Laizet
, “
Vortex line density in counterflowing He II with laminar and turbulent normal fluid velocity profiles
,”
Phys. Fluids
25
,
115101
(
2013
).
22.
S.
Yui
,
M.
Tsubota
, and
H.
Kobayashi
, “
Three-dimensional coupled dynamics of the two-fluid model in superfluid 4He: Deformed velocity profile of normal fluid in thermal counterflow
,”
Phys. Rev. Lett.
120
,
155301
(
2018
).
23.
C. L.
Horner
and
R. A.
Van Gorder
, “
Dynamics of quantized vortex filaments under a local induction approximation with second-order correction
,”
Phys. Fluids
31
,
065103
(
2019
).
24.
S.
Yui
,
H.
Kobayashi
,
M.
Tsubota
, and
W.
Guo
, “
Fully coupled two-fluid dynamics in superfluid 4 He: Anomalous anisotropic velocity fluctuations in counterflow
,”
Phys. Rev. Lett.
124
,
155301
(
2020
).
25.
R. A.
Gingold
and
J. J.
Monaghan
, “
Smoothed particle hydrodynamics: Theory and application to non-spherical stars
,”
Mon. Not. R. Astron. Soc.
181
,
375
(
1977
).
26.
M. B.
Liu
and
G. R.
Liu
, “
Smoothed particle hydrodynamics (SPH): An overview and recent developments
,”
Arch. Comput. Methods Eng.
17
,
25
(
2010
).
27.
Y.
Imoto
,
S.
Tsuzuki
, and
D.
Nishiura
, “
Convergence study and optimal weight functions of an explicit particle method for the incompressible Navier–Stokes equations
,”
Comput. Part. Mech.
6
,
671
(
2019
).
28.
T.
Ye
,
D.
Pan
,
C.
Huang
, and
M.
Liu
, “
Smoothed particle hydrodynamics (SPH) for complex fluid flows: Recent developments in methodology and applications
,”
Phys. Fluids
31
,
011301
(
2019
).
29.
S.
Tsuzuki
, “
Reproduction of vortex lattices in the simulations of rotating liquid helium-4 by numerically solving the two-fluid model using smoothed-particle hydrodynamics incorporating vortex dynamics
,”
Phys. Fluids
33
,
087117
(
2021
).
30.
K.
Kasamatsu
,
M.
Tsubota
, and
M.
Ueda
, “
Nonlinear dynamics of vortex lattice formation in a rotating Bose-Einstein condensate
,”
Phys. Rev. A
67
,
033610
(
2003
).
31.
M.
Tsubota
,
K.
Kasamatsu
, and
M.
Ueda
, “
Vortex lattice formation in a rotating Bose-Einstein condensate
,”
Phys. Rev. A
65
,
023603
(
2002
).
32.
P. R.
Berman
,
Mathematical Preliminaries
(
Springer International Publishing
,
Cham
,
2018
), pp.
33
51
.
33.
J. J.
Monaghan
, “
Smoothed particle hydrodynamics
,”
Annu. Rev. Astron. Astrophys.
30
,
543
(
1992
).
34.
S.
Koshizuka
and
Y.
Oka
, “
Moving-particle semi-implicit method for fragmentation of incompressible fluid
,”
Nucl. Sci. Eng.
123
,
421
(
1996
).
35.
S.
Koshizuka
,
A.
Nobe
, and
Y.
Oka
, “
Numerical analysis of breaking waves using the moving particle semi-implicit method
,”
Int. J. Numer. Methods Fluids
26
,
751
(
1998
).
36.
S.
Tsuzuki
, “
Particle approximation of the two-fluid model for superfluid 4He using smoothed particle hydrodynamics
,”
J. Phys. Commun.
5
,
035001
(
2021
).
37.
K. W.
Madison
,
F.
Chevy
,
V.
Bretin
, and
J.
Dalibard
, “
Stationary states of a rotating Bose-Einstein condensate: Routes to vortex nucleation
,”
Phys. Rev. Lett.
86
,
4443
(
2001
).
38.
D.
Luebke
, “
CUDA: Scalable parallel programming for high-performance scientific computing
,” in
5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro
(
IEEE
,
2008
), pp.
836
838
.
39.
G. S.
Grest
,
B.
Dünweg
, and
K.
Kremer
, “
Vectorized link cell Fortran code for molecular dynamics simulations for a large number of particles
,”
Comput. Phys. Commun.
55
,
269
(
1989
).
40.
M.
Gomez-Gesteira
et al, “
Sphysics—Development of a free-surface fluid solver—Part 1: Theory and formulations
,”
Comput. Geosci.
48
,
289
(
2012
).
41.
P.
Ramachandran
and
K.
Puri
, “
PySPH: A framework for parallel particle simulations
,” in
Proceedings of the 3rd International Conference on Particle-Based Methods (Particles 2013)
, Stuttgart, Germany,
2013
.
42.
R.
Kishor Kumar
,
V.
Loncar
,
P.
Muruganandam
,
S. K.
Adhikari
, and
A.
Balaz
, “
C and Fortran OpenMP programs for rotating Bose–Einstein condensates
,”
Comput. Phys. Commun.
240
,
74
(
2019
).