Introduction

Many quantum fluids, such as superfluid helium-4 (He II), can be considered as a mixture of two miscible fluid components: an inviscid superfluid and a viscous normal fluid consisting of thermal quasiparticles1. A conspicuous feature of the superfluid is the existence of topological defects in the form of quantized vortices2. In 3D space, these vortices appear as density-depleted thin tubes, each carrying a circulating flow with a fixed circulation κ = h/m, where h is Planck’s constant and m is the mass of the bosons constituting the superfluid2. The motion of quantized vortices is responsible for a wide range of phenomena in quantum-fluid systems, such as the emergence of quantum turbulence in He II and atomic Bose-Einstein condensates3,4,5, the initiation of dissipation in type-II superconductors6, the appearance of glitches in neutron star rotation7,8, and the formation of possible cosmic-string network9. As the vortices move through the normal fluid, a mutual friction between the two fluids can arise due to the scattering of the thermal quasiparticles off the vortex cores10,11,12. Understanding the dynamics of quantized vortices in the presence of the normal fluid is therefore of broad significance. However, despite decades of research13,14,15,16,17,18, the lack of experimental data has left the issue of how to accurately model the mutual friction unsettled.

In the pioneering work of Schwarz13,14, a vortex filament model was developed for studying turbulence in He II. In this model, the quantized vortices are described by zero-thickness filaments that are divided into small segments. A vortex segment with a length Δξ located at s would experience a Magnus force \({{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{M}}}}}}}}}={\rho }_{s}\kappa {{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{L}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{s}}}}}}}}})\Delta \xi\) when its velocity uL differs from the local superfluid velocity us. Here \({{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\) is the unit tangent vector along the filament, and ρs is the superfluid density. Besides, any relative motion between the vortex segment and the normal fluid can result in a mutual friction force as proposed by Schwarz \({{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{sn}}}}}}}}}=[-{\gamma }_{0}{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{L}}}}}}}}}))+{\gamma }_{0}^{{\prime} }{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{L}}}}}}}}})]\Delta \xi\), where γ0 and \({\gamma }_{0}^{{\prime} }\) are temperature-dependent empirical coefficients14. By balancing the two forces, Schwarz obtained the vortex equation of motion (see Methods), which has been extensively employed in past vortex research19,20,21,22.

However, a known limitation of the Schwarz model is that the normal-fluid velocity un is prescribed and there is no back action from the vortices to the normal fluid. To fix this issue, a two-way (2W) model was later developed, where un is solved using the Navier-Stokes equation with an added mutual-friction term that couples to the vortices. This model has allowed researchers to explain puzzling observations in He II turbulence18,23. Nonetheless, it was postulated that the coefficients γ0 and \({\gamma }_{{{{{{{{\rm{0}}}}}}}}}^{{\prime} }\) may not be applicable to individual vortices since they were deduced from measurements where un was averaged over an array of vortices15. Over the past two decades, researchers have strived to calculate the friction coefficients in a self-consistent manner15,16,17. These efforts led to the striking prediction of the triple-vortex-ring structure in He II15,17,24. Recently, Galantucci et al. derived the most refined version of the self-consistent two-way (S2W) model where the mutual friction coefficient can be calculated directly from un without relying on any empirical experimental input17.

These different models render distinct normal-fluid flow structures around the quantized vortices, which affect the vortex motion. As an example, we show in Fig. 1a the calculated normal-fluid velocity field around a quantized vortex ring in quiescent He II using all three models (see Methods for details). Unlike the Schwarz model where un = 0, both the 2W model and the S2W model reveal the existence of two oppositely polarized normal-fluid vortex rings sandwiching the quantized vortex ring. These normal-fluid rings affect the local un experienced by the quantized ring and hence can alter the mutual friction dissipation. However, is this triple-ring structure real? If so, which model better describes the true vortex dynamics? These questions are important but have remained open due to the lack of experimental information. In this work, we provide the long-awaited data to show that the S2W model can better reproduce experimental observations. These findings may have the potential to improve our comprehension of various vortex-involved phenomena in quantum two-fluid systems.

Fig. 1: Modeling and imaging quantized vortex rings in He II.
figure 1

a Calculated normal-fluid velocity field un around a quantized vortex ring in He II. Due to the axial symmetry, we only show un in the y-z plane and the vortex ring above the plane (i.e., the green curve). The normal-fluid vortex rings (reddish half circles) are rendered in the same way as in ref. 18. b Schematic diagram of the experimental setup. c Images showing the D2 particles (white dots) trapped on a moving vortex ring in quiescent He II. The dashed ellipse is a fit to the trapped particles' positions. d Obtained vortex-ring profile with the trapped particles (red dots) at different times. e Comparison of the observed ring radius R(t) evolution with model predictions. The error bars on the R(t) data denote the standard deviation derived from the positional uncertainties of the particles in the ellipse fit. The narrow color-shaded areas depict the variation range of the simulation curves for the particle-doped ring, as the radius a of each trapped particle is adjusted from a − Δa to a + Δa.

Full size image

Results

Visualizing quantized vortex rings

To study the vortex motion, we visualize quantized vortices in He II by decorating them with solidified deuterium (D2) tracer particles25,26. This method has already allowed researchers to gain valuable insights into the properties of tangled vortices27,28,29,30,31. In an earlier work by Bewley and Sreenivasan32, they reported on the visualization of a single vortex ring event. However, as these authors pointed out, many tracers had condensed onto the vortex core, which could alter the core size and hence the ring dynamics. The low vortex-line density in that experiment also limited the chance to observe more vortex rings. Based on the insights gained from this valuable work, we have designed our experimental setup to facilitate improved vortex-ring imaging. As shown in Fig. 1b, we control the vortex-line generation by towing a mesh grid in a plexiglass channel (1.6 × 1.6 × 33 cm3) immersed in a He II bath. The bath temperature can be controlled with an accuracy of 1 mK by regulating the vapor pressure. Following the grid motion, a mixture of D2 gas and 4He gas is injected into the channel at about 30 s delay so that the background flow is weak but vortices with a line density of the order 102 cm−2 still remain33,34. The higher line density can increase the likelihood of observing vortex rings, since the reconnections of intersecting vortex lines in the tangle can give rise to them. The D2 gas forms ice particles (density35: ρp = 202.8 g/cm3) with a mean radius of 1.1 μm as determined from their settling velocities (see Methods and Supplementary Fig. 1). When the D2 particles are close to the vortex cores, they get trapped on the vortices due to a Bernoulli pressure caused by the circulating superfluid36. Through extensive trials, we have figured out the optimal injection conditions to achieve the desired particle number density on the vortices. The particles are then illuminated by a laser sheet (thickness 0.8 mm, see Supplementary Fig. 4) and their positions are recorded at 200 Hz by a video camera placed perpendicular to the laser plane. From time to time, we can see vortex rings propagating within the laser sheet. A collection of representative ring events obtained at different temperatures is included in Supplementary Movie 1. We have also captured videos showing how vortex rings are created by vortex-line reconnections (see Supplementary Movie 2). Note that as temperature rises, the lifetime of vortex rings reduces due to increased dissipation, which ultimately decreases the likelihood of capturing high-quality ring events. In the following sections, we present our best data obtained at 1.65 K for comparison with model simulations.

Data analysis and model comparison

To extract useful information on vortex-ring propagation, we focus on analyzing selected events where the rings are decorated by discrete D2 particles and move in He II with negligible background flows. A good example is shown in Fig. 1c where the ring moves downward carrying nine D2 particles (see Supplementary Movie 3). We first use a feature-point tracking routine37 to determine the positions of the trapped particles in each image. Then, the particle positions are fitted with an ellipse. This fitting, which requires at least 5 particles on the ring, allows us to determine both the ring radius R and the orientation of the ring plane (see Methods). Figure 1d shows the extracted ring profile with the trapped particles at different times. The ring shrinks due to the mutual friction dissipation, which leads to an acceleration of its self-induced motion36. By analyzing the angular positions of the particles on the ring over time, we confirm that these particles exhibit little movement along the circumference of the ring (see Methods and Supplementary Fig. 2), a puzzle that has been actively discussed in literature38. In Fig. 1e, we show the obtained ring radius R(t). These data are shown at every five image frames instead of at 200 Hz for better visibility of the error bars. For comparison, we also include the simulated R(t) for a bare vortex ring in quiescent He II with the same initial radius using all three models. The result suggests that the S2W model renders the best agreement with the data.

However, it is well known that the trapped particles can result in additional forces on the vortex core39,40,41, which needs to be considered when simulating the ring’s motion. Following Mineda et al.39 (see Methods), we consider the Stokes drag1 fD = − 6πaμn(uL − un), the gravitational force, and the inertial effect of each trapped particle on the ring. Here μn is the He II dynamic viscosity and a is the particle radius. To evaluate a, we first develop a correlation between the particle’s brightness and its radius by comparing the distributions of these two quantities (see Methods and Supplementary Fig. 3). We then examine the time-averaged brightness of each trapped particle and calculate its radius using the correlation. The obtained radiuses are listed in the Supplementary Table 1. With this information, we can re-calculated R(t) using the three models (see Fig. 1e). Due to the additional Stokes drag, the ring shrinks faster in all three models. Obviously, the Schwarz model overestimates the dissipation and can be rejected. But it becomes less clear whether the S2W model still describes the data better than the 2W model. To make a reliable assessment on these two models, it is imperative to analyze rings with a minimal number of trapped particles, since the uncertainties in the particle positions and sizes could affect both the R(t) data and the simulation curves.

Luckily, we have recorded several unique events where the rings are decorated by only two D2 particles (see Supplementary Movie 4). For these events, the estimated Stokes drag and the gravitational force are only a few percent of the mutual friction (see Methods). Figure 2a shows our best example, where two particles P1 and P2 move in sync while approaching each other due to the shrinkage of the vortex ring. We can measure the separation distance dp(t) between the two particles and their centroid velocity \({u}_{{{{{{{{\rm{p}}}}}}}}}(t)=|\frac{1}{2}({{{{{{{{\bf{u}}}}}}}}}_{1}+{{{{{{{{\bf{u}}}}}}}}}_{2})|\). However, as illustrated in Fig. 2b, dp(t) in general does not equal the vortex-ring diameter 2R(t), and up(t) can differ from the actual ring velocity u(t) since a projection angle θ may exist between the ring’s propagation direction and the laser plane. In order to utilize the experimental data for model comparison, we adopt the following procedures. First, we assume an initial ring radius R(0) and calculate the evolution of the ring’s radius R(t) and velocity u(t) using both the 2W and the S2W models. Next, we evaluate two projection parameters c1 = dp(0)/2R(0) and \({c}_{2}=\cos (\theta )={u}_{{{{{{{{\rm{p}}}}}}}}}(0)/u(0)\). These two parameters remain nearly constant because: 1) the particles do not slide along the vortex core as we learned from the study of rings with 5 or more trapped particles; and 2) the centroid of P1 and P2 moves in a straight path, suggesting a constant projection angle. Using c1 and c2, we can then compare c1R(t) and c2u(t) directly with the experimental data dp(t)/2 and up(t). Finally, we vary R(0) to see which model can render results that simultaneously match dp(t)/2 and up(t).

Fig. 2: Analysis of a vortex ring with two trapped particles.
figure 2

a Images showing two trapped particles (circled) moving in He II at 1.65 K. b A schematic explaining the concept of the projection parameters c1 and c2. c Variation of the corrected brightness I(t) of the two trapped particles. The solid curve shows the prediction of the S2W model that gives the best fit to the ring’s radius and velocity data in d. d Comparison of the projected ring radius and velocity data with model simulations. The tuning range of the simulation curve is set by the constraints on R(0) and u(0), as discussed in the text. The error bars denote the standard deviation associated with the data due to the position uncertainties of the two trapped particles (see Methods).

Full size image

In this analysis, there are a few constraints on the range of R(0) that we can explore. First, R(0) ≥ dp(0)/2 since the two particles cannot be separated by more than the diameter of the ring. Second, u(0) ≥ up(0) due to the projection, which sets an upper limit of R(0) because u(0) drops as R(0) increases. The last constraint comes from the observed particle brightness I after correcting the z-dependance of the laser intensity (see Methods). As shown in Fig. 2c, I for either particles only drops by about 20% during the ring’s propagation. Based on the cross-sectional profile of the laser sheet (see Supplementary Fig. 4), we estimate that the ring can move by at most 0.2 mm perpendicular to the laser plane. This sets an upper limit of the projection angle θ, which constrains u(0) and hence R(0). In Fig. 2d, we show the calculated c1R(t) and c2u(t) using the 2W and the S2W models while R(0) is varied in the range set by all the constraints. Clearly, the experimental data are outside the tuning range of the 2W model. On the other hand, we find that the S2W model can nicely reproduce both dp(t)/2 and up(t) data at R(0) = 140.8 μm. This optimal R(0) happens to be close to dp(0)/2 = 140.6 μm, suggesting that the two particles were located nearly across the diameter of the vortex ring. In conclusion, our 2-particle ring analysis indicates that the S2W model has better fidelity than the 2W model.

Other intriguing observations

Besides model testing, we would also like to report some intriguing phenomena that have not been reported in literature. The first phenomenon is that sometimes we see vortex rings heavily doped with D2 particles can spontaneously flip to the downward direction. A collection of such events is included in Supplementary Movie 5. In Fig. 3a, we superimpose the images of a representative ring taken at different t to show how the ring changes its direction while it shrinks. This phenomenon can be understood by noting that the vortex ring carries a momentum36 \({{{{{{{\bf{P}}}}}}}}(t)={\rho }_{{{{{{{{\rm{s}}}}}}}}}\kappa \pi R{(t)}^{2}\hat{{{{{{{{\bf{n}}}}}}}}}\), where \(\hat{{{{{{{{\bf{n}}}}}}}}}\) is the unit vector normal to the ring plane pointing in the direction of the ring’s motion. The mutual friction and the Stokes drag constantly reduce the ring’s momentum, resulting in the shrinkage of the ring. On the other hand, the gravitational force from the trapped particles continuously generates momentum in the downward direction, which forces the ring to flip downward. To test this understanding, we have conducted simulations using the S2W model. For a heavily doped ring, the exact number N and the radiuses of the trapped particles are hard to determine. Instead, we assume the same radius a for all the trapped particles and adjust both a and N in our trials. For the ring trajectory presented in Fig. 3a, we find that it can be qualitatively reproduced with N = 36 and a = 4.9 μm, as shown in Fig. 3b. Similar flipping can be also achieved with fewer but heavier trapped particles.

Fig. 3: Other intriguing observations of the vortex rings.
figure 3

a A superimposed image showing a heavily doped vortex ring gradually flipping downward while it shrinks. b Simulated motion of a vortex ring with the same initial profile carrying 36 D2 particles (4.9 μm in radius) using the S2W model. c Images showing how the D2 particles trapped on a vortex ring eventually form a cluster that falls freely in He II. d Measured centroid velocity of the trapped D2 particles shown in c. The error bars denote the standard deviation associated with the data due to the uncertainties of the centroid position (see Methods).

Full size image

Another intriguing observation is related to the destiny of the particles on the vortex rings. As a ring shrinks, we often see that the trapped particles form a cluster whose high-speed motion changes abruptly to slowly falling in He II (see Supplementary Movie 6). Figure 3c shows an event where the ring plane is nearly perpendicular to the laser plane. Nonetheless, we can measure the centroid velocity up(t) of the particles. As shown in Fig. 3d, up(t) increases drastically as the ring shrinks. At t = 3.375 s, the trapped particles aggregate to a single cluster and up(t) suddenly drops to the expected settling velocity of about 0.1 mm/s. Our interpretation of this phenomenon is that as the ring shrinks, its velocity relative to the normal fluid becomes so large such that the Stokes drag can pull the trapped particles off the vortex core. Subsequently, the bare ring moves away and diminishes, while the left-behind particles form a cluster that decelerates rapidly to the settling velocity due to the Stokes drag. To test this hypothesis, we estimate the maximum trapping force on a particle from the vortex core as42 fv ≃ ρsκ2/3π and compare it with the Stokes drag fD ≃ 6 πaμnup. For particles with a mean radius a ≃ 1 μm, fD becomes greater than fv when up reaches a threshold value of 5.1 mm/s. This threshold up is close to the observed maximum up in Fig. 3d, which provides a strong support to our interpretation.

Discussion

The results that we have presented provide the first evidence suggesting that the S2W model is more consistent with observed vortex dynamics in He II in comparison to alternative models. This study may stimulate future research in two possible directions. First, the S2W model does not rely on empirical experimental inputs and therefore can be readily adapted for other quantum two-fluid systems, such as BECs43,44, superfluid neutron stars7,8,45, and gravity-mapped holographic superfluid46,47. An accurate evaluation of the mutual friction is particularly important for processes that involve rapid motion of the quantized vortices, such as vortex reconnections, and pinning and depinning of vortices on solid boundaries. The latter process is the key for understanding glitches in neutron star rotation45. The second direction is to examine how the implementation of the S2W model may alter our existing knowledge on quantum turbulence (QT) induced by a chaotic tangle of quantized vortices. For instance, an important topic in QT research is counterflow turbulence where the mutual friction exists at all length scales48,49. Our knowledge on the vortex-tangle properties14,19, disturbances in the normal fluid18,50, and the effect of the mutual friction on the mean-velocity profile23,51,52 may subject to change with future S2W-model simulations.

Methods

Numerical models

Schwarz model

In the framework of Schwarz’s vortex filament model14, all the quantized vortex lines are represented by zero-thickness filaments. The position vector of a filament can be written in the parametric form s = s(ξ, t), where ξ denotes the arc length along the filament. In the presence of the viscous normal fluid, a short segment Δξ of a vortex filament located at s would experience two forces, i.e., the Magnus force \({{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{M}}}}}}}}}={\rho }_{{{{{{{{\rm{s}}}}}}}}}\kappa {{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{L}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{s}}}}}}}}})\Delta \xi\) and the mutual friction force \({{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{sn}}}}}}}}}=[-{\gamma }_{0}{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{L}}}}}}}}}))+{\gamma }_{0}^{{\prime} }{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{L}}}}}}}}})]\Delta \xi\). By balancing these two forces, the velocity of this segment uL = ds/dt can be derived as:

$$d{{{{{{{\bf{s}}}}}}}}/dt={{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{s}}}}}}}}}+\alpha {{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{s}}}}}}}}})-{\alpha }^{{\prime} }{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times [{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{s}}}}}}}}})],$$
(1)

where the coefficients α and \({\alpha }^{{\prime} }\) depend on the empirical mutual friction coefficients γ0 and \({\gamma }_{0}^{{\prime} }\), whose values have been tabulated36. While the normal-fluid velocity un is prescribed, the local superfluid velocity us(s, t) is evaluated as the sum of the background flow velocity us0 and the velocity uin induced at s by all the vortices, which can be calculated using the full Boit-Savart integral19:

$${{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{in}}}}}}}}}({{{{{{{\bf{s}}}}}}}},\, t)=\frac{\kappa }{4\pi }\int\frac{({{{{{{{{\bf{s}}}}}}}}}_{{{{{{{{\bf{1}}}}}}}}}-{{{{{{{\bf{s}}}}}}}})\times d{{{{{{{{\bf{s}}}}}}}}}_{{{{{{{{\bf{1}}}}}}}}}}{|{{{{{{{{\bf{s}}}}}}}}}_{{{{{{{{\bf{1}}}}}}}}}-{{{{{{{\bf{s}}}}}}}}{|}^{3}},$$
(2)

where the integration goes over all the vortex filaments. When we apply the Schwarz model to simulate the motion of a vortex ring in quiescent He II, we set both un and us0 to zero and discretize the initial ring with a resolution Δξ = 0.005 mm. The time evolution of each vortex segment’s position can be obtained through a temporal integration of Eq. (1) using the fourth-order Runge-Kutta method53 with a time step Δt = 10−5 s.

2W model

In the 2W model, the normal-fluid velocity un is no longer prescribed. Instead, it is calculated by solving the classical Navier-Stokes equation with an added mutual friction term18:

$$\frac{\partial {{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}}{\partial t}+({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}\cdot {{{{{{{\boldsymbol{\nabla }}}}}}}}){{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}=-\frac{1}{{\rho }_{{{{{{{{\rm{He}}}}}}}}}}{{{{{{{\boldsymbol{\nabla }}}}}}}}P+{\nu }_{{{{{{{{\rm{n}}}}}}}}}{\nabla }^{2}{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}+\frac{{{{{{{{{\bf{F}}}}}}}}}_{{{{{{{{\rm{ns}}}}}}}}}}{{\rho }_{{{{{{{{\rm{n}}}}}}}}}}$$
(3)

where ρn and ρHe are, respectively, the normal-fluid density and the total density of He II, P is the pressure, νn is the He II kinematic viscosity, and Fns is the mutual friction per unit volume which can be calculated as:

$${{{{{{{{\bf{F}}}}}}}}}_{{{{{{{{\rm{ns}}}}}}}}}({{{{{{{\bf{r}}}}}}}})=\frac{1}{\Delta \Omega ({{{{{{{\bf{r}}}}}}}})}{\int}_{{{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{r}}}}}}}})}(-{{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{sn}}}}}}}}}/\Delta \xi )d\xi$$
(4)

where \({{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{r}}}}}}}})\) denotes that the integration is performed along all the vortex lines in the computational cell ΔΩ(r) = Δx × Δy × Δz located at r. When we simulate the vortex ring dynamics, Eqs. (1) and (3) are solved together to render the positions of the vortex-ring segments s(ξ, t) and the normal-fluid velocity un. The time integration of Eq. (3) is conducted using the second-order Adams-Bashforth method18 with the same time step Δt = 10−5 s, and the spatial differentiation is performed via the second-order finite difference with a spatial resolution Δx = Δy = Δz = 0.0083 mm. The computational domain consists of 1203 grids. We have confirmed that finer spatial and temporal resolutions do not change the simulation results.

S2W model

In the S2W model, the mutual friction force that acts on a vortex segment Δξ is given by17:

$${{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{sn}}}}}}}}}=[-D{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{L}}}}}}}}}))-{\rho }_{{{{{{{{\rm{n}}}}}}}}}\kappa {{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{L}}}}}}}}})]\Delta \xi,$$
(5)

where the only friction coefficient D can be calculated as:

$$D=-4\pi {\rho }_{{{{{{{{\rm{n}}}}}}}}}{\nu }_{{{{{{{{\rm{n}}}}}}}}}/[0.0772+\ln (|{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}^{\perp }-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{L}}}}}}}}}|{a}_{0}/4{\nu }_{{{{{{{{\rm{n}}}}}}}}})].$$
(6)

Here a0 ≃ 1 Å is the vortex-core radius36 and \({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}^{\perp }\) denotes the local normal-fluid velocity at the vortex-segment location that is projected in the plane perpendicular to the segment17. By balancing the Magnus force fM and the revised mutual friction force, the equation of motion for the vortex segment is now given by:

$$d{{{{{{{\bf{s}}}}}}}}/dt={{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{s}}}}}}}}}+\beta {{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{s}}}}}}}}})-{\beta }^{{\prime} }{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times [{{{{{{{{\bf{s}}}}}}}}}^{{\prime} }\times ({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{n}}}}}}}}}-{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{s}}}}}}}}})],$$
(7)

where the coefficients β and \({\beta }^{{\prime} }\) depends on D as derived by Galantucci et al.17. The evolution of the vortex position and un can be obtained by solving Eqs. (3) and (7) with D evaluated self-consistently via Eq. (6). In the S2W simulations, we adopted the same temporal and spatial resolutions as for the 2W simulations.

Note that the hydrodynamic description of the normal fluid is applicable only when the roton mean free path in He II is much smaller than the relevant length scales. This path is about 30 nm at 1 K and decreases with increasing the temperature54. Therefore, as we study a vortex ring of tens to hundreds of microns in diameter at temperatures above 1.6 K, the hydrodynamic description of the normal fluid should be reasonable.

For a quantized vortex ring with a radius R moving in quiescent He II, the self-induced superfluid velocity at the ring’s location is given by36 \({{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{s}}}}}}}}}=\frac{\kappa }{4\pi R}[\ln (8R/{a}_{0})-\frac{1}{2}]\hat{{{{{{{{\bf{n}}}}}}}}}\), which is the same in all three models. However, the local un is different, which leads to the different mutual friction dissipation rate. In the Schwarz model, un = 0 and therefore the highest mutual friction dissipation is expected. In both the 2W model and the S2W model, the back action of the mutual friction in the normal fluid generates two oppositely polarized normal-fluid vortex rings as shown in Fig. 1. In the 2W model, the two normal-fluid rings are concentrically located nearly in the same plane as the quantized vortex ring, whereas in the S2W model the two normal-fluid rings are slightly shifted to above and below the quantized-ring plane. This shift changes the direction of the local un. Nonetheless, the induced local un in both models has a component in the same direction as the local us, which effectively reduces the mutual friction dissipation as compared to that in the Schwarz model.

Effects of the trapped particles

When a vortex segment Δξ carries a trapped particle with a radius a, its equation of motion changes to39:

$$({m}_{{{{{{{{\rm{p}}}}}}}}}+{m}_{{{{{{{{\rm{f}}}}}}}}})\frac{d{{{{{{{{\bf{u}}}}}}}}}_{{{{{{{{\rm{L}}}}}}}}}}{dt}={{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{M}}}}}}}}}+{{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{sn}}}}}}}}}+{{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{D}}}}}}}}}+{{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{g}}}}}}}}}$$
(8)

where the term on the left-hand side represents the inertial effect caused by the trapped particle’s mass \({m}_{{{{{{{{\rm{p}}}}}}}}}={\rho }_{{{{{{{{\rm{p}}}}}}}}}\frac{4}{3}\pi {a}^{3}\) and the fluid’s added mass \({m}_{{{{{{{{\rm{f}}}}}}}}}=\frac{1}{2}{\rho }_{{{{{{{{\rm{He}}}}}}}}}\frac{4}{3}\pi {a}^{3}\). On the right-hand side, besides the Magnus force fM and the mutual friction force \({{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{sn}}}}}}}}}\), two additional forces are included, i.e., the Stokes drag exerted by the normal fluid on the particle fD = − 6πaμn(uL − un) and the gravitational force \({{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{g}}}}}}}}}=({\rho }_{{{{{{{{\rm{p}}}}}}}}}-{\rho }_{{{{{{{{\rm{He}}}}}}}}})\frac{4}{3}\pi {a}^{3}{{{{{{{\bf{g}}}}}}}}\). Other minor effects associated with the acceleration of the superfluid and the normal fluid around the trapped particle are negligible39. This model is accurate when a is much smaller than the separation distance between the particles trapped along the vortex ring, which holds true for the ring events that we selected to analyze.

To get a sense on how large the particle effects are, one may compare the total Stokes drag FD = ∣∑ifD,i∣ and the total gravitational force Fg = ∣∑ifg,i∣ with the total mutual friction force \({F}_{{{{{{{{\rm{sn}}}}}}}}}=|\oint ({{{{{{{{\bf{f}}}}}}}}}_{{{{{{{{\rm{sn}}}}}}}}}/\Delta \xi )d\xi|\), where ∑i means the summation over all the trapped particles and ∮ denotes the integration along the vortex ring. For the 9-particle vortex-ring event shown in Fig. 1, using the particle radiuses obtained through the size analysis (see later discussions in Methods), we estimate that \({F}_{{{{{{{{\rm{D}}}}}}}}}/{F}_{{{{{{{{\rm{sn}}}}}}}}}\) and \({F}_{{{{{{{{\rm{g}}}}}}}}}/{F}_{{{{{{{{\rm{sn}}}}}}}}}\) vary in the range of 10–18% and 4–4.7%, respectively, as R(t) decreases from 312 μm at t = 0 to 150 μm. On the other hand, for the 2-particle vortex-ring event shown in Fig. 2, despite the ring’s smaller initial radius and hence higher propagation speed, the ratios are \({F}_{{{{{{{{\rm{D}}}}}}}}}/{F}_{{{{{{{{\rm{sn}}}}}}}}}\simeq 4.8-10\%\) and \({F}_{{{{{{{{\rm{g}}}}}}}}}/{F}_{{{{{{{{\rm{sn}}}}}}}}}\simeq 0.8-0.83\%\), respectively, as R(t) shrinks from 140.8 μm at t = 0 to 50 μm.

Particle size distribution

We produce solidified D2 tracer particles in He II by slowly injecting a mixture of 5% D2 gas and 95% 4He gas directly into the plexiglass channel immersed in the He II bath. A computer-controlled solenoid valve is used to adjust the injection duration, and a needle valve is adopted to restrict the gas flow rate. Upon the injection, the D2 gas forms solidified ice particles. To evaluate the sizes of the resulted particles, we took images of the particles undergoing freely settling in quiescent He II in each experimental run. As a representative video, Supplementary Movie 7 displays the settling behavior of particles during the experiment run where we recorded the 9-particle ring event and the 2-particle ring event depicted in Fig. 1 and Fig. 2, respectively. By tracking the particles in such videos, we can generate a probability distribution of the particle settling velocity \({u}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{(s)}}}}}}}}}\). The result is shown in the Supplementary Fig. 1a. The \({u}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{(s)}}}}}}}}}\) data can be fitted nicely with a log-normal distribution, from which we can determine that the distribution is peaked at about 0.1 mm/s.

Note that the settling velocity is achieved when the Stokes drag exerting on a D2 particle is balanced by the gravitational force, i.e., \(6\pi a{\mu }_{{{{{{{{\rm{n}}}}}}}}}{u}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{(s)}}}}}}}}}=\frac{4\pi }{3}{a}^{3}({\rho }_{{{{{{{{\rm{p}}}}}}}}}-{\rho }_{{{{{{{{\rm{He}}}}}}}}})g\). This balance leads to \(a={[9{\mu }_{{{{{{{{\rm{n}}}}}}}}}{u}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{(s)}}}}}}}}}/2({\rho }_{{{{{{{{\rm{p}}}}}}}}}-{\rho }_{{{{{{{{\rm{He}}}}}}}}})g]}^{1/2}\). Therefore, knowing the distribution of \({u}_{{{{{{{{\rm{p}}}}}}}}}^{{{{{{{{\rm{(s)}}}}}}}}}\), we can then generate the radius distribution of the D2 particles. As shown in Supplementary Fig. 1b, this distribution is peaked at a ≃ 1.1 μm with a root mean variance Δa of about 0.2 μm.

Positions and radiuses of trapped particles

To evaluate the effects of the trapped D2 particles on the motion of a vortex ring, we need to know the radius and initial position of each individual trapped particle. Using the feature-point tracking routine37, we can determine the coordinates of every particles in the x-z image plane. For particles trapped on the vortex ring, their coordinates (xi, zi) should satisfy the following equation of an ellipse:

$$\frac{{[({x}_{i}-{x}_{0})\cos \phi+({z}_{i}-{z}_{0})\sin \phi ]}^{2}}{{R}_{1}^{2}}+\frac{{[({z}_{i}-{z}_{0})\cos \phi -({x}_{i}-{x}_{0})\sin \phi ]}^{2}}{{R}_{2}^{2}}=1,$$
(9)

where (x0, z0) are the coordinates of the ellipse center, R1 and R2 are, respectively, the semi-major and semi-minor axes of the ellipse, and ϕ is the angle between the ellipse major axis and the x-axis. These five parameters can be uniquely determined through a least squares fit to the positions of the trapped particles when there are at least five particles on the ring. Through this fit, we can determine the vortex ring radius R = R1 and the projection angle θ between the ring’s normal vector \(\hat{{{{{{{{\bf{n}}}}}}}}}\) and the x-z plane (i.e., \(\sin \theta={R}_{2}/{R}_{1}\)). If we set y0 = 0 for the ellipse center at t = 0, the initial yi of each trapped particle can be calculated as \({y}_{i}=[({x}_{i}-{x}_{0})\sin \phi -({z}_{i}-{z}_{0})\cos \phi ]/\tan \theta\). In the Supplementary Table 1, we list the 3D coordinates of all the nine trapped particles for the vortex ring presented in Fig. 1. The coordinate uncertainty comes from the feature-point tracking fit of the particle’s image profile. These coordinates are used in our model simulations.

Knowing the coordinates of these trapped particles, we can also monitor their angular position on the vortex ring over time. We set the ellipse major axis as our reference and denote α as the angle between the position vector OP of a trapped particle with this reference axis, where O is the center of the ring and P denotes the position of the particle. In Supplementary Fig. 2, we show the time variation of α for a few representative particles of the 9-particle ring event. It is clear that these particles do not exhibit significant displacements along the ring.

To evaluate the trapped particle’s radius a, we develop a correlation between a and the particle’s brightness I. For the particles that undergo freely settling in Supplementary Movie 7, we calculate the brightness I of each particle by summing up the counts in the image pixels associated with the particle. A distribution of the particle brightness P(I) can therefore be generated, which is shown in Supplementary Fig. 3. Since I depends on the particle’s surface area and hence a2, we can construct a simple correlation \(I=A{({a}^{2})}^{B}\), where A and B are tuning parameters. For a given pair A and B, we can scale the distribution of a shown in the Supplementary Fig. 1b to generate the distribution of the expected brightness \({I}^{{{{{{{{\rm{(ex)}}}}}}}}}=A{({a}^{2})}^{B}\). We then vary A and B to minimize the difference between the I(ex) distribution and the actual distribution P(I). At the optimal values A* = 1.20 and B* = 1.17, the generated I(ex) distribution agrees nicely with P(I), as shown in Supplementary Fig. 3.

Using the derived correlation \(I={A}^{*}{({a}^{2})}^{{B}^{*}}\), we can calculate the radius ai of a trapped particle i by measuring its brightness Ii. However, we must note that this correlation holds only in a statistical sense. When we apply it to analyze the radiuses of individual particles, there can be intrinsic uncertainties. To improve the reliability, in practice we collect the brightness data of the particle i over the time period that it is clearly observed and then use the time-averaged brightness \({\bar{I}}_{i}\) in the correlation to calculate ai. The uncertainty of the particle radius Δai can be evaluated as the variation range of ai when its brightness is varied from \({\bar{I}}_{i}-\Delta {I}_{i}\) to \({\bar{I}}_{i}+\Delta {I}_{i}\). Both the mean radiuses and the uncertainties are included in the Supplementary Table 1.

Velocity and projection angle for 2P ring

For the 2-particle vortex ring event presented in Fig. 2, the centroid coordinates of the two particles can be easily caculated \({{{{{{{{\bf{x}}}}}}}}}_{{{{{{{{\rm{p}}}}}}}}}=\frac{1}{2}({{{{{{{{\bf{x}}}}}}}}}_{1}+{{{{{{{{\bf{x}}}}}}}}}_{2})\) and \({{{{{{{{\bf{z}}}}}}}}}_{{{{{{{{\rm{p}}}}}}}}}=\frac{1}{2}({{{{{{{{\bf{z}}}}}}}}}_{1}+{{{{{{{{\bf{z}}}}}}}}}_{2})\). To determine the centroid velocity up at time t, we use a method that involves calculating the slopes of linear fits to the xp(t) and zp(t) data collected over a span of 5 consecutive image frames. Specifically, we determine the x-component of up from the linear fit of the xp(t) data, and the z-component from the linear fit of the zp(t) data. The time range for data collection is from t − 2Δt to t + 2Δt. This approach is consistently applied to all velocity calculations.

A constraint on the projection angle θ between the ring’s propagation direction and the x-z image plane can be placed based on the time-variation of the particle’s brightness I(t). This is because \(\tan \theta=\Delta y/\Delta S\), where ΔS = 1.12 mm is the distance traversed by the centroid of the two particles in the x-z plane over the observation time tf = 0.89 s, and Δy is the centroid displacement in the y direction perpendicular to the laser sheet, which can be estimated based on the variation of I(t).

To estimate Δy, we first show the measured brightness Im(t) of each particle in the Supplementary Fig. 4a. The variation of Im(t) is caused by the displacement of the particles in both the y direction and the z direction, since Im(t) is proportional to the laser intensity W which varies primarily in these two directions. To quantify the laser-intensity variations, we then place an optical power meter behind a mask with a narrow slit (20 μm in width) oriented either horizontally or vertically. By moving the horizontal slit in the z direction or by moving the vertical slit in the y direction, we can measure W as a function of y and z. The results are shown in the Supplementary Fig. 4b and c, respectively. The profile of W in each direction can be reasonably fit with a Gaussian function, which renders \(W(y,\, z)\propto {e}^{-2{(y-{y}_{{{{{{{{\rm{c}}}}}}}}})}^{2}/{\sigma }_{{{{{{{{\rm{y}}}}}}}}}^{2}}\cdot {e}^{-2{(z-{z}_{{{{{{{{\rm{c}}}}}}}}})}^{2}/{\sigma }_{{{{{{{{\rm{z}}}}}}}}}^{2}}\), where yc = 0 and zc = − 2.2 mm are the coordinates of the beam’s cross-sectional center, σy = 0.69 mm is the half-thickness of the laser sheet at 1/e2 intensity (i.e., which corresponds to a full thickness at half maximum intensity of 0.82 mm), and σz = 3.5 mm is the sheet’s half-height at 1/e2 intensity.

Finally, we can calculate the corrected brightness \(I(t)={I}_{{{{{{{{\rm{m}}}}}}}}}(t)/{e}^{-2{(z(t)-{z}_{{{{{{{{\rm{c}}}}}}}}})}^{2}/{\sigma }_{{{{{{{{\rm{z}}}}}}}}}^{2}}\). The results are shown in Fig. 2c. The variation of I(t) is entirely due to the particle displacement in the y direction. Since I(t)/I(0) for either particle decreases roughly monotonically by about 20% over the observation time, we can estimate the displacement Δy based on the Supplementary Fig. 4c. For a given initial particle coordinate y(0), we can determine Δy that gives 20% laser-intensity drop. By varying y(0), we find that Δy can reach up to about 0.2 mm. This sets an upper limit \(\tan \theta \le 0.2/1.12=0.18\). Since \(\cos \theta={u}_{{{{{{{{\rm{p}}}}}}}}}(0)/u(0)\), a constraint on u(0) and hence the initial ring radius R(0) can be placed. This constraint together with the other constraints discussed in the paper render the tuning range of the simulation curves shown in Fig. 2d.