Modelling of vorticity, sound and their interaction in two-dimensional superfluids

Superfluidity in two dimensions, first systematically investigated in the 70s in helium thin films [1–4], has sparked major research efforts in recent years, culminating in the 2016 Nobel Prize, awarded for understanding the nature of superfluidity in two dimensions [5–7]. The superfluid phase transition is native to a broad variety of physical systems, such as two-dimensional Bose–Einstein condensates [8, 9], exciton-polariton condensates [10, 11], and topological condensed matter systems [12]. Quantized vortices play a crucial role in these two-dimensional fluids, as their binding into pairs enables the emergence of long-range order and thereby transition to the superfluid phase. Understanding of vortex dynamics provides a pathway for controlling vortices with sound, imaging vortex distributions, understanding quantum turbulence, and engineering dynamical interactions between vortices and sound [13–16].

Vortex dynamics in strongly interacting superfluids is of significance to a range of research fields: in topological condensed matter physics, it is responsible for the superfluid phase transition and the onset of dissipation [5–7]; In astrophysics, the observed glitches in the rotation frequency of neutron stars are thought to result from vortex-unpinning events [17]; evidence for half quantum vortices (HQVs) has been found in superfluid ³He, where the HQVs in the A-phase of superfluid ³He are thought to host Majorana-fermions, bearing promise for fault tolerant topological quantum computing [18–20]; vortices in ³He are interesting as analogues of exotic topological defects [19, 21]—the broken-symmetry-core vortex in superfluid ³He-B corresponding to Witten-strings [22–24]; the HQV in the polar phase of superfluid ³He corresponding Alice-strings [25]; the spin-mass-vortex in ³He-B, which has been proposed as an analogue for composite defects appearing in some grand unified theories of particle physics and even the standard model [26, 27]. Therefore, the ability to determine the flow field induced by an arbitrary configuration of vortices, on an arbitrary and perhaps multiply-connected geometry is of broad importance, as is the ability to predict the strength of vortex-sound interaction.

In the case of Bose–Einstein condensates, pressure- and temperature-waves constitute the relevant sound eigenmodes (first- and second-sound, respectively), while for helium thin films, both these modes of oscillation are suppressed due to the incompressibility of the fluid and the clamping of the normal fluid component. Surface excitations, so-called third-sound waves, become the primary form of sound wave [28]. Recently, both observation of temperature-wave propagation in a two-dimensional Bose–Einstein condensate [29], and real-time measurement and control of third sound on a superfluid helium thin film [30] have been demonstrated.

Quantifying vortex flow fields and their interaction with sound waves has sparked substantial research efforts. Ellis et al [31–33] electrically excited third-sound modes to swirl up an ensemble of vortices in a helium thin film—however, despite elaborate mathematical analysis [34], their modelling of vortex-sound interaction was limited to simple, centred vortex ensembles on a circular resonator. In the case of BECs, a plethora of analyses for vortex-sound interaction has been performed [35–43] and the dispersive interaction has been quantified for centred vortices in simple trap geometries [36, 37].

In this work, we model vortex flow fields and the interaction of sound modes with vortices in a two-dimensional superfluid by exploiting analogies with other areas of physics. We map vortex dynamics onto electrostatics, and superfluid hydrodynamics onto fluid dynamics of an ideal gas. This allows us to draw on technically mature finite-element-modelling (FEM) tools available for these fields. We show how the interactions of sound and the flow field of arbitrary vortex distributions can be computed using these tools on any two-dimensional, not necessarily simply connected, domain. Thus, our work provides a theoretical framework for controlling and imaging vortices, and for engineering a dynamical interaction between sound and vortices.

As an example, we discuss the interaction of sound Bessel modes on a disk-shaped domain with quantized vortices, which is relevant for a number of experiments on superfluid thin films [30, 31, 44, 45] and two-dimensional Bose–Einstein-condensates [29, 46, 47]. The interaction induces splitting between otherwise degenerate sound modes [31, 36, 37]. We show how the vortex number can be extracted from experimental measurement of the splitting. Further, we present a perturbative analytic model, approximating the FEM-simulation for circular geometries, which offers some intuitive insight on the vortex-sound coupling mechanism.

Lastly, we focus our analysis to the prospect of detecting quantized circulation in helium thin films. While vortices in Bose–Einstein condensates can be visualized by optical snapshots [48], and pinned vortices in exciton-polariton condensates can be visualized by optical interferometry, no such direct observation technique exists for helium thin films: the vortex core is an Ångström-size perturbation on an ultra-thin film of transparent liquid, whose flow does not interact dissipatively with the environment. Measurements on helium are important because, unlike exciton-polariton condensates and most Bose–Einstein condensates, the atoms in superfluid helium are strongly interacting, introducing dynamics that can not be modelled through the Gross–Pitaevskii equation, and are not fully understood [49].

We finally propose an experiment where discrete steps due to an increase or decrease in the number of circulation quanta could be observed for the case of a superfluid helium film. We suggest a geometry where vortex pinning around an engineered topological defect leads to experimentally observable quantized steps. Drawing on the finite element model, we discuss how, in this geometry, the interaction with sound can be maximized, so that these steps could be clearly resolved. This would enable the first direct detection of quantized circulation in two-dimensional superfluid helium.

2.1. Sound in two-dimensional superfluids

Superfluid hydrodynamics is generally described by the continuity equation [31]:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\rho }{{\rm{d}}{t}}=-\vec{{\rm{\nabla }}}\cdot (\rho \vec{v}),\end{eqnarray} \tag{ 1 }$$

which derives from mass conservation, and the Euler equation:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\vec{v}}{{\rm{d}}{t}}+(\vec{v}\cdot \vec{{\rm{\nabla }}})\vec{v}=-g\vec{{\rm{\nabla }}}\rho +C,\end{eqnarray} \tag{ 2 }$$

which derives from Newton's second law (momentum conservation). In the case of superfluid helium thin films, $\vec{v}$ is the superfluid flow velocity, $\rho \to h$ is the film height, $g=\tfrac{3{\alpha }_{\mathrm{vdw}}}{{h}^{4}}$ is the linearized Van-der-Waals acceleration [28, 49], and C = 0. In the case of Bose–Einstein-condensate hydrodynamics in the Thomas-Fermi-limit at zero temperature [50, 51], $\vec{v}$ is again the flow velocity, ρ is the density, $g\to {g}_{\mathrm{BEC}}/{M}^{2}$ describes the coupling strength, where M is the mass of an individual atom contributing to the condensate, g_BEC the atom–atom coupling, and $C=-\vec{{\rm{\nabla }}}U/M$ describes the trapping, with U being the extended potential [50] (see appendix, table C2).

We assume small perturbations in the density (BEC) or film height (helium), η, from an equilibrium ρ₀, $\rho (\vec{r},t)={\rho }_{0}+\eta (\vec{r},t)$ with $\eta \ll {\rho }_{0}$. Equations (1) and (2) respectively become:

$$\begin{eqnarray}&&\dot{\eta }=-{\rho }_{0}\vec{{\rm{\nabla }}}\cdot \vec{v}-\vec{v}\,\vec{{\rm{\nabla }}}\eta ,\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&\dot{\vec{v}}+(\vec{v}\cdot \vec{{\rm{\nabla }}})\vec{v}=-g\vec{{\rm{\nabla }}}\eta +C,\end{eqnarray} \tag{ 4 }$$

where we used the chain rule on the rhs in equation (2). Substituting $g\to {c}^{2}=\gamma {RT}/{\rho }_{0}$, where R is the specific gas constant, T the gas temperature and c the speed of sound, we find that the above equations are the linearized Euler and continuity equations, describing small amplitude sound waves in an ideal gas in the isentropic limit [52] (see appendix B for derivation). This allows us to model them using the Aeroacoustics $\to$ Linearized Euler, Frequency Domain (lef) module in COMSOL, with appropriate boundary conditions (see appendix D.1). This provides the sound eigenmodes for an arbitrary bounded geometry. Examples of eigenmodes on a circular and an irregularly shaped two-dimensional resonator with free ('Neumann') boundaries are shown in figure 1(a) [45] (see appendix D.1). In this work, we analyse vortices and sound modes in two-dimensional domains, as this applies to experiments [31, 45, 53] and allows a qualitative discussion of the interaction. However, by choosing the appropriate dimension in the FEM-simulation, the model could be generalized to three dimensions. The finite-element-method allows us to add a background flow field, corresponding to the flow generated by quantized vortices, and find the new sound eigenmodes in the presence of that background flow.

Zoom In Zoom Out Reset image size

2.2. Quantized vortices

A single vortex is described by a quantized circulation around a loop encompassing the vortex core [54, 55] :

$$\begin{eqnarray}&&\oint {\vec{v}}_{v}\cdot {\rm{d}}\vec{l}=\kappa .\end{eqnarray} \tag{ 5 }$$

Here, κ = 2πℏ/M is the circulation quantum. This ensures that the phase acquired by the wave function upon propagation around any loop encompassing the vortex core equals 2π. ${\vec{v}}_{v}$ denotes the vortex-induced velocity field. For the simple case of a point vortex on a plane, the solution is:

$$\begin{eqnarray}&&{\vec{v}}_{v}(r)=\displaystyle \frac{\kappa }{2\pi r}{\hat{e}}_{\theta },\end{eqnarray} \tag{ 6 }$$

where ${\hat{e}}_{\theta }$ is the unit vector in the tangential direction and r the distance from the vortex core. In this work, we describe the quasi-static regime where the motion of vortices during a sound oscillation period is negligible. This is valid in the limit of pinned vortices [31, 33, 56] or low vortex densities where the velocity of the flow field due to neighbouring and image vortices is significantly less than the speed of sound. For Bose–Einstein condensates at zero temperature, the sound velocity c = ℏ/Mξ, where ξ is the healing length, equals the Landau critical velocity [51]. Therefore the quasi-static approximation in Bose–Einstein condensates is valid if the separation between neighbouring vortex cores is significantly larger than their core diameter. This condition is typically fulfilled, and the sound velocity, with typical values of a few mm s⁻¹ [57, 58] is significantly higher than the background flow velocity caused by a typical ensemble of vortices [13, 53]³ .

In comparison to equation (5), Gauss's law of electrostatics in two dimensions reads:

$$\begin{eqnarray}&&\oint \vec{D}\cdot {\rm{d}}\vec{n}=Q,\end{eqnarray} \tag{ 7 }$$

where $\vec{D}$ is the electric displacement field and Q the line charge. By rotating the electric displacement field and replacing it with the flow field ${\vec{v}}_{v}$

$$\begin{eqnarray}&&\left(\begin{array}{c}{D}_{x}\\ {D}_{y}\end{array}\right)\to \left(\begin{array}{c}{v}_{{vy}}\\ -{v}_{{vx}}\end{array}\right),\end{eqnarray} \tag{ 8 }$$

and substituting $Q\to \kappa$, we retrieve the quantized circulation of the vortex defined in equation (5), where ${\vec{v}}_{v}={v}_{v,x}\,{\hat{e}}_{x}+{v}_{v,y}\,{\hat{e}}_{y}$. This provides the analogy between electrostatics and vortex flow [54]. A point charge is a source of divergence (source/sink) for the electric displacement field $\vec{D}$. As is known from potential flow theory, upon the permutation shown in equation (8), a point charge becomes a source of quantized circulation. We model these equations using the Electrostatics(es) module of COMSOL, which allows us to determine the vortex flow field on any two-dimensional geometry. Examples for vortex flow fields on a circular and on an irregular geometry are shown in figure 1(b). Depending on the number of vortices, their positions, and the resonator geometry, the sound eigenmode shape may be significantly altered due to the presence of the vortices. Such an example is shown in figure 1(c). Similarly, quantized circulation around a macroscopic topological defect in a multiply-connected domain can be modelled as shown in section 4.

2.3. Sound-vortex interaction

We can understand the interaction of sound waves and vortices through the change in the kinetic energy of a sound wave caused by addition/subtraction of a vortex. The sound modes are orthogonal to vortex flow fields, which are fully defined by the rotation around the vortex core (see figure 2 (b)). Therefore, the overlap of vortex- and sound velocity fields is zero:

$$\begin{eqnarray}&&{\int }_{A}{\vec{v}}_{v}\cdot {\vec{v}}_{s}\,{\rm{d}}{A}=0,\end{eqnarray} \tag{ 9 }$$

where A is the area of the domain and ${\vec{v}}_{s}$ is the two-dimensional sound velocity distribution. This appears to suggest that there is no coupling between the two flow fields. However, interaction arises due to the change in film height (helium) or density (BEC) associated with the sound wave. The interaction manifests in a splitting of sound resonances due to the presence of vortices, and thus constitutes a dispersive (frequency-shifting) interaction [59]. An example is shown in figure 2(c), where, due to increased density/height on one side of the vortex, the increased kinetic energy due to velocity addition on the right side of the vortex is not fully compensated by the reduction in kinetic energy due to velocity subtraction on the left, resulting in a net increase in energy due to the interaction.

Zoom In Zoom Out Reset image size

We use equations (3) and (4) to calculate the new eigenmodes of sound in the presence of the time-independent background flow ${\vec{v}}_{v}$, assuming a stationary configuration of vortices, with flow field ${\vec{v}}_{v}(\vec{r})$. In this case, the total flow velocity is $\vec{v}(\vec{r},t)={\vec{v}}_{v}(\vec{r})+{\vec{v}}_{s}(\vec{r},t)$, where ${\vec{v}}_{s}$ is the flow field associated with the sound eigenmode. The interaction of persistent currents with sound modes has been quantified for simple, centred vortex distributions on a disk-shaped resonator in superfluid helium thin films [31, 34] and for centred vortices in different trap shapes in Bose–Einstein-condensates [36, 37], but until now there has not been a consistent approach for modelling of a non-trivial vortex distribution in a non-trivial resonator shape, or for multiply connected domains. We confirm the accuracy of this finite element based approach through the comparison with an analytical expression derived for third-sound-vortex interactions on a circular resonator (see appendix A).

To give an experimentally relevant example, in the following we study the interaction of a quantized vortex with sound modes in a disk-shaped resonator with a free ('Neumann') boundary condition (see appendix D.1). This analysis is applicable to geometries used in superfluid helium experiments [30, 44, 45] with a microtoroidal optomechanical resonator of R ∼ 30 μm radius (see figure 3(a)), those of [31, 56, 60, 61], and also experiments with two-dimensional Bose–Einstein-condensates, which are confined by a hard-walled trap [53].

Zoom In Zoom Out Reset image size

Regarding the experimental readout of sound modes, in experiments with helium thin films, the Brownian motion of third sound waves, even at millikelvin temperatures, is high enough to be resolved experimentally in real time [30]. Alternatively, the amplitude of third sound can be tuned by laser heating or cooling [30], amplified by optical absorption heating [44], or electrically excited [31]. For Bose–Einstein condensates, collective modes can be excited by perturbing the condensate, locally exceeding the critical velocity [62].

Solutions of the wave equation on a circular resonator are Bessel modes of the first kind. They are fully quantified by their azimuthal ($m\geqslant 0$) and radial ($n\geqslant 1$) node counts. Modes with $m\ne 0$ can be decomposed into opposite direction travelling waves, that, in the absence of circulation, are degenerate. The presence of a vortex lifts the degeneracy, shifts the mode co-rotating with the vortex flow to a higher frequency, and the counter-rotating mode to a lower frequency (see figure 3(b)). This frequency splitting Δf is experimentally resolved if it is larger than the decay rate Γ of the Bessel mode. As the radial distance of the vortex from the centre of the disk, r_v, increases, the splitting reduces. Firstly, this results from the lower total energy of the vortex flow field as the vortex core approaches the boundary. Secondly, from symmetry arguments the overlap between sound and vortex flow fields is maximized for a centred vortex. The radial dependence of the mode shifts is plotted for the (m,n) = (1, 8) free-boundary-condition Bessel mode in figure 3(c). Each point corresponds to one result of the FEM-simulation, as the vortex is stepped outwards from the centre. As the vortex reaches the outer boundary, the splitting vanishes.

We then compare the results from the FEM-model to the perturbative analytical approach derived in appendix A. We find that the two approaches agree reasonably well, with the discrepancy always less that 10% of the maximal splitting at r_v = 0 (see appendix A). We ascribe the difference to a vortex-induced change in the mode shape due to the nonlinear $(\vec{v}\cdot \vec{{\rm{\nabla }}})\vec{v}$ term in equation (4):

$$\begin{eqnarray}&&(\vec{v}\cdot \vec{{\rm{\nabla }}})\vec{v}=({\vec{v}}_{v}\cdot \vec{{\rm{\nabla }}}){\vec{v}}_{v}+({\vec{v}}_{s}\cdot \vec{{\rm{\nabla }}}){\vec{v}}_{s}+({\vec{v}}_{v}\cdot \vec{{\rm{\nabla }}}){\vec{v}}_{s}+({\vec{v}}_{s}\cdot \vec{{\rm{\nabla }}}){\vec{v}}_{v}.\end{eqnarray} \tag{ 10 }$$

The analytic perturbation theory neglects all of these terms. The FEM model neglects only $({\vec{v}}_{s}\cdot \vec{{\rm{\nabla }}}){\vec{v}}_{s}$, which is a requirement to obtain eigenmodes for the sound waves. This approximation is justified in the small sound amplitude limit $\eta \ll {\rho }_{0}$, where the sound-induced superfluid flow velocity v_s (see equation (A.4)) is small compared to the vortex flow velocity v_v.

Figure 3(d) shows the splitting per vortex as a function of vortex radial position for four different Bessel modes. Critically, the presence of a vortex affects each Bessel mode in a unique fashion. Leveraging this unique fingerprint, the work presented here enabled both the number and the spatial distribution of vortices in a cluster to be extracted independently, by tracking several sound modes simultaneously [45].

One conceptual result of the perturbative analysis is an expression for the frequency splitting that depends only on the profile of the Bessel mode and the radial position of the vortex, and is independent of the details of the vortex flow field. In the case of a single vortex the expression is:

$$\begin{eqnarray}&&{\rm{\Delta }}f=\displaystyle \frac{\kappa \,m}{2{\pi }^{2}}\,\displaystyle \frac{{\int }_{{r}_{v}}^{R}\tfrac{{\rm{d}}r}{r}\,{\eta }^{2}\left(r\right)}{{\int }_{0}^{R}{\rm{d}}r\,r\,{\eta }^{2}\left(r\right)},\end{eqnarray} \tag{ 11 }$$

where m is the Bessel mode azimuthal number, r_v is the radial position of the vortex, and $\eta \left(r\right)={J}_{m}\left({\zeta }_{m,n}\,\tfrac{r}{R}\right)$ is the radial displacement profile of the Bessel mode. J_m is the Bessel function of the first kind of order m, and m and n are respectively the azimuthal and radial mode orders, ζ_m,n is a frequency parameter depending on the mode order and the boundary conditions [63]. As vortex flow fields are subject to linear superposition, we obtain the total splitting simply by adding contributions from single vortices: ${\rm{\Delta }}{f}_{\mathrm{total}}={\sum }_{i}{\rm{\Delta }}{f}_{i}$. An interpretation of this result is that the splitting introduced by a vortex at position r_v is equal to the interaction energy between a centred vortex (r_v = 0) and a sound wave in the region of the disk with radius greater that r_v. So, in some sense, only the fraction of the sound wave at radius larger than the radial position of the vortex core contributes to the splitting. This explains why the splitting per vortex drops rapidly with r_v for low m, high n order sound modes whose kinetic energy is located close to the centre of the disk, while the splitting is sustained at higher radii for higher m, lower n order modes which are more radially extended.

In this section, we investigate the feasibility of observing the quantization of circulation in two-dimensional superfluid helium due to the shift in sound frequencies induced by the addition/subtraction of a single vortex. Remarkably, while quantized vortices are central to the behaviour of two-dimensional superfluids, they have yet to be directly observed in two-dimensional helium. The experimental challenge is significant: the normal-fluid core of a vortex in superfluid helium-4 is roughly one Ångström in diameter [49], the thickness of a superfluid helium film is typically less than 20 nm, and the refractive index of liquid helium is close to that of vacuum (n_He ≈ 1.029). Combined, these characteristics prevent direct optical imaging, as can be performed in Bose–Einstein condensates [47, 64]. In bulk helium, many imaging techniques have relied on the use of some kind of tracer particle [65–67], such as, for instance, micrometer-sized frozen hydrogen crystals. These scatter light and are pulled in to the vortex core, enabling, for example, the recent observation of Kelvin waves [67] in bulk. Naturally, such an approach is significantly more difficult in two-dimensional films due to their few-nanometre thickness.

In order for the vortex-induced frequency splitting Δf (equation (11)) experienced by a third-sound wave [31, 60] to reveal quantized steps, several challenges must be addressed.

First, in order to be resolvable, the splitting should be larger than the linewidth of the third-sound resonances, Δf ≳ Γ, as shown in figure 3(a). Second, any motion of vortices on the resonator surface, as we experimentally observe elsewhere [45], will lead to a continuous evolution of the splitting due to the continuous nature of the splitting function Δf(r_v), see figure 3(c), which may mask the quantized nature of the circulation.

The first challenge can be met by engineering devices that are sufficiently small, which maximizes Δf by increasing the vortex-sound coupling, and by controlling dissipation in these devices in order to reduce Γ. This can be achieved, for instance, by engineering a smooth resonator that is decoupled from its environment by a small connection point [30], or through careful choice of the resonator substrate material [68]. Indeed, Δf ∼ Γ has been recently reported experimentally using microscale on chip optical cavities [45]. A solution to the second challenge is to constrain the position of the circulation around a macroscopic topological defect engineered on the surface of the resonator. For instance, if we replace the topological defect naturally formed by the normal fluid core (of radius a₀) of a superfluid vortex by a microfabricated hole of radius $R\gg {a}_{0}$, the maximal velocity due to the quantized circulation becomes $\tfrac{\kappa }{2\pi R}\ll \tfrac{\kappa }{2\pi a0}$ (see equation (5)). This effectively clips the high velocity region of the flow and is thus energetically favourable. The circulation will then preferentially accumulate around this manufactured defect, up to large values of κ = h/m_He, as has been observed in the spinning up of bulk helium in an annular container [69, 70]. The quantization of the circulation then manifests as quantized values of the splitting experienced by third-sound modes confined to the surface of the resonator.

This approach is in essence a two-dimensional analogue of Vinen's experimental technique for the first observation of circulation quanta in bulk helium [71, 72], where circulation trapped around a vibrating wire lifted the degeneracy between the wire's normal modes of vibration. In figure 4, we propose a practical realization of such a device based on a circular whispering-gallery-mode geometry, as used in our previous work [30, 44]. We utilise the FEM-simulation to design a domain that maximises the splitting Δf, consisting of a single-spoked annular geometry [73] (see figure 4(c)).

Zoom In Zoom Out Reset image size

Next we calculate the superfluid flow field resulting from quantized circulation about the central topological defect, as shown in figure 4(c). For clarity, we assume that this pinned circulation is the only source of circulation on the structure, i.e. there are no vortex cores on the domain. This is calculated through FEM simulation using the 'floating potential' boundary condition for the inner boundary, which enforces both the prescribed circulation strength and parallelism of the flow to the boundary (see appendix D.1 for more details). The results are shown in figure 5(a). The flow is essentially confined to the outer annulus, with negligible flow up and down the spoke and in the central disk. This can be understood by considering the closed contours 1 and 2 which both enclose the central hole. The circulation around both contours must therefore be equal (see equation (5)), implying negligible additional circulation along the extra path contained in contour 1.

Zoom In Zoom Out Reset image size

Figure 5(b) shows an example of a third-sound mode of this spoked resonator (which becomes the (1, 2) eigenmode of a circular resonator as the central hole gets vanishingly small). The presence of the spoke, which connects the annular outer ring to the device pedestal, lifts the degeneracy between the two normal modes, even in the absence of circulation. The mode that has a stronger interaction with the spoke (bottom) experiences an effectively larger resonator and therefore has a lower resonance frequency.

In figure 5(c), we show how this native geometric splitting [60] affects the third-sound mode splitting as a function of the circulation around the central hole. Each black dot represents a finite element simulation of the splitting between the high and low frequency eigenmodes shown in figure 5(b), as a function of the number of circulation quanta around the central defect. This total splitting can be well reproduced through an analytical expression of the form [33]:

$$\begin{eqnarray}&&{s}_{\mathrm{total}}=\sqrt{{s}_{\mathrm{circ}}^{2}+{s}_{\mathrm{geo}}^{2}},\end{eqnarray} \tag{ 12 }$$

where s_geo = 700 Hz is the native geometric splitting for this device (dashed orange line), ${s}_{\mathrm{circ}}=N\times {s}_{\mathrm{circ},0}$ the total circulation-induced splitting with N being the number of circulation quanta and ${s}_{\mathrm{circ},0}$ is the splitting per quantum. The solid red line represents s_total, as given by equation (12). In figure 5(d), we plot the experimentally relevant parameter, which is the splitting increment due to each additional circulation quantum, as a function of the number of circulation quanta already present around the defect. This shows that the geometric splitting due to the spoke (or any unwanted deviation from circularity) will mask the influence of the circulation-induced splitting for small values of the circulation quanta, and reduce the visibility of the steps. For larger values of the circulation, the size of the steps will asymptote towards the value ${s}_{\mathrm{circ},0}=54\,\mathrm{Hz}$, as the appropriate normal mode basis gradually shifts from orthogonal standing waves to counter-propagating waves. Such a large quantized circulation can be experimentally achieved by creating a strong superfluid flow, locally exceeding the critical velocity—for instance by local evaporation of superfluid or by strong driving of third-sound modes [30, 44, 45]. Alternatively, a sub-critical flow can be used to create a high persistent current by reorganizing pre-existing vortex pairs through the Magnus force [31]. In Bose–Einstein condensates, vortices can be created in high numbers by laser stirring [74].

We show that the proposed device would yield quantized steps in the third-sound mode splitting on the order of 50 Hz, a value within reach of current experimental resolution [30, 45]. Note also that since the strength of the vortex-phonon interaction scales inversely with resonator area (see appendix A.3) the magnitude of the splitting can be greatly enhanced by going towards miniature third-sound resonators. This is illustrated in figure 7, which shows the splitting per centred vortex on the (1, 2) Bessel mode as a function of resonator radius. While only of the order of milli-Hertz for early cm-scale capacitively detected third-sound resonators [60, 61] (red and orange dots), it reaches tens to hundreds of Hz with microtoroidal resonators [30, 45] (blue dot), and would attain tens of kHz with micron-radius resonators [75–77] (black dot).

We have developed finite-element modelling tools to compute the interaction between any vortex flow and any sound wave, in arbitrary and potentially multiply-connected two-dimensional domains. This capability offers great versatility, applicable to both BEC superfluids and to thin-film superfluid helium. There is a need for numerical techniques to determine vortex and sound velocity fields and their interactions. In both cases, analytical solutions for the vortex flow field only exist if the domain exhibits a high degree of symmetry. Even if such solutions exist, when departing from simple geometries like a disk, the implementation of the method of images in order to cancel the normal component of the vortex flow on the resonator boundary [78] becomes challenging, and one needs to rely on conformal mapping techniques [53, 79]. For multiply-connected domains, solutions often require an infinite series of images as the domain possesses two or more boundaries, and analytical solutions are only available for simple limit-cases such as a centred annular domain [70].

We verify the validity of our approach by comparing its results to a perturbation theory analysis which we derive in the analytically tractable case of a circular resonator geometry. We derive in this case a useful simple analytical formula, which can be used to compute the vortex-sound coupling for arbitrary configurations of vortices on a disk, without requiring the vortex flow field. Understanding precisely how superfluid vortices and persistent currents couple to sound waves—at the level of a single vortex or circulation quantum—is a crucial capability to shed light on the physics of strongly interacting superfluids, and perform continuous non-destructive measurements of vortex dynamics in these systems [45].

The modelling techniques presented here may help shed light on the validity of phenomenological models such as the point-vortex model [80, 81] in superfluids, as well as further our understanding of quantum turbulence [13–16] and energy dissipation in superfluids [56, 82, 83].

Authors would like to thank Nick Wyatt for his help in setting up the COMSOL simulation in its early stages, Prahlad Warszawski, Andrew Doherty, Matt Davis, Rachpon Kalra, Yasmine Sfendla and Andreas Sawadsky for discussions on vortex dynamics, and Gian-Marco Schnüriger for discussions on wave dynamics. The authors in particular would like to thank Matt Reeves and David Colas for helpful discussions on vortex dynamics and sound in Bose–Einstein condensates. This work was funded by the U.S. Army Research Office through grant number W911NF-17-1-0310. It was also supported by the Australian Research Council through the Centre of Excellence for Engineered Quantum Systems (EQuS, CE110001013); WPB acknowledges the Australian Research Council Future Fellowship FT140100650. CGB acknowledges a Fellowship from the University of Queensland (UQFEL1833877).

A.1. Analytical description of the point vortex flow field ${\vec{v}}_{v}$

The streamfunction Ψ of a point vortex of strength κ in a 2D plane is ${\rm{\Psi }}=-\tfrac{\kappa }{2\pi }\mathrm{ln}\left(r\right)$. The streamfunction Ψ for the the well-known problem of a point vortex inside a circular domain [78], is given in cartesian coordinates by:

$$\begin{eqnarray}&&{\rm{\Psi }}=-\displaystyle \frac{\kappa }{2\pi }\left(\mathrm{ln}\left(\sqrt{{\left(x-{X}_{1}\right)}^{2}+{y}^{2}}\right)-\mathrm{ln}\left(\sqrt{{\left(x-{X}_{2}\right)}^{2}+{y}^{2}}\right)\right).\end{eqnarray} \tag{ A.1 }$$

Here X₁ is the radial offset of the vortex (along the x axis), and ${X}_{2}=\tfrac{{R}^{2}}{{X}_{1}}$ is the radial coordinate of the opposite circulation image-vortex required to enforce no flow accross the resonator boundary [78]. From the streamfunction Ψ, the vortex velocity components are given by:

$$\begin{eqnarray}&&{v}_{{vx}}=\displaystyle \frac{\partial {\rm{\Psi }}}{\partial y};\qquad \mathrm{and}\qquad {v}_{{vy}}=-\displaystyle \frac{\partial {\rm{\Psi }}}{\partial x}.\end{eqnarray} \tag{ A.2 }$$

Using equation (A.2), we plot in figure 6(a) the flow streamlines of a vortex offset by R/2 inside a circular resonator of radius R.

Zoom In Zoom Out Reset image size

A.2. Analytical description of the sound flow field ${\vec{v}}_{3}$

In the following, for simplicity we derive the analytic expressions for helium thin films (i.e. third-sound waves). However, the analysis can be applied to BECs with a straightforward replacement of variables (see appendix C).

The complex surface displacement amplitude η of a travelling superfluid third-sound wave (or alternatively the sound-induced density fluctuations for a BEC) in a circular domain are given by [63]:

$$\begin{eqnarray}&&\eta \left(r,\theta ,t\right)={\eta }_{0}\,{J}_{m}\left({\zeta }_{m,n}\displaystyle \frac{r}{R}\right){{\rm{e}}}^{{\rm{i}}\left(m\theta \pm {\rm{\Omega }}t\right)}.\end{eqnarray} \tag{ A.3 }$$

The associated superfluid flow speed ${\vec{v}}_{3}$ is given by:

$$\begin{eqnarray}&&{\vec{v}}_{3}=\pm \displaystyle \frac{i\,{c}_{3}^{2}}{{\rm{\Omega }}\,{h}_{0}}\vec{{\rm{\nabla }}}\eta \end{eqnarray} \tag{ A.4 }$$

with '+' and '−' signs respectively corresponding to the CW and CCW travelling cases, as in equation (A.3). Note that while the motion of a solid circular membrane would also be given by equation (A.3), its velocity would be different to equation (A.4), leading to dramatically different effective mass scalings [63]. The surface displacement profile ${\mathfrak{R}}\left(\eta \right)$ and instantaneous velocity field ${\mathfrak{R}}\left({\vec{v}}_{3}\right)$ are plotted in figure 6(b), for a CW (m = 1, n = 2) Bessel mode with free boundary conditions. While such a third-sound mode flow is irrotational and therefore not associated with any circulation ($\oint \vec{v}\cdot \vec{{\rm{d}}l}=0$ for any closed loop inside the superfluid), it is associated with a net mass flow (in the CW case, there is more fluid moving clockwise under the wave peak (red) than counter-clockwise under the trough (blue), and similarly there is net CCW fluid motion for the CCW mode). It is this net mass flow which couples to the vortex field, and results in a higher kinetic energy for the sound mode travelling in the same direction as the vortex flow. This argument is developed in the analytical splitting calculation detailed below.

A.3. Analytical description of vortex-sound coupling

Here we derive an analytical expression for the frequency splitting experienced by a third-sound mode due to a vortex inside a circular resonator, and show good agreement with the results of the FEM simulations shown in figure 3. This is a valid approximation to the FEM-model if the change in mode shape due to vortices is small.

The kinetic energy difference ${\rm{\Delta }}E\left(t\right)$ between a sound wave moving with- or against the flow of a quantized vortex is given by:

$$\begin{eqnarray}&&{\rm{\Delta }}E\left(t\right)=\displaystyle \frac{1}{2}\rho {\int }_{\theta =0}^{2\pi }{\int }_{r=0}^{R}{\int }_{z=0}^{{h}_{0}+\eta \left(r,\theta ,t\right)}\left({\left|\left|{\vec{v}}_{3}\left(\vec{r},t\right)+{\vec{v}}_{v}\left(\vec{r}\right)\right|\right|}^{2}-{\left|\left|{\vec{v}}_{3}\left(\vec{r},t\right)-{\vec{v}}_{v}\left(\vec{r}\right)\right|\right|}^{2}\right)r{\rm{d}}r\,{\rm{d}}\theta {\rm{d}}z.\end{eqnarray} \tag{ A.5 }$$

This general expression works for any sound mode and any vortex position. Making the reasonable assumption that ${\vec{v}}_{3}$ and ${\vec{v}}_{v}$ are independent of z, as the inviscid nature of the superfluid precludes any in-plane vorticity and does not require cancellation of the horizontal velocity at z = 0 (no-slip boundary), equation (A.5) becomes:

$$\begin{eqnarray}&&{\rm{\Delta }}E\left(t\right)=2\rho {\int }_{\theta =0}^{2\pi }{\int }_{r=0}^{R}{\vec{v}}_{3}\left(r,\theta ,t\right)\cdot {\vec{v}}_{v}\left(r,\theta \right)\left({h}_{0}+\eta \left(r,\theta ,t\right)\right)r{\rm{d}}r{\rm{d}}\theta .\end{eqnarray} \tag{ A.6 }$$

Since both v_3x and v_vx as well v_3y and v_vy are functions of θ of different parity (see figure 6), $\iint {\vec{v}}_{3}\cdot {\vec{v}}_{v}=0$, and equation (A.6) becomes:

$$\begin{eqnarray}&&{\rm{\Delta }}E\left(t\right)=2\rho {\int }_{\theta =0}^{2\pi }{\int }_{r=0}^{R}{\vec{v}}_{3}\left(r,\theta ,t\right)\cdot {\vec{v}}_{v}\left(r,\theta \right)\eta \left(r,\theta ,t\right)r{\rm{d}}r{\rm{d}}\theta .\end{eqnarray} \tag{ A.7 }$$

This is essentially a form of surface-averaged Doppler shift, weighted by the displacement amplitude η of the mode. Next, we consider the time-averaged energy difference $\left\langle {\rm{\Delta }}E\right\rangle$, averaged over a sound oscillation period T:

$$\begin{eqnarray}&&\left\langle {\rm{\Delta }}E\right\rangle =\displaystyle \frac{1}{T}{\int }_{0}^{T}{\rm{\Delta }}E\left(t\right){\rm{d}}t=2\,\rho {\int }_{r}{\int }_{\theta }r{\rm{d}}r{\rm{d}}\theta \left({v}_{{vr}}\displaystyle \frac{1}{T}{\int }_{0}^{T}{v}_{3r}\eta {\rm{d}}t\right)+\left({v}_{v\theta }\displaystyle \frac{1}{T}{\int }_{0}^{T}{v}_{3\theta }\eta {\rm{d}}t\right),\end{eqnarray} \tag{ A.8 }$$

where we have broken down ${\vec{v}}_{3}$ and ${\vec{v}}_{v}$ into their radial and angular components, respectively v_3r and v_3θ, and v_vr and v_vθ. From equations (A.3) and (A.4), we note that v_3r and η are out-of-phase, while v_3θ and η are in phase. The first integral over time in equation (A.8) reduces therefore to zero, while the second integrates to $\tfrac{1}{2}\left|{v}_{3\theta }\right|\left|\eta \right|$. We therefore get from equations (A.3) and (A.4):

$$\begin{eqnarray}&&\left\langle {\rm{\Delta }}E\right\rangle =\displaystyle \frac{\rho \,m\,{c}_{3}^{2}}{{\rm{\Omega }}\,{h}_{0}}{\int }_{r=0}^{R}r{\rm{d}}r\,{\eta }_{0}^{2}\,\displaystyle \frac{{J}_{m}^{2}\left({\zeta }_{m,n}\tfrac{r}{R}\right)}{r}{\int }_{\theta =0}^{2\pi }\,{v}_{v\theta }{\rm{d}}\theta \end{eqnarray} \tag{ A.9 }$$

which we rewrite, with $\eta \left(r\right)={\eta }_{0}\,{J}_{m}\left({\zeta }_{{\rm{m}},{\rm{n}}}\,\tfrac{r}{R}\right)$, as:

$$\begin{eqnarray}&&\left\langle {\rm{\Delta }}E\right\rangle =\displaystyle \frac{\rho \,m\,{c}_{3}^{2}}{{\rm{\Omega }}\,{h}_{0}}{\int }_{r=0}^{R}\displaystyle \frac{{\rm{d}}r}{r}\,{\eta }^{2}\left(r\right){\int }_{\theta =0}^{2\pi }\,{v}_{v\theta }r\,{\rm{d}}\theta .\end{eqnarray} \tag{ A.10 }$$

We notice here that the integral over θ corresponds to a closed contour integral $\oint \,{\vec{v}}_{v}\cdot {\rm{d}}\vec{l}$, where the contour is a circle of radius r centred at the origin. From equation (5), we know that the value of this contour integral is zero if it does not enclose the vortex core, and κ if it does. The transition occurs for r = offset, the radial offset of the point vortex. We can therefore rewrite equation (A.10) with a modified radial integration lower bound:

$$\begin{eqnarray}&&\left\langle {\rm{\Delta }}E\right\rangle =\displaystyle \frac{\rho \,m\,{c}_{3}^{2}\,\kappa }{{\rm{\Omega }}\,{h}_{0}}{\int }_{r=\mathrm{offset}}^{R}\displaystyle \frac{{\rm{d}}r}{r}\,{\eta }^{2}\left(r\right).\end{eqnarray} \tag{ A.11 }$$

Since for a harmonic oscillator E is proportional to Ω², $\tfrac{{\rm{\Delta }}E}{E}=2\tfrac{{\rm{\Delta }}{\rm{\Omega }}}{{\rm{\Omega }}}$ and the splitting Δf (in Hz) equals:

$$\begin{eqnarray}&&{\rm{\Delta }}f=\displaystyle \frac{{\rm{\Omega }}}{4\pi }\displaystyle \frac{{\rm{\Delta }}E}{E},\end{eqnarray} \tag{ A.12 }$$

with the kinetic energy E of the third-sound mode, for m > 0, given by [63]:

$$\begin{eqnarray}&&E=\displaystyle \frac{1}{2}\int \rho \,{v}^{2}\left(\vec{r}\right){{\rm{d}}}^{3}\left(\vec{r}\right)=\displaystyle \frac{\pi \,\rho \,{c}_{3}^{2}}{2\,{h}_{0}}{\int }_{0}^{R}{\eta }^{2}\left(r\right)r{\rm{d}}r.\end{eqnarray} \tag{ A.13 }$$

Combining equations (A.10) and (A.13), we recover the result shown in equation (11) of the main text:

$$\begin{eqnarray}&&{\rm{\Delta }}f=\displaystyle \frac{\kappa \,m}{2\,{\pi }^{2}}\,\displaystyle \frac{{\int }_{\mathrm{offset}}^{R}\tfrac{{\rm{d}}r}{r}\,{\eta }^{2}\left(r\right)}{{\int }_{0}^{R}{\rm{d}}r\,r\,{\eta }^{2}\left(r\right)}.\end{eqnarray} \tag{ A.14 }$$

We note that, as expected, the splitting does not depend on the superfluid parameters (film thickness, density), and that it is linear in vortex flow field (see equation (A.6)), such that the splitting obeys the superposition principle, whereby the splitting due to an ensemble of vortices is equal to the sum of the splittings per vortex calculated individually. The result equation (A.14) holds for both superfluid helium thin films and Bose–Einstein condensates, with η being the film thickness perturbation/density perturbation, respectively. We numerically verify this result in the FEM simulations, where linearity is generally maintained up to large vortex charges on the order of ∼10² κ, as shown in figures 7(a) and (b). Note that equation (A.14) does not diverge as the vortex offset tends to 0, as $\eta \left(0\right)=0$ for all m > 0 Bessel modes. Interestingly, due to the contour-integral identity used in equation (A.10), the final result does not require any knowledge of the vortex flow field ${\vec{v}}_{v}$. Figure 8 shows a comparison of the analytical and the FEM approaches, in the case of a circular resonator geometry.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the ideal gas, mass conservation and momentum conservation read, respectively [52]:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\rho }{{\rm{d}}{t}}=-\vec{{\rm{\nabla }}}\cdot (\rho \,\vec{u})\end{eqnarray} \tag{ B.1 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\vec{u}}{{\rm{d}}{t}}+(\vec{u}\cdot \vec{{\rm{\nabla }}})\vec{u}=-\displaystyle \frac{1}{\rho }\vec{{\rm{\nabla }}}p,\end{eqnarray} \tag{ B.2 }$$

where $\vec{u}(\vec{r},t)$ is the flow velocity, ρ is the gas density and p the gas pressure. Isentropic flow (i.e. the gas is in thermal equilibrium at all times) for an ideal gas implies [52]:

$$\begin{eqnarray}&&p=\gamma {RT}\rho \quad \mathrm{and}\quad {c}^{2}=\gamma R\,T,\end{eqnarray} \tag{ B.3 }$$

where R is the specific gas constant, T is the gas temperature and c is the speed of sound. We insert equation (B.3) in (B.2) and linearize for small density fluctuations, $\rho (\vec{r})={\rho }_{0}+\alpha (\vec{r})$ with $\alpha \ll {\rho }_{0}$, and recover equations (3) and (4) in the main text.

Here we show how quantities and equations from 2D-electrostatics can be mapped to vortex-induced flow fields, and how acoustics of an ideal gas is mapped to third-sound dynamics.

Table C1.  Electrostatics and vortex flow field. The system is invariant under z-translation, hence we use units and equations in two dimensions.

2D-electrostatics	Vortices
Electric displacement field	Velocity field
$\vec{D}(\vec{r})$ (C m−2)	${\vec{v}}_{v}(\vec{r})$ (m s−1)
Electric line charge	Circulation quantum
Q (C m−1)	κ (m2 s−1)
Gauss's law	Vortex flow equation
$\oint \vec{D}\cdot {\rm{d}}\vec{n}=Q$	$\oint {\vec{v}}_{v}\cdot {\rm{d}}\vec{l}=\kappa$
Perfect electric conductor (ground)	Tangential flow boundary
$\vec{D}\times \vec{n}=0$	${\vec{v}}_{v}\cdot \vec{n}=0$

Table C2.  Acoustics, sound dynamics in a Bose–Einstein condensate in the Thomas-Fermi limit at zero temperature, and third-sound dynamics on a helium thin film. As in table C1, a two-dimensional system is described.

2D-acoustics	Sound in 2D-BEC	Third-sound dynamics
Density perturbation	Density perturbation	Third-sound amplitude
$\alpha (\vec{r},t)$ (kg m−2)	$\eta (\vec{r},t)$ (kg m−2)	$\eta (\vec{r},t)$ (m)
Static density	Static density	Unperturbed film height
ρ0 (kg m−2)	ρ0 (kg m−2)	h0 (m)
Background flow	Irrotational vortex flow	Irrotational vortex flow
${u}_{0}(\vec{r})$ (m s−1)	${v}_{v}(\vec{r})$ (m s−1)	${v}_{v}(\vec{r})$ (m s−1)
Irrotational flow velocity	Sound flow velocity	Third-sound flow velocity
$\delta \vec{u}(\vec{r},t)$ (m s−1)	${\vec{v}}_{s}(\vec{r},t)$ (m s−1)	${\vec{v}}_{3}(\vec{r},t)$ (m s−1)
Static pressure	Atom–atom coupling	Linearized VdW coefficent
p0 (J m−2)	gBEC (J m2)	$g=\tfrac{3{\alpha }_{\mathrm{vdw}}}{{h}_{0}^{4}}$(m s−2)
Speed of sound (acoustics)	Bogoliubov sound velocity	Speed of sound (thin film)
$c=\sqrt{\gamma {RT}}$ (m s−1)	$c=\sqrt{{g}_{\mathrm{BEC}}\cdot {\rho }_{0}/{M}^{2}}$ (m s−1)	${c}_{3}=\sqrt{g\cdot {h}_{0}}$ (m s−1)
Fixed wall boundary	Fixed wall boundary	Free boundary
$\vec{u}\cdot \vec{n}=0$	$\vec{v}\cdot \vec{n}=0$	$\vec{v}\cdot \vec{n}=0$
Fixed pressure boundary	Fixed density boundary	Fixed boundary
p = p0	η = 0	η = 0
Continuity equation (acoustics)	Continuity equation (BEC)	Continuity equation (thin film)
$\dot{\alpha }=-{\rho }_{0}\vec{{\rm{\nabla }}}\cdot \vec{u}-\vec{u}\vec{{\rm{\nabla }}}\alpha$	$\dot{\eta }=-{\rho }_{0}\vec{{\rm{\nabla }}}\cdot \vec{v}-\vec{v}\vec{{\rm{\nabla }}}\eta$	$\dot{\eta }=-{h}_{0}\vec{{\rm{\nabla }}}\cdot \vec{v}-\vec{v}\vec{{\rm{\nabla }}}\eta$
Linearized Euler (acoustics)	Linearized Euler (BEC)	Linearized Euler (thin film)
$\dot{\vec{u}}+(\vec{u}\cdot \vec{{\rm{\nabla }}})\vec{u}=-\tfrac{\gamma {RT}}{{\rho }_{0}}\vec{{\rm{\nabla }}}\alpha$	$\dot{\vec{v}}+(\vec{v}\cdot \vec{{\rm{\nabla }}})\vec{v}=-\tfrac{1}{M}\vec{{\rm{\nabla }}}\left(U+\tfrac{\eta }{M}{g}_{\mathrm{BEC}}\right)$	$\dot{\vec{v}}+(\vec{v}\cdot \vec{{\rm{\nabla }}})\vec{v}=-g\vec{{\rm{\nabla }}}\eta$

D.1. Boundary conditions

In order to solve differential equations on the surface of a two-dimensional resonator, constraints at the boundary have to be specified. Depending on the type of confinement provided by the resonator, the boundary for third sound can be described either by a fixed ('Dirichlet') or a free ('Neumann') boundary condition. A fixed boundary condition η = 0 allows flow in and out of the resonator and the film height at the boundary is fixed to the equilibrium film height. The free boundary condition ${\vec{v}}_{3}\cdot \vec{n}=0$, where $\vec{n}$ is the normal vector on the boundary, allows film height fluctuations at the boundary and prohibits flow in or out of the resonator. In COMSOL, for an ideal gas, the free boundary condition corresponds to a rigid wall, where volume is conserved and the gas pressure can oscillate freely at the boundary. The fixed boundary condition corresponds to fixed pressure, where the gas pressure is fixed at the boundary and the gas can freely flow in and out of the domain [84]. The vortex flow is tangential to the boundary, ${\vec{v}}_{v}\cdot \vec{n}=0$. In the electrostatics analogue, this translates to an electric field which is exactly perpendicular to the boundary, with no tangential component. This corresponds to a perfect electric conductor at the boundary and can be realized by choosing the ground—boundary condition in COMSOL [85]. In order to model a quantized circulation n × κ around a topological defect in the structure, the floating potential boundary condition with built in charge Q = n × κ must be chosen. This boundary condition enforces an electric field orthogonal to the boundary everywhere (due to the equal potential on the boundary), as well as the condition:

$$\begin{eqnarray}&&\oint \vec{D}\cdot \vec{n}\,{dl}=Q.\end{eqnarray} \tag{ D.1 }$$

Upon the substitution of equation (8), this corresponds to a superfluid flow always parallel to the topological defect boundary (i.e. no fluid inflow or outflow), and the quantized circulation condition:

$$\begin{eqnarray}&&\oint \vec{v}\cdot {\rm{d}}\vec{l}=n\,\kappa .\end{eqnarray} \tag{ D.2 }$$

D.2. Notes on implementation in COMSOL^® multiphysics

In the following we describe how superfluid helium thin film can be modelled using the FEM solver COMSOL^® multiphysics 5.0.

A 2D model is set up. The Electrostatics(es) module is used to simulate vortices and a stationary study is created. The resonator outer boundary is set to ground. The circulation quantum κ is defined with adjusted SI-units (${\rm{m}}\,{{\rm{s}}}^{-2}\to {\rm{C}}\,{{\rm{m}}}^{-1}$). At each position where a clockwise vortex is to be modelled, a line charge (out-of-plane) of Q_L = κ is inserted. A counter-clockwise vortex can be modelled by replacing $\kappa \to -\kappa$.

To model third sound, the Aeroacoustics $\to$ Linearized Euler, Frequency Domain(lef) module is added and a Eigenfrequency step is included in the study. In the first, stationary, study step, only the electrostatics interface is solved for, whereas in the second step the Eigenfrequency solver is applied to the acoustics interface. Parameters ρ_sf, A_vdw, h₀ and $g\,=\,3{\alpha }_{\mathrm{vdw}}{h}_{0}^{-4}$ [49] are defined to set the superfluid density ($145\,{\rm{k}}{\rm{g}}\,{{\rm{m}}}^{-3}$ for superfluid helium [86]), the Hamacker constant of the substrate ($2.6\times {10}^{-24}\,{{\rm{m}}}^{5}\,{{\rm{s}}}^{-2}$ for silica [63]), the film thickness and the linearized Van-der-Waals acceleration, respectively. The product RTγ = c² is set to g · h₀ for a uniform superfluid film. Alternatively, a spatially varying function can be defined to reflect a non-uniform film thickness. The boundary condition is set to either rigid wall or fixed pressure (see appendix D.1).

In order to include vortices defined in the Electrostatics(es) module, a critical velocity is defined (v_crit ≈ 60 m s⁻¹ for superfluid helium [87]), and the acoustic background flow field ${\vec{u}}_{0}$ is set to:

$$\begin{eqnarray}&&\left(\begin{array}{c}{u}_{0,x}\\ {u}_{0,y}\end{array}\right)=\left(\begin{array}{c}{D}_{y}\\ -{D}_{x}\end{array}\right),\end{eqnarray} \tag{ D.3 }$$

(where D is the electric displacement field solved for in the first step stationary solver), and truncated at u_0,max = v_crit. This vortex background flow field is treated as constant in time. The FEM-solver computes the velocity perturbation $\delta \vec{u}(t)$, corresponding to the third sound mode, to the total flow: $\vec{u}(\vec{r},t)={\vec{u}}_{0}(\vec{r})+\delta \vec{u}(\vec{r},t)$.

The gas density perturbation $\alpha (\vec{r},t)$ calculated in COMSOL^® can be converted to third sound amplitude by a normalization factor: $N=\alpha (\vec{r},t)/\eta (\vec{r},t)$. It is extracted from

$$\begin{eqnarray}&&{E}_{\mathrm{pot},3}(\eta )={k}_{{\rm{B}}}{T}_{\mathrm{mode}},\end{eqnarray} \tag{ D.4 }$$

where T_mode is the mode temperature. An analytical expression for the potential energy of a third sound mode E_pot,3 is given in [63]. For the simple case of a uniform film thickness and free boundary conditions the conversion is given by

$$\begin{eqnarray}&&N=\displaystyle \frac{3{A}_{\mathrm{vdw}}{\rho }_{\mathrm{sf}}}{2{{\rm{d}}}^{4}{k}_{{\rm{B}}}{T}_{\mathrm{mode}}}\cdot \int {{\rm{d}}}^{2}\vec{r}\,\rho (\vec{r}).\end{eqnarray} \tag{ D.5 }$$

T_mode is the effective temperature of the sound mode, which when thermalized with its environment corresponds to the fridge temperature. It can also be tuned through optomechanical laser heating/cooling [30], increased through laser absorption heating [44] or electrical excitation [31].