Towards synthetic magnetic turbulence with coherent structures

Synthetic turbulence is a relevant tool to study complex astrophysical and space plasma environments inaccessible by direct simulation. However, conventional models lack intermittent coherent structures, which are essential in realistic turbulence. We present a novel method featuring coherent structures, conditional structure function scaling and fieldline curvature statistics comparable to magnetohydrodynamic turbulence. Enhanced transport of charged particles is investigated as well. This method presents significant progress towards physically faithful synthetic turbulence.

Introduction.-Turbulence plays a key role in astrophysical and space plasma environments [1][2][3][4][5][6].However, due to the high computational cost of direct approaches, the effect of turbulence in such environments is difficult to study.This obstacle is often mitigated by splitting the magnetic field in a large-scale coherent component with an analytic description and a small-scale turbulent component, modelled as a Gaussian random field [e.g., 7,8].Such Gaussian random fields can be easily synthesized as a superposition of plane waves with random phases and a prescribed energy spectrum.The transport of energetic charged particles through such fields has been extensively studied [e.g., [9][10][11][12][13].
However, Gaussian random fields can only provide a low-order approximation of magnetic turbulence, neglecting any structure beyond two-point correlations captured by the energy spectrum.They do not exhibit intermittency as observed in first-principles turbulence [e.g., [14][15][16][17].Intermittency was studied in the context of hydrodynamic synthetic turbulence models already by Juneja et al. [18], and its impact on charged particle transport more recently by Pucci et al. [19] and Shukurov et al. [20], finding faster diffusion in structured magnetic fields.
Up until today, there have been several models for synthetic hydrodynamical and magnetohydrodynamical (MHD) turbulence published, such as the p-model on a discrete three-dimensional wavelet space by Malara et al. [21], the minimal multiscale Lagrangian map by Rosales and Meneveau [22], which was applied to MHD turbu- (a) Email: jeremiah.luebke@rub.delence by Subedi et al. [23], a stochastic integral based on the lognormal-model including asymmetric velocity increment statistics by Pereira et al. [24], and its application to MHD turbulence by Durrive et al. [25].All of these models produce three-dimensional, divergence-free vector fields with intermittent statistics, but without coherent geometric features.Recently, Durrive et al. [26] presented a model which embeds Archimedean spirals into a random lognormal vector field.The continuous wavelet cascade by Muzy [27] addresses broken stationarity of discrete wavelet cascades.Related works were recently published by Li et al. [28] and Robitaille et al. [29].
Standard tools of validating synthetic models are the energy spectrum and statistics of field increments.Taken without further decomposition, these quantities provide a global picture of the vector field, hiding the intricate local geometry of magnetic turbulence [see also 30,31].A useful quantity in this regard is the fieldline curvature, which has recently been shown by Kempski et al. [32] and Lemoine [33] to play a key role in the transport of charged particles in magnetic turbulence.Fieldline curvature has previously been discussed in the context of turbulent dynamos [34], as well as hydrodynamic [35][36][37] and magnetohydrodynamic turbulence [38,39].
In this letter we present progress towards a model for synthetic magnetic turbulence featuring intermittent coherent structures.We implement the model as a fast algorithm, which produces a random three-dimensional divergence-free vector field, resembling a turbulent magp-1 arXiv:2401.10573v3[physics.space-ph]6 May 2024 netic field b(x) 1 .The model is a combination of a continuous cascade [27] and the minimal multiscale Lagrangian map [22,23].Additionally, we propose a set of quantities to assess the physical fidelity of synthetic turbulence models, consisting of the energy spectrum, conditional structure function scaling, the fieldline curvature distribution and running diffusion coefficients of charged test particles.Based on these quantities, we compare the proposed synthetic turbulence model with an incompressible resistive MHD turbulence simulation and an intermittent synthetic turbulence model without coherent structures.We also consider the phase-randomized counterparts of the three turbulence models to account for differences in the energy spectra.We conclude by explaining the shape of the MHD fieldline curvature distribution by means of a weighted sum of Gaussian components.
Methods.-We start by extending the continuous cascade in wavelet space [27] to three-dimensional divergence-free vector fields.The continuous cascade at scale l and position x is represented by a log-infinitely divisible process e ω l (x) , which gives the scale-and positiondependent intensity of a vector field v(x).This field is obtained by a vector-valued wavelet transform of l H e ω l over the inertial range scales l min < l < l 0 as The slope of the energy spectrum is given by −2H − 1, torodial wavelets ∇ × lψ l (x)ẑ with ψ(k) = −k 2 e −k 2 and ψ l (x) = ψ(x/l) ensure the zero-divergence condition ∇ • v = 0, a random rotation field R l (x) with correlation length l ensures proper isotropization of these wavelets, and the numerically determined constant A normalizes the field to ⟨v 2 ⟩ = 1.Further, the curl can be moved in front of the spatial convolution operator f * ∇ × g (x) = ∇ × R 3 f (y)g(x − y) dy and the scale integral • • • dl, thus allowing us to express the field in terms of a vector potential v = ∇ × a.
The infinitely divisible process ω l (x) is defined on conelike subsets of the position-scale half-space R d × R >0 equipped with the measure l −d−1 dx dl and has a cumulant generating function ϕ(q).Thus, the moments of the intensity process can be computed as ⟨e qω l ⟩ = (l 0 /l) c d ϕ(q) ∀q ∈ N, where c d = 2 −d π d/2 /Γ(d/2 + 1) comes from the scalespace cone volume.We consider a Gaussian distribution with ϕ(q) = µ/2(p 2 − p) with intermittency parameter µ, in which case e ω l corresponds to the Gaussian multiplicative chaos employed in [24] and related works.
As shown below, the vector field given by Equation ( 1) exhibits (isotropic) anomalous scaling, but lacks coherent features, so we introduce advective structures by adapting the minimal multiscale Lagrangian map (MMLM) for 1 The code is publicly available at https:// github.com/jerluebke/synth-mag-turband archived at doi:10.5281/zenodo.10515965.A pseudocode listing is provided as supplementary material.
MHD [23] to our framework.In short, the MMLM procedure considers only the advective part of the magnetic field evolution equation, i.e. (∂ t + u • ∇) b = 0, which can be solved at time τ with the Lagrangian ansatz b(x τ , τ ) = b(x 0 , 0) and the linearized solution x τ = x 0 + τ u(x 0 , 0).Usually, two initially Gaussian random vector fields representing u(x) and b(x) are deformed on successively finer scales l i by applying the linearized Lagrangian solution with τ ∝ l i to the underlying regular grid, as described in the references.
In our framework, we generate two independent random fields according to Equation ( 1), which are represented as vector potentials and represent the velocity field u(x) = ∇ × a u (x) and the magnetic field b(x) = ∇ × a b (x).When generating u(x) scale-by-scale, we make use of the discretization of the scale integral in Equation ( 1) as a sum which goes from large to small scales.Specifically, we accumulate the deformations of the grid x by the intermediate results u li = ∇ × a u li over all scales l i .We then interpolate the random magnetic vector potential a b (x) at the final deformed grid.Note that this interpolation is not done intertwined with the scale-by-scale generation of a b (x), but only once at the end of the procedure.Thus, we are solving the advection equation for the magnetic vector potential which is exact in two dimensions making this approach especially suitable for strong guide field situations [40].
Herein lies the key difference with the previous MMLM procedure, which applied the interpolation directly to the field b(x) in a scale-by-scale fashion.It should further be noted, that the curl of the intermediate results u li = ∇ × a u li is computed on a uniform grid, so while u li serves as an indicator for the grid deformation, it is not affected by it.
The deformation timescale τ = c l i / max(∥u li ∥) is normalized to the maximal value of the current velocity magnitude and governed by the constant c.This constant is a free parameter, which must be chosen carefully, as it should be large enough for coherent structures to emerge, but not too large to avoid decorrelation.For higher values of c, energy accumulates on smaller scales, which is corrected by reweighting a b lmin (x) after interpolation with k δ/2 , where δ is the deviation of the energy spectrum scaling from the expected −5/3rd scaling.Additionally, we mimic the effect of dissipation by applying a low-pass filter exp(−k 2 /2k 2 0 ) controlled by the artificial dissipation wavenumber k 0 .
For comparison, we perform a direct numerical simulation of incompressible resistive MHD turbulence in a three-dimensional periodic box with resolution N 3 = 1024 3 , no background field and equal diffusivity and resistivity ν = η = 1.2 × 10 −3 [41].We obtain Taylorscale Reynolds numbers R λ = 439 for the velocity field and R λ,m = 94 for the magnetic field.The velocity is driven on Fourier modes 1 ≤ k ≤ 3 with the random forcing proposed in [42], which exhibits very low mean cross-helicity and low noise in the time evolution of turbulence bulk quantities.For this purpose we employed the pseudo-spectral code SpecDyn, which was developed in the context of magnetic dynamo action and is tailored for use on modern HPC systems [43,44].
Results.-We generate sample fields of the continuous cascade process (CC) given by Equation ( 1) and its minimal multiscale Lagrangian mapping extension (LM) corresponding to Equation (2) on a three-dimensional periodic grid with resolution N 3 = 1024 3 .The parameters of both processes are H = 1/3, l 0 = 0.5 and l min ≈ 1.5 dx.The intermittency parameter of the CC process is µ = 0.4, and the additional parameters of the LM process are µ = 0.1, c = 0.2, δ = −0.45,and k 0 = 256.
For visual inspection, slices of field strength and fieldline curvature are plotted in Figure 1, together with a plot of the radially averaged energy spectra.MHD turbulence is characterized by elongated and intricately intertwined coherent structures, while the CC field consists of incoherent intermittent "clouds" of large field strength values.The LM field exhibits thin and intense coherent structures, which are more spatially isolated.The CC spectrum matches the expected −5/3rd scaling well, the LM spectrum exhibits a bottleneck effect at high wavenumbers, and the MHD spectrum is affected by strong resistivity.The bottleneck effect of the LM spectrum is due to the simple advective nature of the Lagrangian mapping, which causes accumulation of energy on small scales.This is counteracted by reweighting the vector potential with k δ/2 and applying artificial dissipation, however, when balancing the relevant parameters c, δ and k 0 , a residual bottleneck effect remains.In order to take potential influences of these different spectra into account, the analysis of curvature and particle transport is also performed with the respective phase-randomized fields.We employ conditional structure functions [45] to study intermittency with respect to the local structure of the field.See also [46] for a discussion of local anisotropy statistics.Given a displacement vector d at a point X, we consider increments δb  Figure 2 shows the normalized scaling exponents ζ p /ζ 3 as well as the raw values of ζ 3 for the three models and the three directions.The further ζ p /ζ 3 deviates from the linear case, the more intermittent is the process, and smaller values of ζ 3 correspond to a rougher process.In the MHD case, the field is the most smooth and non-intermittent parallel to the local mean field direction blocal , while it is equally intermittent in the direction of fluctuations δ b⊥,N and the binormal direction blocal × δ b⊥,N , and being the most rough in the binormal direction.In contrast to this, the CC field shows no significant difference between the three directions, which are all equally rough and intermittent, which is expected from the lack of coherent structures.Lastly, the LM field exhibits again directional dependency, with the parallel direction blocal being the most non-intermittent and the binormal direction blocal ×δ b⊥,N being the most intermittent.However, the field is in the fluctuation direction δ b⊥,N less intermittent compared to the MHD case.
While the conditional structure functions provide an extensive statistical picture of magnetic turbulence, additional geometric insight can be gained from the fieldline curvature (3) Figure 3 shows the distributions p(κ) for the three cases and their random-phase counterparts.The MHD case, in agreement with the literature [38,39], behaves asymptotically as p(κ) ∼ κ→∞ κ −2.5 .Note, as shown by [38], the scaling becomes κ −2 in 2D turbulence, while similar values are expected for general collisionless plasmas.The CC case agrees, apart from being slightly wider, with the Gaussian random-phase fields, which scale distinctly as p(κ) ∼ intermediate case between the previous two cases; its distribution function p(κ) is slightly more narrow than the MHD case and around the slightly right-shifted peak, p(κ) faintly resembles the κ −2.5 scaling, before adjusting to the Gaussian κ −4 scaling.
The extended flat tail for large κ in the MHD case is caused by a significant amount of intermittent sharp fieldline bends scattered throughout the domain [see also 32,33] and the low value of κ peak comes from coherent structures extended on scales comparable to the box size.In contrast to this, fieldline bends in Gaussian fields are distributed in a self-similar way in the domain, thus leading to a much stepper decay of the distribution.
Since the random-phase fields are Gaussian random variables, we expect a universal normalized fieldline curvature distribution p(κ/κ peak ), which is independent of the energy spectrum.This behaviour is observed in numerical experiments and also indicated by [38].However, κ peak does depend on the energy spectrum, e.g.via the slope s in case of a power-law spectrum ∼ k −s .Flatter spectra implicate more energy on small scales, resulting in more contributions from high curvatures and consequently larger κ peak .This connection is illustrated in Figures 4, where p MHD (κ) is modeled as a weighted sum of Gaussian components, similar to the description of the turbulent velocity increment distribution function as a Gaussian scale mixture [47,48].
In addition to the insight into the structure of the fields gained by the previous steps, we also study the transport of charged particles by numerically solving test particle trajectories X(t) according to the Newton-Lorentz equation with a Boris integrator [49] and trilinear interpolation [50].
The magnetic field is normalized to ⟨b 2 ⟩ = 1 and the  particles are parameterized by their normalized gyro radius r g /l 0 = γmcv 0 /qb 0 l 0 , where l 0 denotes the outer scale, b 0 denotes the root mean square strength of the magnetic field, and v 0 , γ = 1/ 1 − v 2 /c 2 , m and q denotes the particle's velocity magnitude, Lorentz factor, mass and charge.
We record mean squared displacements ⟨∆X 2 (t)⟩ = ⟨∥X(t) − X(0)∥ 2 ⟩ and diffusion coefficients D(r g ) = lim t→∞ ⟨∆X 2 (t)⟩/t, once the trajectories have reached diffusive behaviour.The diffusion coefficients obtained from the three models are plotted in Figure 5a, and compared with their random-phase counterparts and quasilinear predictions [10].Figure 5b shows the exemplary time evolution of ⟨∆X 2 (t)⟩ at r g /l 0 = 0.0039, consisting of an initial super-diffusive phase and short sub-diffusive phase, before arriving at stable diffusive behaviour.Particles have the largest diffusion coefficients in MHD turbulence on all scales, which stands in striking difference to the random-phase MHD field, where we find the smallest diffusion coefficients.This behaviour can be explained by the strong deviation of the MHD energy spectrum at high wavenumbers from the −5/3rd spectral slope, and highlights impressively the effectiveness of coherent structures in regard to charged particle transport.The CC case achieves only a very minor increase compared to the random-phase case, and while the LM case performs bet-ter, it is still outperformed by MHD.
Conclusion. -We have presented a novel algorithm for synthetic magnetic turbulence based on a combination of the continuous cascade model, generalized to threedimensional divergence-free vector fields, and the minimal multiscale Lagrangian map.The most important differences to previous works on the MMLM procedure are the explicit cascade structure of the underlying noise, primarily working with the vector potential, and interpolating from deformed to uniform grid once at the end of the procedure instead of interpolating at each scale.These changes considerably strengthen the emergent advective structures.We compare this algorithm with an incompressible resistive MHD turbulence simulation and the pure three-dimensional continuous cascade model, which is intermittent but lacks coherent structures.This comparison is done by means of visual inspection, the energy spectrum, conditional structure function scaling, the fieldline curvature distribution and running diffusion coefficients of charged test particles.
We observe that our algorithm produces turbulence exhibiting pronounced coherent structures, albeit not as densely and intricately organized as MHD coherent structures.This is accompanied by non-trivial conditional structure function scaling, revealing local anisotropy, i.e. relatively low roughness (ζ LM 3,∥ > ζ LM 3,⊥ ) and weak intermittency parallel to the local mean magnetic field, and strong intermittency in the perpendicular direction in agreement with MHD turbulence.However, when directly compared to the MHD case, the field in the parallel direction is clearly rougher (ζ LM 3,∥ < ζ MHD 3,∥ ), and the fluctuation direction is not as intermittent as required.Further, the fieldline curvature distribution resembles the MHD case at small and intermediate curvatures, but exhibits Gaussian scaling at large curvatures.Finally, while charged particle transport is enhanced, it remains outpaced by the MHD case, as expected due to the simpler geometry and smaller length scales of the synthetic coherent structures.
In conclusion, our algorithm presents significant progress towards simulating realistic turbulence.Remaining issues are clearly identified and will guide further improvements in designing synthetic turbulence models.For instance, a feedback mechanism during the algorithm would be highly relevant, which acts on the velocity field u and takes the current state of the deformed grid and the advected magnetic field b into account.An approach based on the Elsässer formulation of the MHD equations appears promising as well.Alternatively, one could aim to design a synthetic scalar curvature field, instead of a full vector field, and make use of recent results linking fieldline curvature and charged particle transport [32,33].Such an approach could build on the description of the non-trivial MHD fieldline curvature distribution as a weighted sum of Gaussian components, as presented in this letter.

Fig. 1 :
Fig. 1: Slice plots of magnetic field strength b/brms and logarithm of fieldline curvature κ/κ peak for three models of magnetic turbulence -Magnetohydrodynamics (MHD), Continuous Cascade (CC) and Lagrangian Mapping (LM).Energy spectra of the models are plotted in the right-most panel and compared with a −5/3rd scaling.The dissipation wavenumber of the MHD simulation kη is indicated on the abscissa.The fieldline curvature is defined as κ = ∥ b • ∇ b∥, with b = b/∥b∥, and normalized by the most frequent value κ peak .

3 Fig. 2 :
Fig. 2: Scaling exponents of conditional structure functions of the magnetic field δb ⊥ p | d ∝ d ζp for the MHD, CC and LM case.The averages of the increments δb ⊥ are conditional on the displacement vector d being aligned with the local mean field b local (∥), the normal fluctuation direction δ b⊥,N ⊥ blocal (δ b), or with the binormal direction δ b⊥,N × blocal (⊥).The scaling exponents are normalized as ζp/ζ3 and the raw values of ζ3 are shown in the top panels.
⊥ = b ⊥ (X+d)−b ⊥ (X) perpendicular to the local mean field b local = b(X + d) + b(X) /2.Averages over these increments are taken conditionally on the direction of d in an instantaneous local basis given by blocal , the normal direction of the fluctuations δ b⊥,N with δb ⊥,N = δb ⊥ − δb ⊥ • blocal blocal and the binormal direction blocal × δ b⊥,N .Based on this, the scaling exponents of the conditional averages ∥δb ⊥ ∥ p | d ∝ d ζp are denoted as ζ p,∥ , ζ p,δ b or ζ p,⊥ depending on d being aligned 10

Fig. 3 :
Fig. 3:Compensated distributions of fieldline curvature κ p(κ/κ peak ) for the MHD, CC and LM case, as well as their respective random-phase cases (dotted lines).The lowcurvature κ 1 scaling, the MHD κ −2.5 scaling and the Gaussian κ −4 scaling are indicated for comparison.The inset shows the respective values of κ peak

κ→∞ κ − 4 .Fig. 4 :
Fig.4: Modelling of the compensated MHD fieldline curvature distribution κ p(κ) as a weighted sum of shifted Gaussian fieldline curvature distributions.The Gaussian fields have power-law spectra S(k) ∼ k −s , where the spectral slope s determines the value of the most dominant curvature κ peak .The weights (not shown) scale analogously to the compensated high-curvature tail with κ −1.5 .

Fig. 5 :
Fig. 5: (a) Diffusion coefficients for gyro radii rg/l0 = 0.001, • • • , 0.25 for the MHD, CC and LM case, as well as their respective random-phase cases (dotted lines).High-and lowenergy predictions from quasilinear theory are given for comparison.The dissipative length scale lη and the gyro radius of the example in the lower panel are indicated on the abscissa.(b) Running diffusion coefficients of particles with gyro radius rg/l0 = 0.0039 as an example to illustrate the temporal evolution of the transport.