2.1. Vortex filament representation

The vortex filament representation characterises vorticity as discrete line singularities with circulation confined to the filament axis, effectively capturing the fundamental features of turbulent flows, including vortex generation, interaction and dissipation. This approach provides a computationally efficient framework for analysing vortical dynamics and energy transfer, particularly in high Reynolds number regimes, where turbulence is governed by the intricate interactions of vortex filaments.

In this representation, vorticity $ \boldsymbol{\omega } = \boldsymbol{\nabla } \times \boldsymbol{u}$ , where $ \boldsymbol{u}$ is the velocity field, is concentrated along one-dimensional curves, known as vortex filaments. The physical domain is an open region $ \varOmega \subset \mathbb{R}^3$ , and the set of $ m$ closed, time-dependent vortex filaments is described by(see Ishida et al. Reference Ishida, Wojtan and Chern2022)

(2.1) \begin{equation} \boldsymbol{\gamma }(:, t) : \bigsqcup _{i=1}^m S^1 \to \varOmega , \end{equation}

where $ \bigsqcup$ denotes a disjoint union, $ S^1$ represents the topological circle, and $ \boldsymbol{\gamma }$ maps each vortex filament into the domain $ \varOmega$ . This framework provides a means to model the essential dynamics of turbulence. Specifically, Shen et al. (Reference Shen, Yao and Yang2024) use quantum vortex tangles to generate synthetic turbulence by combining vortex filaments, effectively capturing the key characteristics of turbulent flows. The ability to vary vortex filament thickness further facilitates the investigation of energy cascades and extreme events within turbulence dynamics.

Accurate representation of vortex structures is crucial for understanding vortex morphology (Xiong & Yang Reference Xiong and Yang2019b ). Vortex filaments, as the limiting case of vortex surfaces (Xiong & Yang Reference Xiong and Yang2017), form a network of elongated, intertwined structures. In contrast, iso-surfaces of the vorticity magnitude $ |\boldsymbol{\omega }|$ typically exhibit shorter, tube-like shapes, often referred to as ‘vortex worms’ (Jiménez et al. Reference Jiménez, Wray, Saffman and Rogallo1993). This distinction becomes particularly pronounced during the energy cascade, where $ |\boldsymbol{\omega }|$ iso-surfaces fragment into smaller structures. However, such fragmentation provides only a partial view of vortex dynamics. Key processes in the energy cascade, including stretching, twisting and reconnection, are more accurately represented by the continuous deformation of vortex filaments.

The motion of vortex filaments is governed by the kinematic equation

(2.2) \begin{equation} \frac {\partial \boldsymbol \gamma }{\partial t}(s, t) = \boldsymbol{u}_\gamma (s, t), \quad s \in \bigsqcup S^1, \; t \in \mathbb{R}, \end{equation}

where $ \boldsymbol{u}_\gamma (s, t)$ is the velocity of the filament. The BS law describes the self-induced motion of the vortex,

(2.3) \begin{equation} \boldsymbol{u}^{BS}_\gamma (s) = \frac {\Gamma }{4 } \int _\gamma \frac {\boldsymbol \gamma '(s') \times [\boldsymbol \gamma (s) - \boldsymbol \gamma (s')]}{|\boldsymbol \gamma (s) - \boldsymbol \gamma (s')|^3} \, \mathrm{d}s', \end{equation}

where $ \Gamma$ is the vortex strength, and $ \mathrm{d}s'$ is the arclength element. To regularise the singularity in the BS law, the RM model (Saffman Reference Saffman1995) introduces a core size $ a$ , modifying the velocity to

(2.4) \begin{equation} \boldsymbol{u}^{RM}_\gamma (s) = \frac {\Gamma }{4 } \int _\gamma \frac {\boldsymbol \gamma '(s') \times [\boldsymbol \gamma (s) - \boldsymbol \gamma (s')]}{(\sqrt {[{\rm e}^{-3/2}\,a^2 + |\boldsymbol \gamma (s) - \boldsymbol \gamma (s')|^2}]^3} \, \mathrm{d}s', \end{equation}

which reduces to the BS law in the limit as $ a \to 0$ . These equations describe vortex transport and interactions, providing a link between local filament dynamics and the large-scale behaviour of turbulence.

Despite its advantages, the vortex filament representation faces challenges in handling high filament densities and complex topological changes. These issues arise from the simplicity of the mathematical formulation and the limitations of existing numerical techniques. For example, reconnection algorithms integral to capturing filament dynamics often rely on subjective criteria, lacking the precision and rigour needed for high-fidelity simulations.

2.2. Implicit vortex filament representation

In explicit methods, vortex filaments are represented by discrete points or lines, requiring continuous updates at each time step as their position and properties are recalculated. In contrast, implicit vortex filaments are defined as zeroiso-surfaces of a scalar field, with their evolution inherently linked to changes in the scalar field itself. As the scalar field evolves, the zeroiso-surface representing the filament adjusts automatically, eliminating the need for explicit updates of each individual point. The computational complexity inherent in explicit methods naturally suggests the potential of quantum computing, particularly given that a complex scalar field can be interpreted as a quantum state, aligning seamlessly with the foundational principles of quantum computation.

Before delving into our proposed quantum implicit vortex filament theoretical framework, we present a representative case to demonstrate its application in capturing the complex dynamics of chaotic turbulence. Figure 1 shows vortex filament extraction from a three-dimensional turbulent flow, where the input is a randomly generated velocity field. Figure 1(a) illustrates the vortex filaments, coloured by helicity, with a yellow box highlighting a region of interest for further analysis. Figure 1(b) zooms in on this highlighted region, providing a detailed view of the filaments coloured by helicity, capturing the intricate entanglement of the vortices. Figure 1(c) focuses on the iso-surfaces of $ \text{Im}(\psi ) = 0$ within the same region, with the filaments coloured by the real part of $ \psi$ , offering a further visualisation of the complex interactions between vortices in the flow.

Let $ \psi$ be a complex scalar field in $ \varOmega$ , with $ \partial \Sigma$ being a closed curve. The winding number of $\psi$ around $\partial \Sigma$ is a topological invariant that quantifies the total phase rotation of $\psi$ along the curve. It is defined as

(2.5) \begin{equation} n_{w} = \frac {1}{2 } \int _{\partial \Sigma } \boldsymbol{\nabla } (\arg \psi ) \boldsymbol{\cdot} \mathrm{d}\boldsymbol{l}, \end{equation}

where $\arg \psi$ represents the phase of the complex scalar field, such that $ \psi = |\psi |\,{\rm e}^{\mathrm{i}\theta }$ and $ \theta = \arg \psi$ , and $\boldsymbol{\nabla } (\arg \psi )$ is the phase gradient, capturing its spatial variation along the curve. The term $\mathrm{d}\boldsymbol{l}$ is an infinitesimal vector element along the curve $\partial \Sigma$ , specifying the direction of integration. This integral measures the total phase change normalised by $ 2$ , counting the number of complete rotations around the curve. Since the phase of $ \psi$ returns to its original value after one full loop unless there is a singularity inside the curve, a non-zero winding number occurs only where the zero-level set of $ \psi$ exists, indicating that phase rotation does not return to zero. A non-zero winding number thus indicates the presence of one or more vortex filaments enclosed by $ \partial \Sigma$ , with its magnitude reflecting the total circulation strength.

Figure 1.

Vortex filaments and the complex scalar field $\psi$ in a three-dimensional turbulent flow.Here, $h_{{helicity}}$ denotes the helicity density. (a) Vortex filaments extracted from a three-dimensional turbulent flow, coloured by helicity. (b,c) Detailed views of the filaments and iso-surfaces of $\mathrm{Im}(\psi ) = 0$ within the highlighted region (bright yellow box) in (a), coloured by helicity and $\mathrm{Re}(\psi )$ , respectively. The detailed views illustrate the complex entanglement of vortices in turbulent flows.

Alternatively, the winding number can also be expressed in a mathematically equivalent form that directly incorporates the complex scalar field $ \psi$ and its gradient, given by

(2.6) \begin{equation} n_{w} = \frac {1}{2 } \int _{\partial \Sigma } \frac {\langle \boldsymbol \nabla \psi , \mathrm{i} \psi \rangle }{|\psi |^2} \boldsymbol{\cdot} \mathrm{d}\boldsymbol l, \end{equation}

where $\langle \boldsymbol{\nabla } \psi , \mathrm{i} \psi \rangle$ represents the inner product (or projection) of the gradient of $\psi$ onto the imaginary axis, scaled by $\psi$ itself. This isolates the change in the phase of $\psi$ , discarding changes in magnitude. Dividing by $|\psi |^2$ normalises $\psi$ , ensuringthat only the phase rotation contributes.

The winding number, as a topological invariant, provides a useful tool for characterising vortex fields in wave-function-based representations. It serves as an analogue to the circulation of a velocity field in classical fluid mechanics, establishing a connection between topological properties and dynamical behaviour. For a closed curve $ \partial \Sigma$ , the circulation of the velocity field $ \boldsymbol{u}$ is defined as

(2.7) \begin{equation} \Gamma _u = \int _{\partial \Sigma } \boldsymbol u \boldsymbol{\cdot} \mathrm{d}\boldsymbol l, \end{equation}

where $ \boldsymbol{u}$ represents the velocity field. By Stokes’ theorem, this circulation corresponds to the flux of vorticity through any surface $ \Sigma$ that is bounded by the curve $ \partial \Sigma$ .

Figure 2.

Schematic of the process for computing the complex scalar field from a given velocity field. The optimisation problem is formulated using an arbitrary closed curve in $\varOmega$ , where the winding number and circulation are related through (2.8). The resulting complex scalar field approximates the velocity field, with its zero-level set corresponding to vortex filaments of strength $h$ .

From (2.6) and (2.7), the winding number of the complex scalar field and the circulation of the velocity field are both defined for any closed curve, thereby linking these two quantities. Following the approach of Weißmann et al. (Reference Weißmann, Pinkall and Schröder2014), we approximate the winding number and circulation as proportional, with the relation

(2.8) \begin{equation} \Gamma _{u} \approx h n_{w}, \end{equation}

where $ h$ is a proportionality constant that connects the two.

Equation (2.8)is not directly suitable for numerical computation, so we reformulate it as an optimisation problem to facilitate implementation. We substitute the expressions for $ n_w$ and $ \Gamma _u$ from (2.6) and (2.7) into (2.8), respectively. Then we rearrange the terms, moving all non-zero terms to the left-hand side, yielding the equation

(2.9) \begin{equation} \frac {1}{2 } \int _{\partial \Sigma } \frac {\langle \boldsymbol \nabla \psi - \mathrm{i} \hat {\boldsymbol u} \psi , \mathrm{i} \psi \rangle }{|\psi |^2} \boldsymbol{\cdot} \mathrm{d}\boldsymbol l \approx 0, \end{equation}

where $ \hat {\boldsymbol u} = 2 \boldsymbol u / h$ . To make this equation hold, we pose the optimisation problem

(2.10) \begin{equation} \min _{\psi } \tfrac {1}{2} \|\boldsymbol \nabla \psi - \mathrm{i} \hat {\boldsymbol u} \psi \|^2 \quad \text{s.t.} \ \|\psi \|^2 = 1. \end{equation}

Equation (2.10) represents a linear functional problem with a normalisation constraint, which can be solved numerically on a discrete grid.

Figure 2 illustrates the procedure for computing the complex scalar field from a given velocity field. The process begins with a closed curve in $\varOmega$ , where the winding number, defined in (2.6), is derived from the complex scalar field, and the circulation, defined in (2.7), is computed from the velocity field. Although these quantities are related conceptually, they are not typically directly proportional. To address this, the problem is reformulated as an optimisation problem, as described in (2.10), which aims to align the winding number and circulation as closely as possible. By solving this optimisation problem for a given velocity field, a complex scalar function is obtained, which serves as an approximate representation of the velocity field. The zero-level set of this scalar field is interpreted as vortex filaments with strength $h$ , and the reconstructed vorticity flux across any surface is approximated to match the vorticity flux of the original velocity field.

Once the complex scalar field is constructed, we define the implicit representation of vortex filaments as follows.

We consider a two-dimensional example to illustrate the representation of the flow field by the complex scalar field. Figure 3 shows the distribution of the magnitude and phase of the complex scalar field, which represents a velocity field induced by four vortex elements, as described by (2.8). In figure 3(b), a closed curve $ \partial \Sigma$ is depicted in $ \mathbb{R}^2$ , with the red and blue circular arrows indicating regions where the winding number around the zero-level set of $ \psi$ is $\pm 1$ , and the black circle highlighting areas where the winding number is zero.

Figure 3.

The distribution of (a) the magnitude and (b) the phase of a two-dimensional complex scalar field, representing a velocity field induced by four vortex elements as described by (2.8), with dark blue diamonds indicating the intersections of the real and imaginary parts of the complex scalar field. In (b), the red and blue circular arrows indicate regions where the winding number around the zero-level set of the complex scalar field is $\pm 1$ , while the black circle highlights areas with winding number zero.

If the zero-level set of $ \psi$ does not intersect the surface $ \Sigma$ enclosed by $ \partial \Sigma$ , then the winding number, computed using (2.6), is zero, implying the absence of net circulation or vorticity flux through the surface. Conversely, a non-zero vorticity flux, as given by (2.8), indicates that the zero-level set of $ \psi$ intersects the surface $ \Sigma$ . In this case, the flow field is characterised by the presence of vortex filaments, which correspond to the zero-level set of $ \psi$ and represent coherent structures within the flow, where the complex scalar field $ \psi$ undergoes a significant phase change.

On the other hand, interpreting the winding number as the number of intersections between the vortex filaments and the surface $ \Sigma$ suggests that these filaments possess a uniform strength. This strength is represented by the constant $ h$ , which links the circulation of the velocity field to the winding number. As a result, each vortex filament contributes equally to the total circulation, providing a simplified yet insightful characterisation of the flow’s topological structure.

For a time-dependent, complex-valued level set function $\psi$ , assuming that the function satisfies the necessary regularity conditions, its time evolution near the zero-level set is governed by the equation

(2.11) \begin{equation} \frac {\partial \psi }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla } \psi = 0, \end{equation}

where $\boldsymbol{u}$ is a velocity field that governs the motion of the zero-level set of $\psi$ . This equation describes the advection of the implicit vortex filament, specifying how it is transported by the flow of the velocity field $\boldsymbol{u}$ . In particular, the evolution of the zero-level set, where $\psi = 0$ , follows the flow determined by $\boldsymbol{u}$ , thereby capturing the dynamics of the vortex filament. This leads to the following theorem.

Theorem 1 (Relation between curve dynamics and implicit representation). Let $\psi$ be a time-dependent complex level set function that satisfies the regularity assumptions, namely, there exists a neighbourhood $\mathcal{U}$ of the zero-level set of $\psi$ such that the gradient matrix of $\psi$ is of full rank in $\mathcal{U}$ . Denote by $\boldsymbol \gamma$ the set of vortex filaments advected by the velocity field $\boldsymbol{u}_{\gamma }$ , and let $\psi$ represent the complex scalar field associated with these vortex filaments, as defined in definition 1 . The evolution of the zero-level set of $\psi$ is governed by ( 2.11 ), where the velocity field $\boldsymbol{u}$ and the relative velocity $\boldsymbol{u}_{\gamma }$ are aligned in the direction normal to the filaments. This is expressed by the condition

(2.12) \begin{equation} (\boldsymbol{u} - \boldsymbol{u}_{\gamma }) \times \frac {\partial }{\partial s} \boldsymbol \gamma (s,:) = \boldsymbol{0}. \end{equation}

The velocity components that primarily influence the evolution of the curves are those defined on the zero-level set. Additionally, the positions of the curves remain unaffected when they are translated along their tangent direction, analogous to rotating a circle about its axis of symmetry. Thus it is sufficient to constrain the normal and binormal components of the velocity field.

3.1. Discretisation of the optimisation formulation

We solve the optimisationprobelm (2.10) on a staggered grid, where the complex scalar field is discretised at the cell centres and the velocity components are staggered at the face centres. The staggered grid arrangement, first introduced by Harlow & Welch (Reference Harlow and Welch1965) in the marker-and-cell method, places velocity components at cell faces while locating scalar quantities at cell centres. This geometric configuration provides dual numerical advantages: it inherently maintains discrete conservation laws for improved solution stability, while simultaneously enabling precise calculations of vorticity flux and velocity circulation that are crucial for vortex dynamics analysis.

This discretisation method follows the approach outlined by Xiong et al. (Reference Xiong, Wang, Wang and Zhu2022). As shown in figure 5, for two adjacent cells, labelled $i$ and $j$ , the complex values $ \psi _i$ and $ \psi _j$ are placed at the respective cell centres, while the shared face between the cells contains the velocity component $ u_{ij}$ . On the $ij$ th face, the term $ \boldsymbol \nabla \psi - \mathrm{i} \hat {\boldsymbol u} \psi$ is discretised as

(3.1) \begin{equation} (\psi _{j} - \psi _{i}) - \tfrac {1}{2} \textrm {i} \hat \eta _{ij} ( \psi _{i} + \psi _{j})+\mathcal{O}(\Delta x^2), \end{equation}

where $ \hat \eta _{ij} = 2 u_{ij}\, \Delta x / h$ , and $ \Delta x$ is the grid cell spacing.

Figure 5.

Discretisation of the optimisationproblem (2.10) on a staggered grid, where the complex scalar field $ \psi$ is defined at the cell centres, and the velocity $ \boldsymbol u$ is defined at the face centres. For adjacent cells $ i$ and $ j$ , the velocity component $ u_{ij}$ at the shared face leads to the discretisation given by (3.1). This formulation allows the optimisation problem to be recast as the minimisation of (3.4) over the sets of cells $ \mathcal{V}$ and faces $ \mathcal{E}$ .

We apply a Taylor expansion to obtain

(3.2) \begin{equation} 1 \pm \frac {1}{2} \textrm {i} \hat \eta _{ij} = \exp \left (\pm \frac {\textrm {i}}{2} \hat \eta _{ij}\right ) + \mathcal{O}(\Delta x^2). \end{equation}

This enables the symmetric representation of (3.1) as

(3.3) \begin{equation} \overline {r}_{ij}\psi _j - r_{ij}\psi _i + \mathcal{O}(\Delta x^2), \end{equation}

where $ r_{ij} = \exp (\textrm {i} \hat {\eta }_{ij}/2)$ , and $ \overline {r}_{ij}$ denotes the complex conjugate of $ r_{ij}$ . The optimisation problem can then be reformulated as

(3.4) \begin{equation} \min _{\psi _k,k\in \mathcal{V}} \tfrac {1}{2}\sum _{ij\in \mathcal{E}} |\overline {r}_{ij}\psi _j - r_{ij}\psi _i|^2\quad \text{s.t.}\ \sum _{k\in \mathcal{V}}|\psi _k|^2 = 1, \end{equation}

where $ \mathcal{V}$ is the set of all cell indices, and $ \mathcal{E}$ is the set of all face indices.

The term $ |\overline {r}_{ij}\psi _j - r_{ij}\psi _i|^2$ can be expanded as

(3.5) \begin{equation} \left ( \overline { \overline {r}_{ij}\psi _j - r_{ij}\psi _i} \right ) \left (\overline {r}_{ij}\psi _j - r_{ij}\psi _i\right )=|\psi _i|^2+|\psi _j|^2-(r_{ij})^2\,\overline {\psi }_j\psi _i-(\overline {r}_{ij})^2\,\overline {\psi }_i\psi _j. \end{equation}

Using this expansion, the optimisationproblem (3.4) can be reformulated as a quadratic form of complex variables under a normalisation constraint

(3.6) \begin{equation} \min _{\|\psi \|^2=1} \psi ^\dagger H \psi, \end{equation}

where the elements of the Hermitian matrix $ H$ are given by

(3.7) \begin{equation} H_{ij} = \begin{cases} - (\overline {r}_{ij})^2 & \text{if the } j\text{th cell is a neighbour of the } i\text{th cell}, \\ d_i & \text{if } i = j, \\ 0 & \text{otherwise}. \end{cases} \end{equation}

Here, $ d_i$ denotes the number of neighbours of the $ i$ th cell and corresponds to the term $ |\psi _i|^2$ . The term $ -(\overline {r}_{ij})^2$ represents the cross terms, specifically $ (r_{ij})^2\, \overline {\psi }_j \psi _i$ and $ (\overline {r}_{ij})^2\, \overline {\psi }_i \psi _j$ . When $ i$ and $ j$ are not neighbours, $ H_{ij} = 0$ . Since $ r_{ij} = \overline {r}_{ji}$ , the matrix $ H$ is Hermitian, transforming the filament extraction problem into a Hermitian eigenvalue problem.

We remark that in (2.10), $ \psi$ represents a continuous complex scalar field defined over the fluid domain. In the discretised formulation (3.6), $ \psi$ is treated as a vector, with each component corresponding to the complex value of the scalar field at the centre of a grid cell. This transition to a discrete representation is necessary for numerical computation.

3.2. Quantum algorithm

We represent the discrete scalar field $\psi$ in (3.6) as the quantum state $|\psi \rangle$ , and interpret the Hermitian matrix as the Hamiltonian operator. Thus (3.6) is reformulated as a quantum eigenvalue problem. To solve this, we employ a VQE.

The quantum circuit $ U(\boldsymbol{\theta })$ is parametrised by the vector $\boldsymbol{\theta }$ , representing the circuit parameters. The initial state is prepared as the zero state, $ |\psi _0\rangle = |0\rangle$ . The structure of the circuit is shown in figure 6, where it consists of $N_{{qubit}}$ qubits and $N_d$ basic blocks. The number of qubits $N_{{qubit}}$ depends on the discretisation of the computational domain. In our case, itis defined as

(3.8) \begin{equation} N_{{qubit}} \sim \log _2 N_{{cell}}. \end{equation}

The blocks share identical structure but with different parameter values. Each block contains four layers of rotational gates, defined as

(3.9) \begin{equation} \begin{aligned} \textrm {R}_y(\theta ) = \left [\begin{matrix}\cos \left (\theta /2\right ) & -\sin \left (\theta /2\right ) \\ \sin \left (\theta /2\right ) & \cos \left (\theta /2\right ) \end{matrix}\right ],\quad \textrm {R}_z(\theta ) = \left [\begin{matrix}\exp \left (-\mathrm{i} \theta /2\right ) & 0 \\0 & \exp \left (\mathrm{i} \theta /2\right ) \end{matrix}\right ], \end{aligned} \end{equation}

followed by a layer of CNOT gates connecting adjacent qubits after every two layers of rotational gates. While rotational gates play a crucial role in quantum circuits (Xiao et al. Reference Xiao, Yang, Shu, Chew, Khoo, Cui and Liu2024), we restrict the use of rotational gates to $\textrm {R}_y$ and $\textrm {R}_z$ for two reasons. First,an $\textrm {R}_x$ gate can be expressed by $\textrm {R}_y$ and $\textrm {R}_z$ , which meansthat including all three does not fundamentally expand the expressibility of the circuit. Second, our circuit structure is inspired by the hardware-efficient ansatz, where minimising the variety of quantum gates enhances practical implementability on near-term quantum devices (Kandala et al. Reference Kandala, Mezzacapo, Temme, Takita, Brink, Chow and Gambetta2017; Huang et al. Reference Huang, Li, Hou, Wu, Yung, Bayat and Wang2022). To expand the search space, the block is further repeated $ N_d$ times. The effect of $ N_d$ on circuit performance is discussed in the next subsection.

Figure 6.

Schematic of the parametrised quantum circuit $ U(\theta )$ for a system of $N_{{qubit}}$ qubits, with $N_{{qubit}} = 5$ used as an illustrative example. The circuit consists of $ N_d$ blocks, each containing four layers of rotational gates with distinct parameters. After every two layers of rotational gates, a layer of CNOT gates is applied to connect adjacent qubits.

The parametrised circuit $ U(\boldsymbol{\theta })$ is applied to the initial state $ |\psi _0\rangle$ , leading to the evolution of the system to the target state $ |\psi \rangle = U(\boldsymbol{\theta })\, |\psi _0\rangle$ . The expectation value of the Hamiltonian operator $ H$ is then computed as the objective function for optimisation $\langle \psi | H | \psi \rangle$ , where the parameters $ \theta$ are adjusted to minimise this function. The optimisation procedure yields the solution

(3.10) \begin{equation} \boldsymbol{\theta }_{{opt}} = \mathop {\arg \min } \limits _{\theta } \langle \psi _0\, | U(\boldsymbol{\theta })^\dagger HU(\boldsymbol{\theta }) |\, \psi _0 \rangle , \end{equation}

which results in the optimal state $ |\psi _{{opt}} \rangle = U(\boldsymbol{\theta }_{{opt}})\, | \psi _0 \rangle$ .

When measuring the objective function $\langle \psi | H | \psi \rangle$ on a quantum computer, the Hamiltonian $ H$ is typically decomposed into a sum of Pauli operators as

(3.11) \begin{equation} H = \sum _{i} c_i P_i, \end{equation}

where $ P_i$ are tensor products of Pauli matrices, and $ c_i$ are real coefficients. The objective function then becomes

(3.12) \begin{equation} \langle \psi\, | H |\, \psi \rangle = \sum _i c_i\, \langle \psi\, | P_i |\, \psi \rangle . \end{equation}

This decomposition allows for the Hamiltonian measurement to be simplified to the evaluation of individual Pauli terms (McClean et al. Reference McClean, Romero, Babbush and Aspuru-Guzik2016).

As the system size increases, the number of Pauli operators grows, leading to a significant increase in measurement costs. To address this, the ‘Pauli operator cut-off’ method is adopted, which reduces the number of measurements while maintaining an appropriate balance between efficiency and accuracy. Further details of this approach are provided in Appendix A.

3.3. Extraction of vortex filaments

The computed quantum state obtained using VQE serves as the basis for identifying the locations of vortex filaments. To determine whether a filament intersects a given surface, the winding number around a closed curve or surface element is calculated by evaluating (2.6). A non-zero winding number indicates the presence of a vortex filament, whereas a zero value suggests no intersection. If a filament is detected, then interpolation is employed to refine the location, accurately pinpointing the crossing point within the discrete grid.

On a discrete grid, we assume that the centres of four adjacent cells, denoted $ jklm$ , form a square with side length $ \Delta x$ . The winding number over this square is approximated byWeißmann et al. (Reference Weißmann, Pinkall and Schröder2014) as

(3.13) \begin{equation} n_{w,jklm} = \frac {1}{2 }\left [\arg \left (\frac {\psi _k}{\psi _j} \right ) + \arg \left (\frac {\psi _l}{\psi _k} \right ) + \arg \left (\frac {\psi _m}{\psi _l} \right ) + \arg \left (\frac {\psi _j}{\psi _m} \right ) \right ], \end{equation}

where $\psi _i$ (with $i = j, k, l, m$ ) denotes the discretised complex scalar field at the grid points corresponding to $j$ , $k$ , $l$ and $m$ . If $ n_{w,jklm} = 0$ , then no vortex filament is considered to intersect the square formed by the four adjacent cells $ jklm$ . If $ n_{w,jklm} \neq 0$ , then it indicates the presence of a vortex filament intersecting the square. In such cases, the exact crossing point (or zero point) is determined using linear interpolation, given by

(3.14) \begin{equation} (1 - v)(1 - u) \psi _j + u(1 - v) \psi _k + (1 - u) v \psi _m + u v \psi _l = 0, \end{equation}

where $ (u, v)$ are local coordinates within the unit square $(0, 1) \times (0, 1)$ , representing the relative position of the intersection of the vortex filament within the square. Solving this equation yields the precise location of the intersection.

Although theoretically effective, quantum computations face practical challenges due to noise and errors. Quantum decoherence, caused by interactions with the environment, can drive a quantum state towards relaxation to the ground state $|0\rangle$ , and quantum gates may introduce gate errors during operations (Bharadwaj & Sreenivasan Reference Bharadwaj and Sreenivasan2025). Other sources, such as limited grid resolution, measurement inaccuracies, and incomplete convergence from a finite number of iterations in the quantum algorithm, also introduce noise into the computed quantum state. This noise complicates the reliable identification of vortex filaments.

To mitigate these challenges, we propose a data-driven machine learning approach for extracting vortex filaments from noisy quantum states. Our neural network comprises two branches: one branch predicts whether a vortex filament crosses the grid surface, analogous to the discretised winding number discussed in (3.13), while the other branch forecasts the local coordinates of the crossing location. The two branches are trained collaboratively, with the loss function incorporating factors such as flux conservation and symmetry, thereby enhancing precision and preserving the physical properties of the vortex filaments. The circuit depth in the training and testing dataset is intentionally reduced to simulate scenarios with noise and error. This method augments the classical interpolation technique, improving resilience to noise and enhancing detection accuracy. By learning patterns in noisy data, the machine learning model can identify filament locations even when classical methods are compromised by noise artefacts. This data-driven approach leverages the robustness of classical algorithms while exploiting the adaptability of machine learning, ensuring accurate and efficient vortex filament extraction in complex, noisy environments. Further details of the machine learning method and its implementation are provided in Appendix B.

The effectiveness of the proposed method is demonstrated in figure 7, where vortex filaments reconstructed using the network approach exhibit near-perfect alignment with the true filament locations. In contrast, filaments extracted using the classical method show significant deviations from the ground truth. This comparison underscores the superior accuracy of the machine-learning-based method, particularly in the presence of noise.

Figure 7.

(a) The vortex filaments reconstructed using the proposed NN method exhibit near-perfect agreement with the target filament locations, where NN stands for neural network. (b) In contrast, the extracted vortex filaments using the classical method show significant deviation from the ground truth.

4.1. Parameter configuration

We analyse the influence of different parameters, and identify the most suitable configurations. One challenge lies in the direct comparison of the extracted vortex filaments with the target ones, especially in complex flow fields where defining positional discrepancies is non-trivial. Alternatively, considering our goal of establishing a connection between the vortex filaments and the velocity field, we reconstruct the velocity field using the RM model referredto in (2.4). We then quantify the discrepancy between the recovered velocity field and the original velocity field using the mean squared error (MSE) metric:

(4.1) \begin{equation} \mathrm{MSE} = \frac {1}{N_{{cell}}} \sum _{{cells}} \left |\boldsymbol{u} - \boldsymbol{u}_{{target}}\right |^2, \end{equation}

where $N_{{cell}}$ denotes the total number of grid cells in the flow field. This metric provides a quantitative assessment of the accuracy of the reconstructed flow field, offering insights into the algorithm’s effectiveness in capturing the dynamics of the system.

We inventory the parameters that require analysis. Figure 6 illustrates the structure of the quantum circuit, where a single block is repeated $N_d$ times. Additionally, the VQE requires iterative optimisation, and the number of iterations $N_i$ is another parameter of interest. Since both parameters are closely related to the optimisation process, we refer to them collectively as optimisation parameters. Meanwhile, the vortex strength parameter $h$ , introduced in (2.8), and the vortex core size $a$ in the RMmodel (2.4), are grouped as vortex parameters, as they illustrate the vortex characteristics.

The choice of optimiser and hyperparameter settings also significantly impacts the performance. In this study, we use the Adam optimiser with initial learning rate 0.01, which decays by a factor 0.99 every 10 steps. These parameters ensure stability and fast convergence during optimisation. While a systematic exploration of these hyperparameters’ effects is beyond the scope of this study, preliminary experiments suggest that they are reasonable choices under the proposed method.

We perform a detailed analysis of the two parameter categories. The optimisation parameters include the circuit depth $N_d$ and the iteration number $N_i$ . The circuit depth $N_d$ determines the scale of the quantum circuit, and in the context of the VQE, it defines the search space for $\psi$ . A deeper circuit can lead to better results but may become less effective due to increased noise and slower convergence (Cerezo et al. Reference Cerezo2021), while also demanding more parameters, further hindering optimisation efficiency. Similarly, the iteration number $N_i$ controls the optimisation steps: more iterations tend to improve results and robustness, but at the cost of increased computational time (Lavrijsen et al. Reference Lavrijsen, Tudor, Muller, Iancu and De Jong2020).

The experiments were conducted using the ‘leapfrogging’ example, shown in figure 3. Initially, we discretised the case into a $16 \times 16$ grid to investigate the relationship between $N_i$ , $N_d$ and the MSE. As shown in figure 8(a), the product of $N_i$ and $N_d$ dominates performance. We refer to this product as the complexity of the VQE, as the number of quantum gates is proportional to $N_d$ , and each gate must be updated $N_i$ times. Along the horizontal axis, increasing $N_d$ reduces the MSE significantly at first, reaching optimal performance within a specific range ( $N_d = 100{-}200$ ), after which performance stabilises or deteriorates. Along the vertical axis, higher complexity ( $N_i N_d = 60\,000$ or $70\,000$ ) achieves lower and more stable MSE. The interaction between $N_d$ and $N_i N_d$ indicates that while complexity dominates performance, selecting an appropriate $N_d$ within the optimal range is crucial for efficient optimisation, highlighting the trade-off between resource allocation and performance.

To further investigate the role of VQE complexity, we estimate the minimum required for convergence across different grid resolutions, ensuring that $N_d$ and $N_i$ remain within an appropriate range for effective optimisation. Figure 8(b) illustrates a linear relationship between the minimum complexity $(N_d N_i)_{min }$ and the grid dimension, given by $\sqrt {N_{{cell}}}$ . The problem is discretised into grids $8 \times 8$ , $16 \times 16$ , $32 \times 32$ and $64 \times 64$ , with optimal $N_d$ values 150, 200 and 250 for the first three smaller grids, and 2000, 2200 and 2400 for the $64 \times 64$ grid due to increased scale of the problem. The results are averaged over three different values of $ N_d$ for each grid to obtain representative $(N_d N_i)_{min }$ values, as shown in the zoomed-in views. The observed linear trend in the plot suggests that the minimum complexity required for convergence is proportional to $\sqrt {N_{{cell}}}$ , offering a practical guideline for estimating computational costs as problem size increases. We found that the product $N_d N_i$ is the primary determinant for VQE performance. Moreover, we identified a linear relationship between the minimum complexity and the grid resolution, providing valuable insights for optimising VQE parameters.

Figure 8.

Results of the optimisation parameter study. (a) Performance on a $16 \times 16$ grid for varying complexities and circuit depths, demonstrating that complexity has a dominant influence, with optimal depths in the range 100–200. (b) Estimated minimum complexity $(N_d N_i )_{min }$ required for convergence as a function of $\sqrt {N_{{cell}}}$ , showing a linear trend.

The vortex parameters considered in this study include the vortex strength $h$ introduced in (2.8), and the vortex core size $a$ introduced in (2.4). In our mathematical model, $h$ represents an estimate of the actual vortex strength. This parameter is straightforward to determine in simplified cases, such as the ‘leapfrogging’ example described earlier, where all vortices share identical strength, and the velocity field is computed using the BS model or RM model. However, in more complex cases such as turbulence, accurately estimating vortex strengths becomes significantly more challenging, as the vortices exhibit varying intensities. For instance, velocity and energy spectra are commonly employed to characterise turbulence (Kolmogorov Reference Kolmogorov1941; Lundgren Reference Lundgren1982). Moreover, a strong vortex can be approximated as the superposition of weaker vortices when $h$ is reduced, but our grid-based approach constrains each grid cell to a single vortex filament. The trade-off between these factors merits further exploration.

The core size $a$ not only prevents singularities but also captures variations in vortex structure, which is essential for accurately modelling more complex vortex behaviours. For instance, in the ‘leapfrogging’ experiments by Lim (Reference Lim1997), Jerrard & Smets (Reference Jerrard and Smets2018) and Wang et al. (Reference Wang, Chen, Chang, Li and Gao2024), even idealised point vortices exhibit deviations from expected behaviour. Consequently, selecting appropriate values for $h$ and $a$ is critical for accurate velocity reconstruction.

To evaluate the influence of these parameters in complex cases, we conducted experiments using a turbulent velocity field, with detailed descriptions provided in § 5.2. The case was discretised on $64 \times 64$ and $128 \times 128$ grids to assess performance under varying resolutions. As shown in figure 9(a), with $a$ fixed at 0.3, the optimal $ h$ value is relatively larger in coarser grids, while the optimal value decreases in finer grids. This indicates that higher grid resolutions tend to favour smaller vortex strengths, aligning with our hypothesis regarding vortex superposition. Figure 9(b) shows the relationship between the core size $ a$ and the MSE for different grid resolutions and vortex strengths. The MSE generally decreases as $ a$ increases, reaching a minimum before slightly increasing at larger $ a$ . Based on these experiments, we observe that the optimal values of $h$ and $a$ in finer grids are generally smaller. On the one hand, this result supports our previous assertion that strong vortices can be effectively represented by the superposition of weaker vortices: finer grids contain more cells, allowing for more vortices and better supporting the superposition. On the other hand, regarding the decrease in $a$ , we speculate that in finer grids, the core size of each vortex becomes smaller, and the presence of more vortices inherently introduces a smoothing effect.

Figure 9.

Results of the vortex parameter investigation. (a) The impact of $ h$ across different grid resolutions. In the $ 64 \times 64$ grid, the optimal value of $ h$ is 0.8, while in the $ 128 \times 128$ grid, it decreases to 0.3.Here, $a$ is fixed at 0.3. (b) The effect of the core size $ a$ . The optimal $a$ forthe $ 64 \times 64$ case $h = 0.8$ is 0.35, and for the $ 128 \times 128$ case $h = 0.3$ it is 0.25, showing a consistent reduction in the optimal core size $a$ as grid cell size decreases across the cases.

4.2. Complexity analysis

We analyse the computational complexity of our proposed approach. As outlined in § 3, the pipeline comprises three main stages: formulation of the optimisation problem, the VQE, and post-processing. The latter two stages are computationally intensive and operate in consecutive order. We conduct a detailed complexity analysis of these two components.

The VQE employs a parametrised quantum circuit to prepare a trial scalar field $\psi$ . The objective function $\langle \psi | H | \psi \rangle$ is then measured, according to which the parameters of each controlled quantum gate are optimised using classical optimisers. After $N_i$ iterations, an optimised set of parameters and the corresponding $\psi$ are obtained.

The primary computational cost of the VQE arises from back-propagation during each optimisation step (Cerezo et al. Reference Cerezo2021). In each iteration, every gate undergoes a gradient-based update. For a grid with $N_{{cell}}$ cells, the problem involves a Hermitian matrix of size $N_{{cell}} \times N_{{cell}}$ , corresponding to $\log (N_{{cell}})$ qubits. According to our circuit design, a basic block in figure 6 consists of $4 \log (N_{{cell}})$ controlled gates (e.g. $R_y$ and $R_z$ ). In a quantum circuit with depth $N_d$ , the total number of controlled gates becomes $4 N_d \log (N_{{cell}})$ . The execution time of the quantum circuit is proportional to the circuit depth $N_d$ . Thus the time complexity of one optimisation step can be expressed as

(4.2) \begin{equation} T(N_{{cell}}, N_d, N_i)_{{one\ step}} \sim \mathcal{O}\big(N_d^2 \log N_{{cell}}\big). \end{equation}

For $N_i$ iterations, the overall complexity of the VQE is

(4.3) \begin{equation} T(N_{{cell}}, N_d, N_i)_{\textit{VQE}} \sim \mathcal{O}\big(N_d^2 N_i \log N_{{cell}}\big). \end{equation}

The term $N_d N_i$ corresponds to the complexity discussed in§ 4.1. Our analysis indicates that the minimum complexity $N_d N_i$ required for convergence scales linearly with the square root of the total grid size, which is approximately $\sqrt {N_{{cell}}}$ . Accordingly, we conjecture that the overall complexity of the VQE can be approximated as

(4.4) \begin{equation} T(N_d, N_{{cell}})_{\textit{VQE}} \sim \mathcal{O}(N_d \sqrt {N_{{cell}}} \log N_{{cell}}). \end{equation}

Since $\mathcal{O}(\sqrt {N_{{cell}}} \log N_{{cell}})$ is smaller than $\mathcal{O}(N_{{cell}})$ , the VQE complexity can be further bounded as

(4.5) \begin{equation} T(N_d, N_{{cell}})_{\textit{VQE}} \sim \mathcal{O}(N_d N_{{cell}}). \end{equation}

The classical post-processing method described in § 3.3 involves two main steps: computation of the winding number and interpolation. For the winding number computation, in an $m \times j \times k$ grid, the discretised winding number must be evaluated for a total of $(m - 1)(j - 1)k + (m - 1)j(k - 1) + m(j - 1)(k - 1)$ faces. Assumingthat the grid dimensions $ m, j, k$ scale as $ N_{{cell}}^{1/3}$ , if the winding number must be computed for every face in the grid, and each face requires a constant amount of computation, then the time complexity of the winding number calculation can be approximated as

(4.6) \begin{equation} T(N_{{cell}})_{{winding}} \sim \mathcal{O}(N_{{cell}}). \end{equation}

The interpolation step is required only for faces where the winding number is non-zero. In the worst-case scenario, interpolation is performed on every face as described in (3.14), yielding a total complexity

(4.7) \begin{equation} T(N_{{cell}})_{{interpolation}} \sim \mathcal{O}(N_{{cell}}). \end{equation}

Consequently, the overall complexity of the classical post-processing method is given by

(4.8) \begin{equation} T(N_{{cell}})_{{post}} \sim \mathcal{O}(N_{{cell}}) + \mathcal{O}(N_{{cell}}) \sim \mathcal{O}(N_{{cell}}). \end{equation}

In addition, the data-driven post-processing approach exhibits a comparable complexity of $\mathcal{O}(N_{{cell}})$ , with further discussion provided in Appendix B.

In summary, the overall complexity of the proposed pipeline is

(4.9) \begin{equation} T(N_d, N_{{cell}}) \sim \mathcal{O}(N_d \sqrt { N_{{cell}}} \log (N_{{cell}})) + \mathcal{O}(N_{{cell}}) \sim \mathcal{O}(N_d N_{{cell}}). \end{equation}

This complexity is considered advantageous. We compare our method with several classical solvers commonly used for linear systems, including QR (Q for the orthogonal or unitary matrix, and R for the upper triangular matrix ) decomposition, Gaussian elimination, and the Lanczos algorithm. The QR decomposition and Gaussian elimination are direct methods, typically requiring $\mathcal{O}(N_{{cell}}^3)$ for an $N_{{cell}} \times N_{{cell}}$ matrix, which makes them computationally expensive for large systems or fine grids (Francis Reference Francis1961; Higham Reference Higham2011). Iterative methods, such as Lanczos methods, can reduce complexity to $\mathcal{O}(m N_{{cell}}^2)$ , where $m$ denotes the number of desired eigenvalues (Lanczos Reference Lanczos1950). And although only the minimum eigenvalue is required in our case, the computational cost remains quadratic at $ \mathcal{O}(N_{{cell}}^2)$ . In contrast, our quantum approach achieves a complexity of $\mathcal{O}(N_d N_{{cell}})$ under ideal conditions, with $N_d$ much smaller than $N_{{cell}}$ , offering a potential speedup.

However, it is important to note that the measurement of both the objective function $ \langle \psi | H | \psi \rangle$ and the final quantum state $ | \psi \rangle$ can be time-consuming, which was not initially considered in our analysis as the simulations were conducted on a virtual quantum machine with full access to state information. As described in § 3.2, the objective function is measured via decomposition of the Hamiltonian into Pauli operators. For an $ N_{{cell}} \times N_{{cell}}$ Hermitian matrix, this results in $ \mathcal{O}(N_{{cell}}^2)$ distinct Pauli terms, leading to a measurement complexity of $ \mathcal{O}(N_{{cell}}^2)$ . To mitigate this overhead, we introduce a Pauli operator cut-off strategy in Appendix A that does not require additional qubits. This strategy selectively removes Pauli terms with negligible contributions, thereby offering a simple and effective way to reduce the measurement burden while preserving noise-robust implementation. Another promising alternative is block-encoding, which embeds a non-unitary Hermitian matrix into a higher-dimensional unitary operator (Gilyén et al. Reference Gilyén, Su, Low and Wiebe2019; Martyn et al. Reference Martyn, Rossi, Tan and Chuang2021). This approach circumvents explicit Pauli decomposition and can be used in conjunction with the Hadamard test for efficient expectation value estimation. Nevertheless, it requires additional ancilla qubits and more complex quantum circuits, increasing susceptibility to noise, particularly in large-scale fluid simulations. In future work, block-encoding could be explored as a complementary technique within our framework.

In addition, the extraction of filaments from $\psi$ requires full knowledge of the final quantum state $ | \psi \rangle$ . On actual quantum hardware, this necessitates quantum state tomography (O’Donnell & Wright Reference O’Donnell and Wright2016), which reconstructs the quantum state from multiple measurement outcomes. The corresponding measurement complexity scales as $ \mathcal{O}(N_{{cell}} / \varepsilon ^2 )$ , where $ \varepsilon$ denotes the target precision. To alleviate this challenge, we discuss two potential strategies. First, rather than reconstructing the entire quantum state, one may extract only essential information using techniques such as classical shadow tomography (Huang, Kueng & Preskill Reference Huang, Kueng and Preskill2020) and sparse tomography (Chen et al. Reference Chen2024), reducing the measurement complexity to $ \mathcal{O}(\log N_{{cell}} / \varepsilon ^2)$ under suitable assumptions. Second, as introduced in Appendix B, the machine-learning-based post-processing framework can incorporate quantum noise and measurement uncertainty into the training process. This data-driven approach may relax the required measurement precision $ \varepsilon$ , thereby further improving overall efficiency.

5.1. Two-dimensional point vortices

We design three different two-dimensional vortex cases to test our method’s ability to recover vortex dynamics. For all configurations, the circuit depth ( $N_d$ ) ranges from 50to 400, and each vortex configuration is discretised on a $16 \times 16$ grid to accurately represent the flow field. Optimisation is performed using the Adam optimiser with $N_i = 300$ iterations to ensure model convergence. As shown in figure 10(a), the MSE for all configurations converges below 0.01, demonstrating that the quantum circuits accurately capture the vortex dynamics. Specifically, the firstconfiguration, labelled‘leapfrogging’, as shown in figure 10(b), features a symmetric arrangement of two positive and two negative vortices, chosen to test the circuit’s ability to capture well-defined, symmetric vortex patterns. The recovered streamlines closely match the ground truth, with maximum distance error $d_{{max}}=0.054$ , indicating that the method performs well in handling simple and symmetric vortex systems. The second and third configurations, ‘Random example A’ and ‘Random example B’, are randomly generated systems with four and five vortices, respectively, as shown in figures 10(c) and 10(d). These set-ups challenge the method’s ability to capture more complex vortex distributions with varying strengths and spatial arrangements. In both cases, the recovered streamlines closely follow the general flow structure with slight deviations, and the recovered vortices align well with the ground truth, demonstrating the method’s accuracy in tracking vortex positions despite minor discrepancies.

Figure 10.

Evaluation of quantum circuit expressibility with three vortex configurations. (a) The MSE converges below 0.01 for all cases. (b) The ‘leapfrogging’ case features a symmetric arrangement of two positive and two negative vortices. (c) ‘Random example A’ and (d) ‘Random example B’ are randomly generated systems with four and five vortices, respectively.

We also assess our method’s ability to capture the time-dependent ‘leapfrogging’ evolution, where four vortices in two pairs exhibit a periodic leapfrog-like motion. Driven by mutual advection, each vortex alternately leaps over and pursues the others, resulting in structured yet dynamic motion with periodic accelerations and decelerations. In numerical implementation, we divide the evolution into 50 time steps over a 10 s period, capturing the vortex behaviour at regular intervals. The calculations are performed on a $32 \times 32$ grid within a quantum circuit of depth 250, optimised over 1200 iterations to ensure computational precision. The results, shown in figure 11, demonstrate the method’s capability to accurately track the evolving vortex motion. With average MSE 0.003 and maximum error 0.046, the results reflect the method’s robustness in characterising the vortex dynamics. By accurately modelling the ‘leapfrogging’ motion, our method proves to be a valuable tool for studying and understanding vortex evolution in various fluid dynamics scenarios.

Figure 11.

Evolution of ‘leapfrogging’ vortices. (a) Target and reconstructed vortices at selected time steps. The inset marks the maximum error region. (b,c) Time histories of maximum error ( $d_{max }$ ) and MSE, showcasing the method’s performance throughout the simulation.

5.2. Two-dimensional turbulent flows

To evaluate the applicability of the method to more complex scenarios, it is applied to two-dimensional turbulent flows. The velocity field is generated using the energy spectrum $E(k) = \exp (-\alpha k^2)$ , where $k = \sqrt {k_x^2 + k_y^2}$ is the wavenumber, representing the magnitude of the wavevector in Fourier space, and $\alpha$ is the decay factor that governs the energy distribution across different wavenumbers. A higher $\alpha$ leads to a faster decay of energy, favouring smaller-scale structures. The Fourier coefficients for the velocity components are given by

(5.1) \begin{equation} \hat {u}(k) = -\frac {k_y}{k}\, E(k)\,{\rm e}^{{\rm i}\,\theta (k)}, \quad \hat {v}(k) = \frac {k_x}{k}\, E(k)\,{\rm e}^{{\rm i}\,\theta (k)}, \end{equation}

where $\hat {u}(k)$ and $\hat {v}(k)$ are the Fourier components of the velocity field in the $x$ - and $y$ -directions, respectively, with $\theta (k)$ representing random phases for isotropy. The inverse Fourier transform reconstructs the physical velocity field, converting the spectral representation into real-space velocity components and generating a spatially coherent turbulent structure. The flow intensity is controlled by scaling the velocity components with a factor $A_{{scale}}$ . In the present computations, $\alpha = 0.2$ , $A_{{scale}} = 8000$ , and the random seed is set to 7 for reproducibility. The flow is discretised on a $128 \times 128$ grid with grid spacing $\Delta x = 0.6$ . Optimisation is performed with circuit depth $N_d = 3000$ and $N_i = 1333$ iterations.

The simulation results presented in figure 12 illustrate the complex turbulent structure and demonstrate the effectiveness of the method in dynamic flow scenarios. Figure 12(a) shows the recovered velocity field $\boldsymbol{u}$ , with colour representing velocity magnitude, and grey streamlines indicating the flow direction, while figure 12(b) displays the corresponding velocity field $\boldsymbol{u}_{{target}}$ for comparison. It is evident that the method effectively recovers the turbulent flow, with the recovered velocity field resembling the target field and capturing key features such as vortex structures and flow direction, demonstrating the algorithm’s capability in reconstructing the target velocity field.

Figure 12.

Results for two-dimensional turbulent flows. (a) The recovered velocity field $\boldsymbol{u}$ , with colour representing velocity magnitude, and grey streamlines indicating the flow direction. (b) Similarly for the target velocity field $\boldsymbol{u}_{{target}}$ .

Figure 13.

Foil extraction results. (a,b) The positions of the extracted vortex filaments and target vortex filaments for the septafoil and cinquefoil cases, respectively. The light-blue dots and lines indicate the positions of extracted points and filaments, and the orange line indicates the ground truth. (c,e) The shapes of the two target vortex filaments. (d,f) Iso-surfaces of $\mathrm{Im}(\psi ) = 0$ , coloured by the real part of the complex scalar field $\mathrm{Re}(\psi )$ . The zero-intersect lines, representing the vortex filament positions, are marked in orange.

5.3. Three-dimensional knotted vortex filaments

We further extend our method to three-dimensional knotted vortex filaments, a representative flow characterised by intricate topological structures and rich dynamical behaviour (Barenghi Reference Barenghi2007; Kleckner & Irvine Reference Kleckner and Irvine2013; Tao et al. Reference Tao, Ren, Tong and Xiong2021). We construct our knotted vortex filaments according to Xiong & Yang (Reference Xiong and Yang2019a , Reference Xiong and Yang2020). The parametric equations of the knot $\boldsymbol{c}(\theta ) = (c_x (\theta ), c_y(\theta ), c_z(\theta ) )$ are

(5.2) \begin{equation} \begin{cases} c_x(\theta ) = \left (R_t + r_t \cos (q \theta ) \right ) \cos (p\theta ), \\ c_y(\theta ) = \left (R_t + r_t \cos (q \theta ) \right ) \sin (p\theta ), \\ c_z(\theta ) = -r_t \sin (q\theta ), \end{cases} \end{equation}

where $R_t$ and $r_t$ are the major radii of the knot.Here, $p$ and $q$ represent the shape of the knotted foil. In our case, $R_t = 1$ and $r_t = 0.5$ , so the foils are neither too big nor too small that the details would be lost during discretisation. We would like to test a septafoil and a cinquefoil, therefore $(p, q)$ is selected to be $(2, 7)$ and $(2, 5)$ . The shapes of the two foils are shown in figures 13(c) and 13(e), respectively. The computational domain is $(-2,-2, -2) \times (2, 2, 2)$ , and is discretised into a $32 \times 32 \times 32$ grid. Referring to the discussion in § 4.1, two cases are optimised with a $N_d = 3500$ circuit depth and 1200 iterations.

Figures 13(a,b) show the comparision between target and extracted filament positions, which shows a good match. Figures 13(d,f) show the iso-surface where $\mathrm{Im}(\psi ) = 0$ , coloured by the values of $\mathrm{Re}(\psi )$ . It can be seen that the extraction results are located along zero sets of $\psi$ . The final MSEs are $0.057$ for septafoils and $0.035$ for cinquefoils , indicating the method’s effective performance in these three-dimensional cases.

5.4. Three-dimensional turbulent flows

We present results for three-dimensional turbulent flows obtained from numerical simulations. A random initial velocity field is constructed within the domain $[- , ]^3$ by assigning Fourier phases to the velocity components. This initial field serves as input to the Navier–Stokes equations, solved using a spectral method to ensure isotropic turbulence properties. The velocity field is then discretised on a uniform $32 \times 32 \times 32$ three-dimensional grid. For vortex extraction, we employ circuit depth 4000 and perform 1300 optimisation steps. The results for $h = 10$ are shown in figure 1. The enlarged view reveals entangled vortex filaments, characteristic of turbulent flows.

Figure 14.

Effect of filament strength $ h$ on convergence and extraction. (a) Convergence of the objective function $\langle \psi | H | \psi \rangle$ , with larger $ h$ yielding smaller final values, indicating improved approximations. Extracted filaments for(b) $ h = 10$ ,(c) $ h = 20$ ,(d) $ h = 30$ , and(e) $ h = 40$ . Increasing $ h$ reduces the number of extracted filaments, emphasising the trade-off discussed in§ 4.1.

The influence of vortex strength, characterised by the parameter $h$ , on the results is further analysed. Figure 14(a) demonstrates the convergence of the objective function $ \langle \psi | H | \psi \rangle$ over iterations for different filament strengths $ h$ . While the objective function for all configurations converges to stable small values, the final convergence value decreases with increasing $h$ . As illustrated in figures 14(b–e), the number of vortex filaments diminishes significantly as $h$ increases. This suggests that a balance between vortex strength and filament density must be achieved, corroborating the trade-off discussed in § 4.1.

Figure 15.

Contours of the scalar field $\psi$ obtained from three-dimensional turbulences. (a) Contours of the magnitude $|\psi |$ at $y = 0$ . (b) Contours of the phase $\arg (\psi )$ at $y = 0$ .

Finally, the distributions of $|\psi |$ and the phase $\arg (\psi )$ are examined across selected cross-sections. For $h = 10$ , displayed in figure 15, contour plots of $|\psi |$ and $\arg (\psi )$ are shown for the planes $y = 0$ and $z = 0$ . The results highlight the spatial structure of $\psi$ , where the scalar field aligns with underlying coherent features. These observations provide further insights into the relationship between $\psi$ and the flow dynamics.

In summary, the method proposed in this study demonstrates effectiveness in three-dimensional turbulence cases. It provides a novel approach for vortex representation and analysis, contributing to the structural understanding of turbulence.