# Fast and Accurate Evaluation of Biot-Savart Integrals over Spatial Curves in Periodic Domains

Juan Ignacio Polanco

► To cite this version:

Juan Ignacio Polanco. Fast and Accurate Evaluation of Biot-Savart Integrals over Spatial Curves in Periodic Domains. SIAM Journal on Scientific Computing, 2025, 47 (4), pp. B822-B845. 10.1137/24M1684682. hal-04783309v2

HAL Id: hal-04783309

https://hal.science/hal-04783309v2

Submitted on 1 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
applied to problems as varied as the simulation of self-propelled swimmers in viscous fluids [36]. Note that, in viscous fluids, vortex methods need to account for effects including vorticity diffusion and energy dissipation due to viscosity. To achieve this, these methods typically deal with vortex particles, which should not be interpreted as physical objects but as an ensemble of point charges generating a fluid flow. In particular, the connectivity of vortex particles is not relevant to vortex particle methods [51].

Here we focus on the conceptually simpler application of vortex methods to superfluid helium-4, where viscous effects are absent near absolute zero and vortex filaments are in fact the main physical object of interest. In this context, the Biot-Savart law is the basic ingredient of the vortex filament model (VFM), which is one of the most common approaches for describing superfluid flows [77, 41, 13]. This model is valid at scales much larger than the atomic vortex thickness, and is therefore well adapted for describing macroscopic vortex motion. Numerically, the standard approach for representing vortex filaments consists in discretising them as a series of connected vortex points. The connectivity is required in order to obtain derived geometrical information such as local tangents and curvatures. Each vortex point evolves in time according to the Biot-Savart law. In fact, the integral (1.1) diverges when evaluated on a vortex point $x = s_0$. The singularity can be avoided by taking into account the finite (but small) vortex core radius $a$ [77], which introduces a logarithmic dependence of the velocity on $a$ (see section 2). It is worth mentioning that this model shares many similarities with slender body theory [16, 54, 63], which describes the dynamics of thin fibers in Stokes flows. Similarly to the VFM, the forces induced by a 3D fiber on the fluid are approximated by a distribution of singularities along the fiber centreline. Furthermore, the velocity induced by a fiber on itself also includes a logarithmic dependence on its small thickness $a$.

This work is motivated by the study of quantum turbulence [13], which is a state of superfluid flows characterised by a wide range of energetically active length scales. Numerically, to investigate such a turbulent state, one needs to compute the non-local interactions between large numbers of vortices, requiring in particular a large number $N$ of discrete vortex points. If one explicitly accounts for all pair interactions, obtaining the velocities at the $N$ points has a $O(N^2)$ cost, which quickly limits the size of the systems which can be numerically studied. In the context of particle simulations, this problem has been solved for a long time, using techniques such as Barnes-Hut (BH) trees [14] or the fast multipole method (FMM) [38], which reduce the complexity to $O(N \log N)$ and $O(N)$ respectively. Such techniques have also been applied to vortex methods for classical [25] and quantum [6] fluids.

Here we are interested in periodic infinite systems, in which a finite set of $N_f$ vortices is replicated an infinite number of times in a spatially periodic fashion. Periodic boundary conditions are commonly used to model a variety of physical systems when one wants to describe phenomena far from boundaries. A natural way of dealing with periodicity in Cartesian domains is via a Fourier series representation, as done for example in Fourier pseudo-spectral methods [22, 19]. These take advantage of the fast Fourier transform (FFT) to efficiently evaluate non-local operations. This suggests the idea of using FFTs to evaluate costly far-field interactions between particles or vortices in periodic systems. However, since particles and vortices are respectively represented as 0D and 1D singularities, one cannot directly describe the associated source fields (e.g. electric charge density or vorticity) using a truncated Fourier series, which in practice precludes the use of FFTs. In particle simulations, one way around this issue is provided by *fast Ewald summation* methods, which are commonly used in molecular
dynamics simulations to speed-up the evaluation of electrostatic interactions between charged particles in periodic systems [32, 46, 26, 29, 4]. There, the basic strategy is to additively split the singular source field onto a smooth field responsible for far-field interactions and a correction field accounting for interactions between nearby particles. Fast Ewald summation methods are characterised by a $O(N \log N)$ complexity in the number of particles $N$. In fluid dynamics contexts, fast Ewald methods have been used to describe 2D point vortex dynamics [61], fluid motions induced by localised forces in the Stokes regime [73, 56, 59], and Darcy flows at the interface between two immiscible fluids [1]. In molecular dynamics simulation benchmarks [3], it has been observed that FFT-based Ewald methods can perform slightly better than the periodic FMM at the same accuracy, actually displaying near linear complexity over a wide range of problem sizes $N$. Moreover, when running at low accuracy levels, the latter can display a slight energy drift over time which is not observed in the former [3].

The aim of this paper is to adapt Ewald methods to the evaluation of the Biot-Savart integral (1.1) in three-dimensional periodic systems and evaluate the relevance of this approach. The paper begins in section 2 with an introduction to the VFM in the context of quantum vortex dynamics. We also propose an approach for accurately estimating the kinetic energy in periodic systems. Section 3 describes the Ewald-based method used to evaluate Biot-Savart integrals. Our approach takes advantage of the nonuniform fast Fourier transform (NUFFT) algorithm to speed-up computations. In section 4, we provide analytical estimates of the approximation errors incurred by the proposed approach in terms of tunable parameters. The relevance of these estimates is then verified in section 5 using numerical experiments of different test cases. That section also showcases the accuracy of energy and impulse conservation by the method, and finishes with numerical evidence of near linear $O(N \log N)$ complexity over a wide range of problem sizes $N$. Finally, section 6 is devoted to conclusions.

**2. The vortex filament model.** The VFM introduced by Schwarz [77] is one of the main approaches used to describe theoretically and numerically the three-dimensional hydrodynamics of quantum fluids such as low-temperature liquid helium [41, 13]. Unlike classical fluids, quantum fluids near absolute zero are characterised by having zero viscosity and being irrotational almost everywhere. In fact, rotational motions are confined to very thin vortex filaments carrying a quantised circulation $\kappa = h/m$, where $h$ is Planck's constant and $m$ the mass of one atom in the case of helium-4. In other words, a straight vortex filament induces a circular motion of the fluid around it with velocity $v = \kappa/(2\pi r)$, where $r$ is the distance to the vortex.

We restrict our attention to the zero temperature limit. Indeed, in finite temperature superfluid helium, quantum vortices can be seen as coexisting and interacting with a viscous normal fluid [12]. This is described by more complex models (some of them based on the VFM) which are still the subject of active research [13]. A second very important aspect which is not discussed here is the reconnection of vortices when they are close to collision. While this phenomenon is crucial for describing energy dissipation and quantum turbulence, its modelling is orthogonal to the subject of this work, and can be disregarded as long as vortex elements stay sufficiently far from each other.

**2.1. Filaments as oriented curves.** In the VFM, one is interested in describing the motion induced by a collection of vortex filaments on themselves. Each filament is represented as an oriented curve {$s(\xi), 0 \le \xi \le L$} where $s(\xi)$ is a location on the curve, $\xi$ is the arc length, and $L$ is the total length of the curve. Here the curve is represented using the natural arc length parametrisation, such that the unit tangent
to the curve is $s'(\xi) = ds/d\xi$ with $|s'(\xi)| = 1$ for all $\xi$. In practice, it is often more convenient to deal with arbitrary parametrisations which will be denoted $s(\tau)$ for $0 \le \tau \le T$, such that $|ds/d\tau| \ne 1$ in general.

We will consider the vortex filaments to be embedded in a triply periodic domain, so that every filament is replicated an infinite number of times along each Cartesian direction. For simplicity, throughout this paper, the domain period is set to $L$ in all directions, but all definitions and results can be readily generalised to different periodicities in each direction. In fluid dynamics, Helmholtz' theorems [45] state that vortex lines cannot end in the fluid. Therefore, in the absence of solid boundaries, they must either be closed curves or extend to infinity. Hence we only consider these two cases, with the additional restriction that infinite curves must be described by a periodic function matching the domain periodicity. Both cases are defined by the property $s(T + \tau) = s(\tau) + nL$ with $n \in \mathbb{Z}^3$. In particular, closed curves satisfy this property with $n = 0$. An example of an infinite periodic curve is $s(\tau) = [\cos(2\pi\tau/L), 0, \tau]$ with $T = L$. This curve satisfies the above property with $n = [0, 0, 1]$, extending infinitely along the third Cartesian direction.

**2.2. The Biot–Savart law.** We now consider a set of $N_f$ vortex filaments and their periodic images, and denote $C = \bigcup_{j=1}^{N_f} C_j$ the set of spatial curves $C_j$ representing their locations. The vorticity field $\omega = \nabla \times v$ associated to these filaments is then

$$ (2.1) \qquad \omega(x) = \kappa \sum_{n \in \mathbb{Z}^3} \oint_C \delta(s - x + nL) \, ds $$

where $\delta(x)$ is the Dirac delta function, and the infinite sum over $n$ accounts for the periodic vortex images. Here $\kappa$ is the circulation of the vortex filaments, related to the magnitude of the velocity induced by them. Throughout this work $\kappa$ will be taken to be constant (as is actually the case of quantum vortices), but in principle one could also consider a variable circulation $\kappa(s)$.

Inverting the curl operator leads to the Biot–Savart law (1.1) describing the velocity induced by the set of vortex filaments on a point $x$. Including periodicity and finite core size effects, this law writes

$$ (2.2) \qquad v(x) = \frac{\kappa}{4\pi} \sum_{n \in \mathbb{Z}^3} \oint_C \frac{(s - x + nL) \times ds}{|s - x + nL|^3} \varphi_v(s - x + nL). $$

Here $0 \le \varphi_v(r) \le 1$ accounts for the finite radius $a$ of the vortex core. As detailed in subsection 2.3, it allows to desingularise the Biot–Savart integral when evaluating the velocity on a vortex. For all practical purposes, it is constant and equal to 1 whenever $x$ is not located on a vortex.

Note that periodicity requires the total vorticity to be zero within a periodic cell, or otherwise the curl operator cannot be inverted (the velocity diverges). This condition corresponds to $\oint_C ds = 0$. This is trivially satisfied for closed filaments, while it requires special care when dealing with infinite unclosed filaments.

**2.3. Desingularisation of the Biot–Savart integral.** In the VFM, one is generally interested in the velocity induced by the vortex filaments on the filaments themselves. In other words, one wants to evaluate (2.2) on positions $x = s_0 \in C$, where the Biot–Savart integral clearly diverges [21]. More precisely, performing a Taylor expansion of $s(\xi)$, one can show that close to $s_0 = s(\xi_0)$ the integrand behaves as [2]

$$ (2.3) \qquad \frac{[s(\xi) - s_0] \times s'(\xi)}{|s(\xi) - s_0|^3} d\xi = \frac{s'_0 \times s''_0}{2|\xi - \xi_0|} d\xi + O(1), $$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Discretisation of a closed vortex filament, and distinction between local and non-local contributions to the Biot-Savart law. Filled grey circles represent vortex discretisation points. The velocity induced on the vortex point $s_0$ includes the local contribution (2.4) of its adjacent segments ($C_0$, in orange), and the non-local contribution of the rest ($\bar{C}_0$, in grey). Periodic images and other filaments (not depicted) also contribute to the non-local velocity. Near the evaluation point $s_0$, the lengths $\ell^-$ and $\ell^+$ of the adjacent segments are represented, as well as the local unit tangent $s_0'$ and curvature vector $s_0''$. Numerically, non-local contributions are computed by approximating line integrals by quadrature sums over each non-local segment. As an example, the open red circles between points $s_j$ and $s_{j+1}$ are the quadrature points associated to that segment (here, using $Q=3$ point Gauss-Legendre quadratures on each segment). These are the points where the integrand is actually evaluated when computing non-local contributions.</figcaption></figure>

where primes denote derivatives with respect to the arc length $\xi$. In particular, $s'_0 \equiv s'(\xi_0)$ and $s''_0 \equiv s''(\xi_0)$ are respectively the unit tangent and the curvature vectors at $s_0$.

The divergence of the Biot-Savart integral is unphysical since (2.1) does not account for the finite (but very small) radius $a$ of the vortex core. This issue can be circumvented in a physically consistent manner by defining the $\varphi_v$ term appearing in (2.2) such that $\varphi_v(\mathbf{r}) = 0$ for $|\mathbf{r}| < \epsilon_v a$, and 1 otherwise, where $\epsilon_v$ is a cut-off coefficient. In other words, $\varphi_v$ is defined such that locations $s$ close to the singularity at $s_0$ are omitted from the integral<fnref n="1" /> [72, Ch. 11]. It is convenient to write $\epsilon_v = e^\Delta/2$, where $\Delta$ is a constant which depends on the actual vorticity profile within the vortex core. This constant can be analytically derived for commonly used vorticity profiles, under the assumption that the local curvature radius $1/|s''|$ stays much larger than the core size $a$. In particular, $\Delta = 1/4$ for a circular vortex core with uniform vorticity [72, Ch. 10].

**2.4. Discretisation of the Biot-Savart integral.** In practice, the VFM is used to describe length scales which are several orders of magnitude larger than $a$. In numerical simulations, vortex filaments are discretised with a typical line resolution $\ell \gg a$. Therefore, directly computing the Biot-Savart integral with a cut-off of order $a$ as described above can lead to large numerical error, especially when quadrature rules are used to estimate the integrals. For this reason, the usual approach is to split the evaluation of the Biot-Savart integral on a discretised filament onto *local* and *non-local* contributions [77], as represented in Figure 1 for a single filament. The local contribution is estimated by analytically integrating the leading-order term of the Taylor expansion (2.3) on the local portion $C_0$ of the discretised curve and excluding

<footnote n="1">¹Interestingly, in the context of Stokes flows, slender body theory has been recently reformulated using a similar cut-off procedure [59].</footnote>
the cut-off region of length $2\epsilon_v a$. This leads to

$$ (2.4) \qquad \mathbf{v}_{\text{local}}(\mathbf{s}_0) = \frac{\kappa}{4\pi} \mathbf{s}_0' \times \mathbf{s}_0'' \left[ \ln\left(\frac{2\sqrt{\ell^-\ell^+}}{a}\right) - \Delta \right], $$

where $\ell^-$ and $\ell^+$ are the lengths of the adjacent segments composing $C_0$. The local velocity is thus oriented along the local binormal direction and is proportional to the local curvature $|s_0''|$. As illustrated in Figure 1, when filaments are discretised by connecting vortex points and $\mathbf{s}_0$ is one of these points, the adjacent segments are usually defined as those linking $\mathbf{s}_0$ to its neighbours. Finally, accounting for the local term, the total induced velocity of a point $\mathbf{s}_0$ on a discretised vortex can be written as

$$ (2.5) \qquad \mathbf{v}(\mathbf{s}_0) = \mathbf{v}_{\text{local}}(\mathbf{s}_0) + \frac{\kappa}{4\pi} \sum_{n \in \mathbb{Z}^3} \int_C' \frac{(\mathbf{s} - \mathbf{s}_0 + n\mathbf{L}) \times d\mathbf{s}}{|\mathbf{s} - \mathbf{s}_0 + n\mathbf{L}|^3}, $$

where the prime over the integral indicates that integration is omitted on the local portion $C_0$ of the filament containing $\mathbf{s}_0$.

**2.5. Energy estimation in periodic systems.** The VFM described above includes no dissipation mechanisms, and energy is expected to be a conserved quantity. It is crucial to be able to accurately estimate the total energy of the system, as this may allow (i) to verify the accuracy of energy conservation at the numerical level or (ii) to estimate energy decay rates when the model is extended with dissipative mechanisms.

In a periodic domain, the kinetic energy per unit mass is defined as

$$ (2.6) \qquad E = \frac{1}{2V} \int_{\Omega} |\mathbf{v}(\mathbf{x})|^2 d^3\mathbf{x} $$

where $\Omega = [0, L]^3$ is the main periodic cell and $V = L^3$ is its volume. In principle, this expression requires knowing the velocity field $\mathbf{v}(\mathbf{x})$ at every point $\mathbf{x}$ in space. In the VFM this is not only impractical – as it would require many evaluations of the Biot-Savart law (2.2) – but is also delicate since the velocity field presents strong gradients near vortex filaments. For these reasons, in non-periodic VFM simulations [75, 5, 40] the above expression is commonly replaced by<fnref n="2" />

$$ (2.7) \qquad E_{\text{non-periodic}} = \kappa \oint \mathbf{v}(\mathbf{s}) \cdot (\mathbf{s} \times d\mathbf{s}), $$

which only requires knowing the vortex geometry and the velocity of the vortex filaments. It is based on a similar expression initially derived as a volume integral [53, p. 218] (see also [72, p. 68]), which is expected to be valid for spatially smooth vorticity fields. It assumes that the velocity and vorticity fields decay to zero at infinity, which is not the case in periodic domains. More importantly, as we shall see below, it does not correctly account for the vortex core size $a$ and thus it is not properly conserved.

An alternative expression for the kinetic energy per unit mass can be obtained from (2.6) using integration by parts [72, p. 68],

$$ (2.8) \qquad E = \frac{1}{2V} \int_{\Omega} \psi(\mathbf{x}) \cdot \omega(\mathbf{x}) d^3\mathbf{x} = \frac{\kappa}{2V} \oint_C \psi(\mathbf{s}) \cdot d\mathbf{s}, $$

<footnote n="2">²Here we express the energy is per unit *density* and not *mass*, since the domain volume $V$ (and thus the fluid mass) is infinite.</footnote>
where the last equality is obtained using (2.1). Note that the boundary terms of the integration by parts vanish in periodic domains. In (2.8), we have introduced the streamfunction vector $\psi(x)$ (or vector potential [33, 72]), which is related to the velocity by $\nabla \times \psi = v$, and thus to the vorticity by Poisson's equation $-\nabla^2 \psi = \omega$. In other words, the kinetic energy can be obtained from a line integral requiring knowledge of streamfunction values on vortex filaments. Interestingly, (2.8) also allows to interpret the tangential streamfunction, $\psi_t = \psi \cdot s'$, as the linear energy density of a vortex point (up to a multiplicative constant $\kappa/2V$), and thus as the contribution of a vortex element to the total kinetic energy. Finally, note that (2.8) is also valid in non-periodic unbounded domains under the same assumptions leading to (2.7). To our knowledge, this expression has never been used before in the context of the VFM. Perhaps one of the reasons is that $\psi(s)$ is a priori not available in VFM simulations, and computing it comes at an additional cost.

**2.6. Obtaining the streamfunction vector.** As mentioned above, the streamfunction is the solution of the Poisson equation $-\nabla^2 \psi = \omega$. In three dimensions, the solution can be explicitly written as a convolution of $\omega$ with the Green's function $G(x, y) = G(x - y) = 1/(4\pi|x - y|)$,

$$ (2.9) \qquad \psi(x) = (G * \omega)(x) = \frac{\kappa}{4\pi} \sum_{n \in \mathbb{Z}^3} \oint_C \frac{ds}{|x - s + nL|} \varphi_\psi(x - s + nL), $$

where once again we have used (2.1) to express the vorticity field. In analogy with (2.2), here $\varphi_\psi$ accounts for the finite radius $a$ of the vortex core. Noting that $\nabla G(r) = -r/(4\pi r^3)$, where $r = |r|$, one can readily show that taking the curl of (2.9) leads to the Biot-Savart law (2.2) when $x \notin C$ (in which case $\varphi_\psi = \varphi_v = 1$ everywhere).

Similarly to the velocity (subsection 2.3), one must properly define the regularisation term $\varphi_\psi$ in (2.9) to avoid the integral from diverging when evaluating it on a filament location $x = s_0 \in C$. As for the velocity, we define it as $\varphi_\psi(r) = 0$ for $|r| < \epsilon_\psi a$, and 1 otherwise. Crucially, the cut-off coefficient $\epsilon_\psi$ need not be equal to the coefficient $\epsilon_v = e^\Delta/2$ used for the velocity. In fact, setting it to $\epsilon_\psi = e^{\Delta-1}/2$ allows to identify the resulting energy (2.8) with a Hamiltonian for the vortex system as argued at the end of this section. As verified numerically in subsection 5.2, this Hamiltonian is accurately conserved in non-dissipative VFM simulations.

Using the same notation as in (2.3), the integrand of (2.9) behaves close to $s_0 = s(\xi_0)$ as

$$ (2.10) \qquad \frac{s'(\xi)}{|s(\xi) - s_0|} d\xi = \frac{s'_0}{|\xi - \xi_0|} d\xi + O(1). $$

Then, in analogy with the local velocity (2.4), integrating over $\xi \in [\xi_0 - \ell^-, \xi_0 - \epsilon_\psi a] \cup [\xi_0 + \epsilon_\psi a, \xi_0 + \ell^+]$ ($C_0$ segments in Figure 1) leads to the local contribution

$$ (2.11) \qquad \psi_{\text{local}}(s_0) = \frac{\kappa s'_0}{2\pi} \left[ \ln\left(\frac{2\sqrt{\ell^--\ell^+}}{a}\right) - (\Delta - 1) \right]. $$

Interestingly, this term is tangent to the filament, and therefore it fully contributes to the kinetic energy (2.8). In the context of vortex dynamics, this local contribution is commonly referred to as vortex tension [62, 13]. Finally, in analogy with (2.5), the streamfunction vector on a discretised filament can be written as

$$ (2.12) \qquad \psi(s_0) = \psi_{\text{local}}(s_0) + \frac{\kappa}{4\pi} \sum_{n \in \mathbb{Z}^3} \int_C' \frac{ds}{|s - s_0 + nL|}, $$
where once again the prime over the integral denotes the omission of local vortex elements around $s_0$ ($C_0$ segments in Figure 1).

To justify the above choice of $\epsilon_\psi$, we consider a circular vortex ring of radius $R \gg a$. From (2.2) (see also [77]), its self-induced translational velocity is

$$ (2.13) \qquad V = \frac{\kappa}{4\pi R} \left[ \ln \left( \frac{8R}{a} \right) - \Delta \right]. $$

The total energy per unit mass $E$, satisfying Hamilton's equation $V = \partial E / \partial p$ where $p = (\kappa/V)\pi R^2$ is the vortex ring impulse, is then given by [71, 80]

$$ (2.14) \qquad E = \frac{\kappa^2 R}{2V} \left[ \ln \left( \frac{8R}{a} \right) - (\Delta + 1) \right]. $$

One can analytically show that $\epsilon_\psi = e^{\Delta-1/2}$ is the only possible choice allowing to obtain (2.14) from expression (2.8). Later in subsection 5.2, we verify that this choice leads to proper energy conservation up to numerical accuracy. This appears to be the case not only for circular vortex rings, but also for complex configurations containing multiple closed vortices of arbitrary shape, and even unclosed vortices extending to infinity (as in Figure 2, centre), in which case energy computations have been deemed impractical in the past [5]. Finally, note that the alternative energy expression (2.7) commonly used in open non-periodic systems is not consistent with Hamilton's equation for a circular ring, as it corresponds to replacing the $\Delta + 1$ term in the ring energy (2.14) with $\Delta$. This suggests that (2.7) does not properly account for the structure of the vortex core.

**2.7. Discretisation of spatial curves and line integrals.** We finish this section by discussing the approach we adopt to represent vortex filaments and estimate line integrals in numerical simulations. In the context of vortex filament simulations for describing quantum fluids, filaments have been almost invariably described by a set of points in space connected by straight segments [77, 74, 6, 40, 87]. Integration of the Biot-Savart law over lines is then straightforward and can even be done analytically over each straight segment. Curve derivatives (local tangents, curvatures) are estimated on the discrete points using e.g. finite difference approximations [6]. While this approach can be implemented with relative ease, it is a low-order discretisation, often requiring very small distances between discretisation points to achieve accurate results. It also presents an inconsistency between the $C^0$ continuity of the lines and the $C^2$ continuity required to estimate curvatures.

Here we take a different route and consider the filaments as smooth curves passing through a set of points, as illustrated in Figure 1. As introduced in subsection 2.1, each curve is parametrised as $s(\tau)$ for some scalar parameter $\tau \in [0, T]$. The numerical degrees of freedom are the positions of discrete vortex points or nodes. Some interpolation method is then applied to evaluate vortex positions in-between nodes. The interpolation also gives direct access to curve derivatives along vortex lines, allowing to estimate local tangent and curvature vectors required for Biot-Savart computations. In the numerical experiments of section 5 we use quintic spline interpolations which have global continuity $C^4$ (see Appendix B for details).

In the context of the Biot-Savart problem, a second difficulty is the need to estimate integrals on vortex lines. Here we estimate line integrals using Gauss-Legendre quadratures on each segment connecting two neighbouring nodes. Concretely, if one considers a closed curve $C$ parametrised as $s(\tau)$ with nodes $\{s_j = s(\tau_j); j = 1, \dots, N\}$,
then an integral over such curve is approximated as

$$ (2.15) \qquad \int_0^T F[s(\tau)] d\tau \approx \sum_{j=1}^N \sum_{i=1}^Q w_i F[s_j^{(i)}] \Delta\tau_j, \quad s_j^{(i)} \equiv s(\tau_j + h_i \Delta\tau_j), $$

where $Q$ is the number of quadrature points per segment (typically one chooses $Q \le 4$ in simulations). Here, $F[s]$ is a functional which may depend on curve locations $s(\tau)$ as well as its derivatives. For convenience, the quadrature rule above is defined in the domain $[0, 1]$, such that $h_i \in [0, 1]$ are the quadrature locations, and $w_i$ the associated weights satisfying $\sum_{i=1}^Q w_i = 1$. Moreover, $\Delta\tau_j = \tau_{j+1} - \tau_j$ is the parameter increment associated to a single segment, with the convention that $\tau_{N+1} = \tau_1 + T$ for a closed or an infinite periodic curve. As an example, the open red circles in Figure 1 are the quadrature points associated to a single filament segment (using $Q = 3$).

### 3. Fast Ewald summation for the Biot-Savart problem.
In this section we adapt fast Ewald summation methods to the evaluation of velocity and streamfunction fields induced by a set of vortex filaments. This requires some adjustments since these methods generally deal with a scalar-valued source term (e.g. electrostatic charge) supported on discrete points, while here the source term (vorticity) is vector-valued and defined on spatial curves. In the context of the VFM, the basic idea of these methods is to split the line integrals in (2.5) and (2.12) onto short- and long-range components. Concretely, including the local terms appearing in those expressions, the streamfunction and the velocity at a vortex location $s_0$ are decomposed as

$$ (3.1) \qquad \psi(s_0) = \psi_{\text{local}}(s_0) + \psi^{(n)}(s_0) + \psi^{(f)}(s_0) - \psi_{\text{self}}^{(f)}(s_0), $$

$$ (3.2) \qquad \mathbf{v}(s_0) = \mathbf{v}_{\text{local}}(s_0) + \mathbf{v}^{(n)}(s_0) + \mathbf{v}^{(f)}(s_0) - \mathbf{v}_{\text{self}}^{(f)}(s_0), $$

where the (n) and (f) superscripts denote near- and far-field components, and $\psi_{\text{self}}^{(f)}$ and $\mathbf{v}_{\text{self}}^{(f)}$ are corrections to the latter which are discussed further below.

#### 3.1. Ewald splitting.
Ewald summation methods split the singular Green's function $G(r)$ associated to Poisson's equation into (i) a fast decaying part $G^{(n)}(r)$ accounting for short-range interactions and (ii) a slowly decaying part $G^{(f)}(r)$, which is well-behaved at $r = 0$ and describes long-range interactions. In three dimensions, this is usually achieved via the identity $\operatorname{erf}(x) + \operatorname{erfc}(x) = 1$,

$$ (3.3) \qquad G(r) = \frac{1}{4\pi r} = \frac{\operatorname{erfc}(\alpha r)}{4\pi r} + \frac{\operatorname{erf}(\alpha r)}{4\pi r} = G^{(n)}(r) + G^{(f)}(r), $$

where $\operatorname{erf}(x) = (2/\sqrt{\pi}) \int_0^x e^{-u^2} du$ and $\operatorname{erfc}(x) = (2/\sqrt{\pi}) \int_x^\infty e^{-u^2} du$ are respectively the error function and the complementary error function. Due to linearity, the solution $\psi = G*\omega$ to the Poisson equation $-\nabla^2\psi = \omega$ is accordingly split into $\psi = \psi^{(n)} + \psi^{(f)}$. Similarly, the velocity field can be decomposed as $\mathbf{v} = \nabla \times \psi = \mathbf{v}^{(n)} + \mathbf{v}^{(f)}$. Importantly, $\alpha$ in (3.3) is the Ewald splitting parameter, which defines an inverse length scale setting the transition between short- and long-range interactions. This is a purely numerical parameter, as the physical fields $\psi$ and $\mathbf{v}$ obtained by adding both contributions are in theory independent of $\alpha$. To facilitate the interpretation of the splitting (3.3), it is helpful to consider the “far-field” vorticity field defined by $\omega^{(f)} = -\nabla^2\psi^{(f)}$, which can be written as the convolution of $\omega$ with $-\nabla^2 G^{(f)}$. In fact this is a Gaussian-filtering operation, since $-\nabla^2 G^{(f)}(\mathbf{r}) = (\alpha/\sqrt{\pi})^3 e^{-(\alpha r)^2}$. In other words, $\psi^{(f)}$ is
the streamfunction associated to a coarse-grained (smoothed) version of the original (singular) vorticity $\omega$. Similarly, $v^{(f)} = \nabla \times \psi^{(f)}$ is the corresponding coarse-grained velocity field.

Taking the gradient of (3.3), one arrives to the corresponding splitting for the Biot-Savart kernel $\nabla G$,

$$ (3.4) \qquad \nabla G(\mathbf{r}) = -\frac{\mathbf{r}}{4\pi r^3} = -[g^{(n)}(\mathbf{r}) + g^{(f)}(\mathbf{r})] \frac{\mathbf{r}}{4\pi r^3} = \nabla G^{(n)}(\mathbf{r}) + \nabla G^{(f)}(\mathbf{r}) $$

with the weight functions

$$ (3.5) \qquad g^{(n)}(\mathbf{r}) = \operatorname{erfc}(\alpha r) + \frac{2\alpha r}{\sqrt{\pi}} e^{-(\alpha r)^2}, \quad g^{(f)}(\mathbf{r}) = \operatorname{erf}(\alpha r) - \frac{2\alpha r}{\sqrt{\pi}} e^{-(\alpha r)^2}. $$

It can be verified that the long-range kernels are non-singular, behaving as $G^{(f)}(\mathbf{r}) = \alpha/(2\pi^{3/2}) + O(r^2)$ and $\nabla G^{(f)}(\mathbf{r}) = -\alpha^3 \mathbf{r}/(3\pi^{3/2}) + O(r^3)$ near $r=0$. On the other hand, the short-range kernels asymptotically decay as $G^{(n)}(\mathbf{r}) \sim e^{-(\alpha r)^2}/[4\pi r(\alpha r)^{3/2}]$ and $\nabla G^{(n)}(\mathbf{r}) \sim \alpha r e^{-(\alpha r)^2}/(2\pi^{3/2}r^2)$ for large $r$.

**3.2. Estimation of short-range interactions.** The near-field streamfunction $\psi^{(n)}$ is obtained by replacing $G(\mathbf{r})$ in (2.9) with $G^{(n)}(\mathbf{r})$ defined in (3.3), resulting in

$$ (3.6) \qquad \psi^{(n)}(s_0) = \frac{\kappa}{4\pi} \sum_{n \in \mathbb{Z}^3} \int_C' f^{(n)}(s - s_0 + nL) \frac{ds}{|s - s_0 + nL|}, $$

where $f^{(n)}(\mathbf{r}) = \operatorname{erfc}(\alpha r)$. We recall that the prime over the integral symbol denotes the omission of local vortex elements around $s_0$ ($C_0$ segments in Figure 1). Similarly, following the decomposition (3.4), the near-field velocity field is explicitly given by the modified Biot-Savart integral

$$ (3.7) \qquad \mathbf{v}^{(n)}(s_0) = \frac{\kappa}{4\pi} \sum_{n \in \mathbb{Z}^3} \int_C' g^{(n)}(s - s_0 + nL) \frac{(s - s_0 + nL) \times ds}{|s - s_0 + nL|^3}. $$

The original integrals in (2.5) and (2.12) are recovered by setting $\alpha = 0$ in the above expressions. In practice, we approximate the above line integrals using quadrature sums over discrete vortex line locations $s_j^{(i)}$ according to (2.15). Moreover, since the short-range kernels $G^{(n)}$ and $\nabla G^{(n)}$ decay exponentially with $(\alpha r)^2$, one can safely define a cut-off distance $r_c$ beyond which short-range interactions can be neglected. The truncation errors associated to the choice of $r_c$ are discussed in section 4.

**3.3. Estimation of long-range interactions.** Since the long-range Green's function $G^{(f)}$ is non-singular, one can use a truncated Fourier series representation to indirectly solve $-\nabla^2 \psi^{(f)} = \omega^{(f)}$. We start by writing the periodic vorticity field as

$$ (3.8) \qquad \omega(\mathbf{x}) \approx \sum_{k \in \mathbb{K}^3} \hat{\omega}_k e^{i\mathbf{k} \cdot \mathbf{x}} \quad \text{where} \quad \hat{\omega}_k = \frac{1}{V} \int_{\Omega} \omega(\mathbf{x}) e^{-i\mathbf{k} \cdot \mathbf{x}} d^3\mathbf{x}. $$

Here $\mathbb{K} = \{\frac{2\pi m}{L}; m = -\lfloor \frac{M}{2} \rfloor, \dots, \lfloor \frac{M-1}{2} \rfloor\}$ is the set of $M$ resolved Fourier wavenumbers in each Cartesian direction and $\lfloor \cdot \rfloor$ denotes the floor operation (for simplicity, we take $M$ to be the same in all directions). We recall that $\Omega = [0, L]^3$ represents the
main periodic cell and $V = L^3$ is its volume. Since the vorticity (2.1) is singular and supported on spatial curves, its Fourier coefficients are then

$$ (3.9) \qquad \hat{\omega}_k = \frac{\kappa}{V} \oint_C e^{-i\mathbf{k}\cdot\mathbf{s}} \mathrm{d}\mathbf{s} \approx \frac{\kappa}{V} \sum_{l=1}^{N_f} \sum_{j=1}^{N_l} \sum_{i=1}^{Q} w_i \Delta\tau_{jl} s'_{ijl} e^{-i\mathbf{k}\cdot\mathbf{s}_{ijl}}, $$

where the last expression is obtained from the quadrature approximation (2.15). The outermost sum is over the $N_f$ vortex filaments of the system, each being discretised by a possibly different number of nodes $N_l$. Moreover, $s_{ijl}$ is a shorthand for $s_l(\tau_j + h_i \Delta\tau_j)$, where $s_l(\tau)$ is the parametrisation of the $l$-th filament, while $s'_{ijl}$ is the derivative with respect to $\tau$ at that location (aligned with the local tangent vector). The triple sum in (3.9) can be interpreted as a sum of vector charges $q_{ijl} = w_i \Delta\tau_{jl} s'_{ijl}$ on locations $s_{ijl}$. Note that using the same quadrature nodes for short- and long-range computations means that interpolated values of $\mathbf{s}$ and $\mathbf{s}'$ can be shared among both components, reducing the computational cost associated to interpolations.

Now, since $\omega^{(f)}$ is the convolution of $\omega$ with a Gaussian kernel, its Fourier coefficients are $\hat{\omega}_k^{(f)} = \hat{\omega}_k e^{-k^2/4\alpha^2}$ where $k = |\mathbf{k}|$. Unlike $|\hat{\omega}_k|$, the amplitudes $|\hat{\omega}_k^{(f)}|$ can be expected to decay quickly for $k \gg \alpha$, which justifies the truncation of the Fourier series. Moreover, solving Poisson's equation amounts to division by $k^2$ in Fourier space. This ultimately allows us to express the far-field velocity in physical space as

$$ (3.10) \qquad \mathbf{v}^{(f)}(\mathbf{x}) \approx \sum_{\substack{\mathbf{k} \in \mathbb{K}^3 \\ |\mathbf{k}| \neq 0}} \hat{v}_{\mathbf{k}}^{(f)} e^{i\mathbf{k}\cdot\mathbf{x}} \quad \text{with} \quad \hat{v}_{\mathbf{k}}^{(f)} = i\mathbf{k} \times \hat{\psi}_{\mathbf{k}}^{(f)} = i\mathbf{k} \times \frac{\hat{\omega}_{\mathbf{k}}^{(f)}}{k^2}. $$

Note that this requires $\hat{\omega}_0 = 0$, i.e. the mean vorticity within a periodic cell must be zero. The far-field streamfunction field $\psi^{(f)}(\mathbf{x})$ can be written similarly to (3.10) from its Fourier coefficients $\hat{\psi}_{\mathbf{k}}^{(f)} = \hat{\omega}_{\mathbf{k}}^{(f)}/k^2$.

When evaluated on a vortex location $s_0$, the obtained far-field velocity and streamfunction vectors include contributions of the local segments adjacent to $s_0$. However, integration over these segments must be excluded according to (2.5) and (2.12). This issue is analogous to the spurious self-interaction electrostatic potential appearing in standard Ewald methods [3]. In the present case, the spurious local integrals are explicitly given by

$$ (3.11) \qquad \psi_{\text{self}}^{(f)}(s_0) = \frac{\kappa}{4\pi} \int_{C_0} f^{(f)}(s - s_0) \frac{\mathrm{d}s}{|s - s_0|}, $$

$$ (3.12) \qquad \mathbf{v}_{\text{self}}^{(f)}(s_0) = \frac{\kappa}{4\pi} \int_{C_0} g^{(f)}(s - s_0) \frac{(s - s_0) \times \mathrm{d}s}{|s - s_0|^3}, $$

where $f^{(f)}(\mathbf{r}) = \text{erf}(\alpha r)$, $g^{(f)}(\mathbf{r})$ is defined in (3.5), and $C_0$ consists of the two segments adjacent to $s_0$ (see Figure 1). They can be estimated using a variant of the quadrature procedure (2.15), thus requiring $2Q$ evaluations of the integrand for each vortex location $s_0$. As indicated by (3.1) and (3.2), these two “self-interaction” terms must be subtracted from the long-range estimations to avoid an unphysical dependence of the results on the Ewald parameter $\alpha$.

As in previous works [44, 3, 66], here we adopt the NUFFT algorithm [31] to efficiently approximate the Fourier sums (3.9) and (3.10). These can be respectively
obtained using type-1 (non-uniform to uniform) and type-2 (uniform to non-uniform) NUFFTs [37]. When relying on the NUFFT, the accuracy of long-range computations is then controlled by the internal NUFFT parameters and the truncation wavenumber $k_{max} = (2\pi/L)[(M-1)/2]$. These two sources of error are respectively discussed in Appendix A and section 4.

**4. Truncation error estimates.** We now provide estimates of the root-mean-square errors associated to the short- and long-range cut-offs $r_c$ and $k_{max}$. Similar estimates have already been provided for the electrostatic problem [50, 30], but these do not directly apply to the VFM since (i) the quantities of interest are not the same, and (ii) the singular sources in the present case are spatial curves and not points. We also show that the effect of the parameters $\alpha$, $r_c$ and $k_{max}$ on accuracy can be reduced to a unique non-dimensional coefficient $\beta$, thus greatly simplifying the parameter selection procedure.

The derivation of the error estimates is detailed in Appendix C. In summary, the estimates make the simplifying assumptions that (i) vortex filaments are homogeneously distributed in the spatial domain, (ii) the relative orientation of two vortex elements at a distance $r > r_c$ is completely decorrelated, and (iii) $k_{max}$ is much larger than the typical curvature $|s''|$ of the vortices. Introducing the non-dimensional cut-off coefficients $\beta^{(n)} = \alpha r_c$ and $\beta^{(f)} = k_{max}/2\alpha$, we show that both short- and long-range errors decay exponentially with the square of these respective coefficients, which justifies reducing them into a unique coefficient $\beta$. This finally leads to the truncation error estimates

$$ (4.1) \qquad \varepsilon_v \approx \kappa \left[ \frac{\mathcal{L}}{\sqrt{\pi \alpha V}} + \frac{1}{\beta} \left( \frac{\mathcal{L}}{8\pi V} \right)^{1/2} \right] e^{-\beta^2}, $$

$$ (4.2) \qquad \varepsilon_\psi \approx \frac{\kappa}{2\alpha} \left[ \frac{\mathcal{L}}{\alpha V} + \frac{1}{\beta^{1/2}} \left( \frac{\mathcal{L}}{8\pi V} \right)^{1/2} \right] \frac{e^{-\beta^2}}{\beta^{3/2}}, $$

where the short- and long-range contributions to the errors correspond to the first and second terms within each square bracket. The accuracy of the method is thus mainly driven by a unique non-dimensional parameter $\beta$. This sets both physical- and Fourier-space cut-offs $r_c = \beta/\alpha$ and $k_{max} = 2\beta\alpha$ for a given value of the inverse splitting distance $\alpha$, which is left as a free parameter that can be adjusted to optimise performance (subsection 5.3).

**5. Numerical experiments.** We now perform different numerical experiments to evaluate the accuracy, conservation properties and computational complexity of the proposed FFT-based method. All computations are performed in double precision arithmetic. The methods are implemented in the open-source VortexPasta.jl package<fnref n="3" /> [68] written in Julia [17]. The online documentation includes installation instructions, usage examples and tutorials, a description of the model and numerical methods, hints for running on computing clusters and graphical processing units (GPUs), and an extensive API reference. This should allow a new user to quickly get started with running simulations even without any prior knowledge of the Julia language.

In all numerical tests below, the domain is periodic with period $L = 2\pi$ in all directions. The vortex core size and the vorticity profile parameter – both appearing in

<footnote n="3">³Available at https://github.com/jipolanco/VortexPasta.jl under the MPL-2.0 license.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Visualisation of considered test cases within a periodic cell of size $L^3 = (2\pi)^3$. Left: trefoil knot with characteristic size $R = \pi/3$. Centre: infinitely extended straight lines. Arrows indicate the line orientation. Right: 40 ellipses with random locations, orientations, sizes and aspect ratios. For visualisation purposes, curves crossing the domain boundaries are broken onto multiple parts, which are periodically wrapped to fit within the periodic cell. Colours allow to distinguish between individual vortex lines.</figcaption></figure>

the locally induced velocity (2.4) – are respectively set to $a = 10^{-8}$ and $\Delta = 1/4$, and the vortex circulation is set to $\kappa = 1$. Vortex lines are represented as smooth curves using quintic splines, which allows to estimate derivatives and evaluate quantities in-between discretisation points. Moreover, line integrals are estimated using $Q = 3$ quadrature points per vortex segment (see subsection 2.7).

**5.1. Accuracy.** To numerically assess the accuracy of the proposed method, we consider three different test cases of varying complexity. The test cases are detailed below and illustrated in Figure 2.

1. *Trefoil knot.* The trefoil knot consists in a single knotted curve (Figure 2, left), interacting with itself and its periodic images. The trefoil knot curve is parametrically defined by

$$ (5.1) \qquad \mathbf{X}(\theta) = R \begin{bmatrix} \sin(\theta) + 2 \sin(2\theta) \\ \cos(\theta) - 2 \cos(2\theta) \\ -\sin(3\theta) \end{bmatrix} \quad \text{for } \theta \in [0, 2\pi], $$

where $R$ determines the size of the trefoil knot. We take $R = L/6$ so that periodicity effects are non-negligible. The curve is discretised by evaluating (5.1) on $N = 512$ equally spaced values of $\theta$.

2. *Infinite lines.* This test case consists of a set of straight infinitely extended lines inducing a three-dimensional velocity field. Each periodic cell is crossed by three pairs of lines aligned positively and negatively with each Cartesian direction (Figure 2, centre). The filament orientations are such that the velocity induced on each line is non-zero and positively aligned with the line ($v \cdot s' > 0$). Since the lines have zero curvature, the velocity induced by a line on itself is zero, and therefore the velocity of a vortex point is completely due to non-local interactions. Therefore, this test case emphasises the accuracy of long-range computations. We discretise each filament with $N_l = 128$ independent nodes (for $l = 1, \dots, 6$).

3. *Random ellipses.* The third test case is closer to a disordered configuration relevant to turbulent vortex flows. It consists of 40 ellipses which are randomly positioned and oriented within the spatial domain (Figure 2, right). The size and aspect ratio of each ellipse is also random. Concretely, the minor and major radii are
<figure><img src="image_4.png" /><figcaption>FIG. 3. Accuracy of full Biot-Savart computations for three different test cases. Relative root-mean-square error of filament velocity (+) and streamfunction values (×) at varying non-dimensional cut-off parameter $\beta$. The reference solution is obtained using $\beta = 8$. The Ewald splitting parameter is kept constant, $\alpha = 24/L$. The black dashed lines represent the total velocity error estimate (4.1). The grey dashed lines represent the total streamfunction error estimate (4.2).</figcaption></figure>

random values uniformly and independently distributed in $[L/16, L/4]$. We discretise each ellipse with $N_l = 128$ nodes (for $l = 1, \dots, 40$).

**Empirical relative errors.** To numerically quantify the truncation errors in each test case, we start by computing a reference solution with very high accuracy. Concretely, we set the cut-off coefficient to $\beta = 8$ and compute NUFFTs with a nominal relative tolerance $\sim 10^{-14}$. Here, the reference solution corresponds to a set of velocity and streamfunction values, $v_i^{\text{ref}}$ and $\psi_i^{\text{ref}}$, on each discrete filament node $s_i$. We then compare this reference solution with the results $v_i$ and $\psi_i$ obtained by varying $\beta$. The splitting parameter is set to $\alpha = 24/L$ in all numerical experiments. Accuracy estimates are obtained by evaluating the relative $\ell_2$ errors

$$ (5.2) \qquad \varepsilon_v^{\text{rel}} = \frac{\|\mathbf{v} - \mathbf{v}^{\text{ref}}\|}{\|\mathbf{v}^{\text{ref}}\|} \quad \text{and} \quad \varepsilon_{\psi}^{\text{rel}} = \frac{\|\psi - \psi^{\text{ref}}\|}{\|\psi^{\text{ref}}\|}, $$

where the $\ell_2$ norm of a vector field $\mathbf{u}$ evaluated at $N$ discrete locations is defined by $\|\mathbf{u}\|^2 = \frac{1}{N} \sum_{i=1}^N |\mathbf{u}_i|^2$. Here $|\cdot|$ denotes the vector magnitude.

The obtained relative errors are plotted in Figure 3 (lines with symbols). In almost all cases, the relative error decays exponentially with the squared cut-off parameter $\beta^2$, until roundoff errors become dominant. The decays are in remarkably good agreement with the rms error estimates (4.1) and (4.2) (dashed lines), which in almost all cases provide upper bounds to the actual errors. This confirms the relevance of the truncation estimates in a variety of configurations.

**5.2. Timestepping and conservation properties.** The VFM as described in section 2 is expected to conserve the total energy (2.8) and the hydrodynamic impulse $\mathbf{p} = \frac{\kappa}{2V} \oint_C \mathbf{s} \times d\mathbf{s}$. To verify this, we perform a VFM simulation in which the position of each vortex point $\mathbf{s}$ evolves in time according to $\frac{d\mathbf{s}(t)}{dt} = \mathbf{v}(\mathbf{s}(t), t)$. Here the right hand side is the Biot-Savart velocity (2.2) evaluated at location $\mathbf{s}(t)$ and at time $t$. The test case we consider here is inspired by the classical problem of two coaxial circular vortex rings which travel in the same direction. Due to their mutual interaction, the rings are expected to continuously change their sizes and therefore their translational velocities. As a result, they pass through one another in a cyclic
fashion in “leapfrogging” motion. This phenomenon has received much attention in
classical viscous fluids [79] and has also been investigated in VFM simulations [84].

Here we consider a slightly more complex variant of the classical problem: instead of perfectly circular vortex rings we simulate elliptical ones, so that not only their size changes in time but also their shape. The two elliptical vortex rings evolve in a cubic periodic domain of size L. The vortices are initially located on two parallel planes – both orthogonal to the x axis – at a distance L/4 from each other. Their minor and major radii are respectively set to R_a = L/4 and R_b = R_a√2. Their size comparable to L means that the effect of periodic images cannot be fully neglected. Besides, the ellipses are rotated relative to each other so that their major axes are respectively aligned with the y and z axes. This initial configuration is visible in Figure 4(a).

Aiming at a relative accuracy of ~10⁻⁶, each vortex is discretised with a relatively small number of points N = 32. Such a coarse discretisation can lead to errors in the local terms (2.4) and (2.11) due to the use of Taylor expansions. For this reason, we subdivide the local filament segments (C₀ in Figure 1) so that Taylor expansions are applied on a smaller central portion (about 10% of the local segment length), while 3-point Gauss–Legendre quadratures are used in the remaining portions not in contact with the singularity. Each discrete vortex point sᵢ is evolved in time according to dsᵢ/dt = v(sᵢ). For the temporal discretisation we adopt the standard fourth-order explicit Runge–Kutta (RK4) scheme. Spline coefficients are updated each time the points sᵢ are advanced. While the vortices can stretch and shrink over time, we do not perform any remeshing of the filaments for the sake of simplicity (but this is indeed needed in more complex problems). To guarantee near 6-digit accuracy, we set β = 3.5 (Figure 3), while the splitting parameter is α = 7/L so that the short-range cut-off distance is r_c = β/α = L/2.<fnref n="4" />

To preserve stability, the maximum allowed timestep $\Delta t$ in VFM simulations is roughly proportional to $\ell_{\min}^2$ [77, 41] where $\ell_{\min}$ is the smallest spacing between vortex points – in this test case $\ell_{\min} \approx 2\pi R_a/N$. Physically, this condition can be understood by considering a small-amplitude sinusoidal perturbation of wavelength $\lambda$ on an otherwise straight vortex. The perturbation is physically expected to rotate about the vortex axis with a period $T_{kw}(\lambda) \approx \frac{2\lambda^2}{\kappa} [\ln(\frac{\lambda}{\pi a}) + \frac{1}{2} - (\Delta + \gamma)]^{-1}$, where $\gamma \approx 0.57721$ is the Euler–Mascheroni constant [81, 7]. Such periodic motion is known as a Kelvin wave. Numerically, the timestep $\Delta t$ must therefore be small enough to capture the Kelvin wave oscillations associated to the smallest resolved scale $\lambda \sim \ell_{\min}$. In practice, we find that setting $\Delta t = C T_{kw}(\ell_{\min})$ with $C = 1$ is enough for the RK4 scheme to remain stable. We finish by mentioning that the above restriction is actually imposed by the local self-interaction term (2.4) [20, 28]. This suggests the use of splitting methods [60, 18] – perhaps in combination with the Hasimoto transform [43, 10] to deal with the local term – to relax such strong restriction in future works.

Different temporal snapshots of the simulation are visualised in Figure 4(a). The two vortices travel in the x direction while passing through one another as their size and shape change. In particular, at time t = 2.5 the first vortex is going through the second one, while the opposite has just occurred at time t = 7.5. This cyclic motion is also visible in Figure 4(b), which shows the temporal variation of the total vortex length L. As expected, this quantity is not conserved, displaying relatively slow variations of

<footnote n="4">⁴The condition $r_c \le L/2$ ensures that, in short-range computations, a pair of vortex points sees each other at most once. This simplifies the implementation as one does not need to explicitly deal with periodic images. In molecular dynamics this is called the minimum image convention [4].</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Energy and impulse conservation of two leapfrogging elliptical rings. (a) Vortex geometry at different time instants. Vortices mainly advance in the x direction, visiting 4 periodic boxes (light grey cubes) in the shown simulation time. Also shown is the symmetry axis aligned with the x direction (dash-dotted line). (b) Relative change of the total vortex length $L$ over time. (c) Absolute value of the relative change of kinetic energy $E$ and axial impulse $p_x$ over time. The dashed line corresponds to an incorrect estimation of the energy $E' = E + \delta E$ discussed in the main text. Circles in (b) and (c) correspond to the snapshots in (a). Times are normalised by $R_a^2/\kappa$.</figcaption></figure>

up to 2% which are clearly correlated with the physical leapfrogging motion. These slow variations coexist with higher-frequency fluctuations which are explained by the non-circular shape of each vortex. On the other hand, as seen in Figure 4(c), the total kinetic energy (2.8) and the total hydrodynamic impulse are approximately conserved. Indeed, they display temporal variations uncorrelated with the leapfrogging motion which are well within the $10^{-6}$ relative accuracy of the calculations.

We stress that, as discussed in subsection 2.6, observing a high-accuracy energy conservation rests on the proper choice of near-singularity cut-off $\epsilon_\psi = e^{\Delta^{-1/2}}$ for the streamfunction integral (2.12). Remarkably, while this choice was justified in subsection 2.6 for circular vortex rings, it appears to stay valid for more complex vortex geometries. If instead one used the same cut-off $e^\Delta/2$ as for the velocity integral (exchanging $\Delta - 1$ with $\Delta$ in (2.11)), this would result according to (2.8) in an extra (negative) energy contribution $\delta E = -\frac{\kappa^2}{4\pi V}L$. The resulting energy $E' = E + \delta E$ would not be conserved since the vortex length $L$ is allowed to fluctuate in time. This is verified in our simulations, where $E'$ (dashed line in Figure 4(c)) displays fluctuations which are correlated with the fluctuations of $L$ and are about 3 orders of magnitude larger than those of $E$.

**5.3. Performance.** We finally investigate the computational complexity of the proposed method in terms of the total number of vortex points $N$. We start by discussing how to choose the splitting parameter $\alpha$ to achieve optimal complexity for a fixed value of the accuracy parameter $\beta$. Intuitively, $\alpha$ controls the relative amount of work done by short- and long-range computations. From the point of view of short-range computations, one may attempt to choose $\alpha$ so that it reduces their expected complexity from $O(N^2)$ to $O(N)$. Assuming a spatially homogeneous vortex point distribution, the runtime associated to the short-range part can be estimated as $T^{(n)} \propto QN^2(r_c/L)^3$. In terms of $\alpha$ and $\beta$, this becomes $T^{(n)} \propto QN^2\beta^3/(\alpha L)^3$.
<caption>TABLE 1</caption>

Parameters used in performance tests and summary of benchmark results. $\epsilon_{tol}$, nominal relative tolerance; $\beta$, cut-off parameter; $C_{\beta}$, constant prefactor in (5.3); $\sigma$ and $w$, NUFFT parameters (see Appendix A). The last two columns summarise the obtained runtimes: (1) prefactor $C_{fit}$ of the least-squares fit $T_{fit}(N) = C_{fit}N \log_{10}N$ and (2) relative fitting error $\epsilon_{fit} = \|T_{fit} - T\|/\|T\|$.

<table><thead><tr><th>ε<sub>tol</sub></th><th>β</th><th>C<sub>β</sub></th><th>σ</th><th>w</th><th>C<sub>fit</sub> (μs)</th><th>ε<sub>fit</sub> (%)</th></tr></thead><tbody><tr><td>10<sup>-3</sup></td><td>2.0</td><td>2.06</td><td>1.5</td><td>4</td><td>1.41</td><td>6.07</td></tr><tr><td>10<sup>-4</sup></td><td>2.5</td><td>1.92</td><td>1.5</td><td>6</td><td>2.07</td><td>2.73</td></tr><tr><td>10<sup>-6</sup></td><td>3.5</td><td>1.71</td><td>1.5</td><td>8</td><td>4.52</td><td>2.25</td></tr><tr><td>10<sup>-8</sup></td><td>4.0</td><td>1.64</td><td>1.5</td><td>10</td><td>6.56</td><td>4.01</td></tr><tr><td>10<sup>-10</sup></td><td>4.5</td><td>1.57</td><td>1.5</td><td>12</td><td>8.95</td><td>1.81</td></tr><tr><td>10<sup>-12</sup></td><td>5.0</td><td>1.52</td><td>1.5</td><td>14</td><td>12.17</td><td>0.70</td></tr><tr><td>10<sup>-14</sup></td><td>5.5</td><td>1.47</td><td>1.5</td><td>16</td><td>15.63</td><td>0.99</td></tr></tbody></table>

Therefore, choosing

$$ (5.3) \qquad \alpha(N) = C_{\beta}N^{1/3}/L $$

can be expected to lead to linear complexity, namely $T^{(n)} \propto QN(\beta/C_{\beta})^3$. While $C_{\beta}$ should be regarded as a constant when $N$ is varied, its optimal value may vary with the accuracy parameter $\beta$, hence the notation. The next question is how does the choice (5.3) impact the complexity of long-range computations. These are dominated by the cost of the NUFFTs, which basically consist of $\propto QN$ spreading and interpolation operations and $O(M^3 \log M)$ 3D FFTs (see Appendix A). The former are already linear in $N$ and do not depend on $\alpha$. As for the latter, the size of each 1D FFT is $M \approx k_{max}L/\pi = 2\beta\alpha L/\pi$, so that choosing $\alpha$ according to (5.3) leads to $M \propto \beta C_{\beta} N^{1/3}$ and therefore $O(N \log N)$ complexity. In conclusion, choosing $\alpha$ as in (5.3) can be expected to lead, at worst, to $O(N \log N)$ complexity of the full method.

To numerically assess the performance of the method, we adapt the random ellipses test case introduced in subsection 5.1, illustrated in Figure 2 (right). Each ellipse is now discretised with $N_l = 32$ vortex points, and the Biot-Savart velocity is evaluated on all vortex points for systems with up to 8192 vortices ($N = 2^{18}$ points). The benchmarks are executed on a single core of an Intel Core i7-12700H laptop processor using Julia 1.10.3. To speed up short-range computations, a modified cell list algorithm [58] is used to find neighbouring points, in which the spatial domain is divided into cells of size $r_c/3$ in each direction. We consider 7 different accuracy levels, determined by $\beta$ and the NUFFT parameters as detailed in Table 1. In each case, we choose $C_{\beta}$ so that the time spent in short- and long-range computations is approximately equal. Empirically, we find that $C_{\beta} \approx 2.6\beta^{-1/3}$ roughly fulfills this criterion, but the proper choice this prefactor is likely dependent on the actual implementation and on the machine where the tests are run.

In Figure 5 we present the results of the numerical experiments. Each timing is obtained as the average over 10 runs. Panel (a) shows the runtime associated to evaluating the Biot-Savart velocity on all vortex points for a wide range of problem sizes $N$. The proposed method displays near linear $O(N \log N)$ complexity for all different accuracy levels (least-square fits are detailed in Table 1). For comparison, we also show the result of a naïve computation of the Biot-Savart velocity (blue crosses), which displays the expected $O(N^2)$ scaling. Note that the naïve computation does not account for periodic boundary conditions, and thus the results are not directly
<figure><img src="image_6.png" /><figcaption>FIG. 5. Performance of Biot-Savart velocity computations. (a) Runtime as a function of the number of vortex points N. Circles: proposed method for relative tolerances $10^{-3}$, $10^{-6}$, $10^{-10}$ and $10^{-14}$ (from darker to lighter). Dashed lines represent least-squares fits $T_{fit}(N) = C_{fit}N \log_{10}N$; see Table 1 for values of $C_{fit}$. Blue crosses: naïve Biot-Savart computation in non-periodic domain. (b) Influence of requested accuracy on runtimes, for different values of N. Except for the rightmost circles, errors are relative to the highest accuracy case ($10^{-14}$ relative tolerance). Dotted vertical lines represent the nominal relative tolerances (Table 1). Dashed black lines correspond to the empirical scaling $T \propto p_{\epsilon}^2 N \log_{10} N$ where $p_{\epsilon} = -\log_{10} \epsilon_{tol}$ is the estimated number of precision digits.</figcaption></figure>

comparable to those obtained by the proposed method.

From Figure 5(a), it is clear that increasing the accuracy level can have an important impact on runtimes. This is quantified in more detail in panel (b) for a few different problem sizes $N$ and for all different tolerances listed in Table 1. Empirically, runtimes are seen to roughly increase as $T \propto p_{\epsilon}^2$ (black dashed lines), where $p_{\epsilon} = -\log_{10} \epsilon_{tol}$ is the nominal number of accuracy digits. The figure also confirms that the nominal tolerances $\epsilon_{tol}$ (vertical dotted lines) are reasonable (and in fact pessimist) estimates of the actual errors $\epsilon_v^{rel}$ (markers). To be clear, except for the $\epsilon_{tol} = 10^{-14}$ case, the plotted errors are computed via (5.2) by taking the $\epsilon_{tol} = 10^{-14}$ case as reference.

**6. Conclusions.** We have introduced an efficient technique for the numerical evaluation of the Biot-Savart law in periodic systems when the source field (e.g. current or vorticity) is supported on three-dimensional spatial curves. The approach is adapted from fast Ewald summation methods commonly used to accelerate particle simulations. The present work is relevant to applications in electromagnetism and fluid dynamics, where the source fields respectively induce magnetic or velocity fields around the curves. As an application, we have considered the vortex filament model (VFM) describing superfluid helium flows near absolute zero. In this context, the Fourier-space representation of the induced fields is a very attractive feature of the method, as it allows the efficient evaluation of physically relevant quantities such as energy spectra [11], coarse-grained velocity fields [8, 55] and velocity circulation [64, 69].

Several tunable parameters are introduced by the method, including an inverse “splitting” length scale $\alpha$ and two cut-off parameters $r_c$ and $k_{max}$ in physical and Fourier spaces. We have shown that the choice of the latter two can be reduced to that of a single non-dimensional coefficient $\beta$ setting the accuracy of the method, while $\alpha$ can be tuned to maximise performance. For a fixed accuracy level, choosing $\alpha \propto N^{1/3}$ leads to near-linear $O(N \log N)$ complexity of the method, where $N$ is the
number of vortex discretisation points. To our knowledge, this parameter selection procedure contrasts with the seemingly common approach in molecular dynamics which consists in first specifying reasonable cut-offs $r_c$ and $k_{max}$ in dimensional units, and then choosing a value of $\alpha$ which maximises accuracy given these cut-offs [30, 9, 3].

In terms of performance, our implementation is already capable of simulating turbulent systems with a reasonably large number of vortices using thread-based parallelism. Moreover, our NUFFT implementation (Appendix A) enables the use of GPUs to accelerate the computation of long-range interactions. The short-range part may also benefit in the near future from a GPU implementation of the cell lists algorithm for neighbour finding. Besides, these core algorithms are very amenable to distributed-memory parallelism (via MPI) using standard domain decomposition approaches, which could enable computations on multiple CPU or GPU nodes. Such strategies will be investigated in the future.

A second challenge concerns the timestep stability requirement $\Delta t \propto \ell_{min}^2$, which becomes very restrictive as the smallest resolved scale $\ell_{min}$ is decreased. Within this work, this issue is partially remediated by representing spatial curves using highly continuous splines, which allows to use relatively large values of $\ell_{min}$. Potential strategies for further improving on these limitations can include the use of splitting [18] or multirate [76] timestepping methods, in which the evolution equations would be split onto a fast term responsible for oscillations (e.g. the local self-induced velocity) and a slower term (e.g. all non-local contributions).

The proposed method can be readily generalised to cases where the vortex circulation $\kappa$ varies along the filaments. In particular, it may be used in the implementation of vortex filament methods for classical flows [25]. Another possible extension concerns the simpler two-dimensional (2D) case, in which vortex filaments are replaced by point vortices. Periodic point vortex systems have received much attention [56, 81, 52, 59], in part due to their link with statistical mechanics and with 2D turbulence. In that case, Ewald summation can be achieved using a function splitting similar to (3.3) [42, 24]. Finally, and more directly relevant to the motivation of this work, the proposed method may be applied to the study of *finite-temperature* liquid helium-4, in which the superfluid flow coexists and interacts with a normal viscous fluid. A very relevant model for describing this regime consists in coupling the VFM for the superfluid with the incompressible Navier-Stokes equations governing the normal fluid [48, 88, 35]. When the latter are solved using standard pseudo-spectral methods [22], the Fourier-space representation included in the present method can enable an accurate and efficient two-way transfer of interaction forces between the vortices and the normal fluid.

**Acknowledgments.** The author acknowledges support from the French Agence Nationale de la Recherche through the QuantumVIW project (Grant No. ANR-23-CE30-0024-04).

**Appendix A. NUFFT details and influence on Biot-Savart accuracy.**

Here we consider the error induced by the NUFFT tolerance on the long-range component of Biot-Savart computations. Briefly, a type-1 NUFFT can be summarised in three steps: (i) spreading values from a set of non-uniform locations onto a uniform grid using some smoothing kernel, (ii) performing a regular FFT on that grid, and (iii) undoing the spreading operation by deconvolution in Fourier space [37]. The type-2 NUFFT is the adjoint operation, roughly consisting of the same steps in opposite order. To reduce aliasing errors, the first two steps are done in a grid of $\sigma M$ points in each direction, where $M$ is the number of desired Fourier modes and $\sigma > 1$ is an oversampling factor. The accuracy of the NUFFT is then determined by (i) the choice
<figure><img src="image_7.png" /><figcaption>FIG. 6. Velocity and streamfunction errors associated to the NUFFT precision, for the test cases illustrated in Figure 2. Relative root-mean-square error of long-range filament velocity (+) and streamfunction values (×) at varying NUFFT spreading widths w. NUFFTs are computed using the backwards Kaiser-Bessel spreading kernel (A.1) with an oversampling factor σ = 1.5. Dashed lines represent an exponential decay ~ 10⁻⁰.⁹w. In all cases, the long-range cut-off parameter is set to β(f) = 8, ensuring subdominant cut-off errors. The Ewald splitting parameter is kept constant, α = 24/L.</figcaption></figure>

of spreading kernel, (ii) the width of the kernel (usually in number of grid points w), and (iii) the oversampling factor σ [31, 70, 15]. The computational complexity of a 3D NUFFT on N non-uniform points can be estimated to be O(w³N + (σM)³ log M).

To evaluate NUFFTs, here we use the (backwards) Kaiser-Bessel function as spreading kernel [47, 70]. For a support in x ∈ [-1, 1], the kernel and its Fourier transform are

$$ (A.1) \qquad \phi(x) = \frac{\sinh(\zeta\sqrt{1-x^2})}{\pi\sqrt{1-x^2}} \quad \leftrightarrow \quad \hat{\phi}(k) = I_0(\sqrt{\zeta^2 - k^2}), $$

where $I_0$ is the zero-th order modified Bessel function of the first kind. The shape parameter is chosen as $ζ = γwπ(1 - 1/2σ)$ to minimise NUFFT errors [70]. Here $γ = 0.995$ is an empirical “safety factor” [15] which slightly improves accuracy. For performance reasons, the kernel $φ(x)$ is evaluated via a fast and accurate piecewise polynomial approximation [15, 78] which takes advantage of SIMD vectorisation on modern CPUs. Our NUFFT implementation [67] is written in Julia and is available on https://github.com/jipolanco/NonuniformFFTs.jl (MIT License). Compared to other existing libraries available in Julia [15, 49], our implementation allows real-valued non-uniform data, in which case it takes advantage of real-to-complex FFT routines of the FFTW library [34]. Moreover, it uses the KernelAbstractions.jl [23] package to provide a platform-agnostic GPU implementation, enabling the acceleration of vortex simulations using different kinds of GPUs.

In Figure 6, we evaluate the effect of varying the NUFFT accuracy on the relative error of long-range Biot-Savart computations, for the cases illustrated in Figure 2. Concretely, we set the oversampling factor to $σ = 1.5$ and vary the spreading width $w$. For comparison, we compute a reference solution with $σ = 2$ and $w = 16$, which ensures that its accuracy is dominated by roundoff error. As seen in Figure 6, the relative error associated to the long-range streamfunction and velocity decay exponentially as ~ 10⁻⁰.⁹w when using an oversampling factor $σ = 1.5$. For $σ = 2$, the same errors decay as ~ 10⁻⁵w (not shown), but at the expense of costlier FFTs.
## Appendix B. Discretisation of spatial curves.

The evaluation of the Biot-Savart law (2.2) and its local regularisation (2.4) on vortex filaments require the estimation of local tangent and curvature vectors from a set of discrete points {$s_j$}$_{j=1}^N$ representing a spatial curve. This assumes that curves can be parametrised by some sufficiently regular function $s(\tau)$ whose derivatives on discrete locations $\tau_j$ (such that $s(\tau_j) = s_j$) are well defined. Furthermore, the use of quadratures to estimate line integrals (see (2.15)) requires being able to evaluate (interpolate) $s(\tau)$ and its derivatives at any arbitrary location $\tau$. To achieve these requirements, we choose to discretise curves using parametric quintic splines. First, note that in general the nodes $s_j$ are not equispaced, and that the choice of the parametrisation (the $\tau_j$ values) is arbitrary. We adopt the common choice of setting $\tau_j = \tau_{j-1} + |s_j - s_{j-1}|$ for $j \in \{2, \dots, N\}$, $\tau_1 = 0$, where $|\cdot|$ represents the Euclidean distance. In this case, $\tau$ is a rough approximation of the curve arc length $\xi$, tending to $\xi$ as the number of points increases.

We consider parametric splines of order $k$ [27]. Here by convention the polynomial order is $k-1$, so that $k=4$ and 6 respectively correspond to cubic and quintic splines. The basic idea is to parametrise a spatial curve in terms of a B-spline basis,

$$ (B.1) \qquad s(\tau) = \sum_{j=1}^{N} c_j b_j(\tau), \quad \text{for } \tau \in [0, T], $$

where $b_j(\tau)$ is a B-spline basis function. The basis functions are fully defined by the choice of spline order $k$ and spline knots $t_j$, which are a set of locations in $[0, T]$. Because the curves considered here are periodic, we can simply set $t_j = \tau_j$ for all $j$ (free-ended curves would require some extra care). An important and convenient property of the B-splines is their compact support. Namely, for even $k$, the B-spline $b_j(\tau)$ is zero outside of the interval $[t_{j-k/2}, t_{j+k/2}]$.

Spline interpolation consists in determining the coefficients $c_j$ given the equalities $s(\tau_i) = s_i$ for $i \in \{1, \dots, N\}$. This leads to a linear system of the form $A_{ij}c_j = s_i$ where $A_{ij} = b_j(\tau_i)$. Thanks to the compact support of the B-splines, $A$ is a cyclic banded matrix, which enables the use of specialised algorithms to efficiently solve the system. More precisely, with our choice of knots, $A$ has $k-1$ bands, as well as a few out-of-bands entries in the top right and bottom left corners due to periodicity. In particular, for cubic splines ($k=4$), $A$ is almost tridiagonal with two out-of-bands entries (one on each corner). A standard solution is to convert the system to a fully tridiagonal one using the Sherman-Morrison formula and then use the Thomas algorithm to solve the system in $O(N)$ time [86]. A similar procedure [57] can be applied to efficiently solve the cyclic pentadiagonal system arising from periodic quintic spline interpolation ($k=6$).

The derivative of a spline $s(\tau)$ of order $k$ is itself a spline of order $k-1$, whose coefficients $c'_j$ can be directly obtained from the coefficients $c_j$ using a local differentiation formula. Moreover, splines are efficiently evaluated at locations $\tau$ using de Boor's algorithm [27]. Therefore, knowing the coefficients $c_j$ one can evaluate curve positions and derivatives at any location $\tau$.

## Appendix C. Derivation of truncation error estimates.

### C.1. Short-range errors.
The absolute error associated to the truncation of the short-range velocity integral (3.7) at a single location $s_0$ is given by

$$ (C.1) \qquad v_{\text{err}}^{(n)}(s_0) = \frac{\kappa}{4\pi} \sum_{n \in \mathbb{Z}^3} \int_{r>r_c} g^{(n)}(r) \frac{r \times ds}{|r|^3} \quad \text{with } r = s - s_0 + nL, $$
where the integral is over all vortex elements at a distance $|r| = r > r_c$ from the point $s_0$. We first note that, for large $r$, the weight function $g^{(n)}(r)$ defined in (3.5) behaves as $g^{(n)}(r) = \frac{2\alpha r}{\sqrt{\pi}} e^{-(\alpha r)^2} [1 + O(r^{-5/2})]$. Secondly, we assume that the vortices are homogeneously distributed in the spatial domain, so that the total vortex length within a spherical shell of radii $[r, r+dr]$ with infinitesimal thickness $dr$ is the product between the mean vortex line density and the volume of the shell, namely $d\mathcal{L} = (\mathcal{L}/V)(4\pi r^2 dr)$. Here $\mathcal{L}$ is the total vortex length within a periodic cell of volume $V = L^3$. Thirdly, we assume that the orientations of the vortex elements $ds$ within the shell are random and independent of the separation vector $r$. Under these assumptions, we can formally estimate the root mean square (rms) short-range velocity error to be

$$ (C.2) \qquad \varepsilon_v^{(n)} \approx \frac{\kappa}{4\pi} \int_{r_c}^{\infty} \frac{2\alpha r}{\sqrt{\pi}} e^{-(\alpha r)^2} \frac{r}{r^3} \frac{d\mathcal{L}}{r^3} = \frac{2\kappa \mathcal{L}}{\sqrt{\pi} \alpha V} \int_{\beta^{(n)}}^{\infty} x e^{-x^2} dx, $$

where $\beta^{(n)} = \alpha r_c$. This integrates exactly to

$$ (C.3) \qquad \varepsilon_v^{(n)} \approx \frac{\kappa \mathcal{L}}{\sqrt{\pi} \alpha V} e^{-(\beta^{(n)})^2}. $$

The same procedure can be applied to the short-range streamfunction (3.6). Considering that $\text{erfc}(x) = \frac{e^{-x}}{x^{3/2}} [1 + O(x^{-2})]$ for large $x$, one obtains the associated rms error

$$ (C.4) \qquad \varepsilon_{\psi}^{(n)} \approx \frac{\kappa}{4\pi} \int_{r_c}^{\infty} \frac{e^{-(\alpha r)^2}}{(\alpha r)^{3/2}} \frac{d\mathcal{L}}{r} = \frac{\kappa \mathcal{L}}{\alpha^2 V} \int_{\beta^{(n)}}^{\infty} \frac{e^{-x^2}}{\sqrt{x}} dx = \frac{\kappa \mathcal{L}}{2\alpha^2 V} \Gamma\left(\frac{1}{4}, (\beta^{(n)})^2\right). $$

Here $\Gamma(s, x) = \int_x^\infty u^{s-1} e^{-u} du$ is the upper incomplete gamma function, which behaves asymptotically as $\Gamma(s, x) = x^{s-1} e^{-x} [1 + O(x^{-1})]$. Replacing this in (C.4) finally leads to the estimate

$$ (C.5) \qquad \varepsilon_{\psi}^{(n)} \approx \frac{\kappa \mathcal{L}}{2\alpha^2 V} \frac{e^{-(\beta^{(n)})^2}}{(\beta^{(n)})^{3/2}}. $$

In both cases, the truncation error is dominated by the exponential decay with the square of the non-dimensional cut-off parameter $\beta^{(n)}$.

**C.2. Long-range errors.** The rms error associated to the long-range velocity field is defined by

$$ (C.6) \qquad \left(\varepsilon_v^{(f)}\right)^2 = \frac{1}{V} \int_{\Omega} \left|\mathbf{v}^{(f)}(\mathbf{x}) - \tilde{\mathbf{v}}^{(f)}(\mathbf{x})\right|^2 d^3\mathbf{x}, $$

where $\mathbf{v}^{(f)}(\mathbf{x})$ is the Fourier-truncated velocity field defined in (3.10), and $\tilde{\mathbf{v}}^{(f)}(\mathbf{x})$ is the non-truncated velocity field – obtained by replacing $\mathbb{K}$ with $\mathbb{K}_\infty = \{\frac{2\pi m}{L}; m \in \mathbb{Z}\}$ in (3.10). Using Parseval’s theorem, the error can be equivalently written as

$$ (C.7) \qquad \left(\varepsilon_v^{(f)}\right)^2 = \sum_{\mathbf{k} \in (\mathbb{K}_\infty \setminus \mathbb{K})^3} |\hat{\mathbf{v}}_{\mathbf{k}}^{(f)}|^2 \le \sum_{|\mathbf{k}| > k_{\max}} |\hat{\mathbf{v}}_{\mathbf{k}}^{(f)}|^2, $$

where $\hat{\mathbf{v}}_{\mathbf{k}}^{(f)}$ are the Fourier coefficients of the long-range velocity. The last inequality results from converting from a “cubic” truncation ($k_i > k_{\max}$ for each Cartesian component $i$) to a spherical truncation $|\mathbf{k}| > k_{\max}$.
In physical terms, the velocity error (C.7) can be directly related to the kinetic energy spectrum $E(k)$, defined in a periodic domain as

$$ (C.8) \qquad E(k) = \frac{1}{2\Delta k} \sum_{|k| \in \mathcal{I}_k} |\hat{v}_k|^2 \quad \text{for } k = m\Delta k, m \in \mathbb{N}, $$

where $\hat{v}_k$ is a Fourier coefficient of the (unsmoothed) velocity field, $\Delta k = 2\pi/L$ is the distance between two successive discrete wavenumbers, and the sum is over all 3D wavenumbers $k$ within the spherical shell $\mathcal{I}_k \equiv [k - \frac{\Delta k}{2}, k + \frac{\Delta k}{2})$. Formally, the total kinetic energy (per unit mass) is then $E = \int_0^\infty E(k) dk$. In the present context, (C.8) will allow us to quantify the effect of truncating Fourier-space computations at a certain wavenumber $k_{\max}$.

At sufficiently small scales – smaller than the typical distance between vortices and the typical curvature radius of the vortices – vortex filaments can be considered as isolated straight lines. Under that assumption, one can show that the energy spectrum is $E(k) = \frac{\kappa^2 \mathcal{L}}{4\pi L^3} k^{-1}$ for sufficiently large $k$. To see this, one can consider the energy spectrum associated to a single straight vortex. Without loss of generality, the vortex is aligned with the $z$ direction and passes through the origin, so that its length within the cubic box is $\mathcal{L} = L$ and its induced velocity field is invariant in $z$. We can thus consider the simpler problem of a point vortex in two-dimensional space. By (3.9), the Fourier coefficients of its associated vorticity field are $\hat{\omega}_k = \frac{\kappa}{L^2} e^{-ik \cdot x} e_z = \frac{\kappa}{L^2} e_z$, where $x = (x, y) = (0, 0)$ is the vortex location, $k = (k_x, k_y)$, and $e_z$ is a unitary vector. Then, the associated velocity coefficients are $\hat{v}_k = ik \times \frac{\hat{\omega}_k}{k^2}$, whose squared magnitudes are $|\hat{v}_k|^2 = \frac{\kappa^2}{L^4 k^2}$ (since $\hat{\omega}_k$ is orthogonal to $k$). Besides, one can estimate the number of 2D Fourier modes intersecting the shell $\mathcal{I}_k$ to be on average $N_k = 2\pi m = kL$. Since $|\hat{v}_k|^2$ only depends on the amplitude of $k$, one can assume that all non-zero $\hat{v}_k$ within the shell $\mathcal{I}_k$ are approximately equal. Replacing in (C.8), we finally obtain the energy spectrum $E(k) \approx \frac{L}{4\pi} N_k |\hat{v}_k|^2 \approx \frac{\kappa^2 \mathcal{L}}{4\pi L^3} k^{-1} = \frac{\kappa^2 \mathcal{L}}{4\pi L^3} k^{-1}$.

This prediction is accurately verified in Figure 7 (blue curves) for the three different vortex configurations illustrated in Figure 2. Note that the relatively slow decay as $k^{-1}$ suggests that the integral of this spectrum (and thus the total energy) presents a logarithmic divergence. In fact, physically this is not the case, since the spectrum is regularised near the wavenumber $k_a = 2\pi/a$ associated to the vortex core size $a$, which is orders of magnitude smaller than the scales described by the VFM.

To provide an estimate of the error in the computation of the coarse-grained velocity $v^{(f)}$, we now consider its associated kinetic energy spectrum $E^{(f)}(k) = E(k) e^{-k^2/2\alpha^2} = \frac{\kappa^2 \mathcal{L}}{4\pi V} k^{-1} e^{-k^2/2\alpha^2}$ for large $k$. Note that the integral of $E^{(f)}(k)$ is the contribution of the long-range velocity field to the total kinetic energy $E$. Examples of $E^{(f)}(k)$ are shown in Figure 7 for two different values of the smoothing parameter $\alpha$. In the figure, the long-range truncation is at $k_{\max} = 64$ and the domain period is $L = 2\pi$. It is clear that, if $\alpha$ is too large (e.g. $\alpha L = 240$ in the figure, corresponding to $\beta^{(f)} = k_{\max}/2\alpha \approx 0.84$), then this level of truncation is insufficient to capture the whole long-range kinetic energy. A smaller value of $\alpha$ ($\alpha L = 24$ in the figure) means that the full coarse-grained energy can be captured when truncating at the same $k_{\max}$.

Assuming the truncation wavenumber $k_{\max}$ to be sufficiently large, we can now estimate an upper bound for the absolute truncation error (C.7) associated to the
<figure><img src="image_8.png" /><figcaption>FIG. 7. Kinetic energy spectra associated to three different test cases. The domain period is $L = 2\pi$ and the spectra are truncated at $k_{max} = 64$ (the number of Fourier modes in each direction is $M = 128$). Each panel shows the truncated energy spectrum associated to the unsmoothed velocity field ($\alpha L = \infty$), as well as the energy spectra associated to two Gaussian-smoothed velocity fields ($\alpha L = 240$ and 24). The spectra are normalised by the total kinetic energy $E_{total}$ obtained via (2.8). At large wavenumbers $k$, the spectra are well approximated by $E(k) = \frac{\kappa^2 \mathcal{L}}{4\pi V} k^{-1} e^{-k^2/2\alpha^2}$ (black dashed lines), where $\mathcal{L}$ is the total vortex length in each test case.</figcaption></figure>

long-range velocity,

$$ (C.9) \qquad \frac{(\varepsilon_v^{(f)})^2}{2} \approx \int_{k_{\max}}^{\infty} E^{(f)}(k) dk = \frac{\kappa^2 \mathcal{L}}{4\pi V} \int_{k_{\max}}^{\infty} \frac{e^{-k^2/2\alpha^2}}{k} dk = \frac{\kappa^2 \mathcal{L}}{8\pi V} E_1[2(\beta^{(f)})^2], $$

where $E_1(x) = \int_x^\infty \frac{e^{-u}}{u} du$ is the exponential integral function. For large $x$, $E_1(x) = \frac{e^{-x}}{x} [1 - \frac{1}{x} + O(x^{-2})]$, and therefore $(\varepsilon_v^{(f)})^2 \approx \frac{\kappa^2 \mathcal{L}}{8\pi V} e^{-2(\beta^{(f)})^2} / (\beta^{(f)})^2$. Finally, the estimated error associated to the long-range velocity field is

$$ (C.10) \qquad \varepsilon_v^{(f)} \approx \kappa \left(\frac{\mathcal{L}}{8\pi V}\right)^{1/2} \frac{e^{-(\beta^{(f)})^2}}{\beta^{(f)}}. $$

A similar estimation can be done for the long-range truncation error $\varepsilon_\psi^{(f)}$ associated to the streamfunction. This is obtained by integrating

$$ (C.11) \qquad \frac{(\varepsilon_\psi^{(f)})^2}{2} \approx \int_{k_{\max}}^{\infty} \frac{E^{(f)}(k)}{k^2} dk \approx \frac{\kappa^2}{4\alpha^2} \frac{\mathcal{L}}{4\pi V} \frac{e^{-2(\beta^{(f)})^2}}{[2(\beta^{(f)})^2]^2}, $$

where we have used $\int_x^\infty \frac{e^{-u}}{u^2} du = \frac{e^{-x}}{x} - E_1(x) = \frac{e^{-x}}{x^2} [1 + O(x^{-1})]$. This results in

$$ (C.12) \qquad \varepsilon_\psi^{(f)} \approx \frac{\kappa}{2\alpha} \left(\frac{\mathcal{L}}{8\pi V}\right)^{1/2} \frac{e^{-(\beta^{(f)})^2}}{(\beta^{(f)})^2}. $$

Similarly to short-range errors, the long-range truncation errors are dominated by the exponential decay with the square of the cut-off parameter $\beta^{(f)}$. We stress that the estimates (C.10) and (C.12) are expected to be valid for sufficiently large values of $k_{max}$ (or $\beta^{(f)}$). Furthermore, they can be expected to overestimate the actual errors due to the inequality in (C.7), which discards the contributions of some large-wavenumber modes which are actually resolved in the simulations (more precisely, a factor of $\sim 6/\pi$ Fourier modes is discarded by the above estimates).
### C.3. Combined truncation errors.
To avoid unnecessary computations and achieve optimal accuracy for a given computational cost, one would like the errors associated to the short- and long-range components to be equivalent. In light of the above estimations, it seems natural to set their respective non-dimensional cut-off parameters to the same value, $\beta^{(n)} = \beta^{(f)} = \beta$, which also helps reducing the size of the parameter space. This finally leads to the following total truncation error estimates:

$$ (C.13) \qquad \varepsilon_v \approx \kappa \left[ \frac{\mathcal{L}}{\sqrt{\pi \alpha V}} + \frac{1}{\beta} \left( \frac{\mathcal{L}}{8\pi V} \right)^{1/2} \right] e^{-\beta^2}, $$

$$ (C.14) \qquad \varepsilon_{\psi} \approx \frac{\kappa}{2\alpha} \left[ \frac{\mathcal{L}}{\alpha V} + \frac{1}{\beta^{1/2}} \left( \frac{\mathcal{L}}{8\pi V} \right)^{1/2} \right] \frac{e^{-\beta^2}}{\beta^{3/2}}. $$

In summary, the accuracy of the method is mainly controlled by a unique non-dimensional parameter $\beta$, from which both physical- and Fourier-space cut-offs $r_c$ and $k_{max}$ can be obtained. This still leaves the inverse splitting distance $\alpha$ as a free parameter that can be adjusted to optimise performance (subsection 5.3).

## REFERENCES

[1] D. M. AMBROSE, M. SIEGEL, AND S. TUPOVA, *A small-scale decomposition for 3D boundary integral computations with surface tension*, J. Comput. Phys., 247 (2013), pp. 168–191, https://doi.org/10.1016/j.jcp.2013.03.045.

[2] R. J. ARMS AND F. R. HAMA, *Localized-Induction Concept on a Curved Vortex and Motion of an Elliptic Vortex Ring*, Phys. Fluids, 8 (1965), pp. 553–559, https://doi.org/10.1063/1.1761268.

[3] A. ARNOLD, F. FAHRENBERGER, C. HOLM, O. LENZ, M. BOLTEN, H. DACHSEL, R. HALVER, I. KABADSHOW, F. GÄHLER, F. HEBER, J. ISERINGHAUSEN, M. HOFMANN, M. PIPPIG, D. POTTS, AND G. SUTMANN, *Comparison of scalable fast methods for long-range interactions*, Phys. Rev. E, 88 (2013), p. 063308, https://doi.org/10.1103/PhysRevE.88.063308.

[4] A. ARNOLD AND C. HOLM, *Efficient Methods to Compute Long-Range Interactions for Soft Matter Systems*, in Advanced Computer Simulation Approaches for Soft Matter Sciences II, C. Holm and K. Kremer, eds., Advances in Polymer Science, Springer, Berlin, Heidelberg, 2005, pp. 59–109, https://doi.org/10.1007/b136793.

[5] A. W. BAGGALEY AND C. F. BARENGHI, *Spectrum of turbulent Kelvin-waves cascade in superfluid helium*, Phys. Rev. B, 83 (2011), p. 134509, https://doi.org/10.1103/PhysRevB.83.134509.

[6] A. W. BAGGALEY AND C. F. BARENGHI, *Tree Method for Quantum Vortex Dynamics*, J. Low Temp. Phys., 166 (2012), pp. 3–20, https://doi.org/10.1007/s10909-011-0405-6.

[7] A. W. BAGGALEY AND J. LAURIE, *Kelvin-wave cascade in the vortex filament model*, Phys. Rev. B, 89 (2014), p. 014504, https://doi.org/10.1103/PhysRevB.89.014504.

[8] A. W. BAGGALEY, J. LAURIE, AND C. F. BARENGHI, *Vortex-Density Fluctuations, Energy Spectra, and Vortical Regions in Superfluid Turbulence*, Phys. Rev. Lett., 109 (2012), p. 205304, https://doi.org/10.1103/PhysRevLett.109.205304.

[9] V. BALLENEGGER, J. J. CERDA, O. LENZ, AND CH. HOLM, *The optimal P3M algorithm for computing electrostatic energies in periodic systems*, J. Chem. Phys., 128 (2008), p. 034109, https://doi.org/10.1063/1.2816570.

[10] V. BANICA, G. MAIERHOFER, AND K. SCHRATZ, *Numerical Integration of Schrödinger Maps via the Hasimoto Transform*, SIAM J. Numer. Anal., 62 (2024), pp. 322–352, https://doi.org/10.1137/22M1531555.

[11] C. F. BARENGHI, V. S. L'VOV, AND P.-E. ROCHE, *Experimental, numerical, and analytical velocity spectra in turbulent quantum fluid*, Proc. Natl. Acad. Sci. USA, 111 (2014), pp. 4683–4690, https://doi.org/10.1073/pnas.1312548111.

[12] C. F. BARENGHI, L. SKRBEK, AND K. R. SREENIVASAN, *Introduction to quantum turbulence*, Proc. Natl. Acad. Sci. USA, 111 (2014), pp. 4647–4652, https://doi.org/10.1073/pnas.1400033111.

[13] C. F. BARENGHI, L. SKRBEK, AND K. R. SREENIVASAN, *Quantum Turbulence*, Cambridge University Press, Cambridge, 2023, https://doi.org/10.1017/9781009345651.

[14] J. BARNES AND P. HUT, *A hierarchical O(N log N) force-calculation algorithm*, Nature, 324 (1986), pp. 446–449, https://doi.org/10.1038/324446a0.
[15] A. H. BARNETT, J. MAGLAND, AND L. AF KLINTEBERG, *A Parallel Nonuniform Fast Fourier Transform Library Based on an "Exponential of Semicircle" Kernel*, SIAM J. Sci. Comput., 41 (2019), pp. C479-C504, https://doi.org/10.1137/18M120885X.

[16] G. K. BATCHELOR, *Slender-body theory for particles of arbitrary cross-section in Stokes flow*, J. Fluid Mech., 44 (1970), pp. 419-440, https://doi.org/10.1017/S002211207000191X.

[17] J. BEZANSON, A. EDELMAN, S. KARPINSKI, AND V. B. SHAH, *Julia: A fresh approach to numerical computing*, SIAM Rev., 59 (2017), pp. 65-98, https://doi.org/10.1137/141000671.

[18] S. BLANES, F. CASAS, AND A. MURUA, *Splitting Methods for differential equations*, May 2024, https://arxiv.org/abs/2401.01722. To appear in Acta Numer. (2024).

[19] J. P. BOYD, *Chebyshev and Fourier Spectral Methods*, Dover Publications, Mineola, N.Y, second edition ed., Dec. 2001.

[20] T. F. BUTTKE, *A numerical study of superfluid turbulence in the self-induction approximation*, Journal of Computational Physics, 76 (1988), pp. 301-326, https://doi.org/10.1016/0021-9991(88)90145-3.

[21] A. J. CALLEGARI AND L. TING, *Motion of a Curved Vortex Filament with Decaying Vortical Core and Axial Velocity*, SIAM J. Appl. Math., 35 (1978), pp. 148-175, https://doi.org/10.1137/0135013.

[22] C. CANUTO, M. Y. HUSSAINI, A. QUARTERONI, AND T. A. ZANG, *Spectral Methods in Fluid Dynamics*, Springer Berlin Heidelberg, Berlin, Heidelberg, 1988, https://doi.org/10.1007/978-3-642-84108-8.

[23] V. CHURAVY, *KernelAbstractions.jl*, Zenodo, Dec. 2024, https://doi.org/10.5281/zenodo.4021259.

[24] B. CICHOCKI AND B. U. FELDERHOF, *Electrostatic interactions in two-dimensional Coulomb systems with periodic boundary conditions*, Physica A, 158 (1989), pp. 706-722, https://doi.org/10.1016/0378-4371(89)90487-1.

[25] G.-H. COTTET AND P. D. KOUMOUTSAKOS, *Vortex Methods: Theory and Practice*, Cambridge University Press, Cambridge, 2000, https://doi.org/10.1017/CBO9780511526442.

[26] T. DARDEN, D. YORK, AND L. PEDERSEN, *Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems*, J. Chem. Phys., 98 (1993), pp. 10089-10092, https://doi.org/10.1063/1.464397.

[27] C. DE BOOR, *A Practical Guide to Splines*, Applied Mathematical Sciences, Springer-Verlag, New York, 1978.

[28] F. DE LA HOZ, C. J. GARCÍA-CERVERA, AND L. VEGA, *A Numerical Study of the Self-Similar Solutions of the Schrödinger Map*, SIAM J. Appl. Math., 70 (2009), pp. 1047-1077, https://arxiv.org/abs/27862547.

[29] M. DESERNO AND C. HOLM, *How to mesh up Ewald sums. I. A theoretical and numerical comparison of various particle mesh routines*, J. Chem. Phys., 109 (1998), pp. 7678-7693, https://doi.org/10.1063/1.477414.

[30] M. DESERNO AND C. HOLM, *How to mesh up Ewald sums. II. An accurate error estimate for the particle-particle-particle-mesh algorithm*, J. Chem. Phys., 109 (1998), pp. 7694-7701, https://doi.org/10.1063/1.477415.

[31] A. DUTT AND V. ROKHLIN, *Fast Fourier Transforms for Nonequispaced Data*, SIAM J. Sci. Comput., 14 (1993), pp. 1368-1393, https://doi.org/10.1137/0914081.

[32] P. P. EWALD, *Die Berechnung optischer und elektrostatischer Gitterpotentiale*, Ann. Phys., 369 (1921), pp. 253-287, https://doi.org/10.1002/andp.19213690304.

[33] R. P. FEYNMAN, R. B. LEIGHTON, AND M. L. SANDS, *The Feynman Lectures on Physics. Volume 2: Mainly Electromagnetism and Matter*, Basic Books, New York, 2011.

[34] M. FRIGO AND S. G. JOHNSON, *The design and implementation of FFTW*, Proc. IEEE, 93 (2005), pp. 216-231, https://doi.org/10.1109/JPROC.2004.840301.

[35] L. GALANTUCCI, A. W. BAGGALEY, C. F. BARENGHI, AND G. KRSTULOVIC, *A new self-consistent approach of quantum turbulence in superfluid helium*, Eur. Phys. J. Plus, 135 (2020), p. 547, https://doi.org/10.1140/epjp/s13360-020-00543-0.

[36] M. GAZZOLA, B. HEJAZIALHOSSEINI, AND P. KOUMOUTSAKOS, *Reinforcement Learning and Wavelet Adapted Vortex Methods for Simulations of Self-propelled Swimmers*, SIAM J. Sci. Comput., 36 (2014), pp. B622-B639, https://doi.org/10.1137/130943078.

[37] L. GREENGARD AND J.-Y. LEE, *Accelerating the Nonuniform Fast Fourier Transform*, SIAM Rev., 46 (2004), pp. 443-454, https://doi.org/10.1137/S003614450343200X.

[38] L. GREENGARD AND V. ROKHLIN, *A fast algorithm for particle simulations*, J. Comput. Phys., 73 (1987), pp. 325-348, https://doi.org/10.1016/0021-9991(87)90140-9.

[39] F. GROTTO AND S. MORLACCHI, *Decay of time correlations in point vortex systems*, Physica D, (2024), p. 134169, https://doi.org/10.1016/j.physd.2024.134169.

[40] R. HÄNNINEN, *Dissipation enhancement from a single vortex reconnection in superfluid helium*,
Phys. Rev. B, 88 (2013), p. 054511, https://doi.org/10.1103/PhysRevB.88.054511.

[41] R. HÄNNINEN AND A. W. BAGGALEY, *Vortex filament method as a tool for computational visualization of quantum turbulence*, Proc. Natl. Acad. Sci. USA, 111 (2014), pp. 4667-4674, https://doi.org/10.1073/pnas.1312535111.

[42] H. HASIMOTO, *On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres*, J. Fluid Mech., 5 (1959), pp. 317-328, https://doi.org/10.1017/S0022112059000222.

[43] H. HASIMOTO, *A soliton on a vortex filament*, J. Fluid Mech., 51 (1972), pp. 477-485, https://doi.org/10.1017/S0022112072002307.

[44] F. HEDMAN AND A. LAAKSONEN, *Ewald summation based on nonuniform fast Fourier transform*, Chem. Phys. Lett., 425 (2006), pp. 142-147, https://doi.org/10.1016/j.cplett.2006.04.106.

[45] H. HELMHOLTZ, *Über Integrale der hydrodynamischen Gleichungen, welche den Wirbelbewegungen entsprechen.*, J. Reine Angew. Math., 55 (1858), pp. 25-55.

[46] R. HOCKNEY AND J. EASTWOOD, *Computer Simulation Using Particles*, CRC Press, 1988, https://doi.org/10.1201/9780367806934.

[47] J. JACKSON, C. MEYER, D. NISHIMURA, AND A. MACOVSKI, *Selection of a convolution function for Fourier inversion using gridding (computerised tomography application)*, IEEE Trans. Med. Imaging, 10 (Sept./1991), pp. 473-478, https://doi.org/10.1109/42.97598.

[48] D. KIVOTIDES, C. F. BARENGHI, AND D. C. SAMUELS, *Triple Vortex Ring Structure in Superfluid Helium II*, Science, 290 (2000), pp. 777-779, https://doi.org/10.1126/science.290.5492.777.

[49] T. KNOPP, M. BOBERG, AND M. GROSSER, *NFFT.jl: Generic and Fast Julia Implementation of the Nonequidistant Fast Fourier Transform*, SIAM J. Sci. Comput., 45 (2023), pp. C179-C205, https://doi.org/10.1137/22M1510935.

[50] J. KOLAFA AND J. W. PERRAM, *Cutoff Errors in the Ewald Summation Formulae for Point Charge Systems*, Mol. Simul., 9 (1992), pp. 351-368, https://doi.org/10.1080/08927029208049126.

[51] P. KOUMOUTSAKOS, *Multiscale Flow Simulations Using Particles*, Annu. Rev. Fluid Mech., 37 (2005), pp. 457-487, https://doi.org/10.1146/annurev.fluid.37.061903.175753.

[52] V. S. KRISHNAMURTHY AND T. SAKAJO, *The N-vortex problem in a doubly periodic rectangular domain with constant background vorticity*, Physica D, 448 (2023), p. 133728, https://doi.org/10.1016/j.physd.2023.133728.

[53] H. LAMB, *Hydrodynamics*, Dover Publications, 6 ed., 1945, https://archive.org/details/hydrodynamics00lamb.

[54] E. LAUGA AND T. R. POWERS, *The hydrodynamics of swimming microorganisms*, Rep. Prog. Phys., 72 (2009), p. 096601, https://doi.org/10.1088/0034-4885/72/9/096601.

[55] J. LAURIE AND A. W. BAGGALEY, *Vorticity locking and pressure dynamics in finite-temperature superfluid turbulence*, Phys. Rev. Fluids, 8 (2023), p. 054604, https://doi.org/10.1103/PhysRevFluids.8.054604.

[56] D. LINDBÖ AND A.-K. TORNBERG, *Spectrally accurate fast summation for periodic Stokes potentials*, J. Comput. Phys., 229 (2010), pp. 8994-9010, https://doi.org/10.1016/j.jcp.2010.08.026.

[57] X.-G. LV AND J. LE, *A note on solving nearly penta-diagonal linear systems*, Appl. Math. Comput., 204 (2008), pp. 707-712, https://doi.org/10.1016/j.amc.2008.07.012.

[58] W. MATTSON AND B. M. RICE, *Near-neighbor calculations using a modified cell-linked list method*, Comput. Phys. Commun., 119 (1999), pp. 135-148, https://doi.org/10.1016/S0010-4655(98)00203-3.

[59] O. MAXIAN, A. MOGILNER, AND A. DONEV, *Integral-based spectral method for inextensible slender fibers in Stokes flow*, Phys. Rev. Fluids, 6 (2021), p. 014102, https://doi.org/10.1103/PhysRevFluids.6.014102.

[60] R. I. MCLACHLAN AND G. R. W. QUISPEL, *Splitting methods*, Acta Numer., 11 (2002), pp. 341-434, https://doi.org/10.1017/S0962492902000053.

[61] J. J. MONAGHAN AND R. J. HUMBLE, *Vortex Particle Methods for Periodic Channel Flow*, J. Comput. Phys., 107 (1993), pp. 152-159, https://doi.org/10.1006/jcph.1993.1132.

[62] D. W. MOORE AND P. G. SAFFMAN, *The motion of a vortex filament with axial flow*, Phil. Trans. R. Soc. Lond. A, 272 (1972), pp. 403-429, https://doi.org/10.1098/rsta.1972.0055.

[63] Y. MORI, L. OHM, AND D. SPIRN, *Theoretical Justification and Error Analysis for Slender Body Theory*, Commun. Pure Appl. Math., 73 (2020), pp. 1245-1314, https://doi.org/10.1002/cpa.21872.

[64] N. P. MÜLLER, J. I. POLANCO, AND G. KRSTULOVIC, *Intermittency of Velocity Circulation in Quantum Turbulence*, Phys. Rev. X, 11 (2021), p. 011053, https://doi.org/10.1103/PhysRevX.11.011053.

[65] P. K. NEWTON AND G. CHAMOUN, *Vortex Lattice Theory: A Particle Interaction Perspective*,
SIAM Rev., 51 (2009), pp. 501–542, https://doi.org/10.1137/07068597x.

[66] M. PIPPIG AND D. POTTS, *Parallel Three-Dimensional Nonequispaced Fast Fourier Transforms and Their Application to Particle Simulation*, SIAM J. Sci. Comput., 35 (2013), pp. C411–C437, https://doi.org/10.1137/120888478.

[67] J. I. POLANCO, *Nonuniform FFTs.jl*. Zenodo, Nov. 2024, https://doi.org/10.5281/zenodo.14637607.

[68] J. I. POLANCO, *VortexPasta.jl*. Zenodo, Jan. 2025, https://doi.org/10.5281/zenodo.14749388.

[69] J. I. POLANCO, N. P. MÜLLER, AND G. KRSTULOVIC, *Vortex clustering, polarisation and circulation intermittency in classical and quantum turbulence*, Nat. Commun., 12 (2021), p. 7090, https://doi.org/10.1038/s41467-021-27382-6.

[70] D. POTTS AND G. STEIDL, *Fast Summation at Nonequispaced Knots by NFFT s*, SIAM J. Sci. Comput., 24 (2003), pp. 2013–2037, https://doi.org/10.1137/S1064827502400984.

[71] P. H. ROBERTS AND R. J. DONNELLY, *Dynamics of vortex rings*, Phys. Lett. A, 31 (1970), pp. 137–138, https://doi.org/10.1016/0375-9601(70)90193-3.

[72] P. G. SAFFMAN, *Vortex Dynamics*, Cambridge University Press, Jan. 1993, https://doi.org/10.1017/cbo9780511624063.

[73] D. SAINTILLAN, E. DARVE, AND E. S. G. SHAQFEH, *A smooth particle-mesh Ewald algorithm for Stokes suspension simulations: The sedimentation of fibers*, Phys. Fluids, 17 (2005), p. 033301, https://doi.org/10.1063/1.1862262.

[74] D. C. SAMUELS, *Velocity matching and Poiseuille pipe flow of superfluid helium*, Phys. Rev. B, 46 (1992), pp. 11714–11724, https://doi.org/10.1103/PhysRevB.46.11714.

[75] D. C. SAMUELS, *Vortex Filament Methods for Superfluids*, in Quantized Vortex Dynamics and Superfluid Turbulence, C. F. Barenghi, R. J. Donnelly, and W. F. Vinen, eds., Lecture Notes in Physics, Springer, Berlin, Heidelberg, 2001, pp. 97–113, https://doi.org/10.1007/3-540-45542-6-9.

[76] A. SANDU, *A Class of Multirate Infinitesimal GARK Methods*, SIAM J. Numer. Anal., 57 (2019), pp. 2300–2327, https://doi.org/10.1137/18m1205492.

[77] K. W. SCHWARZ, *Three-dimensional vortex dynamics in superfluid $^4$He: Line-line and line-boundary interactions*, Phys. Rev. B, 31 (1985), pp. 5782–5804, https://doi.org/10.1103/PhysRevB.31.5782.

[78] D. S. SHAMSHIRGAR, J. BAGGE, AND A.-K. TORNBERG, *Fast Ewald summation for electrostatic potentials with arbitrary periodicity*, J. Chem. Phys., 154 (2021), p. 164109, https://doi.org/10.1063/5.0044895.

[79] K. SHARIFF AND A. LEONARD, *Vortex Rings*, Annu. Rev. Fluid Mech., 24 (1992), pp. 235–279, https://doi.org/10.1146/annurev.fl.24.010192.001315.

[80] I. S. SULLIVAN, J. J. NIEMELA, R. E. HERSHBERGER, D. BOLSTER, AND R. J. DONNELLY, *Dynamics of thin vortex rings*, J. Fluid Mech., 609 (2008), pp. 319–347, https://doi.org/10.1017/S0022112008002292.

[81] W. THOMSON, *3. Vibrations of a Columnar Vortex*, Proc. R. Soc. Edinb., 10 (1880/ed), pp. 443–456, https://doi.org/10.1017/S0370164600044151.

[82] A. VAN KAN, A. ALEXAKIS, AND M.-E. BRACHET, *Intermittency of three-dimensional perturbations in a point-vortex model*, Phys. Rev. E, 103 (2021), p. 053102, https://doi.org/10.1103/PhysRevE.103.053102.

[83] L. J. VERMEER, J. N. SØRENSEN, AND A. CRESPO, *Wind turbine wake aerodynamics*, Prog. Aerosp. Sci., 39 (2003), pp. 467–510, https://doi.org/10.1016/S0376-0421(03)00078-2.

[84] D. H. WACKS, A. W. BAGGALEY, AND C. F. BARENghi, *Coherent laminar and turbulent motion of toroidal vortex bundles*, Phys. Fluids, 26 (2014), p. 027102, https://doi.org/10.1063/1.4864659.

[85] J. B. WEISS AND J. C. MCWILLIAMS, *Nonergodicity of point vortices*, Phys. Fluids, 3 (1991), pp. 835–844, https://doi.org/10.1063/1.858014.

[86] M. YARROW, *Solving periodic block tridiagonal systems using the Sherman-Morrison-Woodbury formula*, in 9th Comput. Fluid Dyn. Conf., Buffalo, NY, U.S.A., June 1989, American Institute of Aeronautics and Astronautics, https://doi.org/10.2514/6.1989-1946.

[87] S. YUI, H. KOBAYASHI, M. TSUBOTA, AND R. YOKOTA, *Quantum turbulence coupled with externally driven normal-fluid turbulence in superfluid $^4$He*, May 2021, https://doi.org/10.48550/arXiv.2105.09499, https://arxiv.org/abs/2105.09499.

[88] S. YUI, M. TSUBOTA, AND H. KOBAYASHI, *Three-Dimensional Coupled Dynamics of the Two-Fluid Model in Superfluid $^4$He: Deformed Velocity Profile of Normal Fluid in Thermal Counterflow*, Phys. Rev. Lett., 120 (2018), p. 155301, https://doi.org/10.1103/PhysRevLett.120.155301.