Towards Highquality Visualization of Superfluid Vortices
Abstract
Superfluidity is a special state of matter exhibiting macroscopic quantum phenomena and acting like a fluid with zero viscosity. In such a state, superfluid vortices exist as phase singularities of the model equation with unique distributions. This paper presents novel techniques to aid the visual understanding of superfluid vortices based on the stateoftheart nonlinear KleinGordon equation, which evolves a complex scalar field, giving rise to special vortex lattice/ring structures with dynamic vortex formation, reconnection, and Kelvin waves, etc. By formulating a numerical model with theoretical physicists in superfluid research, we obtain highquality superfluid flow data sets without noiselike waves, suitable for vortex visualization. By further exploring superfluid vortex properties, we develop a new vortex identification and visualization method: a novel mechanism with velocity circulation to overcome phase singularity and an orthogonalplane strategy to avoid ambiguity. Hence, our visualizations can help reveal various superfluid vortex structures and enable domain experts for related visual analysis, such as the steady vortex lattice/ring structures, dynamic vortex string interactions with reconnections and energy radiations, where the famous Kelvin waves and decaying vortex tangle were clearly observed. These visualizations have assisted physicists to verify the superfluid model, and further explore its dynamic behavior more intuitively.
Superfluid dynamics, vortex structure, visual analysis
1 Introduction
Superfluidity is a special state of matter [1] first discovered in liquid helium (HeliumII) in 1930s [2]. Since then, various superfluid phenomena have been found in different areas of physics, from condensed matter physics, high energy physics, to astrophysics.
Superfluids have many surprising properties, some of which are counterintuitive. For example, superfluids have zero viscosity, so they can move frictionlessly; superfluids can exhibit extremely high thermal conductivity, so heat can be conducted in almost no time. In one particular situation, liquid helium in superfluidity can flow out of the containing cup by itself and move along the wall without any external force [3]. These interesting properties and behaviors result from macroscopic quantum phenomena, which cause strong correlation among the superfluid particles. Therefore, superfluids behave very differently than classical fluids [1].
Studying superfluids has several potential applications.
For example, superfluid liquid helium has been used as a moderator to cool down neutrons.
The BoseEinstein condensates (BEC
Superfluid modeling. Early research on superfluids relied on real experiments [2]. Physicists had to design complicated methods to create a specific physical environment for a matter to reach the state of superfluidity, so that they can measure the physical properties of superfluids, e.g., thermal conductivity, dynamic viscosity, etc. However, it is highly challenging to obtain highprecision experimental measurements [4]. Hence, it is very difficult to find out how superfluids behave subject to varying physical parameters and system configurations purely by real experiments.
This motivates physicists to develop mathematical models [5, 6] to describe the dynamic behavior of superfluids. However, unlike classical fluids, which are governed by the wellknown NavierStokes equations, there is no wellestablished superfluid model that has been thoroughly verified and is general in handling a large variety of situations. So far, only a few simple models [7, 8] have been experimentally affirmed, and they do not involve superfluid vortices. Research on superfluid modeling is still in active development. The research in this paper is a collaborative work with theoretical physicists who study superfluid modeling, in which we aim to develop simulation and visualization methods to help them explore potential properties of superfluid models, as well as verify and compare their models with real observations.
Superfluid vortices. One distinctive characteristic of superfluids is the formation of quantum vortices [9, 10], which, in contrast to vortices in classical fluids, exhibit a quantized flux circulation of some quantity, e.g., phase. This leads to some special phenomena not observed in classical fluids. For example, a steady regularlydistributed vortex lattice can be experimentally observed in rotating superfluids, where the vortices are completely isolated from one to another. Superfluid vortices are basic elements in quantum turbulence. In addition to experimental observations in condensed matter physics, they are also applied to nuclear physics, highenergy physics, and cosmology, e.g., cosmic strings and vortex tangles in early Universe simulations [11, 5, 6].
Vortices are important features in fluids with extensive amount of work devoted to study them [12, 13] and to perform visual analysis [14, 15, 16, 17, 18, 19] on them in classical fluids. Since superfluid vortices differ from their classical counterparts, mainly due to the phase singularity of the model equation, we cannot directly apply existing flow visualization methods to superfluid vortices. Moreover, studying superfluid vortices is an emerging research topic in physics, and much more is to be understood about their physical properties. Hence, we work closely with physicists to design simulation and visualization methods to explore superfluid vortices based on their expertise.
Overview of this work. In this work, our goal is to develop visualization methods to aid physicists to visualize and explore superfluid vortices, with an emphasis on applications to support intuitive scientific study. To visualize superfluid vortices, data sets need to be first prepared typically by means of simulations based on the governing model equations. Here in particular, we employ and solve the model equation derived in [6], which is a very recent superfluid model, theoretically proven to be a more accurate one that can cover the entire velocity range, from zero to the speed of light. Note that while [6] focuses on the model derivation, this work develops methods for highquality superfluid vortex visualization, where we numerically solve the model equation with efficiency (which has not been done in [6]), enhance the quality of the resulting superfluid data sets, effectively identify and extract superfluid vortices, and then visualize different superfluid phenomenon, as well as perform associated visual study with domain experts in superfluid research. In summary, this work has the following key technical contributions:

First, to develop superfluid visualization, we need prepare appropriate data sets. As far as we know, there are no publiclyavailable superfluid flow data sets that could be used for visualization. Thus, we derive methods to generate superfluid flow data sets for visualization by numerically solving the model equation. To obtain a stable solution of superfluid flow in 2D, we propose transforming the underlying equation to polar coordinate domain to eliminate the numerically unstable temporalspatial cross term, and devise a stable hybrid method. In 3D, rather than develop complicated numerical solver, we instead adopt model approximation and develop an efficient but stable highorder finite difference scheme.

Second, to enhance the visualization quality, we devise a temporal filtering method to suppress the high frequency noiselike waves, which significantly improves the quality of the data sets, and thus the quality of the visualizations. Note that this is a very important step, since without such a filtering, noiselike waves may seriously influence the vortex identification, which could lead to erroneous visualizations.

Third, due to the phase singularity of the superfluid model equation at the vortices, existing approaches are difficult to apply for vortex visualization. Thus, we propose to utilize the singularity instead, and develop a new technique based on velocity circulation when identifying superfluid vortices, especially in 3D. In theory, we can compute such a circulation along any closed curve around a grid point, and for simplicity, a closed curve on a plane in 3D. However, ambiguities in numerical computation of circulation may fail to find certain important vortex structures. Hence, we devise an orthogonalplane strategy to obtain complete circulation values over the field, so that we can faithfully identify a full set of superfluid vortices. Furthermore, together with filtering (smoothing), we can easily extract vortex tubes suitable for visualization.

Fourth, as a direct application of our techniques, we perform visual study and analysis of superfluid vortices for both steady and unsteady scenarios, and present various superfluid phenomena as our visualization results. In the steady case, we study superfluids in a rotating frame, and explore the regularlydistributed vortex lattice and vortex ring structures. In addition, we examine and verify the Feynman relation [20] in 2D in order to verify the appropriateness of the physical model. In the unsteady case, we study superfluids in an inertial frame, and evolve and visualize vortex strings with reconnections, Kelvin wave formation and transmission, leapfrogging, as well as vortex tangle (a turbulent state of superfluids). With these visualizations, scientists can obtain higherquality visual analysis of vortex dynamics for superfluid studies as well as scientific experiments.
In this work, we consider different superfluid conditions and produce assorted visualization results, either steady or unsteady; e.g., Fig. 1 shows the dynamic formation process of a 3D superfluid vortex ring structure discovered and visualized by our method, which has never been observed before; note that the similar result presented in [6] was generated by the numerical simulation technique presented in this paper, but visualized with direct volume rendering of the field magnitude (density), which is less accurate; we will explain later that such visualization is also not generally applicable. Moreover, we verify and compare our method with related physical models together with domain physicists, with positive expert feedbacks. Furthermore, this work is the first attempt we are aware of in developing visualization method based on superfluid model simulation data sets; it is useful for visual exploration and analysis of vortex structures in superfluids, and contributes to open a new area of visualization on macroscopic quantum phenomena in fluids.
2 Related work
Superfluid vortices have been observed and studied in various experiments. AboShaeer et al. [9] used laser beams to confine and rotate BoseEinstein condensates (BEC) and produced a vortex lattice structure in superfluids, while Lathrop et al. [10] observed turbulence with vortices in superfluid liquid helium.
Superfluid simulation. The rotating BEC in superfluid simulations often adopts the GrossPitaevskii equation (GPE), which is a lowspeed limit model [21] often used as a tool for studying superfluid vortex dynamics. Tsubota et al. [21] solved the GPE in 2D to simulate vortex lattice in rotating superfluids, while Kasamatsu et al. [22] solved the GPE in 3D and observed vortex tubes with lattice formation. Zuccher et al. [23] analyzed the quantum vortex reconnection, while Caplan et al. [24] explored vortex ring dynamics. This work adopts a novel superfluid model [6], which covers the entire velocity range and incorporates the previous models as its nonrelativistic limit. The work of [6] focuses on deriving the equation for a more general and accurate superfluid model. In contrast, this work focuses instead on methods towards highquality superfluid visualizations, including simulation and vortex identification techniques for superfluid data set preparation and vortex structure visualization.
Vortex core identification. There are a number of methods for identifying vortex cores in classical fluid flow fields. Hunt et al. [25] proposed the Qcriterion, while Jeong and Hussain [26] proposed the criterion. Jiang et al. [27] provided a taxonomy of identification methods according to multiple criteria. While most of the existing techniques are local, Weißmann et al. [28] formulated a global method to identify vortex core lines over a vector field based on quantum mechanics analogy.
Different from the above methods, which are designed for classical fluid flows, we develop a novel method to identify superfluid vortices from a complexvalued flow field. Our method considers the phase singularity of vorticity in the superfluid model, and is able to effectively identify vortices in superfluid flows.
Vortical flow visualization. There are several approaches to visualize vortices in classical fluid flows, one of which is the lineintegralconvolution (LIC) visualization. Wiebel et al. [29] embedded streamlines in an LIC texture to explore boundaryinduced vortices. Krishnan et al. [30] attached streamlines to flow surfaces in the visualizations. Yu et al. [31] created a hierarchy of streamline bundles for multiresolution exploration of flow fields. Chaudhuri et al. [32] enhanced salient features in streamline visualizations by multiscale analysis.
Featurebased methods are another common approaches in vortical flow visualizations. JankunKelly et al. [33] proposed to detect and visualize vortices in engineering environments. Schneider et al. [34] used the criterion with the largest contours to extract isosurfaces, revealing vortical features in a turbine flow. Schafhitzel et al. [14] visualized hairpin vortices with isosurface rendering. Treib et al. [16] developed an interactive visualization system to show extremely detailed terascale turbulence simulations on a desktop PC.
Other vortical flow visualization methods focus on different portions and aspects of vortices in classical flows. Sadlo et al. [35] visualized vorticity transport in an incompressible flow. Laney et al. [36] used MorseSmale complex to study turbulent mixing in RayleighTaylor instability. Weinkauf et al. [37] extracted vortex cores from swirling particle motion in unsteady flows. Helgeland et al. [38] extracted and visualized vortex structures in wallbounded turbulent flows, while Johnson et al. [39] proposed a framework for interactive visualization of transitional flows.
Vortical flow visualization has also led to practical applications in many different disciplines. In medical study, Soni et al. [40] used particle trajectories to visualize flows in small bronchial tubes; Zachow et al. [41] developed visualizations for exploring nasal flow; Köhler et al. [42] used line predicates to extract vortices in cardiac MRI flows, and Born et al. [43] also used line predicates to extract blood flow in MRI data. In turbine study, Shaffi et al. [19] visualized vortices in turbulent wake flows produced in a simulated wind farm, while Koehler et al. [15] studied the flow around the dragonfly wings.
However, very little work has been done on simulating and visualizing superfluids. Zuccher et al. [23] extracted density isosurface to visualize quantum vortex reconnection in nonrelativistic superfluids. The method has limited accuracy and its extraction results are datadependent. Very recently, Guo et al. [44] visualized the vortices in a superconductor simulation, focusing on extracting the topology of vortex lines. They targeted to extract spatiotemporal vortex information, and adopted a sophisticated graphbased searching algorithm based on a global formulation. In contrast, our goal in this work is to visualize and analyze the spatial and temporal structures of vortices in superfluids based on our simulation data sets, thus requiring efficient methods to support extensive physical analysis.
With this in mind, we developed an efficient superfluid simulation and a vortex identification/visualization method based on the circulation field and isosurface extraction, which is local, fast, and parallelizable.
Using our method, vortex tubes
3 Superfluid Model
Unlike classical fluids, superfluids are described by a complexvalued scalar field (see in Eq. 1). Here, we first review the nonlinear KleinGordon (NLKG) equation derived and formulated in [6], which is the fundamental model equation used in our simulation and visualization.
Nonlinear KleinGordon equation. The reason we choose the nonlinear KleinGordon equation as our governing equation for superfluid flows is that such an equation can be used for the entire velocity range, especially the relativistic region in which superfluidity occurs. The equation comes from the quantum field theory and it has the following form:
(1) 
where is a complex scalar field defined as
(2) 
with magnitude and phase ; includes a nonlinear term and a mass term to describe the selfinteraction among the superfluid particles; and is time. Note that is parameterized by the nonlinearity constant , the field magnitude at infinity , and the field mass . Here, can be used to derive the hydrodynamic density and velocity as: and , which will be used later in our superfluid vortex visualization.
Note that in superfluids, since the flow velocity results from the gradient of phase ( is a scalar field), it is curl/vorticityfree, i.e., if is nonsingular. Thus, the existence of superfluid vorticity is encoded into the phase singularity of [20, 6]. The formation of superfluid vortices requires certain conditions, e.g., the rotation of the system or the interaction of the vortex rings; when superfluid vortices are formed, they never dissipate locally until they reconnect (where the dissipation comes from the radiating waves at the reconnections), and their topologies can only be changed by the reconnection of vortex lines or the decomposition of a vortex ring into smaller ones.
Superfluids in rotation. The rotation of superfluids not only leads to vorticity, but also changes the Hamiltonian of the system. Hence, superfluid vortices spontaneously arrange into a regular lattice structure in favor of low energy under constant rotation. Such a nontrivial configuration is a unique characteristic of superfluids from the macroscopic quantum phenomena.
To model superfluids (with vortices) in a rotating frame, we need to include a rotation term in Eq. 1 as:
(3) 
There are two known approaches to model :

Coordinate transformation formulation [6]. This approach employs coordinate transformation and considers superfluid flows directly in a rotating frame without approximation, and formulates as:
(4) where is the constant (timeinvariant) angular velocity of the entire system and is the distance vector to the center of rotation. Such a coordinate transformation is a timedependent rotation which involves noninertial frame of reference (the rotating frame). As shown in [6], the rotation term in Eq. 4 can be obtained via such a coordinate transformation. Although this is a physicallyconsistent mathematical formulation, it contains a cross term (i.e., ) consisting of commutative spatial and temporal derivatives of the field , which induces numerical instability in the simulation and requires a sophisticated numerical solver to obtain stable solutions even for steady and lowspeed flows.

Currentcurrent formulation [5]. This alternative approach uses virtual forces to drive the superfluids in order to produce a rotation. It can be considered as a model simulating local rotation systems, i.e., the angular velocity is positiondependent; Although this formulation is not as accurate as the coordinate transformation formulation, stable numerical solutions can be achieved with high computational efficiency. Here is formulated as:
(5) where is a coupling constant; is a Gaussianshaped trapping potential; and is the angular velocity of the constant rotation.
4 Model Simulation
Since there are no publicly available solver and data sets for superfluid flows based on NLKG (with and without rotation), we first need to numerically solve the model equation with enough accuracy to simulate the superfluid flows before we perform vortex visualization. To simulate highquality superfluid flows, we devise new numerical solvers for the superfluid model equation in both 2D and 3D respectively to produce vortex structures in a relatively high resolution grid.
For 2D simulations, we focus on rotating superfluids, where we adopt the coordinate transformation formulation and propose to transform the equation into polar coordinate domain, where a new hybrid numerical method based on Fourier transform and finite difference can be derived. Such a simulation in 2D can accurately demonstrate the essence of vortex lattice in rotating superfluids, e.g. the Feynman relation (discussed later), and can be readily implemented in the numerical computations in comparison with the 3D scenarios.
For 3D simulations, we consider both rotating and freeevolving superfluids. The freeevolving ones, as long as initialized properly with appropriate parameters, can eventually produce the vortex tangle due to the nonlinearity of the model. Here, we adopt the currentcurrent formulation to reduce the numerical complexity and propose a highorder finitedifference scheme to solve the model equation in Cartesian domain. Note that only highorder schemes can correctly preserve the vortex shapes.
It is worth mentioning that the NLKG equation we study is dimensionless, i.e., we enforce the speed of light , so that other physical quantities are scaled accordingly for more convenient and more stable numerical simulations as in many traditional fluid solvers.
4.1 2D Simulation
If we perform 2D superfluid simulations with the coordinate transformation formulation directly in Cartesian coordinate space, the cross derivatives in would induce strong numerical instability no matter how we tune the simulation parameters such as the time step. To overcome this issue, we propose to transform and compute the entire model equation in polar coordinate domain.
Polar coordinate domain formulation. By transforming the model equation to polar coordinate domain, we can reduce the numerical complexity and simplify as:
(6) 
where denotes the azimuthal angle, and the Laplacian operator is transformed to be:
(7) 
To improve the simulation quality, we propose to decouple the cross term in Eq. 6, since it can easily cause instabilities.
Polar spectral discretization. To decouple the remaining cross derivative in , we perform Fourier transform of the scalar field in polar coordinate domain along the dimension of :
(8) 
and insert the expansion into Eq. 3 with rotation term formulated by Eq. 6 to obtain:
(9) 
where the hat sign ( ) indicates the Fourier transform and is the wavenumber. It is important to note that with Fourier transform, we can separate the spatial and temporal derivatives to avoid the cross derivatives and thus enhance the stability of the simulation.
Radial finitedifference discretization. To perform simulations with Eq. 9, we further discretize the domain over radial dimension and time. For temporal discretization, we employ the semiimplicit CrankNicolson scheme[45] with temporal averaging on the linear terms to compute the final discretization in Fourier domain. For spatial discretization, we approximate both the first and second radial derivatives on the righthandside of Eq. 9 in polar coordinate domain by secondorder central finitedifference schemes. The combination of the two discretizations is written as:
(10)  
(11) 
where is the time step index; is the time step size set to be 0.1; and the other configuration parameters are: domain radius is 120; is 1/240; is 1.0; is 6.4; is 1.0; is 5.0; the polar grid resolution is 256 (radial)1024 (angular); and is the spatial derivative operator in the radial dimension, which takes the following form:
(13)  
where is the uniform spacing along the radial dimension and is the number of samples along the radius. The grid points are cellcentered, so that , to avoid the polar singularity like in [46]. Hence, the radial derivatives for the first grid point next to the pole does not require across the pole, since its summed coefficient is zero according to Eq. 13. The outer boundary along the radial dimension (not explicitly stored) is treated by a Neumann condition with finite difference discretization. The above treatments lead to a tridiagonal linear system, which can be solved efficiently in the simulation over time in Fourier domain. To obtain the final simulation result, we further apply an inverse Fourier transform.
Due to the nonlinearity, exact numerical stability analysis is difficult to be performed for the above derivation. However, the system is stable and accurate depending on the choice of . We empirically found that a large may induce instability and sacrifice accuracy, and setting balances between efficiency and accuracy, and improves the stability in the 2D simulations.
Initialization. Traditionally, the study of rotating superfluids often adopts the GPE [21] (which is effectively a lowspeed limit of the NLKG equation), which also verifies the experimental observations that rotating superfluids may converge to a vortex lattice state [9]. Such a state satisfies a stationary ansatz: , where is the global temporal derivative of the phase, and gives the appearance of the vortex lattice structure. To reach such a steady but nonstationary state, we adopt the following setting for the initial condition:
(14) 
where and stand for the initial distributions of the first and second time steps; is the azimuthal angle of polar coordinates; is the field magnitude at infinity; is an arbitrary profile to introduce anisotropy and perturbations, where we choose to use a bowl shape with slight azimuthal deformation; is an integer, such that the initial domain is configured to have a circulation of at the pole but nowhere else; this ensures the final distribution of stable vortex lattice, whose circulation is also in total; specifies the global temporal derivative of the phase at the beginning; and is the time step. Note that is not a physically conserved quantity. Since the final state may not fully converge in practice, the temporal derivative of the phase is not constant in space, and can only be measured by spatial averaging. Hence, , which is initially specified in Eq. 14, may generally not equal to in the final state.
4.2 3D Simulation
We may extend Section 4.1 to 3D simulations of superfluids in rotation by using spherical coordinates and adopting spherical harmonics. However, this may complicate the solution and introduce excessive computational overhead. Hence, we make a compromise between accuracy and computational efficiency and adopt a model approximation approach where the currentcurrent formulation is used in 3D simulations with rotation. In the following, we consider two different scenarios of 3D superfluid simulations:
Scenario 1: 3D rotating superfluids. Since the currentcurrent formulation does not couple the spatial and temporal derivatives, we can efficiently discretize the equation. In detail, we handle the first and second derivatives of both the spatial and temporal terms in Cartesian coordinates by second and fourth order central finitedifference schemes in time and space, respectively, in explicit form and with conditional stability. Note that highorder spatial discretization should be used to reduce the discretization anisotropy, which may distort the shape of the vortex lattice. Hence, the discretization leads to the following linear system:
(15)  
(16) 
where is the spatial component of in the currentcurrent formulation: and is the discrete Laplacian operator. In our formulation, we introduce a global rotation about the zaxis, so has no zcomponent. Moreover, since the linear system is diagonal, we can solve it efficiently with the following parameters: the domain side length ; ; ; ; ; ; and the grid resolution is . The system is conditionally stable depending on the choice of , and large reduce simulation accuracy and increase system instability. In our simulations, we set to be 0.004 to ensure both stability and accuracy. The periodic boundary condition is applied at the domain boundary and arbitrary initialization can be used inside the domain provided that there is an initial perturbation with sufficient energy on .
By adjusting the trapping potential , we can obtain different vortex lattice solutions. For example, by using a steep Gaussianshaped trapping potential (which is a 3D Gaussian with smaller deviation and larger magnitude), we can obtain the vortex ring lattice shown in Fig. 2 (left). However, by using an oblate Gaussianshaped trapping potential (which is a more flat 3D Gaussian with larger deviation and smaller magnitude), some straight vortex lines could emerge in the solution, see Fig. 2 (right).
Scenario 2: 3D superfluid vortex dynamics. On the other hand, we can stop the rotation by giving zero angular velocity and create 3D simulations to show how vortices dynamically evolve and interact over time without a global rotation. In this case, since vortices do not occur spontaneously, we need to prepare some initial vortices in the simulation domain.
Our first approach to prepare initial vortices is to extract the vortices produced from Scenario 1 and take them to create the initial condition before evolution. Specifically, we take a volume of data from the previous simulation result as , and set according to Eq. 14, where is tweaked to be compatible with other parameters. The asymmetry in the initial condition results in complex vortex motion, such as vortex reconnection and Kelvin waves. However, this initial condition is restricted by the number of vortex lines generated from the 3D rotating superfluids, which is relatively sparse. Thus, it cannot produce sufficiently dense vortex interaction dynamics, and then vortex tangle. Moreover, the vortex tangle always decays during interaction, meaning that we will have fewer and fewer vortices as time proceeds.
To generate more interesting and more complex vortex dynamics with much denser vortex tangle, we employ a boundary injection approach, which is to artificially “shoot” closed vortex rings from the boundary of the simulation domain. This is physically feasible to do since we can produce the vortex rings easily with simple mechanism similar to producing a smoke ring. To do so, we first prepare a small volume of containing a single closed ring. Such a small volume of closed ring can be directly shifted to the boundary edge of the whole simulation domain due to the periodic boundary setting. The ring can be added to the simulation by multiplying the prepared field volume with and at some randomly selected time step . We then continuously place randomly shifted rings at the boundary to simulate rings coming from the boundary. Finally, we stop placing rings when the underlying vortex dynamics produces sufficiently dense tangles, and we let them evolve freely to generate chaotic vortex structures.
4.3 Temporal Filtering
Due to nonlinear selfinteraction among superfluid particles, very highfrequency waves could appear in the simulation results; this situation is physically intrinsic but not numerical due to the compressibility nature of the model equation. Since different bands of these waves interfere with one another in different directions, the composed waves look like noise in spatial domain, see Fig. 3 (top), which significantly deteriorate the visual inspection of vortices. In principle, these highfrequency waves are plausible solutions to the NLKG equation. However, they only act as smallscale compression waves and are of little interest to the domain scientists for studying superfluid vortices, since superfluid vortices are much more structured and stationary with much lower temporal frequency than these noiselike waves. To get a clearer physical picture, the highfrequency waves must be removed during the simulation. To address this issue, we may introduce a damping term into the model equation (Eq. 1) as has been applied in BoseEinstein condensates simulations [21], but this method will result in excessive energy loss and hinder the formation of a vortex lattice, which requires finite in the final steady state.
From preliminary simulation results, we observe that such highfrequency waves are temporally incoherent even though they vary significantly over time; in addition, compared to the variation of vortices, they change much faster over time. Hence, we adopt a temporal filtering method to suppress the influence of these noiselike waves. In detail, we construct a temporal lowpass filter [47] in both 2D and 3D simulations, and apply it to the magnitude of only. Mathematically, the filter is expressed as:
(17) 
where controls the filter strength, typically chosen to be close to 1.0. With this temporal filtering, we can significantly enhance the quality of the simulation results (see Fig. 3 (bottom)) for better vortex core identification and clear vortex structure visualization. This filtering technique differs from the traditional filtering methods, since we isolate magnitude and phase as independent components, where it is nontrivial to conduct spectral analysis. However, this technique has proven to be very effective since the dynamics is described by a complexvalued equation, where the magnitude and phase may have independent modes. Furthermore, it is easy to see that our technique does not influence the stationary solution, where the magnitude has only the lowest mode with the phase being a single finite mode.
5 Vortex Identification and Visualization
After simulating and obtaining the superfluid flow, we can identify vortex cores from the data set in order to extract the vortex structures for visualization. Note that we aim to visualize the structures for vortex cores only, but not the entire structure of vortices. Thus, we ignore other structures such as the spiral patterns, which are fluctuating structures around the vortex cores that occur at very small scales. To overcome the difficulty that superfluid vorticity is singular, our key idea is to employ velocity circulation and develop a specific circulationbased numerical method with an orthogonalplane strategy to efficiently identify vortex cores while avoiding ambiguity in the circulation computation with accuracy up to the grid resolution. Here, we want to emphasize that subgrid scale vortices cannot be faithfully differentiated without ambiguity; we will discuss in more detail later.
5.1 Observation: Circulation in Superfluids
As stated in Eq. 2, the superfluid flow field is denoted as , where (phase) is singular at the vortex cores. It is known from both existing theory and our visualization (see Fig. 3(c)&(d)) that superfluid vortex cores coincide with the singularity in phase space, where the phase value is undefined, e.g., the branch point A in Fig. 3(d). This is a unique characteristic of the superfluid model different from the classical model, where vortices do not possess singularities for the related physical quantities over the flow field. In addition, superfluids are potential (irrotational) flows except at the vortex cores, where phase singularities cause the vortices.
From the above analysis, it is easy to think of circulation, which can be employed for vortex core identification. If we compute the circulation along a small closed loop (either 2D or 3D) around a singular point in the phase field, e.g., a closed loop around point A in Fig. 3(d), we will find that such circulation is an integer multiple of ; the integer number indicates the number of vortex cores inside the closed loop. However, if we compute the circulation around any nonsingular point, e.g., point B in Fig. 3(d), even if the phase is discontinuous, we will find that the circulation is always zero. This is the key criterion that we will explore to formulate our vortex core identification method, with which we will employ to visualize and analyze the vortices.
5.2 Vortex Core Identification
In general, the circulation is computed as an integral:
(18) 
where v is the flow velocity, is the phase of the complexvalued scalar field , and is a closed loop in the phase field. Since superfluid vortices are usually very small with infinitesimally small vortex cores, we use a small loop (could be very arbitrary) to enclose a point, so that the resulting circulation is always near the vortices, and zero elsewhere. Hence, we can produce a welldefined circulation field, and then threshold it to identify vortices. Here, we use a threshold value , above which the points are considered as vortices; such a value can allow certain tolerance of circulation fluctuation from numerical evaluation with very coarse grid points. However, can be finetuned manually, but should help to successfully separate vortex points from nonvortex points.
It is worthwhile to note that existing works on identifying quantum vortices usually rely on density isocontour [23] (equivalently the magnitude of ), which identifies vortices simply by thresholding the density value, and it drops to zero near vortex cores but being flat for other regions. However, when vortices undergo reconnections, not only their shapes change topologically, but also energies radiate out, see the red boxes in Fig. 4, which show a 3D simulation of freeevolving superfluid vortex dynamics after a reconnection. Such energy radiation significantly affects the superfluid density, so the density isocontour approach may generate unwanted nonvortex structures and false positive superfluid vortex structures, see Fig. 4 (a). In Fig. 4 (b), our vortex identification method, which makes use of the circulation property, can correctly capture the vortices without producing additional nonvortex structures, see the corresponding red box.
As mentioned in Section 2, Guo et al. [44] proposed a comprehensive approach to identify superconductor vortices by detecting phasesingularityinduced circulation. We also adopt the circulation to identify superfluid (quantum) vortices. However, we propose a more practical framework, which is local and fully parallelizable, enabling efficient computation on the GPU.
5.3 Numerical Evaluation
To compute a circulation field, one may arrange a dense rectangular grid over the spatial domain and compute the circulation integral (Eq. 18) along a small loop around each grid point. However, since is periodic (range: ), discontinuities appear across the period of , resulting in unstable evaluation of and leading to erroneous numerical integration. To overcome this issue, we first note that for a complexvalued field , its phase is defined as:
(19) 
where and are the real and imaginary parts of the field, respectively. If we denote and at location , we can reformulate in Eq. 18 as:
(20) 
which can be used for consistent numerical integration since both and are continuous. In practice, the field magnitude variation around a vortex may be very large, and the integrand may change steeply near the singularity due to division by . It is obvious that direct numerical integration using the trapezoidal rule requires smooth variation of the integrand to preserve accuracy. When sampling at locations very close to the vortex cores, the integral becomes inaccurate due to large variations. On the other hand, the circulation is used only to indicate singularity, which is independent of the field magnitude. Therefore, we normalize the field to unitary magnitude before computing the circulation, and reformulate the circulation as the integral of the gradient of the normalized phase field:
(21) 
where and are the real and imaginary parts of the normalized field. We will discretize in our practical computations.
Closed loop selection. In principle, the closed loop can be selected arbitrarily (note that in 2D, the loop is a closed planar curve, while in 3D, the loop is a closed 3D curve), but arbitrary loops could unnecessarily introduce additional computation and inaccuracy due to data interpolation. Therefore, we align loops with the simulation grid, such that we can simply and exactly fetch the data from the neighborhood grid points. In 3D, we form loops on axisaligned planes (either x, y, or zplane) to construct a circulation loop at each grid point. It is worth to note that when using a small loop, the integral may not be sufficiently accurate due to large variation of the data near the vortex core. However, when using a large loop, more than one vortex core may be included in the loop and will be identified as a single vortex core, which cannot reflect the true structure of the vortices up to the precision of grid resolution. As a compromise, we use a small loop confined to an axisaligned plane, which contains the nearest neighbor of the grid point on that plane as the integral loop, see Fig. 5(a), but we discard the four corners to make the loop closer to a circle in order to avoid artifacts that the identified vortex will have a boxshaped crosssection due to a boxlike path. By dividing the path L into segments, the numerical circulation can be reexpressed as:
(22) 
where is the th segment between two consecutive sample points, and the corresponding integral can be easily approximated numerically by the following trapezoidal rule:
(23)  
(24)  
(25) 
where and are the locations of the starting and end points of the th segment, respectively.
Note that since we use a loop with finite crosssection, we do not expect to locate the vortex cores exactly since they are infinitesimally small. Instead, we confine the vortex core within a finite tube with approximately twice the size of the grid spacing. Although we are not sure exactly the location of the vortex core, we are certain that the crosssection includes vortex cores whose mutual distances are smaller than the radius of the crosssection, implying that such a circulationbased scheme involves the aggregate of all smallscale vortices, and the shape of the tube indicates the spatial structure of the vortex cores; such spatial structure can be determined by the grid resolution. For vortices in subgrid scale, uncertainties may still occur, which can only be rectified by increasing the grid resolution.
Circulation ambiguity. For 2D simulations, the circulation field constructed above can effectively identify the vortex cores. However, for 3D simulations, if we compute the circulation field on planes all with one specific orientation, e.g., all on yz planes, we may miss some core vortex structures, see Figs. 5(a) and 6(a). The reason here is that the vortex core lines, which are parallel to the circulation loop planes, may not be effectively identified by the loops on those parallel planes.
In general, a vortex core line may or may not pass through a local circulation loop plane depending on the orientation of the plane. Since vortex cores may not overlap with the grid points, the numerical evaluation of circulations on planes almost parallel to the orientation of the vortex core lines will give closetozero value; hence, the detection of a portion of the vortex core lines may fail, and the grid point may be misidentified as a nonvortex point. Fig. 5 (a) shows the vortex core line in orange (arrow) and two parallel circulation planes at grid points A and B. At A, the vortex core line passes through the loop on the circulation plane, so we can successfully identify the vortex core. However, at the nearby point B, since the local vortex core line is parallel to the circulation plane, the vortex core identification fails even though the figure suggests the existence of a vortex core at B. Such ambiguous case becomes indistinguishable from the real case of a nonvortex point.
Orthogonalplane strategy. To resolve such ambiguity issue, we propose to use a set of three orthogonal planes at each grid point and compute the circulation values independently on each plane, see Fig. 5(b). Note that these three orthogonal planes form a complete set of basis planes that can be used to construct any planes that pass through a grid point for circulation evaluation. This means that vortex core lines that pass through the vicinity of a grid point with arbitrary orientation can be effectively captured by the circulation values from at least one of the three evaluations.
Thus, taking the maximum of the three circulation values (one for each grid planes) as our final circulation value can effectively and stably identify all vortex cores at the precision of given grid resolution without the circulation ambiguity issue described in the previous subsection. Fig. 6 compares the vortex core identification results, where in Fig. 6(b), a complete vortex identification results can be obtained using our proposed strategy. Also note that such a strategy, which is local in computation, can be efficiently parallelized by our GPU implementation.
5.4 Vortex Visualization
After vortex identification, we obtain a binary field to indicate whether each grid point is a vortex point or not, and we perform vortex visualization based on such a field. Note that since we compute circulation with a loop size of two grid cells for vortex identification, such a binary field only indicates whether there are vortex cores inside the loop, and there is almost no chance for a grid point to be exactly the vortex core. On the other hand, vortex cores are infinitesimally small, but visualizing them requires structures with finite size, usually using a tube structure. Thus, locating exactly the vortex core position is not necessary for visualization, and the binary field already forms a good candidate set for vortex tube generationï¼ which can be used for visualization based on isosurface extraction and rendering.
However, if we directly use this binary field to generate the vortex tube, severe spatial aliasing artifacts (note that such aliasing is due to insufficient spatial sampling and is different from those noiselike waves in Section 4.3) would occur because the binary field is like a rectangle function which is discontinuous across the edges. To remove aliasing, we apply a Laplaciandiffusion filter with several iterations to smooth the binary field, which effectively removes the high frequency components beyond the sampling rate. Hence, we can generate the vortex tubes by extracting isosurfaces from such a smoothed field. In our visualizations, we construct the vortex tube isosurfaces based on an isovalue (ranged ) that relates to the radius of the vortex tube, and we set it as , but it can be varied to make the tube thicker or thinner, but no thicker than several cells depending on the degree of smoothing.
Note that the purpose of such a smoothing is to remove aliasing artifacts arisen from the sampling of high frequency signals. Since the smoothing filter is isotropic, the vortex tube generated from the smoothed field still represents the vortex core structure. However, smoothing may introduce uncertainty to visualization, e.g., very close vortex tubes may be merged into one vortex tube. However, these uncertainties are restricted in scales which are roughly smaller than the size of the smoothing kernel. To maintain faithful vortex structure, the above smoothing cannot be excessive. To control the degree of smoothing, we control the kernel size of the Laplaciandiffusion filter as well as the number of iterations, which can be finetuned empirically.
6 Visual Analysis
Once we identify the vortex cores and extract the vortex tubes around, we can render the vortex tubes in order to visualize the superfluid vortex structures and perform visual analysis to answer some important questions with our physicist collaborators. In addition, this research work also serves to help our collaborators explore various interesting phenomena and properties of superfluid vortices. In summary, our collaborators are particularly interested in the following questions on superfluids:

Whether the underlying superfluid model conforms to the real experiments and describes the true physics.

How the superfluid vortices distribute when they are under constant rotation, in both 2D and 3D.

Whether the Feynman relation is satisfied.

How superfluid vortices reconnect to form new vortices.

How waves (including the helical Kelvin waves) are generated and transmitted through vortexvortex interactions.

How vortex tangles are formed and how they behave, etc.
These questions cannot be directly answered by traditional analysis with the original model equation, even with the help of numerical simulations. In the following, we will visualize superfluid vortices and analyze their properties for these questions. Note that the supplementary video presents animations for the resulting vortex structures in superfluids, for both steady and unsteady (dynamic) scenarios.
6.1 Steady Vortex Dynamics
For superfluids in constant rotation, steady states can be reached with stable vortex distributions with appropriate parameter settings. Hence, we can study both 2D and 3D steady vortex structures in the presence of vortex dynamics.
Vortex lattice formation in 2D. The 2D simulation of superfluids in constant rotation can be directly visualized using the scalar field magnitude. Through the visualization, we can then observe regularlydistributed vortex lattice structures, where the vortices are completely isolated and located repulsively to maintain low energy, see Fig. 7. In detail, we maintain the physical validity of the simulation by imposing a constraint , so that the rotation does not violate the limit of the speed of light. Moreover, it is also important to ensure a sufficiently large in the initial condition, so that we can maintain the vortex lattice when the field converges to a finite .
Verification by the Feynman relation. By increasing the speed of rotation in the system, the number of vortices increases accordingly, see some of the simulation results in Fig. 7. Note that we also present the corresponding phase fields, revealing that the vortices always coincide with the phase singularities. By automatically detecting the isolated vortex cores in our results with progressively increased rotations, we can obtain a linear relation between the number of vortex cores () and , see Fig. 8, showing that our simulation results approximately conform to the Feynman relation [20].
Comparison to real experiments. Similar vortex lattice structures have already been observed in real experiments with BEC [9], which can be considered as an experimental results of lowspeed limit of the superfluid flow. Fig. 9 shows a sidebyside comparison between the real BEC experimental result (Fig. 9 (a)) and our simulation result (Fig. 9 (b)). It is clear that the spatial distribution of the vortex lattice structure is quite similar between the real experiment and our simulation (see the triangle shapes formed by the vortex cores, and the repeated hexagonal structure of nearby vortex cores for each vortex core). However, the two results cannot be made exactly matchable since it is difficult to use the same conditions and system parameters for both experiment and simulation. Nevertheless, the close similarity of the vortex lattice distribution gives strong support for the correctness of the simulation as well as the underlying physical model.
Vortex ring formation in 3D. In 3D rotating superfluids, the trapping potential is an oblate Gaussian function, such that the vortices are completely confined within the simulation grid, where the vortex ring structures are formed, see Fig. 10 for the result, and Fig. 1 for the snapshots of the related simulation process. In Fig. 1, since the whole magnitude of the velocity is almost the same on an isosurface, we use the xcomponent of velocity to colorcode the surface, which is helpful to infer the orientation of the vortex lines, e.g., the pairing of antiparallel vortex lines, as marked by the red box in Fig. 11.
Furthermore, it is interesting to see that the resulting steady vortex ring structures are nearly symmetric about the center of rotation, but bending to form closed loop structures. Note that experimental physicists have been working hard to produce these 3D structures in real experiments, but without success so far. To the best of our knowledge, we are the first in the literature to simulate these vortex ring 3D structures in rotating superfluids.
6.2 Unsteady Vortex Dynamics
Superfluids without constant rotation will not have steady solutions. The vortices will always change their forms over time, resulting in very complicated vortex structures. Here, we study the properties for both global and local vortex dynamics, where we set the system parameters as: domain size , , , , time step , and resolution is . In particular, we study two important local vortex dynamics that form the basis of the global vortex dynamics: the vortexvortex reconnection and the Kelvin wave formation and transmission. In addition, we also study vortex tangle and the infinite leapfrogging phenomenon, which is a very special vortex dynamics.
Local vortex reconnection. Although global vortex dynamics (including vortex tangle) in superfluids is rather complicated, it can be decomposed into several characteristic basic local dynamics. One of them is the local vortex reconnection illustrated in Fig. 12(a). From the figure, it is clear that when two vortex strings are close to each other at certain velocity, they may collide. The result of the collision will be a reorganization of the vortex string topology by combining the old vortex strings and reproducing the new ones. Such reconnection phenomenon has been observed experimentally in superfluid liquid Helium [48], see results in Fig. 13; we can see that our simulated and visualized results can produce similar dynamics, as compared to Fig. 12(a).
Note that like the comparison for vortex lattice formation, it is difficult to have exactly the same sidebyside comparison for the vortex reconnection between experiment and simulation since it is hard to have the same system conditions and parameter settings. However, we can still make some comparison and find out the similarity between their spatial topological structures, which gives strong evidence that the model equation we use and simulation methods we develop conform to the real physics.
Kelvin wave. Another important basic local dynamics is the Kelvin wave formation and transmission, which can be easily triggered and observed when vortex structures reconnect. Here, we demonstrate this phenomenon with Fig. 12(b), where a vortex ring approaches a vortex line. After colliding with the vortex line, the vortex ring will be absorbed and a proportion of its energy will be sent out in the form of helical Kelvin waves from the collision points to the two sides of the vortex line. Finally, due to the periodic boundary condition, the waves will collide with its two wavefronts along the vortex line and a new vortex ring will be “split out” from the vortex line but with smaller ring size, revealing that Kelvin wave is actually a driving force for successive reconnection events in superfluids.
Global vortex dynamics. We produce global vortex dynamics in a nonrotating system using initial vortex lines taken from intermediate solutions of a 3D rotating superfluid simulation. Without rotation, the vortices will not remain steady and will evolve freely over time, see Fig. 11. Long vortex lines will be further split by vortex reconnections into shorter vortex lines, and as time proceeds, merging process will also happen to reconnect shorter vortex lines into longer ones. Such split and merge may exist simultaneously in the domain, and the whole process iterates over time to produce more complicated vortex structures with irregular topology, see Fig. 11(b).
It is important to note that during each vortex reconnection in global vortex dynamics, some energy will be released out from the vortices in the form of waves, see Fig. 15 for a corresponding zoomin view produced by using direct volume rendering. It is such energy release that decays the vortices in the whole system, see again Fig. 11 (a & b). Hence, certain amount of energy is needed to add into the system to maintain the overall vortex structures during the reconnections.
Vortex tangle. Using a few vortex lines taken from a 3D rotating superfluid cannot provide sufficiently dense vortex strings to achieve the state of vortex tangle. As described earlier in Section 4.2, if we inject vortex rings randomly from the domain boundary, we can effectively create dense amount of vortex reconnections and generate the state of vortex tangle, see Fig. 14. Note that the vortex tangle is decaying quantum turbulence, where the vortices are reduced by energy release during the massive reconnection processes.
Leapfrogging phenomena. Lastly, we study one interesting vortex dynamics phenomenon called “leapfrogging.” Leapfrogging phenomena have been observed in classical fluids, and can also be simulated [49]. However, it has not yet been observed in superfluids. By initializing with two coplanar vortex rings with the same orientation, we can simulate stable leapfrogging phenomena, see Fig. 12(c). Please refer to the supplementary video for the corresponding animated simulations and visualizations.
Our results appear similar to the leapfrogging in classical fluids, where the two rings alternatively move forward through the other with varying radius. However, due to vorticity diffusion and turbulence, leapfrogging in classical fluids can be destroyed easily, and cannot last for a long time. In the case of superfluids, due to zero viscosity, vortices do not merge or dissipate, thus leading to more stable and longlasting leapfrogging phenomena. If we consider periodic boundaries in the spatial domain, a twovortexring leapfrogging system can proceed indefinitely.
6.3 Implementation and Discussion
Implementation details. We implemented and ran our superfluid simulation and vortex structure visualization on a workstation with an Intel Xeon E52630 v3 CPU, 32GB system memory and an NVIDIA GTX TITAN X GPU. The solvers and renderers are all written in C++ and CUDA. The simulations and visualizations are accelerated by the GPU due to its nice parallelism feature.
The average time for producing the simulations is 9 to 10 seconds for 10000 iterations in 2D with a grid resolution of in polar coordinates. The average time is around 0.1 second per iteration in 3D with a grid resolution of in Cartesian coordinates. The 2D renderings are straightforward, but it takes even more time than simulation because we use Matlab for visualization, and call our CUDA kernels through the MEX interface of Matlab. For 3D simulation, the vortex identification is very fast; it takes a few seconds to obtain the circulation field and volumetrically render the structure at a grid resolution of , such that the 3D simulation can be previewed without much performance drop. For higher quality, we build our offline volume raytracer for the generated vortex tube isosurfaces, where the rendering cost is as high as a few minutes per frame.
Initialization. In general, we require the initialization to contain large enough circulation, but we can use different distributions of the field, which does not influence the final vortex phenomena, for no matter steady or unsteady superfluid vortex dynamics. For steady solutions on rotating superfluids, the solution always converges to the steady vortex lattice (or vortex ring in 3D) structure, while for unsteady solutions, we are interested in either local dynamics (such as reconnections and Kelvin waves) or vortex tangle. In either case, although some initial conditions are artificial, the underlying motions are strictly under the corresponding physical law by the model, which can result in believable phenomena that could be observed in practice, and our method can produce some phenomena that can be truly observed in real experiments.
Accuracy analysis. Since we use highorder numerical methods in both 2D and 3D, the original discretization of the model equation is accurate enough for visualization purposeï¼ especially if we use high resolution grid. However, due to strong waves which are eventually of high frequency, we employ a temporal filtering method to make the simulation data clear enough for vortex identification. This is based on an assumption that vortices and compression waves are completely separable in time frequency. However, if there are vortices that vary as fast as waves, these vortices will be removed, which introduces a certain degree of uncertainty in our final results. However, by our experiments and observations, such phenomenon is very rare in most cases, as the waves are finally of extremely high frequency everywhere.
On the other hand, during vortex visualization, we employ circulation with a loop size of two grid cells to identify the vortex cores. Here, we also make an assumption that all the vortex cores are of the same scale, which is true for superfluids and in contrast to classical fluids where vortices can have varying scales. With circulationbased approach, we are uncertain about the exact location of the vortex core, where some smallscale (scale smaller than two grid cells) variation of the vortex structure is missing. However, we are sure that the vortex tubes we generated and visualized can truly reflect the large scale (scale larger than two grid cells) structure of the vortex variations since the vortex cores are completely included inside the circulation loops.
Since largescale structures are more meaningful than small fluctuations for most of the cases, the visualization we generated can be useful enough for scientific comparison and visual analysis, as done in Sections 6.1 & 6.2. For more accurate analysis involving smallerscale structures, we only need to increase the grid resolution, and our vortex identification and visualization methods can automatically obtain smaller scale structures since the circulation loop is always of two grid cell size.
Feedback from physicists. Our visualization results have been given to two research groups of superfluid scientists for feedback. One of such groups are our collaborators at Nanyang Technological University, led by Professor Emeritus Kerson Huang from MIT, who are all theoretical physicists. They felt that our visualizations and animations have more intuitively helped them observe the superfluid structure from the physical model, which has given them deeper understanding, especially in 3D, and find out more interesting and useful properties for superfluid vortices, as evidenced by a recent paper published by the group [6]. Another group of scientists are from the research group led by Professor Daniel Lathrop at University of Maryland, who are all experimental scientists on superfluids. They gave us very positive feedback for our simulation and visualization results, with a particular interest on Kelvin wave simulation and visualization. With our results, both groups of scientists can easily analyze the steady and unsteady behavior of superfluid vortices on a computer, and even discover new vortex structures more intuitively by trying different sets of system parameters. This can provide some guidance in experimental observations, which indicates the significance of our simulation and visualization work to help support research in fundamental science.
Limitations. Since the vortex cores exist as singularities in the phase field of the NLKG equation and are infinitesimally small, we cannot exactly locate them and visualize them with thinenough structures at proper grid resolution. In addition, our visualization of superfluid vortices is done by extracting isosurfaces of a smoothed circulation field, and the vortex tube size is larger than the grid spacing. Thus, very dense vortex tangles cannot be visualized clearly. Increasing the grid resolution can alleviate such problem and make the tube thinner to better reflect the vortex core geometry, but at the expense of more computing resource and time. Finally, the Laplacian smoothing on the circulation field may also remove some vortices given the limited grid resolution.
7 Conclusion
In this work, we aim for highquality superfluid visualization and to help our collaborators visually study superfluid vortices. To achieve this goal, we solve the underlying nonlinear KleinGordon equation to prepare proper datasets for visualization, and develop an effective vortex core identification method based on a particular loop selection model and an orthogonalplane strategy, which are local and efficient for parallel implementation on GPU. Visual analysis of superfluid vortex structures is conducted with domain experts by considering both steady and unsteady scenarios, where we observe vortex lattice and ring structures in 2D & 3D, approximate the Feynman relation in 2D, and produce local vortex reconnections, Kelvin waves, vortex tangles and leapfrogging phenomena in 3D. These results helped physicists to better understand the stationary as well as dynamic properties of superfluid vortices. For verification purpose, we compared our visualization with experimental results and the Feynman relation.
As a future work, we would like to explore the use of adaptive grids to produce higher resolution vortex core regions to more accurately locate vortex cores and reveal finerscale structures. Such result will be useful for comparative studies with real experiments. Moreover, we plan to explore multiscale visualization, where we visualize smallscale vortex core structures on top of largescale structures for studying the relation and interaction between vortex core structures across different scales.
Acknowledgment
The authors would like to thank all anonymous reviewers for their constructive comments. This work is supported by National Science Foundation of China (NSFC)  Outstanding Youth Foundation (Grant No. 61502305), ShanghaiTech University startup funding, Guangdong highlevel personnel of special support program (Grant No. 2016TQ03X319), and the CUHK strategic recruitment fund and direct grant (4055061).
Footnotes
 BEC is a state of matter of a dilute gas of bosons cooled to temperatures very close to absolute zero.
 vortex tubes are basic vortex structure in 3D superfluid visualization
References
 D. Tilley and J. Tilley, Superfluidity and Superconductivity, ser. Halsted press book. Wiley, 1974.
 P. Kapitza, “Viscosity of Liquid Helium below the Point,” Nature, vol. 141, p. 74, Jan 1938.
 J. F. Annett, Superconductivity, superfluids and condensates, ser. Oxford master series in physics. Oxford: Oxford Univ. Press, 2004.
 R. Donnelly, Quantized Vortices in Helium II, ser. Cambridge Studies in American Literature and Culture. Cambridge University Press, 1991, no. v. 2.
 K. Huang, C. Xiong, and X. Zhao, “ScalarField Theory of Dark Matter,” International Journal of Modern Physics A, vol. 29, no. 13, p. 1450074, 2014.
 C. Xiong, M. R. R. Good, Y. Guo, X. Liu, and K. Huang, “Relativistic superfluidity and vorticity from the nonlinear KleinGordon equation,” Physical Review D, vol. 90, p. 125019, Dec 2014.
 L. Tisza, “Transport phenomena in helium ii,” Nature, vol. 141, no. 3577, 1938.
 L. Landau, “The theory of superfluidity of helium ii,” Journal of Physics USSR, vol. 5, no. 71, 1941.
 J. R. AboShaeer, C. Raman, J. M. Vogels, and W. Ketterle, “Observation of Vortex Lattices in BoseEinstein Condensates,” Science, vol. 292, no. 5516, pp. 476–479, 2001.
 G. P. Bewley, D. P. Lathrop, and K. R. Sreenivasan, “Superfluid Helium: Visualization of quantized vortices,” Nature, vol. 441, pp. 588 – 588, June 2006.
 K. Huang, H.B. Low, and R.S. Tung, “Scalar field cosmology ii: Superfluidity, quantum turbulence, and inflation,” International Journal of Modern Physics A, vol. 27, no. 26, p. 1250154, 2012.
 R. J. Adrian, “Hairpin Vortex Organization in Wall Turbulence,” Physics of Fluids, vol. 19, no. 4, 2007.
 V. Suponitsky, J. Cohen, and P. Z. Baryoseph, “The Generation of Streaks and Hairpin Vortices from a Localized Vortex Disturbance Embedded in Unbounded Uniform Shear Flow,” Journal of Fluid Mechanics, vol. 535, pp. 65–100, 7 2005.
 T. Schafhitzel, K. Baysal, M. Vaaraniemi, U. Rist, and D. Weiskopf, “Visualizing the Evolution and Interaction of Vortices and Shear Layers in TimeDependent 3D Flow,” IEEE Trans. on Visualization and Computer Graphics, vol. 17, no. 4, pp. 412–425, Apr. 2011.
 C. Koehler, T. Wischgoll, H. Dong, and Z. Gaston, “Vortex Visualization in Ultra Low Reynolds Number Insect Flight,” IEEE Trans. on Visualization and Computer Graphics, vol. 17, no. 12, pp. 2071–2079, Dec. 2011.
 M. Treib, K. Bürger, F. Reichl, C. Meneveau, A. Szalay, and R. Westermann, “Turbulence Visualization at the Terascale on Desktop PCs,” IEEE Trans. on Visualization and Computer Graphics, vol. 18, no. 12, pp. 2169–2177, Dec. 2012.
 M. Hummel, H. Obermaier, C. Garth, and K. I. Joy, “Comparative Visual Analysis of Lagrangian Transport in CFD Ensembles,” IEEE Trans. on Visualization and Computer Graphics, vol. 19, no. 12, pp. 2743–2752, Dec. 2013.
 J. Clyne, P. Mininni, and A. Norton, “PhysicallyBased Feature Tracking for CFD Data,” IEEE Trans. on Visualization and Computer Graphics, vol. 19, no. 6, pp. 1020–1033, June 2013.
 S. Shafii, H. Obermaier, R. Linn, E. Koo, M. Hlawitschka, C. Garth, B. Hamann, and K. I. Joy, “Visualization and Analysis of VortexTurbine Intersections in Wind Farms,” IEEE Trans. on Visualization and Computer Graphics, vol. 19, no. 9, pp. 1579–1591, Sept. 2013.
 R. P. Feynman, “Application of quantum mechanics to liquid helium,” Progress in Low Temperature Physics, ed. C. J. Gorter, vol. I, 1955.
 M. Tsubota, K. Kasamatsu, and M. Ueda, “Vortex lattice formation in a rotating BoseEinstein condensate,” Physical Review A, vol. 65, p. 023603, Jan 2002.
 K. Kasamatsu, M. Machida, N. Sasa, and M. Tsubota, “Threedimensional dynamics of vortexlattice formation in BoseEinstein condensates,” Physical Review A, vol. 71, p. 063616, Jun 2005.
 S. Zuccher, M. Caliari, A. W. Baggaley, and C. F. Barenghi, “Quantum Vortex Reconnections,” Physics of Fluids, vol. 24, no. 12, 2012.
 R. M. Caplan, “Study of Vortex Ring Dynamics in the Nonlinear Schrödinger Equation Utilizing GPUAccelerated HighOrder Compact Numerical Integrators,” Ph.D. dissertation, Claremont Graduate University, 2012.
 J. C. R. Hunt, A. A. Wray, and P. Moin, “Eddies, Stream, and Convergence Zones in Turbulent Flows,” Center For Turbulence Research, vol. Report CTRS88, pp. 193–208, 1988.
 J. Jeong and F. Hussain, “On the identification of a vortex,” Journal of Fluid Mechanics, vol. 285, pp. 69–94, 2 1995.
 M. Jiang, R. Machiraju, and D. Thompson, “Detection and visualization of vortices,” in The Visualization Handbook. Academic Press, 2005, pp. 295–309.
 S. Weimann, U. Pinkall, and P. Schröder, “Smoke Rings from Smoke,” ACM Trans. on Graphics, vol. 33, no. 4, pp. 140:1–140:8, Jul. 2014.
 A. Wiebel, X. Tricoche, D. Schneider, H. Jänicke, and G. Scheuermann, “Generalized Streak Lines: Analysis and Visualization of Boundary Induced Vortices,” IEEE Trans. on Visualization and Computer Graphics, vol. 13, no. 6, pp. 1735–1742, Nov./Dec. 2007.
 H. Krishnan, C. Garth, and K. I. Joy, “Time and Streak Surfaces for Flow Visualization in Large TimeVarying Data Sets,” IEEE Trans. on Visualization and Computer Graphics, vol. 15, no. 6, pp. 1267–1274, Nov./Dec. 2009.
 H. Yu, C. Wang, C.K. Shene, and J. H. Chen, “Hierarchical Streamline Bundles,” IEEE Trans. on Visualization and Computer Graphics, vol. 18, no. 8, pp. 1353–1367, Aug. 2012.
 A. Chaudhuri, T.Y. Lee, H.W. Shen, and R. Wenger, “Exploring Flow Fields Using SpaceFilling Analysis of Streamlines,” IEEE Trans. on Visualization and Computer Graphics, vol. 20, no. 10, pp. 1392–1404, Oct. 2014.
 M. JankunKelly, M. Jiang, D. Thompson, and R. Machiraju, “Vortex visualization for practical engineering applications,” IEEE Trans. on Visualization and Computer Graphics, vol. 12, no. 5, pp. 957–964, Sept 2006.
 D. Schneider, A. Wiebel, H. Carr, M. Hlawitschka, and G. Scheuermann, “Interactive Comparison of Scalar Fields Based on Largest Contours with Applications to Flow Visualization,” IEEE Trans. on Visualization and Computer Graphics, vol. 14, no. 6, pp. 1475–1482, Nov./Dec. 2008.
 F. Sadlo, R. Peikert, and M. Sick, “Visualization tools for vorticity transport analysis in incompressible flow,” IEEE Trans. on Visualization and Computer Graphics, vol. 12, no. 5, pp. 949–956, Sept 2006.
 D. Laney, P.T. Bremer, A. Mascarenhas, P. Miller, and V. Pascucci, “Understanding the Structure of the Turbulent Mixing Layer in Hydrodynamic Instabilities,” IEEE Trans. on Visualization and Computer Graphics, vol. 12, no. 5, pp. 1053–1060, Sept./Oct. 2006.
 T. Weinkauf, J. Sahner, H. Theisel, and H.C. Hege, “Cores of Swirling Particle Motion in Unsteady Flows,” IEEE Trans. on Visualization and Computer Graphics, vol. 13, no. 6, pp. 1759–1766, Nov./Dec. 2007.
 A. Helgeland, B. A. P. Reif, Ø. Andreassen, and C. E. Wasberg, “Visualization of Vorticity and Vortices in WallBounded Turbulent Flows,” IEEE Trans. on Visualization and Computer Graphics, vol. 13, no. 5, pp. 1055–1066, Sept./Oct. 2007.
 G. Johnson, V. Calo, and K. Gaither, “Interactive visualization and analysis of transitional flow,” Visualization and Computer Graphics, IEEE Trans. on, vol. 14, no. 6, pp. 1420–1427, Nov 2008.
 B. Soni, D. Thompson, and R. Machiraju, “Visualizing Particle/Flow Structure Interactions in the Small Bronchial Tubes,” IEEE Trans. on Visualization and Computer Graphics, vol. 14, no. 6, pp. 1412–1419, Nov./Dec. 2008.
 S. Zachow, P. Muigg, T. Hildebrandt, H. Doleisch, and H.C. Hege, “Visual Exploration of Nasal Airflow,” IEEE Trans. on Visualization and Computer Graphics, vol. 15, no. 6, pp. 1407–1414, Nov./Dec. 2009.
 B. Köhler, R. Gasteiger, U. Preim, H. Theisel, M. Gutberlet, and B. Preim, “SemiAutomatic Vortex Extraction in 4D PCMRI Cardiac Blood Flow Data Using Line Predicates,” IEEE Trans. on Visualization and Computer Graphics, vol. 19, no. 12, pp. 2773–2782, Dec. 2013.
 S. Born, M. Pfeifle, M. Markl, M. Gutberlet, and G. Scheuermann, “Visual Analysis of Cardiac 4D MRI Blood Flow Using Line Predicates,” IEEE Trans. on Visualization and Computer Graphics, vol. 19, no. 6, pp. 900–912, June 2013.
 H. Guo, C. L. Phillips, T. Peterka, D. Karpeyev, and A. Glatz, “Extracting, tracking, and visualizing magnetic flux vortices in 3d complexvalued superconductor simulation data,” IEEE Trans. on Visualization and Computer Graphics, vol. 22, no. 1, pp. 827–836, Jan 2016.
 J. Anderson, Computational Fluid Dynamics: The Basics with Applications. McGrawHill Education, 1995.
 K. Mohseni and T. Colonius, “Numerical treatment of polar coordinate singularities,” Journal of Computational Physics, vol. 157, no. 2, pp. 787 – 795, 2000.
 A. V. Oppenheim, R. W. Shafer, and J. R. Buck, DiscreteTime Signal Processing, 2nd ed. Upper Saddle River, New Jersey: Prentice Hall, 1998.
 E. Fonda, D. P. Meichle, N. T. Ouellette, S. Hormoz, and D. P. Lathrop, “Direct observation of Kelvin waves excited by quantized vortex reconnection,” Proceedings of the National Academy of Sciences of the United States of America, vol. 111, no. 1, pp. 4647–4734, March 2014.
 N. Riley and D. Stevens, “A note on leapfrogging vortex rings,” Fluid Dynamics Research, vol. 11, no. 5, pp. 235–244, 1993.