Weaving classical turbulence with quantum skeleton

Weiyu Shen State Key Laboratory for Turbulence and Complex Systems, College of Engineering, Peking University, Beijing 100871, PR China Jie Yao Advanced Research Institute of Multidisciplinary Sciences, Beijing Institute of Technology, Beijing 100081, PR China Yue Yang [email protected] State Key Laboratory for Turbulence and Complex Systems, College of Engineering, Peking University, Beijing 100871, PR China HEDPS-CAPT, Peking University, Beijing 100871, PR China

Abstract

Matter entanglement is a common chaotic structure in both quantum and classical systems. Turbulence can be pictured as a tangle of vortex filaments in superfluids and viscous vortices in classical fluids. However, it is hard to explain how the statistical properties of turbulence arise from elemental structures. Here we use the quantum vortex tangle as a skeleton to generate an instantaneous classical turbulent field with intertwined vortex tubes. Combining the quantum skeleton and tunable vortex thickness makes the synthetic turbulence satisfy key statistical laws and provides valuable insights for elucidating energy cascade and extreme events. By manipulating the elemental structures, we customize turbulence with desired statistical features. This bottom-up approach of “weaving” turbulence provides a testbed for analyzing and modeling turbulence.

Turbulence is recognized as one of the most complex physical systems [1]. Delving into the underlying structure is the foundation for understanding and regulating complex systems. Although bottom-up approaches have achieved significant success in physics [2], chemistry [3], and biology [4, 5], it remains very challenging to decompose a turbulent field into elementary structures and explain the statistical properties of turbulence that arise from the superposition of simple flow fields.

Vortices are the key components of turbulence [6], constantly accommodating and releasing stress within the fluid, and shaping its behavior, akin to the sinews and muscles of a living organism [7, 8]. They are not randomly distributed in turbulence, but rather form coherent patterns [9], such as hairpin vortices [10] near walls, reflecting some underlying order beneath the seemingly random fluid motion. Various structural-based vortex models, with spherical [11], tubular [12], sheet [13], or spiral [14] shapes, have been suggested to represent fine-scale structures of turbulence [15]. On the other hand, massive observations of turbulence [16, 17, 18, 19] reveal that the intense vortices are not a simple accumulation of elemental structures; the chaotic vortices, they always entangle with each other, forming a complex network [20]. Therefore, identifying simple universal structures that can capture essential features and dynamics of turbulence is a formidable challenge.

An inverse approach to understanding turbulence assembles a turbulent field from the bottom-up – using fine structures as building blocks. By comparing synthetic and real turbulence, we can distinguish the effects of different elemental structures or their combinations on turbulence statistics. Existing synthetic methods of turbulence, such as Fourier modes [21], linear-eddy models [22], fractal models [23], and multiscale Lagrangian map [24] can generate random velocity fields with some coherence. Still, they struggle to simultaneously reproduce fine-scale vortices, turbulence statistics, and intermittency in 3D turbulence (table S1).

To shed light on the elemental structure of turbulence, we develop a novel bottom-up method to “weave” an instantaneous turbulent field using intertwined vortices with customizable fine-scale features. This provides a unique way to generate turbulent fields for assessing turbulence theory and developing structure-based models. It eliminates the need for simulating flow evolution and can be employed as tailored flow data for further turbulence simulations and machine learning training.

Skeleton of turbulence. Considering the vortices should be chaotic and entangled, we propose that quantum turbulence [25] can serve as the “skeleton” for classical turbulence. Superfluid turbulence, governed by quantum physics, manifests as a tangle of vortex filaments, exhibiting remarkable statistical similarities with classical turbulence [26, 27]. We simulated the vortex tangle of superfluid helium II using the vortex filament method [28] in a cube of size $\mathcal{L}=0.1$ cm with periodic boundary conditions (see movie S1 and Supplementary materials, Materials and Methods A). After a rapid initial transition, the length of superfluid vortex filaments reaches a statistically stationary, homogeneous, and isotropic state (Figs. 1A and B). In this stage, the frequency spectrum of the vortex-line density fluctuation obeys the $-5/3$ scaling, consistent with low-temperature experiments [29] (inset in Fig. 1B). These quantum vortex filaments are used as centerlines, i.e., the skeleton, of intertwined viscous vortex tubes in weaving turbulence.

Subsequently, we generated the axisymmetric Burgers vortex tubes from the skeleton. Here, the radial vorticity distribution is $\omega_{s}(\rho)=\Gamma\exp\left[-\rho^{2}/(2\sigma^{2})\right]/(2\pi\sigma^{2})$ with the radial distance $\rho$ from the centerline and the normalized circulation $\Gamma=1$ . This simple exact solution of the Navier–Stokes equations is an attractive candidate for modeling fine-scale turbulence and encapsulates the stretching by the local strain and the dissipation by viscosity. In particular, the vortex tube has a variable core size $\sigma(s)$ along arc-length $s$ of the centerline, which characterizes the evolution of vortices with multiscale features.

We have developed a general method to construct the vorticity for vortex tubes with arbitrary centerline topology, differential twist, and variable thickness (see Supplementary materials, Materials and Methods B, and Supplementary Text, Section II for details). Figure 1C shows an example of a constructed 3D turbulent flow field with the quantum skeleton (in Fig. 1A). This “woven turbulence”, a special synthetic turbulent field, consists of interwind, worm-like vortices. Statistically, it is homogeneous and isotropic, as evidenced by the Gaussian probability density functions (PDFs) of fluctuating velocity components (Fig. S5).

Flow statistics and energy cascade. Turbulence statistics are estimated from the instantaneous woven field based on the celebrated Kolmogorov 1941 theory [30], with measuring the turbulent kinetic energy $k_{t}=3.12$ and enstrophy $\left\langle\Omega\right\rangle=2.72\times 10^{3}$ (See Supplementary Text, Section III A for details). We calculated the key statistics, the mean dissipation rate $\left\langle\epsilon\right\rangle=3.67$ , the kinetic viscosity $\nu=6.67\times 10^{-4}$ , and the Taylor Reynolds number $Re_{\lambda}=161.4$ (table S2).

The various structural features introduced in our bottom-up construction lead to a multi-scale nature of the woven turbulent flow field. The vortex tubes are intertwined in chaotic patterns, with varying entanglement degrees and core sizes (Fig. 2A). These features affect the range of scales in turbulence. All critical length scales are calculated and annotated in Fig. 2B. The quantum skeleton encodes large-scale information including the integral length scale $L=k_{t}^{3/2}/\left\langle\varepsilon\right\rangle$ and average distance $l$ between vortices, which characterize the denseness of vortex tubes. The viscous vortices generated on the skeleton expand the range of scales from the zero measure of quantum filaments to the finite Kolmogorov scale $\eta=\left(\nu^{3}/\varepsilon\right)^{1/4}$ , which characterizes the viscous dissipation of kinetic energy. By this means, the integration of the quantum skeleton and vortex tubes endows the woven turbulence with a wide range of length scales as in classical turbulence.

To explore the energy cascade [31] in the woven turbulence, figure 2C compares the rescaled 3D energy spectra of classical, woven, and quantum turbulence. The energy spectrum of woven turbulence agrees well with the model result of classical turbulence [32] and direct numerical simulation (DNS) result and satisfies the scaling of $k^{-5/3}$ in the inertial range. The compensated energy spectrum in the inset of Fig. 2C reveals the bottleneck effect [33] between the inertial and dissipation ranges, which demonstrates that woven turbulence captures the cascade of turbulence at small scales. By contrast, the spectrum of quantum turbulence only shows the scaling of $k^{-5/3}$ , whereas it is distinct from that of classical and woven turbulence at small scales due to the lack of finite-sized tubular geometry.

The cascade picture elucidates why woven turbulence can recover classical turbulence statistics from tangled filaments by introducing viscous vortex tubes. Figure 2D sketches the bifurcation of energy cascades in woven and quantum turbulence. In quantum turbulence, the reconnection of vortex filaments generates high-frequency Kelvin waves, and the energy is dissipated into heat by phonon radiation at the atomic scale [34]. In woven turbulence, the quantum Kelvin-wave cascade is transformed into the classical Richardson cascade by supplementing the viscous vortex scales below the vortex filament spacing $l$ , so that the energy in the inertial range can be transferred and then dissipated by viscosity at the Kolmogorov scale. This cascade is supported by the spectral energy flux and transfer function in Fig. S7.

Longitudinal velocity structure functions further confirm the consistent statistics in woven and classical turbulence. In Fig. 2E, the second-order structure function $D_{LL}(r)$ , rescaled by Kolmogorov’s scaling of $r^{-2/3}$ in the inertial range, is consistent with that of classical homogeneous isotropic turbulence. The woven turbulence also displays the negative third-order structure function $D_{LLL}(r)$ with the scaling of $r^{-1}$ . This unique characteristics of classical turbulence is absent in Gaussian random fields and quantum turbulence. Moreover, the higher-order structure functions of woven turbulence have the anomalous scalings in classical turbulence (Fig. S6).

Intermittency and extreme events. The coexistence of ordered and random natures causes classical turbulence to exhibit intermittency and extreme events, which are highly challenging to present in synthetic turbulence. In woven turbulence, distributions of the local enstrophy $\Omega=|\bm{\omega}|^{2}$ and dissipation (divided by viscosity) $\Sigma=\varepsilon/\nu=2S_{ij}S_{ij}$ , measuring the strength of the velocity gradient tensor $S_{ij}=\partial u_{i}/\partial x_{j}$ [35], are highly intermittent (Figs. 3A and B) as in classical turbulence. In isotropic turbulence, their volume averages $\left\langle\Omega\right\rangle=\left\langle\Sigma\right\rangle$ are the same, but the distributions of their local structures are different. Isosurfaces of $\Omega$ and $\Sigma$ with the same threshold occupy different parts of the vortex tube (Fig. 3C). Dissipation isosurface partly envelops the enstrophy one. Our vortex model with varying core sizes explains that the mismatch between the enstrophy and dissipation is due to vortex stretching. In Fig. 3D, these two quantities are co-located in a vortex tube with a constant core size, while variable core sizes lead to more intense dissipation being concentrated where the vortex is stretched.

The flow structure of extreme events can be classified and identified through $Q(\bm{x})=-\left(S_{ij}S_{ji}\right)/2$ and $R(\bm{x})=-\operatorname{det}\left[S_{ij}\right]$ , the second and third invariants of $S_{ij}$ [36]. The joint PDF of $R$ and $Q$ shows the distribution of the local topology of extreme events (Fig. 3E). The $R$ - $Q$ plane is divided into four regions (vortex stretching, vortex compressing, axial strain, and biaxial strain) by $R=0$ and the Vieillefosse line $27R^{2}+4Q^{3}=0$ [37]. The extreme events are dominated by vortex stretching and compressing. However, unlike the classical teardrop shape, the $R$ - $Q$ PDF in woven turbulence has a butterfly shape that lacks the fine-scale structural transition from vortices (upper part) to strains (lower part). It is mainly attributed to our vortex model, which is compact in a tube region. In contrast, real turbulence diffuses vorticity throughout the entire space, so that extreme events involve low-vorticity regions near $Q=0$ .

Another notable difference is the symmetrical and asymmetrical $R$ - $Q$ PDFs in woven and real turbulence, respectively. We performed a DNS for decaying turbulence with the initial woven field, as described in Section IV of the Supplementary Text. When the vortex tubes start to interact with each other, well before the onset of fully developed turbulence, a part of the axial strains are rapidly converted into biaxial strains, resulting in an asymmetric $R$ - $Q$ shape (Figs. S9 and S10). Subsequently, small-scale vortices produced by the vortex interaction fill the gap between the vortex tubes, and the $R$ - $Q$ shape evolves into a typical teardrop shape, implying that the local topology highly depends on the small-scale flow structures rather than the large-scale vortex skeleton.

Customizing turbulence. Besides unveiling insights into turbulent structures, the method of weaving turbulence also enables customizing a turbulent field with desired properties, such as the Reynolds number, turbulent kinetic energy, and helicity [38, 39], by precisely controlling the geometry and topology of elemental vortical structures (see Fig. 4A for example). In particular, this method can generate turbulent fields that are valuable in theoretical study but difficult or costly to obtain through experiments and numerical simulations.

The quantum skeletons with more vortex filaments and higher entanglement correspond to larger turbulent kinetic energy and Reynolds number. A more valuable question is how the vortex geometry affects turbulence. We manipulated the core size and internal twist [40] of vortex tubes (Fig. 4A) with the same quantum skeleton. Increasing core size $\sigma$ can reduce the width of the inertial range of woven turbulence, and increase the Reynolds number with a power-law scaling of $Re_{\lambda}\sim\sigma^{-3/10}$ (Fig. 4D). In general, the scale $\sigma$ of vortex tubes is smaller than the scale $l$ of the quantum skeleton (see Fig. 2B). Increasing the core size of coherent vortices brings the scales of viscous vortices closer to large scales of the skeleton, leading to an overall contraction of length scales.

The internal twist of vortex lines within vortex tubes is a key topological component of helicity, altering the flow chirality [38]. The local twist rate [40] $\eta_{t}$ of the vortex lines along a periodic vortex centerline $\mathcal{C}$ contributes to the twisting component $T_{w}=\oint_{\mathcal{C}}\eta_{t}/(2\pi)\mathrm{~{}d}s$ in the topological decomposition of total helicity [41, 42]

	$\displaystyle H$	$\displaystyle\equiv\iiint\bm{u}\cdot\bm{\omega}\,\mathrm{d}\mathcal{V}$		(1)
		$\displaystyle=\sum_{i\neq j}\Gamma_{i}\Gamma_{j}L_{k,ij}+\sum_{i}\Gamma_{i}^{2}(W_{r,i}+T_{w,i}),$		(1)

where each index labels a single vortex tube, and the writhe $W_{r}$ and link $L_{k}$ are only related to the knotting and linking of the vortex skeleton. In Fig. 4B, we constructed a helical turbulent field with highly coiled vortex lines as right-handed spirals inside the vortex tube. Figure 4E shows the chirality from the helicity spectrum. Right-handed twist waves attenuate the left-handed polarization of the helicity spectrum at a given wavenumber, resulting in a predominant right-handed chirality. In contrast, woven turbulence without internal twist has a balanced helical wave decomposition with $H\approx 0$ .

Modifying the cross-sectional shape of elemental vortices is effective in customizing turbulence statistics. As an example, we use the vortex sheet, another candidate of the elemental structure, to weave a turbulent field in Fig. 4C. The resulting energy spectrum exhibits the scaling law of $k^{-7/3}$ (Fig. S12A), which is steeper than the $k^{-5/3}$ contributed by vortex tubes. This difference is consistent with the theoretical scaling of an ensemble of elemental vortices, i.e., $E\sim k^{-2}$ for sheet-like structures [12] is steeper than $E\sim k^{-1}$ for tube-like structures. Moreover, the $R$ - $Q$ joint PDF in the turbulence woven by vortex sheets has higher values near $Q=0$ than that woven by vortex tubes (Fig. 4F), because vortex sheets cover larger spatial extent than vortex tubes, resulting in increased regions with low vorticity and high strain rate in turbulence.

Discussion. In summary, our bottom-up approach of “weaving” turbulence, utilizing adjustable fine structures as building blocks, offers a testbed for exploring and understanding elemental structures of turbulence. By utilizing quantum vortex filaments as “skeleton”, we have generated classical turbulent fields consisting of intertwined vortex tubes. The combination of the quantum skeleton and tunable tube thickness ensures that the woven turbulence adheres to key statistical laws, including the -5/3 scaling of the energy spectrum, negative third-order structure function, and positive spectral energy flux in the inertial range. By manipulating the elemental structures, we can tailor turbulent flow fields with different Reynolds numbers, helicities, and local flow geometries.

The customizable woven turbulence can be used to assess turbulence theory [43], develop structure-based turbulence models [15], and offer a playground for studying multi-physics coupled flows [44, 45] and complex fluids [46, 47]. Although the present construction method may not capture all the intricacies of real turbulence, it opens up an avenue for exploring turbulence in the context of matter research and seeking connections between classical/quantum turbulence and other complex systems. Its low cost facilitates generating large training data for machine learning of turbulence [48]. By harnessing this powerful tool, one can incorporate diverse elemental structures, such as Lundgren’s spiral vortex [14], into the woven flow field, and employ quantum wall flows [49] to weave wall turbulence with complex boundary conditions. This approach promises to illuminate valuable structures and their functionalities in understanding and modeling turbulence.

Acknowledgments

The authors are grateful to A. W. Baggaley for providing the code of quantum turbulence simulation, and thank C. F. Barenghi and S. Xiong for helpful comments. Numerical constructions and simulations were carried out on the Tianhe-2A supercomputer in Guangzhou, China. This work has been supported by the National Natural Science Foundation of China (Grant Nos. 11925201 and 11988102), the National Key R&D Program of China (No. 2020YFE0204200), and the Xplore Prize.

Author contributions

W.S. and Y.Y. conceived the idea and designed the research; Y.Y. supervised the project. W.S. constructed woven turbulence, simulated quantum turbulence, and processed data. J.Y. simulated classical turbulence. W.S. and Y.Y. wrote the original draft; all authors discussed the results and reviewed the final manuscript. All the authors have given approval for the manuscript.

Refer to caption — Figure 1: Constructing a turbulent field in classical viscous flow with vortex filaments of quantum turbulence. (A) A tangle of vortex filaments in the statistically stationary stage of quantum turbulence obtained from the simulation of superfluid helium II using the vortex filament method. The simulation is implemented in a cube of size $\mathcal{L}=0.1$ cm with periodic boundary conditions (movie S1). (B) Temporal evolution of the quantum vortex-line density $L_{q}$ , the vortex length per unit volume. After rapid growth at early times, the length of quantum vortex filaments gradually reaches a statistically stationary stage. The inset shows the frequency spectrum (power spectral density) of the fluctuations of $L_{q}$ during the stationary stage. The dashed line denotes the scaling of $f^{-5/3}$ , where $f$ is the frequency. (C) Vortex-surface visualization of an instantaneous woven classical turbulent field with the quantum skeleton in (A). The skeleton size is non-dimensionalized into a periodic box of size $\mathcal{L}=2\pi$ . Each viscous vortex tube is constructed with a quantum filament as the tube centerline, and its vorticity profile is a Gaussian function. The core size $\sigma(s)=0.015(1+0.5\sin s)$ of the vortex tube varies along the arc length parameter $s$ . The isosurface is the normalized vortex-surface field [50] $\phi_{v}=0.1$ and color-coded by the helicity density $h=\bm{u}\cdot\bm{\omega}$ , the dot product of the fluid velocity $\bm{u}$ and the vorticity $\bm{\omega}=\nabla\times\bm{u}$ .

Weaving classical turbulence with quantum skeleton
Supplementary material

Weiyu Shen

State Key Laboratory for Turbulence and Complex Systems, College of Engineering,
Peking University, Beijing 100871, PR China

Jie Yao

Advanced Research Institute of Multidisciplinary Sciences,
Beijing Institute of Technology, Beijing 100081, PR China

Yue Yang

State Key Laboratory for Turbulence and Complex Systems, College of Engineering and HEDPS-CAPT
Peking University, Beijing 100871, PR China

Materials and Methods

.1 Simulation of quantum turbulence

We generate a tangle of vortex filaments in superfluid turbulence without the normal fluid using the vortex filament method (VFM) [28]. The vortex filaments are further used as centerlines in constructing viscous vortex tubes with finite thickness.

In the VFM, the vortex filament moves according to the total induced velocity solved by the de-singularized Biot-Savart integral

\frac{\mathrm{d}\bm{s}}{\mathrm{d}t}=\frac{\kappa_{q}}{4\pi}\bm{s}^{\prime}\times\bm{s}^{\prime\prime}\ln\left(\frac{2\sqrt{l_{+}l_{-}}}{\mathrm{e}^{1/2}a_{0}}\right)+\frac{\kappa_{q}}{4\pi}\int_{\bm{s}_{1}\neq\bm{s}}\frac{(\bm{s}_{1}-\bm{s})\times\mathrm{d}\bm{s}_{1}}{|\bm{s}_{1}-\bm{s}|^{3}}.

(1)

Here, vortex core radius is $a_{0}=10^{-8}$ cm for superfluid ⁴He, circulation $\kappa_{q}=h/m=9.97\times 10^{-4}$ cm²/sec is fixed by quantum constraint of the Planck constant $h$ and bare mass $m$ of a helium atom, $l_{\pm}$ are the lengths of the line segments connected to $\bm{s}$ after discretization, and $\bm{s}^{\prime}$ and $\bm{s}^{\prime\prime}$ are unit tangent and normal vectors at point $\bm{s}$ , respectively.

The computational box is a cube of size $\mathcal{L}=0.1$ cm with periodic boundary conditions. A tree algorithm with opening angle $\theta_{\mathrm{tree}}=0.4$ is used to speed up the evaluation of Biot–Savart integrals. The third-order Adams–Bashforth method is used for time advance with the stepping $\textrm{d}t=5\times 10^{-5}$ sec. The spatial resolution is based on the cutoff scale $\delta$ which is the length of a filament segment. Through vortex reconnection, the size of the smallest vortex filament is cut off at scale $\delta$ [51].

We performed two cases to simulate quantum vortex tangle with different degrees of entanglement by adjusting $\delta=0.01$ cm and $0.005$ cm (Movies S1 and S2). Both initial configurations consist of 50 vortices at random positions and with random orientations. After a rapid transition, the tangles become statistically isotropic and homogeneous (Figs. 1B and S1). We chose the vortex tangle at $t=10$ sec in the statistically stationary stage as the skeleton of woven classical turbulence. The two VFM simulations provide two sets of vortex tangles with different complexities.

.2 Construction of woven classical turbulence

In the construction of woven classical turbulence, first, we used the VFM to simulate the quantum vortex tangle. Then, we reconstructed the control points on the filaments into cubic spline curves defined piecewise by polynomial parametric equations (Supplementary Section II.1). Based on these centerlines $\mathcal{C}$ , the 3D vorticity field for the vortex tubes was constructed in a periodic box of size $\mathcal{L}=2\pi$ on $N^{3}=1024^{3}$ uniform grid points by mapping from the curved cylindrical coordinates centered on $\mathcal{C}$ to Cartesian coordinates (Supplementary Sections II.2 and II.3). Finally, boundary processing was performed so that the vector field satisfies periodic boundary conditions. Two types of boundary problems need to be handled separately. When a centerline crosses the boundary (Type I), we set ghost points for the centerline. When the vorticity for a vortex tube with finite thickness crosses the boundary, while the centerline does not cross the boundary (Type II), we consider the effect of the other $3^{3}-1=26$ periodic boxes surrounding the original domain. The construction process is summarized in Fig. S2.

.3 Direct numerical simulations of classical turbulence

We performed a direct numerical simulation for solving the incompressible Navier–Stokes equations

\frac{\partial{\bm{u}}}{\partial t}+\bm{u}\cdot\nabla\bm{u}=-\nabla p+\nu\nabla^{2}\bm{u}+\bm{f},\quad\nabla\cdot\bm{u}=0

(2)

in a periodic box with the open-source pseudo-spectral code HIT3D [52], where $p$ is the pressure, $\nu$ the kinematic viscosity, and $\bm{f}$ an external forcing for maintaining turbulence [53]. The forcing operates at large scales within a spherical shell of radius $|\bm{k}|\leq 2$ in Fourier space. The domain size is $L^{3}=(2\pi)^{3}$ with $1024^{3}$ grid points, and the Taylor Reynolds number is $Re_{\lambda}=161.4$ , which matches the Reynolds number of the woven turbulence. The turbulence statistics in this simulation are compared with those in the woven turbulence.

I Comparison of different synthetic turbulence methods

We compare different synthetic turbulence methods in Table S1. Most existing methods can capture the scaling laws of low-order statistics of real turbulence, such as energy spectra and second-order spatial correlations. However, only our method of weaving turbulence can simultaneously reproduce low-order statistics, intermittency, and coherent vortices.

Spectral methods [54, 21] construct a turbulent field as a superposition of various Fourier modes. They can successfully achieve low-order statistical properties of the flow, such as the scaling laws of energy spectra or second-order structure functions. However, they cannot render the intermittency and high-order statistics of turbulent flows. Linear eddy models [22] mainly estimate turbulent transport and mixing of scalars.

Fractal models include wavelet-based [55] and multiplicative methods [56, 57, 23, 58, 59], which can reveal fundamental scaling properties and have been applied to model subgrid-scale stress in large eddy simulations of turbulence. These methods artificially prescribe high-order statistics but lack dynamic description. Therefore, the synthetic flows based on fractal models do not contain the basic coherent structures of turbulence.

The multiscale Lagrangian map method [24, 60] can generate a 3D non-Gaussian synthetic velocity field by deforming the superposition of random-phase Fourier modes with a prescribed spectrum. It can reproduce many statistical and geometric properties. However, the main difference from real turbulence is that the regions of dominant enstrophy and low pressure are in sheet-like shapes rather than filaments.

Methods	Low-order statistics	Intermittency	Coherent vortices
Gaussian random field	$\times$	$\times$	$\times$
Spectral methods	$\checkmark$	$\times$	$\times$
Linear-eddy models	$\checkmark$	$\times$	$\times$
Fractal models	$\checkmark$	$\checkmark$	$\times$
Multiscale Lagrangian map	$\checkmark$	$\checkmark$	sheet-like
Weaving turbulence	$\checkmark$	$\checkmark$	customizable

Table S1: Comparison of different synthetic turbulence methods. Low-order statistics include energy spectra and second-order spatial correlations. Intermittency is reproduced by the scaling law of high-order statistics, such as high-order velocity structure functions. The check or cross mark denotes that the feature is or is not reproduced, respectively.

II Numerical methods for weaving a turbulent field

We detail the theoretical and numerical methods for constructing a 3D divergence-free vorticity field from a tangle of quantum vortex filaments. Section II.1 describes how to transform a vortex filament composed of multiple spatial control points into multiple smooth parameter equations based on cubic spline polynomials. Section II.2 describes how to generate a 3D vorticity fields using vortex filaments with their centerlines in a curved cylindrical coordinate system. Section II.3 presents a numerical algorithm for constructing a vorticity field on a Cartesian grid. Section II.4 describes the boundary processing, and discusses the effect of different boundary processing methods on flow statistics.

II.1 Cubic-spline reconstruction of vortex filaments

In the simulation of quantum turbulence with the vortex filament method, a vortex filament is determined by a sequence of discrete control points. To use these vortex filaments as the skeleton of vortices, we need to transform the control points of a spatial curve into an explicit parametric equation.

Under the periodic boundary conditions, all vortex filaments are either closed loops or infinitely extended with periodic repetitions. We can trace each filament from any point along itself (crossing the periodic boundary if necessary), and it will always return to the starting point, so each vortex filament can be considered closed. For a single sequence of discrete control points $\widetilde{\bm{x}}_{i}=(x_{i},y_{i},z_{i}),~{}i=1,2,\ldots,N_{p}$ , we introduce the discrete arc-length parameter

\widetilde{s}_{i}=\left\{\begin{aligned} &-\left|\widetilde{\bm{x}}_{N_{p}}-\widetilde{\bm{x}}_{1}\right|,\quad i=0,\\ &0,\quad i=1,\\ &\sum_{j=2}^{i}\left|\widetilde{\bm{x}}_{j}-\widetilde{\bm{x}}_{j-1}\right|,\quad i=2,3,\ldots,N_{p},\\ &\left|\widetilde{\bm{x}}_{N_{p}}-\widetilde{\bm{x}}_{1}\right|+\sum_{j=2}^{N_{p}+1}\left|\widetilde{\bm{x}}_{j}-\widetilde{\bm{x}}_{j-1}\right|,\quad i=N_{p}+1,\end{aligned}\right.

(3)

and total discrete length $\widetilde{L}=\widetilde{s}_{N_{p}+1}$ . Then we obtain the normalized arc-length parameter

\widetilde{\zeta}_{i}=\frac{2\pi\widetilde{s}_{i}}{\widetilde{L}}\in[0,2\pi],\quad i=1,2,\ldots,N_{p}+1.

(4)

We use the cubic spline parametric equation

\bm{c}(\zeta)=\bm{\alpha}_{i}\zeta^{3}+\bm{\beta}_{i}\zeta^{2}+\bm{\gamma}_{i}\zeta+\bm{\delta}_{i},\quad\zeta\in[\widetilde{\zeta}_{i},\widetilde{\zeta}_{i+1}),\quad i=1,2,\ldots,N_{p},

(5)

to smoothly connect the control points piecewise. In each segment $\zeta\in[\widetilde{\zeta}_{i},\widetilde{\zeta}_{i+1})$ , we determine four vector coefficients $[\bm{\alpha}_{i},\bm{\beta}_{i},\bm{\gamma}_{i},\bm{\delta}_{i}]$ by solving the linear system

\begin{bmatrix}\bm{c}(\widetilde{\zeta}_{i})\\ \bm{c}(\widetilde{\zeta}_{i+1})\\ \dfrac{\textrm{d}\bm{c}}{\textrm{d}\zeta}(\widetilde{\zeta}_{i})\\ \dfrac{\textrm{d}\bm{c}}{\textrm{d}\zeta}(\widetilde{\zeta}_{i+1})\end{bmatrix}=\begin{bmatrix}\widetilde{\zeta}_{i}^{3}&\widetilde{\zeta}_{i}^{2}&\widetilde{\zeta}_{i}&1\\ \widetilde{\zeta}_{i+1}^{3}&\widetilde{\zeta}_{i+1}^{2}&\widetilde{\zeta}_{i+1}&1\\ 3\widetilde{\zeta}_{i}^{2}&2\widetilde{\zeta}_{i}&1&0\\ 3\widetilde{\zeta}_{i+1}^{2}&2\widetilde{\zeta}_{i+1}&1&0\end{bmatrix}\begin{bmatrix}\bm{\alpha}_{i}\\ \bm{\beta}_{i}\\ \bm{\gamma}_{i}\\ \bm{\delta}_{i}\end{bmatrix},

(6)

where the positions

\left\{\begin{aligned} &\bm{c}(\widetilde{\zeta}_{i})=\widetilde{\bm{x}}_{i}\\ &\bm{c}(\widetilde{\zeta}_{i+1})=\widetilde{\bm{x}}_{i+1}\end{aligned}\right.

(7)

and first-order derivatives

\left\{\begin{aligned} &\frac{\textrm{d}\bm{c}}{\textrm{d}\zeta}(\widetilde{\zeta}_{i})=\frac{\widetilde{\bm{x}}_{i+1}-\widetilde{\bm{x}}_{i-1}}{\widetilde{\zeta}_{i+1}-\widetilde{\zeta}_{i-1}}\\ &\frac{\textrm{d}\bm{c}}{\textrm{d}\zeta}(\widetilde{\zeta}_{i+1})=\frac{\widetilde{\bm{x}}_{i+2}-\widetilde{\bm{x}}_{i}}{\widetilde{\zeta}_{i+2}-\widetilde{\zeta}_{i}}\end{aligned}\right.

(8)

of the two endpoints of a segment are calculated from the positions and center differences of the control points, respectively. Substituting Eqs. (7) and (8) into Eq. (6) yields

\begin{bmatrix}\bm{\alpha}_{i}\\ \bm{\beta}_{i}\\ \bm{\gamma}_{i}\\ \bm{\delta}_{i}\end{bmatrix}=\begin{bmatrix}-\dfrac{2}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{3}}&\dfrac{2}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{3}}&\dfrac{1}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{2}}&\dfrac{1}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{2}}\\ \dfrac{3(\widetilde{\zeta}_{j}+\widetilde{\zeta}_{j+1})}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{3}}&-\dfrac{3(\widetilde{\zeta}_{j}+\widetilde{\zeta}_{j+1})}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{3}}&-\dfrac{\widetilde{\zeta}_{j}+2\widetilde{\zeta}_{j+1}}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{2}}&-\dfrac{2\widetilde{\zeta}_{j}+\widetilde{\zeta}_{j+1}}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{2}}\\ -\dfrac{6\widetilde{\zeta}_{j}\widetilde{\zeta}_{j+1}}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{3}}&\dfrac{6\widetilde{\zeta}_{j}\widetilde{\zeta}_{j+1}}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{3}}&\dfrac{\widetilde{\zeta}_{j+1}(2\widetilde{\zeta}_{j}+\widetilde{\zeta}_{j+1})}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{2}}&\dfrac{\widetilde{\zeta}_{j}(\widetilde{\zeta}_{j}+2\widetilde{\zeta}_{j+1})}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{2}}\\ \dfrac{\widetilde{\zeta}_{j+1}^{2}(3\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{3}}&\dfrac{\widetilde{\zeta}_{j}^{2}(\widetilde{\zeta}_{j}-3\widetilde{\zeta}_{j+1})}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{3}}&-\dfrac{\widetilde{\zeta}_{j}\widetilde{\zeta}_{j+1}^{2}}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{2}}&-\dfrac{\widetilde{\zeta}_{j}^{2}\widetilde{\zeta}_{j+1}}{(\widetilde{\zeta}_{j}-\widetilde{\zeta}_{j+1})^{2}}\end{bmatrix}\begin{bmatrix}\widetilde{\bm{x}}_{i}\\ \widetilde{\bm{x}}_{i+1}\\ \dfrac{\widetilde{\bm{x}}_{i+1}-\widetilde{\bm{x}}_{i-1}}{\widetilde{\zeta}_{i+1}-\widetilde{\zeta}_{i-1}}\\ \dfrac{\widetilde{\bm{x}}_{i+2}-\widetilde{\bm{x}}_{i}}{\widetilde{\zeta}_{i+2}-\widetilde{\zeta}_{i}}\end{bmatrix}.

(9)

A special treatment is needed for where two consecutive control points are on opposite sides of a boundary. To make sure that the spline curves cross the boundary correctly and meet the periodic boundary conditions, we set ghost points on both sides of the boundary. Finally, we reconstruct all $N_{v}$ vortex filaments composed of discrete points $\widetilde{\bm{x}}_{i}$ into $N_{v}$ smooth cubic splines $\mathcal{C}^{(k)}:\bm{c}^{(k)}(\zeta),~{}k=1,2,\ldots,N_{v}$ , determined by Eq. (5) with continuous parameters $\zeta\in[0,2\pi)$ and coefficients in Eq. (9).

For any given vortex filament, the total arc-length of the spline is no less than the total discrete length of the initial sequence of control points, i.e.,

L=\oint\left|\frac{\textrm{d}\bm{c}(\zeta)}{\textrm{d}\zeta}\right|d\zeta\geqslant\widetilde{L},

(10)

where the equality only holds when the filament is a straight line. Therefore, the continuous arc-length parameter

s(\zeta)=\int_{0}^{\zeta}\left|\frac{\textrm{d}\bm{c}(\zeta)}{\textrm{d}\zeta}\right|d\zeta\in[0,L)

(11)

of a spline $\mathcal{C}$ is inconsistent with the discrete arc-length parameter $\widetilde{s}_{i}\in[0,\widetilde{L})$ of the control points on $\mathcal{C}$ , and the mapping $\widetilde{s}_{i}\mapsto\widetilde{\zeta}_{i}\mapsto s(\zeta=\widetilde{\zeta}_{i})$ is bijective by Eqs. (4) and (11).

II.2 Theoretical construction of the vorticity field

We construct the vorticity field $\bm{\omega}$ for multiple tube-like vortices with tunable internal structures. These vortices are generated based on their centerlines described Eq. (5). This construction method with its numerical algorithm is an extension of that in Refs. [61, 62, 40] by incorporating various simple vortex models and entanglement of multiple vortices.

For a single vortex centerline $\mathcal{C}$ , we introduce the curved cylindrical coordinate system $(s,\rho,\theta)$ [63, 62] based on the Frenet–Serret frame $(\bm{T},\bm{N},\bm{B})$ in a tubular region surrounding curve $\mathcal{C}$ . The system satisfies

\bm{x}(s,\rho,\theta)=\bm{c}(s)+\rho\cos\theta\bm{N}(s)+\rho\sin\theta\bm{B}(s)

(12)

within a tubular region

\Omega_{C}=\left\{\bm{x}\left|\bm{x}\in\mathbb{R}^{3},\min_{s}\right|\bm{x}-\bm{c}(s)\mid<R_{v}\right\},

(13)

where $R_{v}$ denotes the radius of $\Omega_{C}$ . The Frenet–Serret formulas on $\mathcal{C}$ are

\left\{\begin{aligned} \frac{\mathrm{d}\bm{T}}{\mathrm{d}s}&=\kappa\bm{N},\\ \frac{\mathrm{d}\bm{N}}{\mathrm{d}s}&=-\kappa\bm{T}+\tau\bm{B},\\ \frac{\mathrm{d}\bm{B}}{\mathrm{d}s}&=-\tau\bm{N},\end{aligned}\right.

(14)

where $\bm{T}\equiv\mathrm{d}\bm{c}/\mathrm{d}s$ denotes the unit tangent, $\bm{N}\equiv(\mathrm{d}\bm{T}/\mathrm{d}s)/|\mathrm{d}\bm{T}/\mathrm{d}s|$ the unit normal, $\bm{B}\equiv\bm{T}\times\bm{N}$ the unit binormal, $\kappa$ the curvature, and $\tau$ the torsion of $\mathcal{C}$ . Based on coordinates ( $s,\rho,\theta$ ), we specify the vorticity field

\bm{\omega}(s,\rho,\theta)=\omega_{s}(s,\rho,\theta)\bm{e}_{s}+\omega_{\rho}(s,\rho,\theta)\bm{e}_{\rho}+\omega_{\theta}(s,\rho,\theta)\bm{e}_{\theta}

(15)

of the tubular region, where the local frame is spanned by unit vectors

\left\{\begin{array}[]{l}\bm{e}_{s}=\bm{T},\\ \bm{e}_{\rho}=\cos\theta\bm{N}+\sin\theta\bm{B},\\ \bm{e}_{\theta}=-\sin\theta\bm{N}+\cos\theta\bm{B}.\end{array}\right.

(16)

Figure S3A illustrates a segment of the tubular region $\Omega_{C}$ in the curved cylindrical coordinate system $(s,\rho,\theta)$ .

We construct $\bm{\omega}$ based on the geometric relation between vortex lines and surfaces. The vortex surface field (VSF) $\phi_{v}$ is defined to satisfy the constraint [50]

\bm{\omega}\cdot\nabla\phi_{v}=0,

(17)

so that $\bm{\omega}$ is tangent to the isosurface of $\phi_{v}$ everywhere except at $\bm{\omega}=0$ , and every isosurface of $\phi_{v}$ is a vortex surface consisting of vortex lines.

For vortex tubes with the variable core size $\sigma(s)$ and local twist rate $\eta_{t}(s,\phi_{v})$ [40], the vorticity of vortex tubes with differential twist and variable thickness is specified as

\bm{\omega}(s,\rho,\theta)=\Gamma f_{t}(s,\rho)\left[\underbrace{\bm{e}_{s}}_{\text{flux}}+\underbrace{\frac{\rho}{\sigma(s)(1-\kappa(s)\rho\cos\theta)}\frac{\mathrm{d}\sigma(s)}{\mathrm{d}s}\bm{e}_{\rho}}_{\text{tube thickness}}+\underbrace{\frac{\rho\eta_{t}(s,\phi_{v})}{1-\kappa(s)\rho\cos\theta}\bm{e}_{\theta}}_{\text{twist}}\right],

(18)

with the circulation $\Gamma$ , Gaussian kernel function

\displaystyle f_{t}(s,\rho)=\left\{\begin{aligned} &\frac{1}{2\pi\sigma(s)^{2}}\exp\left[\frac{-\rho^{2}}{2\sigma(s)^{2}}\right],\quad\rho\in[0,R_{v}),\\ &0,\quad\rho\in[R_{v},+\infty)\end{aligned}\right.

(19)

as in the Burgers vortex model, and the normalized VSF

\phi_{v}(s,\rho)=2\pi\sigma(s)^{2}f_{t}(s,\rho)\in[0,1].

(20)

Here, the three terms on the right-hand side of Eq. (18) represent the vorticity flux, tube-thickness, and twist terms of $\bm{\omega}$ , respectively. The vector field constructed by Eq. (18) is proved to be divergence-free [40]. The radius $R_{v}$ of the constructed tubular region $\Omega_{C}$ is large enough to ensure that almost all circulation is included. The cross-sectional profile of the vortex tubes is shown in Fig. S3B.

For constructing sheet-like vortices, we introduce a 2D oblate normal kernel function

\displaystyle f_{s}(s,\rho,\theta)=\left\{\begin{aligned} &\frac{1}{2\pi\sigma_{1}\sigma_{2}}\exp\left[-\frac{1}{2}\left(\frac{\rho^{2}\cos^{2}\theta}{\sigma_{1}^{2}}+\frac{\rho^{2}\sin^{2}\theta}{\sigma_{2}^{2}}\right)\right],\quad\rho\in[0,R_{v}),\\ &0,\quad\rho\in[R_{v},+\infty)\end{aligned}\right.

(21)

with $\sigma_{2}\gg\sigma_{1}$ . We truncate the elliptical distribution at $\rho=R_{v}$ , so that the vorticity on a cross section is concentrated within a sheet-like region, which is centered on $\mathcal{C}$ , and has thickness $2\sigma_{1}$ and width $2R_{v}$ (see Fig. S3C). However, this operation reduces the total circulation, so we specify the vorticity

\bm{\omega}(s,\rho,\theta)=\frac{\Gamma}{\Phi}f_{s}(s,\rho,\theta)\bm{e}_{s}

(22)

with a compensation coefficient

\Phi=\int_{0}^{R_{v}}\int_{0}^{2\pi}\frac{\rho}{2\pi\sigma_{1}\sigma_{2}}\exp\left[-\frac{1}{2}\left(\frac{\rho^{2}\cos^{2}\theta}{\sigma_{1}^{2}}+\frac{\rho^{2}\sin^{2}\theta}{\sigma_{2}^{2}}\right)\right]\mathrm{d}\theta\mathrm{d}\rho

(23)

to make the circulation equal to $\Gamma$ .

Finally, we obtain the vorticity

\bm{\omega}=\sum_{k=1}^{N_{v}}\bm{\omega}^{(k)}

(24)

for a turbulent field consisting of $N_{v}$ intertwined vortex tubes by adding up the vorticities separately constructed from the skeleton $\mathcal{C}^{(k)},~{}k=1,2,\ldots,N_{v}$ .

II.3 Numerical construction of the vorticity field

It is straightforward to extend the numerical algorithm in Ref. [40] to compute Eq. (18) or (22) on a Cartesian grid with $N_{x}\times N_{y}\times N_{z}$ uniform grid points. For a given parametric curve $\mathcal{C}:c(\zeta)$ with $\zeta\in\left[0,2\pi\right)$ , we divide $\mathcal{C}$ into $N_{C}$ segments by $N_{C}$ dividing points

\bm{c}_{i}=\bm{c}\left(\zeta_{i}\right),\quad i=1,2,\ldots,N_{C},

(25)

with $\zeta_{i}=(i-1)\Delta\zeta$ and $\Delta\zeta=2\pi/N_{C}$ . Note that $\zeta$ is not necessary to be an arc-length parameter $s$ because of the one-to-one mapping between $\zeta$ and $s$ . Then the space in the proximity of curve $\mathcal{C}$ can be divided into $N_{C}$ subdomains

\Omega_{i}=\left\{\bm{x}\mid\left(\bm{x}-\bm{c}_{i}\right)\cdot\bm{T}_{i}\geqslant 0\text{ and }\left(\bm{x}-\bm{c}_{i+1}\right)\cdot\bm{T}_{i+1}<0\right\},

(26)

with

\bm{T}_{i}=\frac{\bm{c}_{i+1}-\bm{c}_{i}}{\left|\bm{c}_{i+1}-\bm{c}_{i}\right|},\quad i=1,2,\ldots,N_{C},

(27)

where subscripts $N_{C}+1$ and $1$ are equivalent. For a given $\bm{x}$ , we first use Eq. (26) to determine the subdomains $\Omega_{i}$ containing $\bm{x}$ . The subscripts of all the $\Omega_{i}$ containing $\bm{x}$ are denoted by a set

\widetilde{I}_{\zeta}(\bm{x})=\left\{j\mid\bm{x}\in\Omega_{j}\right\}.

(28)

For each $j\in\widetilde{I}_{\zeta}(\bm{x})$ , the parameter of $\mathcal{C}$ is approximated by

\widetilde{\zeta}_{j}=\frac{\zeta_{j+1}\left(\bm{x}-\bm{c}_{j}\right)\cdot\bm{T}_{j}+\zeta_{j}\left(\bm{c}_{j+1}-\bm{x}\right)\cdot\bm{T}_{j}}{\left|\bm{c}_{j+1}-\bm{c}_{j}\right|}.

(29)

At $\widetilde{\bm{c}}_{j}=\bm{c}\left(\widetilde{\zeta}_{j}\right)$ , we use the second-order finite difference scheme to calculate the Frenet-Serret frame

\left\{\begin{aligned} \widetilde{\bm{T}}_{j}&=\bm{T}\left(\widetilde{\zeta}_{j}\right),\\ \widetilde{\bm{N}}_{j}&=\bm{N}\left(\widetilde{\zeta}_{j}\right),\\ \widetilde{\bm{B}}_{j}&=\bm{B}\left(\widetilde{\zeta}_{j}\right),\end{aligned}\right.

(30)

as well as

\left\{\begin{aligned} &\widetilde{\kappa}_{j}=\kappa\left(\widetilde{\zeta}_{j}\right),\\ &\widetilde{\frac{\textrm{d}\sigma}{\textrm{d}s}}_{j}=\frac{\textrm{d}\sigma}{\textrm{d}s}\left(\widetilde{\zeta}_{j}\right)\end{aligned}\right.

(31)

in Eq. (18). In addition, the distance from $\bm{c}\left(\widetilde{\zeta}_{j}\right)$ is calculated by

\widetilde{\rho}_{j}=\left|\bm{x}-\widetilde{\bm{c}}_{j}\right|,

(32)

and azimuth-related functions are calculated by

\left\{\begin{aligned} \cos\widetilde{\theta}_{j}=\frac{\left(\bm{x}-\widetilde{c}_{j}\right)\cdot\widetilde{\bm{N}}_{j}}{\widetilde{\rho}_{j}},\\ \sin\widetilde{\theta}_{j}=\frac{\left(\bm{x}-\widetilde{c}_{j}\right)\cdot\widetilde{\bm{B}}_{j}}{\widetilde{\rho}_{j}}.\end{aligned}\right.

(33)

We approximate Eq. (18) or (22) as

\bm{\mathcal{\omega}}(\bm{x})=\sum_{j\in\widetilde{I}_{\zeta}(\bm{x})}\bm{\widetilde{\omega}}_{j},

(34)

for vortex tubes

\bm{\widetilde{\omega}}_{j}=\left\{\begin{aligned} &\Gamma f_{t}\left(\widetilde{\zeta}_{j},\widetilde{\rho}_{j}\right)\left[\widetilde{\bm{T}}_{j}+\frac{\widetilde{\rho}_{j}}{\sigma\left(\widetilde{\zeta}_{j}\right)\left(1-\widetilde{\kappa}_{j}\widetilde{\rho}_{j}\cos\widetilde{\theta}_{j}\right)}\widetilde{\frac{\textrm{d}\sigma}{\textrm{d}s}}_{j}\left(\cos\widetilde{\theta}_{j}\widetilde{\bm{N}}_{j}+\sin\widetilde{\theta}_{j}\widetilde{\bm{B}}_{j}\right)\right.\\ &\qquad+\left.\frac{\widetilde{\rho}_{j}\eta_{t}\left(\widetilde{\zeta}_{j},\phi_{v}\left(\widetilde{\zeta}_{j},\widetilde{\rho}_{j}\right)\right)}{1-\widetilde{\kappa}_{j}\widetilde{\rho}_{j}\cos\widetilde{\theta}_{j}}\left(-\sin\widetilde{\theta}_{j}\widetilde{\bm{N}}_{j}+\cos\widetilde{\theta}_{j}\widetilde{\bm{B}}_{j}\right)\right],\quad 1>\widetilde{\kappa}_{j}\widetilde{\rho}_{j}\cos\widetilde{\theta}_{j}\\ &\mathbf{0},\quad 1\leqslant\widetilde{\kappa}_{j}\widetilde{\rho}_{j}\cos\widetilde{\theta}_{j}\end{aligned}\right.

(35)

or vortex sheets

\bm{\widetilde{\omega}}_{j}=\left\{\begin{aligned} &\frac{\Gamma}{\Phi}f_{s}\left(\widetilde{\zeta}_{j},\widetilde{\rho}_{j},\widetilde{\theta}_{j}\right)\widetilde{\bm{T}}_{j},\quad 1>\widetilde{\kappa}_{j}\widetilde{\rho}_{j}\cos\widetilde{\theta}_{j},\\ &\mathbf{0},\quad 1\leqslant\widetilde{\kappa}_{j}\widetilde{\rho}_{j}\cos\widetilde{\theta}_{j}\end{aligned}\right.

(36)

with computed and given variables. Finally, we construct $N_{v}$ intertwined vortex tubes by adding $N_{v}$ vorticity fields calculated by Eq. (34) using Eq. (24). The procedure for the numerical construction of $\bm{\omega}(\bm{x})$ is summarized in Algorithm 1.

Input:

\bm{x}

\widetilde{\bm{x}}_{i}

dx

dy

dz

f_{t}(\zeta,\rho)

f_{s}(\zeta,\rho,\theta)

\phi_{v}(\zeta,\rho)

\eta_{t}(\zeta,\phi_{v})

\sigma(\zeta)

\sigma_{1}

\sigma_{2}

R_{v}

\Gamma

N_{x}

N_{y}

N_{z}

N_{v}

N_{p}

, and

N_{C}

Output:

\bm{\omega}(\bm{x})

1 initialization;

3for $k\leftarrow 1$ to $N_{v}$ do

4 for $i\leftarrow 1$ to $N_{p}$ do

5 Calculate

\widetilde{s}_{i}

by Eq. (3);

6 Calculate

\widetilde{\zeta}_{i}

by Eq. (4) ;

7 Calculate

\bm{\alpha}_{i},\bm{\beta}_{i},\bm{\gamma}_{i},\bm{\delta}_{i}

by Eq. (9);

8 Obtain

\bm{c}(\zeta)

by Eq. (5) ;

10 end for

11 for $i_{x}\leftarrow 1$ to $N_{x}$ do

12 for $i_{y}\leftarrow 1$ to $N_{y}$ do

13 for $i_{z}\leftarrow 1$ to $N_{z}$ do

14 for $i\leftarrow 1$ to $N_{C}$ do

15 Calculate

\widetilde{\bm{T}}_{i}

by Eq. (27);

16 Divide the space in the proximity of curve

\bm{c}(\zeta)

into the subdomain

\Omega_{i}

by Eq. (26);

17 Obtain

\widetilde{I}_{\zeta}

by Eq. (28) at

\bm{x}((i_{x}-1)dx,(i_{y}-1)dy,(i_{z}-1)dz)

;

19 end for

20 for $i\leftarrow 1$ to $N_{C}$ do

21 if $i\in\widetilde{I}_{\zeta}$ then

22 Calculate

\widetilde{\zeta}_{j}

by Eq. (29);

23 Calculate

\widetilde{\bm{T}}_{j},\widetilde{\bm{N}}_{j}

, and

\widetilde{\bm{B}}_{j}

by Eq. (30);

24 Calculate

\widetilde{\kappa}_{j}

and

\widetilde{\frac{\textrm{d}\sigma}{\textrm{d}s}}_{j}

by Eq. (31);

25 Calculate

\widetilde{\rho}_{j}

by Eq. (32);

26 Calculate

\cos\tilde{\theta}_{j}

and

\sin\tilde{\theta}_{j}

by Eq. (33);

27 if vortex tube then

28 Calculate

\bm{\omega}^{(k)}(\bm{x})

by Eqs. (34) and (35) with computed and given variables;

30 end if

31 if vortex sheet then

32 Calculate

\Phi

by Eq. (23);

33 Calculate

\bm{\omega}^{(k)}(\bm{x})

by Eqs. (34) and (36) with computed and given variables;

35 end if

37 end if

39 end for

41 end for

43 end for

45 end for

\bm{\omega}\leftarrow\bm{\omega}+\bm{\omega}^{(k)}

47 end for

Algorithm 1 Calculation of

\bm{\omega}(\bm{x})

for intertwined vortex tubes/sheets

II.4 Boundary processing

We present two types of boundary problems at the bottom of Fig. S2. For centerlines that cross the boundary (Type I), we add ghost points to extend the centerlines. For vorticity that crosses the boundary due to the tube thickness (Type II), we account for the effect of the other $3^{3}-1=26$ periodic boxes surrounding the computational domain. These treatments increase the computational cost.

To examine the influence of boundary processing on flow statistics, we compare the results from the two processing methods in Fig. S4. We used six coplanar boxes around the flow field for face processing and all 26 boxes for complete processing. The boundary vorticity is negligible compared to the whole vorticity field, so it has a negligible effect on flow statistics. Note that the boundary processing should be applied to the initial field to recover periodic boundary conditions in numerical simulations of flow evolution.

III Comparison of turbulence statistics

We compare important statistics between woven turbulence and real classical turbulence. First, we use the second-order structure function to estimate the mean dissipation rate in Section III.1. Then, we apply the statistical theory for homogeneous isotropic turbulence (HIT) to obtain other flow statistics and relevant length scales in Sections III.1 and III.2. We also compute the kinetic energy spectra and nth-order longitudinal velocity structure functions and compare them in woven and classical turbulence in Section III.3. Finally, we evaluate and compare the kinetic energy flux and transfer function in Section III.4.

III.1 Flow statistics

The velocity in woven turbulence nearly satisfies the Gaussian distribution, as shown by the probability distribution functions (PDFs) in Fig. S5, implying that the flow field with intertwined vortex tubes is statistically isotropic.

In the statistical theory of HIT, the second-order longitudinal velocity structure function

D_{LL}(r)=\left\langle\left[u_{1}\left(\bm{x}+\bm{e}_{1}r\right)-u_{1}(\bm{x})\right]^{2}\right\rangle

(37)

measures the covariance of the velocity difference between two points, $\bm{x}+\bm{e}_{1}r$ and $\bm{x}$ , along the same direction $\bm{e}_{1}$ . The celebrated Kolmogorov 1941 (K41) theory [30] suggests

D_{LL}(r)=C_{2}\left\langle\varepsilon\right\rangle^{2/3}r^{2/3}

(38)

in the inertial range, where $\left\langle\varepsilon\right\rangle$ is the mean dissipation rate and $C_{2}$ is the Kolmogorov constant. The scaling exponent of $D_{LL}$ for HIT reaches $2/3$ in the inertial range but drops below $2/3$ in the energy-containing and dissipation ranges. Therefore, we can estimate the mean dissipation rate from the inertial range by

\left\langle\varepsilon\right\rangle=\left[\frac{C_{2}}{\max(r^{-2/3}D_{LL})}\right]^{-3/2}

(39)

where we adopt $C_{2}=2.2$ as a widely accepted value [32].

We use the mean enstrophy $\left\langle\Omega\right\rangle=\langle\left|\bm{\omega}\right|^{2}\rangle$ to estimate the effective kinematic viscosity

\nu=\frac{\left\langle\varepsilon\right\rangle}{\left\langle\Omega\right\rangle}.

(40)

Finally, we calculate the Taylor microscale

\lambda=\sqrt{\frac{10k_{t}}{\left\langle\Omega\right\rangle}}

(41)

and Taylor Reynolds number

Re_{\lambda}\equiv\frac{u^{\prime}\lambda}{\nu}=\frac{k_{t}}{\nu}\sqrt{\frac{20}{3\left\langle\Omega\right\rangle}},

(42)

where $u^{\prime}=\sqrt{\left\langle u_{i}u_{i}\right\rangle/3}$ is the turbulent root-mean-square (RMS) velocity and $k_{t}=3{u^{\prime}}^{2}/2$ is the turbulent kinetic energy. Table S2 summarizes the computed values of statistics for woven turbulence as shown in Fig. 1C.

Flow statistic	$k_{t}$	$\left\langle\Omega\right\rangle$	$\left\langle\varepsilon\right\rangle$	$\nu$	$Re_{\lambda}$
Value	3.12	5434.3	3.67	$6.76\times 10^{-4}$	161.4

Table S2: Flow statistics of woven turbulence: turbulent kinetic energy

k_{t}

, mean enstrophy

\left\langle\Omega\right\rangle

, mean dissipation rate

\left\langle\varepsilon\right\rangle

, kinematic viscosity

\nu

, and Taylor Reynolds number

Re_{\lambda}

III.2 Length scales

We calculate the length scales in turbulent flow fields. The Kolmogorov length scale

\eta\equiv\left(\frac{\nu^{3}}{\left\langle\varepsilon\right\rangle}\right)^{\frac{1}{4}}

(43)

is the characteristic length scale of the smallest turbulent motion in viscous flow. The integral length scale

L\equiv\frac{k_{t}^{3/2}}{\left\langle\varepsilon\right\rangle}

(44)

characterizes the large eddies in viscous flow. The core size

\sigma(s)=0.015(1+0.5\sin s)

(45)

characterizes the length scale of a viscous vortex tube.

For quantum turbulence, the average distance between quantum vortices is estimated by

l_{q}=L_{q}^{-1/2},

(46)

where $L_{q}$ is the quantum vortex line density (i.e., the vortex length per unit volume). Since the quantum skeleton is scaled dimensionlessly into a periodic box of size $\mathcal{L}=2\pi$ , the average distance between vortices in woven turbulence is

l=l_{q}\frac{\mathcal{L}}{\mathcal{L}_{q}},

(47)

where the box size of quantum turbulence is $\mathcal{L}_{q}=0.1$ cm. Table S3 lists the length scales of woven turbulence.

Length scale	$\eta$	$\sigma$	$\lambda$	$l$	$L$	$\mathcal{L}$
Value	$3.03\times 10^{-3}$	$(0.75\sim 2.25)\times 10^{-2}$	$7.57\times 10^{-2}$	$3.14\times 10^{-1}$	1.50	$2\pi$

Table S3: Length scales of woven turbulence: Kolmogorov length scale

\eta

, core size

\sigma

, Taylor microscale

\lambda

, average distance between vortices

l

, integral length scale

L

, and box size

\mathcal{L}

III.3 Energy spectrum and structure functions

We calculate the kinetic energy spectra and nth-order longitudinal velocity structure functions of woven turbulence and compare them with direct numerical simulation (DNS) and model results of HIT with identical flow parameters. The 3D kinetic energy spectrum

E(k)=\iint\hat{\bm{u}}^{*}(\bm{k})\cdot\hat{\bm{u}}(\bm{k})~{}\mathrm{d}A(k)

(48)

represents the contribution to turbulent kinetic energy by wavenumbers from $k$ to $k+\mathrm{d}k$ , where $\hat{\bm{u}}(\bm{k})$ denotes the Fourier modes of velocity field $\bm{u}$ at wavenumber $\bm{k}$ , the superscript “*” denotes the complex conjugate, and $\mathrm{d}A(k)$ is the area of a surface element on a sphere with radius $k=\left|\bm{k}\right|$ in the Fourier space. We plot the rescaled 3D energy spectrum of woven turbulence in Fig. S6A and the compensated energy spectrum in Fig. S6B.

The nth-order longitudinal velocity structure function

D_{n}(r)=\left\langle\left[u_{1}\left(\bm{x}+\bm{e}_{1}r\right)-u_{1}(\bm{x})\right]^{n}\right\rangle

(49)

measures the statistical correlation of velocity differences along the direction of a unit vector $\bm{e}_{1}$ . The K41 theory [30] suggests

D_{LL}(r)\sim\left(\left\langle\varepsilon\right\rangle r\right)^{2/3}

(50)

and

D_{LLL}(r)\sim\left\langle\varepsilon\right\rangle r.

(51)

To account for the intermittency effects into higher-order structure functions, we adopt the She and Leveque 1994 (SL94) model [64]

D_{n}(r)\sim\left(\left\langle\varepsilon\right\rangle r\right)^{\zeta(n)}

(52)

with

\zeta(n)=\frac{1}{9}n+2\left[1-\left(\frac{2}{3}\right)^{\frac{n}{3}}\right].

(53)

We plot the rescaled longitudinal 2nd- to 5th-order velocity structure functions of woven turbulence in Figs. S6C-F. Considering the intermittency corrections, the 4th- and 5th-order structure functions are rescaled with the SL94 scaling in Eq. (53).

To compare statistics in woven and real classical turbulence, we performed a DNS of the incompressible Navier–Stokes equation (see Materials and Methods Section C for more details). Moreover, we used a simple model [32]

D_{LL}(r)=2{u^{\prime}}^{2}\left(\frac{r^{2}}{r^{2}+\left(15C_{2}\right)^{3/2}\eta^{2}}\right)^{2/3-\mu/18}\left(\frac{r^{2}}{r^{2}+\left[4/\left(3C_{2}\right)\right]^{3}L^{2}}\right)^{1/3+\mu/18}

(54)

with the intermittency exponent $\mu=0.25$ . The corresponding $E(k)$ can be calculated from Eq. (54) by

E(k)=\frac{1}{2}k^{3}\frac{d}{dk}\left(\frac{1}{k}\frac{dE_{11}(k)}{dk}\right),

(55)

where

E_{11}(k)=\frac{2{u^{\prime}}^{2}}{\pi}\int_{0}^{\infty}\left(1-\frac{1}{2{u^{\prime}}^{2}}D_{LL}(r)\right)\cos(kr)dr

(56)

is the one-dimensional energy spectrum.

Figure S6A shows that the energy spectrum of woven turbulence agrees well with those from the model and DNS. It not only shows the $-5/3$ scaling law in the inertial range but also exhibits the bottleneck effect of energy transfer between the inertial and dissipation ranges in Fig. S6B.

Figures S6B-F show that the rescaled 2nd- to 5th-order longitudinal structure functions have the same scaling law in Eq. (52) in the inertial range as in classical turbulence. Note that the plateaus for the inertial range are short due to the low Reynolds number in the woven turbulence. The negative 3rd-order structure function in Fig. S6D suggests the negative skewness of velocity fluctuations in woven turbulence as in real classical turbulence. This consistency reveals the woven and real turbulent fields have similar fine-scale vortices, which is an essential difference from the Gaussian random field. Moreover, the amplitudes of higher-order structure functions in woven turbulence have significant differences from those observed in DNS. These results may be improved by using more complex fine-scale elemental vortices.

III.4 Energy flux and transfer

The spectral kinetic energy flux and transfer support the classical Richardson cascade in woven turbulence in Fig. 2D. We calculate the spectral kinetic energy flux by

\Pi_{K}(k)=\sum_{\left|\bm{k}^{\prime}\right|\leqslant k}\operatorname{Re}\left\{\bm{u}^{*}\left(\bm{k}^{\prime}\right)\cdot\mathcal{F}[(\bm{u}\cdot\bm{\nabla})\bm{u}]\left(\bm{k}^{\prime}\right)\right\}.

(57)

This represents the net transfer of energy from all eddies with wavenumber smaller than $k$ to those with wavenumber larger than $k$ . Then, we obtain

\Pi_{K}(k)=-\int_{0}^{k}T(k)dk=\int_{k}^{\infty}T(k)dk,

(58)

where $T(k)$ is the spectral kinetic energy transfer function, which encodes the energy gain or loss at $k$ .

In the inertial subrange, we expect

\Pi_{K}=-\frac{\partial}{\partial t}\int_{0}^{k}Edk=\left\langle\varepsilon\right\rangle,

(59)

since $\left\langle\varepsilon\right\rangle$ is the rate of loss of energy from large eddies. Thus we can also estimate the mean dissipation rate using $\langle\hat{\varepsilon}\rangle=\max(\Pi_{K})=2.67$ . Note that $\langle\hat{\varepsilon}\rangle\approx 73\%\langle\varepsilon\rangle$ is lower than $\langle\varepsilon\rangle$ which is estimated by the second-order structure function in Eq. (39).

In the classical turbulence theory, the kinetic energy cascade removes energy from large scales, transfers it to small scales through the inertial range, and dissipates it by viscosity at the smallest scales. In Fig. S7A, the energy flux increases with $k$ in the large-scale, energy containing range. It is transferred from large to small scales with a constant energy flux in the inertial range and decreases with $k$ in the small-scale, energy dissipation range. Note that the inertial range is very short due to the low Reynolds number in the woven turbulence. Hence, the energy transfer function is negative in the large-scale energy-containing range, positive in the small-scale dissipation range, and around zero in the inertial range. We compare the spectral kinetic energy flux and transfer function of woven turbulence and DNS results of classical turbulence in Fig. S7. The spectral energy transfer of woven turbulence agrees with that of classical turbulence. This suggests that the viscous vortices introduced in the woven turbulence cause the shift from the quantum cascade to the classical one.

III.5 Helicity decomposition

The total helicity

H=\iiint h\,\mathrm{d}\mathcal{V}

(60)

is the volume integral of the helicity density $h=\bm{u}\cdot\bm{\omega}$ .

For multiple closed vortex tubes, the helicity can be topologically decomposed into

H=\sum_{i\neq j}\Gamma_{i}\Gamma_{j}L_{k,ij}+\sum_{i}\Gamma_{i}^{2}(W_{r,i}+T_{w,i}),

(61)

where $\Gamma_{i}$ is the circulation of vortex tube $i$ , $L_{k,ij}$ is the Gauss linking number between vortex tubes $i$ and $j$ , and $W_{r,i}$ and $T_{w,i}$ are the writhing and twisting numbers of vortex tube $i$ , respectively [65, 41, 42]. The writhing number measures the bending or knotting of the vortex centerline $\mathcal{C}$ , and the twisting number measures the twisting of the vorticity field around the centerline.

The writhing number is computed by

W_{r}=\frac{1}{4\pi}\oint_{\mathcal{C}}\oint_{\mathcal{C}}\frac{\left(\bm{x}-\bm{x}^{*}\right)\cdot\mathrm{d}\bm{x}\times\mathrm{d}\bm{x}^{*}}{\left|\bm{x}-\bm{x}^{*}\right|^{3}},

(62)

where $\bm{x}$ and $\bm{x}^{*}$ are two points on $\mathcal{C}$ . The twisting number

T_{w}=\frac{1}{2\pi}\oint_{\mathcal{C}}\eta_{t}\mathrm{~{}d}s

(63)

is defined by the local twist rate

\eta_{t}\left(s\right)=\left(\bm{N}_{s}\times\bm{N}_{s}^{\prime}\right)\cdot\bm{T}

(64)

of vortex lines $\mathcal{C}^{*}$ along the centerline $\mathcal{C}$ , where $\bm{N}_{s}$ is a radial unit vector from $\mathcal{C}$ to $\mathcal{C}^{*}$ in plane $S_{C}$ , $\bm{N}_{s}^{\prime}=\mathrm{d}\bm{N}_{s}/\mathrm{d}s$ , and $\bm{T}$ is the unit tangent vector of $\mathcal{C}$ [40, 63]. The Gauss linking number is given by

L_{k,ij}=L_{k}(\mathcal{C}_{i},\mathcal{C}_{j})=\frac{1}{4\pi}\oint_{\mathcal{C}_{i}}\oint_{\mathcal{C}_{j}}\frac{\left(\bm{x}_{i}-\bm{x}_{j}\right)\cdot\mathrm{d}\bm{x}_{i}\times\mathrm{d}\bm{x}_{j}}{\left|\bm{x}_{i}-\bm{x}_{j}\right|^{3}},

(65)

where $\bm{x}_{i}$ and $\bm{x}_{j}$ are points on $\mathcal{C}_{i}$ and $\mathcal{C}_{j}$ , respectively.

The helical wave decomposition [66] can separate a flow field into purely right- or left-handed polarized parts. In Fourier space, each wave vector has two helical wave modes for the velocity field

\hat{\bm{u}}(\bm{k})=\hat{\bm{u}}^{+}(\bm{k})+\hat{\bm{u}}^{-}(\bm{k})=\hat{u}^{+}(\bm{k})\bm{h}^{+}(\bm{k})+\hat{u}^{-}(\bm{k})\bm{h}^{-}(\bm{k}),

(66)

where $\bm{h}^{\pm}$ are the eigenvectors of the curl operator that satisfy $\mathrm{i}\bm{k}\times\bm{h}^{\pm}=\pm|\bm{k}|\bm{h}^{\pm}$ , and $u^{\pm}$ are the corresponding Fourier coefficients. In physical space, the velocity field

\bm{u}(x)=\bm{u}^{R}(x)+\bm{u}^{L}(x)+\nabla\phi

(67)

can be written as a sum of complex helical waves and the gradient of a harmonic potential $\phi$ , with $\bm{u}^{R}(\bm{x})=\int\hat{\bm{u}}^{+}(\bm{k})\mathrm{d}\bm{k}$ and $\bm{u}^{L}(\bm{x})=\int\hat{\bm{u}}^{-}(\bm{k})\mathrm{d}\bm{k}$ . Here, $\bm{u}^{R}$ is the projection of $\bm{u}$ onto the vector space spanned by all eigenfunctions with positive eigenvalues $(+1|\bm{k}|)$ of the curl operator, and is the right-handed component of $\bm{u}$ ; $\bm{u}^{L}$ is a linear combination of eigenfunctions with negative eigenvalues $(-1|\bm{k}|)$ , and is the left-handed component of $\bm{u}$ . For a flow that is at rest at infinity, $\nabla\phi$ is a constant vector field and can be eliminated by choosing an appropriate inertial frame. The vorticity field can be decomposed similarly as

\bm{\omega}(x)=\bm{\omega}^{R}(x)+\bm{\omega}^{L}(x),

(68)

with $\bm{\omega}^{R}=\nabla\times\bm{u}^{R}$ and $\bm{\omega}^{L}=\nabla\times\bm{u}^{L}$ .

Then, the helicity can be decomposed as

H=H^{R}+H^{L}=\int_{V}\bm{\omega}^{R}\cdot\bm{u}^{R}\mathrm{~{}d}V+\int_{V}\bm{\omega}^{L}\cdot\bm{u}^{L}\mathrm{~{}d}V.

(69)

The cross-terms $\bm{\omega}^{L}\cdot\bm{u}^{R}$ and $\bm{\omega}^{R}\cdot\bm{u}^{L}$ vanish due to the orthogonality of the eigenfunctions of the curl operator, and $H^{R}$ and $H^{L}$ are positive and negative definite, respectively.

IV Evolution from woven turbulence

We study the evolution of an initial woven turbulent field to examine the dynamics of woven turbulence. We find that woven turbulence gradually evolves into the regular HIT in a decaying process. This evolution is quantified by the turbulent kinetic energy, mean dissipation rate, and joint PDF of the second and third invariants ( $Q$ and $R$ ) of the velocity-gradient tensor.

IV.1 DNS with the initial woven turbulence

We first simulated a quantum vortex tangle with the cutoff scale of $0.005$ cm (Movie S2 and Materials and Methods, Section C), and chose the vortex tangle at ending time $t=10$ sec as the skeleton of woven turbulence. The initial vorticity field is constructed with core size $\sigma=0.03(1+0.5\sin s)$ varying along the arc length parameter $s$ (Fig. S8, upper left). This field is equivalent to the woven turbulent field in Fig. 1C in a smaller subdomain as $(1/2)^{3}$ of the original one, which facilitates further DNS and analysis.

We solve the 3D incompressible Navier–Stokes equations in the vorticity–velocity form using the pseudo-spectral method [39, 40] in a periodic box of size $L=2\pi$ on uniform grid points $N^{3}=512^{3}$ . The numerical solver removes aliasing errors using the two-third truncation method with the maximum wavenumber $k_{\max}\approx N/3$ . The time integration is treated by the explicit second-order Runge–Kutta scheme in physical space, with the adaptive time step ensuring the small enough Courant–Friedrichs–Lewy number for numerical stability and accuracy. In this DNS, we set the kinematic viscosity to $\nu=7.2\times 10^{-4}$ to match the effective viscosity calculated from woven turbulence.

IV.2 Transition from woven turbulence to real turbulence

Figure S8 illustrates the evolution of vortices in decaying of woven turbulence. The entangled vortex tubes approach each other, collide and reconnect under their self-induced motion. This interaction produces a variety of small-scale structures that fill the gaps between the initial vortex tubes. With the intense vortex dynamics, the vortex tubes break down into smaller fragments and then dissipate by viscous effects.

In terms of the turbulent kinetic energy and mean dissipation rate in Fig. S9, the woven turbulence first evolves into HIT in a rapid transition. After the non-dimensionalized time $t^{*}=t/(\langle\sigma\rangle^{2}/\nu)\approx 1$ , the turbulent flow decays slowly towards a quiescent state, with decay laws of $K\sim t^{-10/7}$ and $\left\langle\varepsilon\right\rangle\sim t^{-17/7}$ as in decaying HIT. The shape of the energy spectra remains almost the same through the temporal evolution (Fig. S9C). This suggests that the initial woven turbulence resembles real turbulence.

We examine how the symmetric butterfly-shaped joint PDF of $R$ and $Q$ evolves into a classic asymmetric teardrop shape during the transition from woven turbulence to real turbulence. Figure S10 depicts the teardrop formation. When the vortex tubes contact, the $R$ - $Q$ PDF starts to migrate along the right branch of the Vieillefosse line. Around $t^{*}=0.8$ , the joint PDF recovers the classic teardrop shape. As illustrated in Figs. S10A and D, the major difference in vortex morphology is the emergence of small-scale threads among vortex tubes during the rapid transition, while the large-scale vortex skeleton has very minor changes. Thus, the symmetry breaking of the $R$ - $Q$ PDF seems to be due to the generation of small-scale vortices.

V Customization of turbulence

We describe how to customize a turbulent field by weaving different elemental vortices, which serves as a testbed for various turbulence studies. The turbulence statistics can be manipulated by changing the core size (Section V.1), internal structure (Section V.2), and cross-sectional shape (Section V.3) of vortex tubes, with the same quantum vortex skeleton.

V.1 Core size

The core size of the vortex tube affects the Reynolds number and the width of the inertial range (Fig. 4 D). Using the quantum skeleton with a cutoff scale of $\delta=0.01$ cm (Movie S2 and Materials and Methods, Section C), we generate different woven turbulent fields with constant core sizes ranging from $\sigma=0.015$ to $\sigma=0.09$ (Fig. S11), corresponding to the Reynolds number from $Re=69$ to 118.

V.2 Internal structure

The internal structure of vortices influences small-scale statistics of classical turbulence. The vortex lines inside the vortex tube may not be parallel to the vortex centerline and may coil within the tube [20, 40]. The internal twist of a vortex is closely related to its helicity [38, 40]. We can construct vorticity fields with highly twisted vortices by precisely controlling the local twist rate $\eta_{t}$ . In Fig. 4B, we created helical woven turbulence with a quantum skeleton cutoff scale of $\delta=0.01$ cm, and constant core size of $\sigma=0.03$ , and local twist rate of $\eta_{t}=50$ . The right-handed internal twist wave fills the flow field with positive helicity density and right-handed chirality.

V.3 Cross-sectional shape

With the same quantum skeleton, the cross-section of the elemental vortex can be set in various shapes. We generate the turbulent field woven by sheet-like vortices using Eqs. (22), (21), and (23) with sheet thickness of $\sigma_{1}=0.02$ and width of $R_{v}=\sigma_{2}=0.4$ (Fig. 4C). Compared with the classical or vortex-tube woven turbulence, the statistics of the vortex-sheet woven turbulence in Fig. S12 show the steeper scaling law of $k^{-7/3}$ and shorter inertial range, and remain the negative third-order velocity structure function.

VI Supplementary movies

Movie S1: Evolution of quantum vortex tangle of superfluid Helium II using the vortex filament method with cutoff scale $\delta=0.005$ cm.

Movie S2: Evolution of quantum vortex tangle of superfluid Helium II using the vortex filament method with cutoff scale $\delta=0.01$ cm.

References

[1] Sreenivasan, K. R. Turbulence and the tube. Nature 344, 192–193 (1990).
[2] Coulais, C., Teomy, E., de Reus, K., Shokef, Y. & van Hecke, M. Combinatorial design of textured mechanical metamaterials. Nature 535, 529–32 (2016).
[3] Ding, L. et al. Polymer semiconductors: Synthesis, processing, and applications. Chem. Rev. 123, 7421–7497 (2023).
[4] Henderson, R. & Unwin, P. N. Three-dimensional model of purple membrane obtained by electron microscopy. Nature 257, 28–32 (1975).
[5] Guido, N. J. et al. A bottom-up approach to gene regulation. Nature 439, 856–60 (2006).
[6] Cardesa, J. I., Vela-Martín, A. & Jiménez, J. The turbulent cascade in five dimensions. Science 357, 782–784 (2017).
[7] Küchemann, D. Report on the I.U.T.A.M. symposium on concentrated vortex motions in fluids. J. Fluid Mech. 21, 1–20 (1965).
[8] Moffatt, H. K., Kida, S. & Ohkitani, K. Stretched vortices – the sinews of turbulence; large-Reynolds-number asymptotics. J. Fluid Mech. 259, 241–264 (1994).
[9] Hussain, A. K. M. F. Coherent structures and turbulence. J. Fluid Mech. 173, 303–356 (1986).
[10] Adrian, R. J. Hairpin vortex organization in wall turbulence. Phys. Fluids 19, 041301 (2007).
[11] Synge, J. L. & C., L. C. On a statistical model of isotropic turbulence. Trans. R. Soc. Canada 37, 45–79 (1943).
[12] Townsend, A. A. On the fine-scale structure of turbulence. Proc. R. Soc. Lond. A 208, 534–542 (1951).
[13] Corrsin, S. Turbulent dissipation fluctuations. Phys. Fluids 5, 1301–1302 (1962).
[14] Lundgren, T. S. Strained spiral vortex model for turbulent fine structure. Phys. Fluids 25, 2193–2203 (1982).
[15] Pullin, D. I. & Saffman, P. G. Vortex dynamics in turbulence. Annu. Rev. Fluid Mech. 30, 31–51 (1998).
[16] She, Z.-S., Jackson, E. & Orszag, S. A. Intermittent vortex structures in homogeneous isotropic turbulence. Nature 344, 226–228 (1990).
[17] Jiménez, J., Wray, A. A., Saffman, P. G. & Rogallo, R. S. The structure of intense vorticity in isotropic turbulence. J. Fluid Mech. 255, 65–90 (1993).
[18] Ishihara, T., Kaneda, Y. & Hunt, J. C. R. Thin shear layers in high Reynolds number turbulence—DNS results. Flow Turbulence Combust. 91, 895–929 (2013).
[19] Matsuzawa, T., Mitchell, N. P., Perrard, S. & Irvine, W. T. M. Creation of an isolated turbulent blob fed by vortex rings. Nat. Phys. 19, 1193–1200 (2023).
[20] Xiong, S. & Yang, Y. Identifying the tangle of vortex tubes in homogeneous isotropic turbulence. J. Fluid Mech. 874, 952–978 (2019).
[21] Fung, J. C. H., Hunt, J. C. R., Malik, N. A. & Perkins, R. J. Kinematic simulation of homogeneous turbulence by unsteady random Fourier modes. J. Fluid Mech. 236, 281–318 (1992).
[22] Kerstein, A. R. A linear- eddy model of turbulent scalar transport and mixing. Combust. Sci. Technol. 60, 391–421 (1988).
[23] Juneja, A., Lathrop, D. P., Sreenivasan, K. R. & Stolovitzky, G. Synthetic turbulence. Phys. Rev. E 49, 5179–5194 (1994).
[24] Rosales, C. & Meneveau, C. A minimal multiscale Lagrangian map approach to synthesize non-Gaussian turbulent vector fields. Phys. Fluids 18 (2006).
[25] Barenghi, C. F., Skrbek, L. & Sreenivasan, K. R. Introduction to quantum turbulence. Proc. Natl. Acad. Sci. U. S. A. 111, 4647–52 (2014).
[26] Kivotides, D. & Leonard, A. Quantized turbulence physics. Phys. Rev. Lett. 90 (2003).
[27] Polanco, J. I., Muller, N. P. & Krstulovic, G. Vortex clustering, polarisation and circulation intermittency in classical and quantum turbulence. Nat. Commun. 12, 7090 (2021).
[28] Hanninen, R. & Baggaley, A. W. Vortex filament method as a tool for computational visualization of quantum turbulence. Proc. Natl. Acad. Sci. U. S. A. 111 Suppl 1, 4667–74 (2014).
[29] Bradley, D. I. et al. Fluctuations and correlations of pure quantum turbulence in Superfluid 3He-B. Phys. Rev. Lett. 101, 065302 (2008).
[30] Kolmogorov, A. N. The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. Proc. R. Soc. Lond. A 434, 9–13 (1991).
[31] Yao, J. & Hussain, F. Vortex reconnection and turbulence cascade. Annu. Rev. Fluid Mech. 54, 317–347 (2022).
[32] Su, H., Yang, Y. & Pope, S. B. Simple model for the bottleneck effect in isotropic turbulence based on Kolmogorov’s hypotheses. Phys. Rev. Fluids 8 (2023).
[33] Lohse, D. & Muller-Groeling, A. Bottleneck effects in turbulence: Scaling phenomena in r versus p space. Phys. Rev. Lett. 74, 1747–1750 (1995).
[34] Vinen, W. F. Decay of superfluid turbulence at a very low temperature: The radiation of sound from a Kelvin wave on a quantized vortex. Phys. Rev. B 64 (2001).
[35] Buaria, D., Pumir, A., Bodenschatz, E. & Yeung, P. K. Extreme velocity gradients in turbulent flows. New J. Phys. 21, 043004 (2019).
[36] Meneveau, C. Lagrangian dynamics and models of the velocity gradient tensor in turbulent flows. Annu. Rev. Fluid Mech. 43, 219–245 (2011).
[37] Moffatt, H. K. Extreme events in turbulent flow. J. Fluid Mech. 914, F1 (2021).
[38] Scheeler, M. W., van Rees, W. M., Kedia, H., Kleckner, D. & Irvine, W. T. M. Complete measurement of helicity and its dynamics in vortex tubes. Science 357, 487–491 (2017).
[39] Meng, Z., Shen, W. & Yang, Y. Evolution of dissipative fluid flows with imposed helicity conservation. J. Fluid Mech. 954, A36 (2023).
[40] Shen, W., Yao, J., Hussain, F. & Yang, Y. Role of internal structures within a vortex in helicity dynamics. J. Fluid Mech. 970, A26 (2023).
[41] Scheeler, M. W., Kleckner, D., Proment, D., Kindlmann, G. L. & Irvine, W. T. M. Helicity conservation by flow across scales in reconnecting vortex links and knots. Proc. Natl. Acad. Sci. U. S. A. 111, 15350–15355 (2014).
[42] Shen, W., Yao, J., Hussain, F. & Yang, Y. Topological transition and helicity conversion of vortex knots and links. J. Fluid Mech. 943, A41 (2022).
[43] Davidson, P. A. Turbulence, an introduction for scientists and engineers (Oxford University Press, 2015).
[44] Aharon-Steinberg, A. et al. Direct observation of vortices in an electron fluid. Nature 607, 74–80 (2022).
[45] Dong, C. et al. Reconnection-driven energy cascade in magnetohydrodynamic turbulence. Sci. Adv. 8, eabn7627 (2022).
[46] Wu, K. T. et al. Transition from turbulent to coherent flows in confined three-dimensional active fluids. Science 355 (2017).
[47] Duclos, G. et al. Topological structure and dynamics of three-dimensional active nematics. Science 367, 1120–1124 (2020).
[48] Novati, G., de Laroussilhe, H. L. & Koumoutsakos, P. Automating turbulence modelling by multi-agent reinforcement learning. Nat. Mach. Intell. 3, 87–96 (2021).
[49] Baggaley, A. W. & Laizet, S. Vortex line density in counterflowing He II with laminar and turbulent normal fluid velocity profiles. Phys. Fluids 25 (2013).
[50] Yang, Y. & Pullin, D. I. On Lagrangian and vortex-surface fields for flows with Taylor–Green and Kida–Pelz initial conditions. J. Fluid Mech. 661, 446–481 (2010).
[51] Araki, T., Tsubota, M. & Nemirovskii, S. K. Energy spectrum of superfluid turbulence with no normal-fluid component. Phys. Rev. Lett. 89, 145301 (2002).
[52] Chumakov, S. G. A priori study of subgrid-scale flux of a passive scalar in isotropic homogeneous turbulence. Phys. Rev. E 78, 036313 (2008).
[53] Machiels, L. Predictability of small-scale motion in isotropic fluid turbulence. Phys. Rev. Lett. 79, 3411–3414 (1997).
[54] Lee, S., Lele, S. K. & Moin, P. Simulation of spatially evolving turbulence and the applicability of Taylor’s hypothesis in compressible flow. Phys. Fluids A 4, 1521–1530 (1992).
[55] Arneodo, A., Bacry, E. & Muzy, J. F. Random cascades on wavelet dyadic trees. J. Math. Phys. 39, 4142–4164 (1998).
[56] Eggers, J. & Grossmann, S. Effect of dissipation fluctuations on anomalous velocity scaling in turbulence. Phys. Rev. A 45, 2360–2369 (1992).
[57] Benzi, R. et al. Extended self-similarity in turbulent flows. Phys. Rev. E 48, R29–R32 (1993).
[58] Biferale, L., Boffetta, G., Celani, A., Crisanti, A. & Vulpiani, A. Mimicking a turbulent signal: Sequential multiaffine processes. Phys. Rev. E 57, R6261–R6264 (1998).
[59] Scotti, A. & Meneveau, C. A fractal model for large eddy simulation of turbulent flow. Physica D 127, 198–232 (1999).
[60] Rosales, C. & Meneveau, C. Anomalous scaling and intermittency in three-dimensional synthetic turbulence. Phys. Rev. E 78, 016313 (2008).
[61] Xiong, S. & Yang, Y. Construction of knotted vortex tubes with the writhe-dependent helicity. Phys. Fluids 31 (2019).
[62] Xiong, S. & Yang, Y. Effects of twist on the evolution of knotted magnetic flux tubes. J. Fluid Mech. 895, A28 (2020).
[63] Chui, A. Y. K. & Moffatt, H. K. The energy and helicity of knotted magnetic flux tubes. Proc. R. Soc. Lond. A 451, 609–629 (1995).
[64] She, Z.-S. & Leveque, E. Universal scaling laws in fully developed turbulence. Phys. Rev. Lett. 72, 336–339 (1994).
[65] Berger, M. A. & Field, G. B. The topological properties of magnetic helicity. J. Fluid Mech. 147, 133–148 (1984).
[66] Waleffe, F. The nature of triad interactions in homogeneous turbulence. Phys. Fluids A 4, 350–363 (1992).