BxC Toolkit: Generating Tailored Turbulent 3D Magnetic Fields

Turbulent states are ubiquitous in plasmas, and the understanding of turbulence is fundamental in modern astrophysics. Numerical simulations, which are the state-of-the-art approach to the study of turbulence, require substantial computing resources. Recently, attention shifted to methods for generating synthetic turbulent magnetic fields, affordably creating fields with parameter-controlled characteristic features of turbulence. In this context, the BxC toolkit was developed and validated against direct numerical simulations (DNSs) of isotropic turbulent magnetic fields. Here, we demonstrate novel extensions of BxC to generate realistic turbulent magnetic fields in a fast, controlled, geometric approach. First, we perform a parameter study to determine quantitative relations between the BxC input parameters and the desired characteristic features of the turbulent power spectrum, such as the extent of the inertial range, its spectral slope, and the injection and dissipation scale. Second, we introduce in the model a set of structured background magnetic fields, B 0, as a natural and more realistic extension to the purely isotropic turbulent fields. Third, we extend the model to include anisotropic turbulence properties in the generated fields. With all these extensions combined, our tool can quickly generate any desired structured magnetic field with controlled, anisotropic turbulent fluctuations, faster by orders of magnitude with respect to DNSs. These can be used, e.g., to provide initial conditions for DNSs or easily generate synthetic data for many astrophysical settings, all at otherwise unaffordable resolutions.


INTRODUCTION
Fluids, gases, and plasmas are often found in turbulent states of motion, making the study of turbulence a central topic in many research fields, from hydrodynamics to modern astrophysics (Biskamp 2003;Brandenburg & Lazarian 2013;Galtier 2016;Goedbloed et al. 2019).The importance of turbulence largely comes from the ubiquitous presence of such a state of motion.The present work focuses on (but is not restricted to) astrophysical applications, where turbulence is observed in numerous environments: in accretion disks and astrophysical jets, in solar/stellar atmospheres and winds, in molecular clouds or galaxies, etc. (Barnes 1979;Parker 1979;Balbus & Hawley 1998;Schekochihin & Cowley 2007;Beresnyak & Lazarian 2019).Nonetheless, this work builds on an issue that is shared by all turbulencerelated fields of study.The state-of-the-art approach to study turbulent states is to perform direct numerical daniela.maci@kuleuven.besimulations (DNSs), which are effective, yet often extremely expensive in terms of computational resources required.
With the aim of reducing the need of expensive numerical simulations, recent studies (Juneja et al. 1994;Cametti et al. 1998;Zimbardo et al. 2000;Ruffolo et al. 2006;Malara et al. 2016;Lübke et al. 2023Lübke et al. , 2024) ) have been focusing on the development of software that can generate synthetic data of turbulent quantities, e.g.magnetic fields.The general approach used in synthetic turbulence is to avoid solving physical equations numerically, using instead simplified models and algorithms that are able to mimic properties that are characteristic of turbulence.In this paper, we develop an alternative synthetic model based on the previously presented BxC code (Durrive et al. 2022), which is intended as a general, versatile model, not restricted to any specific application.BxC, which stands for magnetic field from multiplicative chaos, is a code fully implemented in Python that can rapidly and inexpensively produce data-cubes of turbulent magnetic fields of order ∼ 1000 3 points, and more.The present work is based on a double intent: having an easily customizable turbulent power spectrum and giving the users the possibility to generate fields that are closer to real astrophysical scenarios, by introducing physical characteristics such as anisotropy, while keeping the model as simple and efficient as possible.
In almost all synthetic models, great importance is given to the power spectrum of turbulent fields, which is typically expected to show a power law decay in a range of values that extends from the injection scale to the dissipation scale (i.e. in the inertial range).This concept was formalized by Kolmogorov's theory (Kolmogorov 1991), and it is now established in the literature as a distinctive characteristic of turbulence.In addition to the shape of the power spectrum, distinctive properties such as anisotropy and intermittency should be input-controlled in any synthetic model.The concept of anisotropy is strictly related to the presence of a guide field, which causes differences in the parallel and perpendicular directions in terms of scaling laws, and hence energy transfer mechanisms.Intermittency is a property related to the statistical distribution of spatial increments of the fields.In practice, in a turbulent field, one would expect the distribution to show heavily non-Gaussian tails for small increments in the magnetic field B, while tending to Gaussianity for large increments.BxC convincingly mimics and reproduces the intermittent and non-Gaussian character of a turbulent magnetic field, as demonstrated in Durrive et al. (2022), where fields generated by BxC have been successfully compared and validated against results from a magnetohydrodynamic (MHD) DNS.The two fields (the one generated by BxC and the one obtained with the DNS) show the same properties in terms of power spectrum, probability distribution function (PDF) of increments (both in the magnetic field B and the current density j), structure functions, and spectrum of exponents, hence proving the tubulent and intermittent character of the fields generated by BxC.In addition to showing power-law decaying power spectra, intermittency, and anisotropy, which are properties that may be more or less important according to the specific application for which the model is designed, synthetic algorithms can differ from each other in terms of spatial dimensionality (e.g.1D, 2D, or 3D) and on whether the fields are time-dependent or not.Although we primarily focus on static fields, Durrive et al. (2022) showed how BxC has also the potential of reproducing a time-dependent magnetic field evolution, by varying continuously the input control-parameters.
Given the many features to take into consideration, a variety of models have been proposed so far.In this context, two classes of synthetic models have attracted widespread attention: wavelet-based algorithms (e.g.Juneja et al. 1994;Cametti et al. 1998;Malara et al. 2016), which are very efficient in reproducing intermittency in the fields, and Fourier-based sampling procedures (e.g.Zimbardo et al. 2000;Ruffolo et al. 2006), in which the model can reproduce anisotropy but intermittency is not present.Recently, also a combination of both approaches has been proposed by Lübke et al. (2023Lübke et al. ( , 2024)).BxC does not belong to either of these categories, but, as suggested by the name itself, uses a completely alternative approach linked to the concept of multiplicative chaos (Kahane 2000a,b;Durrive et al. 2020), which is a wider research field on (hydrodynamic) turbulence (see Rhodes & Vargas 2014 for a review on the topic).Indeed, BxC has been inspired by recent developments in hydrodynamics, which suggest a relation between random fields and incompressible hydrodynamic turbulent fields (Chevillard et al. 2011).At the same time, BxC differs from this approach in terms of fundamental structures and visual aspects of the generated fields.Specifically, Chevillard et al. (2011) succeeds in incorporating in the model mathematical properties of turbulence, but lacks visual turbulent appearance in vortices and current sheets.BxC has been developed with the aim of reproducing also visual characteristics of turbulent magnetic fields, meaning that BxC-generated fields also have an actual "look-and-feel" resemblance to turbulent fields, such as hierarchically structured magnetic eddies and current sheets.This is because the actual non-linear transformation of fractional Gaussian fields in our toolkit builds in spiral shapes and multilayered current sheets, which are typical for all DNS realizations (see Section 2 for more details on the algorithm).So far, this aspect has been underexplored and it fundamentally distinguishes BxC from Chevillard et al. (2011), as well as most other synthetic models.Indeed, to the best of our knowledge, no other synthetic turbulence model allows for the fully customizable production of turbulent fields including both structures and higher order statistics at the same time.
Assuming the turbulent and intermittent character of the fields generated using BxC, in this paper we introduce several novel capabilities into BxC to construct more realistic turbulent magnetic fields.Section 2 briefly describes the important features of the BxC model, including recent modifications that have been been introduced compared to the original algorithm (Durrive et al. 2022).Sections 3 and 4 contain the results of this study.In Section 3 we present the results of a parameter study conducted on the controlparameters, which are varied independently from one another in order to assess and quantify their effects on the power spectrum.In Section 4 we generalize the model to reproduce more realistic scenarios.In particular the isotropic fully turbulent fields generated by BxC are extended to turbulent fields with a background structure (Section 4.1) and to fields featuring anisotropic turbulence (Section 4.2).Finally, the conclusions of this study are presented in Section 5.

METHODOLOGY
The BxC toolkit is a Python implemented code that generates data cubes of turbulent magnetic fields using a combination of analytical formulas and geometric constructions.The work presented here builds on the preexisting algorithm, implementing new features that render the generated fields even more realistic.In order to contextualize the new features that are now included in the model, we first briefly summarize the relevant aspects of the BxC algorithm, as presented in Durrive et al. (2022).
The basic idea behind the BxC model is to start from a white noise vector field and transform it in a strategic way, in order to obtain magnetic and current vector structures with statistical properties that are typical of a turbulent magnetic field.The white noise vector field is generated from a normal distribution with zero mean (µ = 0) and unit standard deviation (σ = 1), but both µ and σ are input parameters that can be user-defined.The same value of mean and standard deviation is used for all three vector components.As a first step, the white noise vector field is used to generate a fractional Gaussian field (FGF), R, which, in turn, is used to generate the final turbulent magnetic field B. In both these steps, the fields are transformed analytically using a modified version of Biot-Savart's law.Before it is used to generate B, R is subject to a non-linear, geometrically inspired transformation.In the original formulation, Biot-Savart's law relates the magnetic field B to the current density field j through the convolution: The modified versions implemented in the BxC model are, instead: (3) A clear difference between the standard Biot-Savart law and the ones implemented in the model is the replacement of the current density j in Eq. ( 2) with a white noise vector s, which generates the FGF, R. In Eq. (3), the current density is replaced by a non-linear transformation of R, which thereby introduces deviations from Guassianity and generates the turbulent magnetic field, B. We may denote this non-linear transformation as c = SR, where S is a geometrically controlled transcendental operator defined in terms of gradient and norm of R. The specific choice of S is inspired by two visual characteristic features of turbulent fields: first, the omnipresence of spiral-shaped current sheets swirling around through space; second, the appearance of intense sheets surrounded by more diffuse ones.The final formula for S is: where T i and T d are top-hat functions that effectively reproduce the intense and diffuse sheets respectively, and: Eq. ( 5) is a 3D generalization of an Archimedean spiral, in which the norm of the FGF |R| plays the role of radius and θ R (in Eq. ( 6)) is defined accordingly in terms of the gradient of |R| (i.e. it is interpreted as an angle).Eq. ( 4) defines a scalar field dependent on |R| which multiplies R component-wise, hence preserving the vectorial nature of c.
In practice, the BxC algorithm does not artificially build the magnetic field directly, but it generates the turbulent magnetic field starting from a carefully constructed "current density" field, c.This is the key aspect of the process that eventually yields structures in the fields.The actual current density vector can be computed afterwards from the relation j = ∇ × B. In addition to the replacement of the current density, other changes are applied in order to replace Eq. (1) with Eq. ( 2)-(3): the domain is limited to cubes of side L R and L B respectively, the power law behaviour is extended to the free exponents h R and h B , and finally regularization is performed through η R and η B to avoid singularities at the origin.This process naturally introduces in the code a set of input parameters that will be the focus of the parameter study in the next section.It is also worth mentioning that building the magnetic field using Eq.(3) has the intrinsic advantage of ensuring a divergence-free field, without additional constraints on c (see Appendix A).
The magnetic fields obtained using the model described above are fully isotropic by construction, since there is no preferential direction imposed at any point.
However, in many realistic astrophysical environments, magnetic turbulence is often found to develop anisotropy along some preferential direction, i.e. with respect to a background or guide magnetic field.With the purpose of reproducing turbulence in such environments, we detail here how the model is modified such that it allows for either structured background magnetic fields and/or the choice of a preferential direction.One of the main concerns in achieving this generalization is to retain the same algorithmic structure in order to maintain similar computational times and not overcomplicate the model itself.The most immediate way to achieve this direction-based variation of the fields is to exploit the input parameters of the model.In principle, almost all the input parameters present in the model can be adapted to become direction-dependent.However, the limited set of direction-dependent parameters discussed below was chosen upon consideration of what would physically render the fields anisotropic.In practice, having parameters that are direction-dependent means that instead of considering the same value in every direction, we now give as input three different values corresponding to the x−, y−, and z−direction respectively.In particular we will show (see Section 4.2) how anisotropy can already be controlled by introducing direction-dependent values of: • white noise standard deviation (σ); • integration region for R (L R ).
To recover isotropic fields with this new generalized approach, it is sufficient to set each parameter value equal in all directions, retrieving the original BxC algorithm (Durrive et al. 2022).

PARAMETER STUDY
One of the main advantages of BxC lies in the possibility to customize the generated fields not only from a visual point of view, but also from a statistical one.In previous work Durrive et al. (2022) already suggested this possibility, showing that variations in the input parameters reflect in variations in the statistical aspects of the fields.It was demonstrated how parameter changes influence multiple aspects of the generated fields in terms of spectral properties, but the original paper (Durrive et al. 2022) did not provide specific relations between the model's input parameters and statistical properties.Therefore, in this work we perform an in-depth, quantitative study on the influence of the input parameters on the final turbulent fields.
We vary each parameter independently from the others, in order to assess its impact on the power spectrum.Each parameter has been limited to values that produce reasonable turbulent fields both visually and statistically.Figure 1 shows a typical 3D turbulent field generated by BxC.The figures in the left panels are all related to the magnetic field B, while the ones on the right show the current density field j.The top panels show the norm of the fields, which display a fluid-like behaviour typical of turbulence.In the magnetic field picture, one can recognize vortex-like structures in which the intensity smoothly alternates from high to low and vice-versa.Similarly, the current density plot exhibits few high intensity thin sheets in the shape of spirals, surrounded by more diffuse ones.The same visual aspects are confirmed in the middle row, in which we show 2D cuts (one for each Cartesian plane) of isocontours corresponding to different intensity levels.In the left panel of the bottom row, we show magnetic field lines, which confirm that the field is isotropic and there is no preferential direction.The right panel of the bottom row, instead, shows 3D isocontours of the current density field, confirming that the current density is indeed distributed in sheet-like structures and that high intensity structures are very sparse and not volume-filling.
Details on the ranges considered for each parameter can be found in Table 1, together with the fixed values used for each parameter when varying the others (i.e.reference values).All parameters with length units are expressed in our 3D box units, which has a total side length of 1 and is centered around the origin.For all the runs, a grid resolution of N = 1024 3 is used.It is worth pointing out, already at this point, that the results of varying two of these parameters, h R and η R , will not be further discussed in this paper as our rigorous analysis has clearly shown that they do not affect the power spectrum in any significant way.The parameter study was performed fully on isotropic field, such as the one shown in Figure 1 where we took all reference values as in Table 1.Limiting the parameter study to isotropic fields allowed us to isolate the effects of the input parameters from variations caused by anisotropic features.
In order to consistently assess the effect of varying each input parameter, we first identify a set of characteristic features of the power spectrum, highlighted schematically in Figure 2.For a generic turbulent power spectrum we focus on three points of interest in the (k, P )-plane (i.e.wavenumber versus spectral power): the point at which the inertial range begins, the point at which it ends, and the maximum wavenumber (i.e. the point at which the power spectrum ends).The subscripts i, d, and max stand for injection, dissipation and maximum respectively, and are used to distinguish the   three different points.Moreover, the spectral slope in the inertial range, which we indicate with ζ, is an important feature that is taken into consideration in this analysis.Among the complete set of features shown in Figure 2, we selected particularly P i , k i , k d , and ζ.The analysis led us to identify a series of relations between input parameters and desired spectral feature.For each of these relations, a fit was performed for the purpose of user-friendly customization of the fields.Explicit equations for each fit can be found at the end of this section in Table 2 together with the relative maximum absolute and relative error.The input parameters L R and L B control the integration region of the FGF R and the magnetic field B respectively, as described by Eq. ( 2) and (3).It is intuitive to expect a relation between the size of these integration regions and the large scale structures of the fields.Focusing on the power spectrum features, the large scale structures are related to the injection scale, i.e. the point at which the power-law behavior of the power spectrum begins.In the context of our analysis, this means that the parameters L R and L B are expected to affect the quantities k i (scale) and P i (energy) as shown in    point is beyond the scope of this work, considering that the maximum relative error on the fit is ≈ 18%, and will be pursued elsewhere.However, we also considered the correlation length k c as an alternative quantity of interest that could be less affected by the domain discretization.The results of this analysis are shown in Figure 4.As expected, there are fewer overlappings when considering the correlation length.However, the relationship is no longer linear, but follows a Gaussian curve.The parameter for the fit, together with the maximum relative error, can be found in Table 2.
Figure 5 shows the results obtained for variations in the L B parameter.Also here, the top panel shows the power spectra colored according to different values of the parameter, while the same-color dashed horizontal line indicates the corresponding energy at the injection scale.As expected, also in this case we have a clear effect on the injection point: while the scale at which the inertial range begins is not changed, the energy related to the injection scale increases as we consider larger values of L B .The bottom panel of Figure 5 shows a quadratic fit for the relation between L B and P i .

Controlling the slope of the power spectrum through h B
The parameter h B appears in the exponent of Eq. (3), which then presents a spectrum with power-law like behaviour.Therefore, we expect such parameter to control the spectral exponent of the power-law.The results of varying the h B parameter confirm our expectations and are shown in Figure 6.The top panel shows the power spectra colored according to the parameter values.Here the dashed lines indicate the minimum and maximum slopes, corresponding to the lowest and highest value of h B respectively.Also in this case, a fit was performed in order to find a specific quantitative relation that describes the dependence of ζ on the parameter h B .Such relation is found to be quadratic, as shown in the bottom panel of Figure 6.Moreover, additional runs with higher values of h B have been performed (reaching h B = 4) which indicate that the slope saturates at a roughly ζ = 0.7, making h B = 3 an effective maximum value to be considered.
In the analysis for the L R and L B parameters, one can notice that the injection scale and spectral power are the only aspects of the power spectrum affected by the parameter change.On the other hand, when varying h B , changes in the slope are accompanied by changes related to the injection and dissipation scale as well.In particular, changing h B influences the spectral power at both the injection and the dissipation points, but the relative wavenumbers do not change.For instance, the injection spectral power decreases as the parameter h B increases, while the injection scale stays constants for all runs.The top panel of Figure 7 shows the fit performed to find the relation between h B and P i .As indicated in the figure as well, the fit for this relation was performed using only values h B ≥ 1.4.Initial fits were performed on the entire set of h B values, considering quadratic and exponential fitting curves.In both cases, the resulting fit was poor, amounting to a maximum relative error of ≈ 77 and ≈ 23, respectively.In addition to finding a proper fitting relation, we recall that visual aspects of the generated fields have also been taken into consideration.For small values of h B (i.e.h B = 1.0, 1.2), the fields lose the characteristic shape and structures that are typical of turbulence.This is clearly visible in the middle panel of Figure 7, in which a 3D box of the norm of B is shown for h B = 1.0.The field lacks the typical visual appearance of turbulent fields and it resembles more a Gaussian field, in contrast to the field shown in Figure 1.In this regard, the visual aspect is supported by a statistical analysis as well.The bottom panel of Figure 7 shows the comparison between the PDF of increments in B for h B = 1 and h B = 2.0 (i.e.our reference value).The distribution related to h B = 2.0 strongly deviates from a Gaussian distribution, which is in contrast to the behaviour exhibited for h B = 1, suggesting that the latter case is not intermittent.In addition, the extremely steep shape of the power spectrum would suggest that the field generated with h B = 1 is not turbulent.All the above considerations led us to exclude lower values of h B and consider only those values that generate mag- netic fields that can be classified as turbulent both from a visual and a statistical perspective.However, it is worth pointing out that a field with the same characteristics (i.e.turbulence mainly present at large scales) has been recently used by Pezzi & Blasi (2024) in the context of cosmic ray transport.Similarly, it can be noted that the analysis conducted on the turbulent character of fields generated by low h B values might also suggest the possibility of having an user-controllable level of intermittency.This aspect is left for future studies.

Controlling the dissipation scale and spectral power through η B
The last parameter considered in this study is η B , which is the regularization parameter in the construction of the magnetic field in Eq. ( 3).From a physical and mathematical point of view, η B determines the relative size of the small scales in our system, i.e. the dissipation scales.
The results of varying η B are shown in Figure 8.The top panel shows the power spectra for each value of the parameter considered.Here, the vertical dashed lines represent the dissipation scale and increasing the value of η B causes a shift toward larger scales.The fit shown in the bottom panel of Figure 8 clearly shows that the relationship between η B and k d is very well described by a quadratic curve.Looking at the power spectra overview, one can easily notice that the dissipation scale is not the only aspect of the power spectrum that is affected since we clearly see that η B affects the spectral slope as well as the injection spectral power, while still leaving the injection scale k i unchanged.Hence, η B effectively controls both dissipation scale k d and energy at once.Although the feature P d is not actively discussed in this work, from the top panel of Figure 8 it is visible that as the dissipation scale k d moves toward larger scales, the dissipation spectral power P d decreases.
As we are mostly interested in spectral slope and injection scale aspects, we performed a fit for both the slope ζ and injection spectral power level P i variations.The results are shown in Figure 9.The top panel shows the relation between η B -ζ, which, with a maximum relative error of 0.064, is in good agreement with a linear fit.The bottom panel shows the relation between η B and P i , which in this case is quadratic.

Interplay between parameters
The analysis conducted in the previous sections clearly shows that there is an interplay between parameters, especially with regard to h B and η B .Concerning the Table 2. List of all fitting curves obtained for each characteristic feature of the power spectrum, with related maximum absolute and relative error.The notation (par; * ), where par=LB, hB, or ηB, is used when a feature does not depend on "par" only, but the other parameters on which it depends are kept constant to the reference value.

Fit
Relative error (max) ki(LR) ≈ −68.5LR + 16.4 0.18 effects of changing η B , it was previously mentioned that one of the objectives of this study is to be able to control the dissipation scale without affecting either the slope or the injection spectral power of a given power spectrum.
The analysis above showed that this is not achievable by only changing the parameter η B .However, it is still possible to obtain the same result using a combination of η B and h B , since we found that the effects on the slope and injection spectral power caused by these two parameters are opposite with respect to each other.This means that changing η B , while also changing h B at the same time so that it forces the slope to be constant to the original value, will indeed result in a power spectrum that varies only in the dissipation scale.For this reason, a series of runs was performed varying both η B and h B .
Figure 10 shows the fitting curves of the relationship between h B and ζ for different values of η B .All the curves are quadratic, but there is a clear pattern dependent on η B : as the latter increases, ζ tends to smaller values.In order to obtain an expression that could relate the slope to both parameters, fits have been performed on the coefficients of the h B -ζ relations.Eventually, we obtain a relation of the form: in which the coefficients A(η B ), B(η B ), and C(η B ) are all quadratic expressions in η B .

Background magnetic fields
One of the main purposes of this work is to extend BxC to reproduce more realistic scenarios, allowing for structures more complicated than purely isotropic turbulent magnetic fields.Indeed, in most actual physical situations turbulence develops on top of an existing, background magnetic field.In this section, a novel method to generate turbulent magnetic fields with arbitrary background structures is presented.
Keeping in mind possible applications of BxC to astrophysics, a set of demonstrative background magnetic field topologies have been implemented in the BxC toolkit.An arbitrary number of additional, user-defined background fields can be easily added to the code.In order of increasing complexity, we now allow for a background uniform magnetic field, a magnetic arcade (with or without shear) (Terradas et al. 2015), and a cylindrical flux tube with a force-free Gold-Hoyle model (Vandas & Romashets 2017).The topology of the implemented background field is shown in the left column of Figure 11.The uniform background magnetic field, shown in the upper left panel of Figure 11, is very simple but represents a variety of physical scenarios that are extremely relevant in the context of astrophysical plasmas, since turbulence is often found and studied in relation to an existing guide field (e.g.Roytershteyn et al. 2015;Dong et al. 2022).The arcade and the flux-rope structure have been chosen as possible models of interest, since the divergence-free nature of magnetic fields implies that flux tubes are cornerstone ingredients in any 3D setting, while arcades occur in many plasma setups from laboratory, over solar, to astrophysical scenarios (e.g.Cheng & Cjoe 2001;Xia et al. 2012;Ryutova 2018).
The general principle on which this method is based is to first generate the background magnetic field B 0 , inde- pendently of the turbulent, BxC-generated field B turb , with the only restriction that both of them have to be defined on the same grid and that the background field must (analytically) be divergence-free.Eventually these two fields, B 0 and B turb , are summed up point-wise, obtaining a turbulent magnetic field B str with a background structure.In the implemented model, both fields are normalized to the maximum value of their norm, leading to fields with maximum value of order unity.At the same time, the purely turbulent magnetic field B turb is built using a filtering function in order to ensure that the background structure does not lose its characteristic geometry due to excessively strong turbulent fluctuations.In our examples, such a filtering function is chosen to be the norm of the background magnetic   field Γ = |B 0 |, hence the intensity of the turbulent field is modulated according to the pointwise value of |B 0 |.This yields turbulent fields with lower intensity at the points in which the background magnetic field is weaker, thus allowing for the field topology to still be relevant.The filtering function Γ is implemented inside the modified Biot-Savart law for the magnetic field, meaning that Eq. (3) becomes: Implementing the filtering function inside the modified Biot-Savart law is an important step, necessary to preserve the divergence-free property of the fields.It can be proven (see Appendix A) that the divergencefree property of the (modified) Biot-Savart law is independent of c, due to properties of the convolution.In practice, this allows us to multiply c by any scalar field and still generate fields that are analytically divergencefree.Eventually, the final formula to obtain a turbulent magnetic field with the structure of a background magnetic field is: where γ here is simply a constant value that gives the user the possibility to regulate the intensity of turbulence with respect to the background magnetic field in addition to (and independently from) the scaling/normalization performed inside the integral.In contrast to the pointwise scaling performed by Γ, which deals with the relative intensity of the turbulent field within the box itself, γ scales the entire turbulent field as to regulate its intensity compared to the background field.Figure 11 shows turbulent magnetic fields with each background structure implemented in the model, for different values of γ.The middle column shows lightly turbulent fields, while the right column shows highly turbulent fields, obtained for low and high values of γ respectively.

Anisotropy
As described in Section 2, our BxC algorithm can use direction-dependent model parameters for generating anisotropic fields.The introduction of anisotropy in the fields is an important feature in order to have realistic synthetic turbulent fields that users can customize according to the application.In numerous astrophysical plasmas, turbulence is found to develop anisotropically, meaning that there exists a preferential direction with respect to a background (or guide) field, modifying energy transfer along parallel and perpendicular directions from large to small scales.Such anisotropy is still heavily researched in magnetic field turbulence (e.g.Montgomery & Turner 1981;Matthaeus et al. 1996;Müller et al. 2003;Horbury et al. 2008) and precise scaling laws have not been determined in general.The most well-studied case is known as "critical balance" (CB) (Sridhar & Goldreich 1994;Goldreich & Sridhar 1995, 1997), for which scaling laws for parallel and perpendicular components of the field have been experimentally determined and validated with physical models (see e.g.Oughton & Matthaeus 2020 for an overview of CB theory and its applications).For this reason, our BxC tool has been extended to provide the freedom of obtaining different scalings in the parallel and perpendicular directions.However, the feature of anisotropy will be discussed here relatively qualitatively, as currently more dedicated work is required to quantify the precise effect of the numerical value of all parameters.
The approach used to obtain anisotropic fields is the following: first, we identified the parameters in the algorithm that could potentially induce differences in the resulting fields based on direction.Second, these parameters were generalized from a single value for all directions to a direction-dependent three-value implementation.Lastly, different runs were made and visualized to test that the generated fields were effectively anisotropic.The analysis on the generated field was conducted both from a visual point of view and a statistical one.Visually, the aim was to reproduce turbulent structures, particularly current sheets, that are aligned with a preferential direction.This is often found in simulations that reproduce turbulent magnetic fields in the presence of a guide field (e.g.Roytershteyn et al. 2015;Dong et al. 2022).For the statistical analysis, we compute the power spectrum in its parallel and perpendicular component, where the parallel direction is assumed to be along ẑ.After computing the 3D power spectrum P (k x , k y , k z ), we integrate in the (k x -k y ) plane summing over concentric shells as follows: where k max is defined as the miminum between max(k x ) and max(k y ).One can then directly visualize the 2D spectrum.The result of this study are shown in Figure 12, for different cases considered, together with the isotropic case (first column) for easier comparison.
The first parameter to be considered for controlling and introducing anisotropy is the standard deviation of the white noise vector.For the isotropic case, this vector field was initialized from a normal distribution of zero mean (µ = 0) and unit standard deviation (σ = 1) in all three directions.Allowing σ to have different values for x, y, and z has the desired effect of inducing anisotropy.In particular, wanting to mimic the presence of a background field in the z-direction, σ z has been increased compared to σ x and σ y , which are both set to unity.For what concerns the resulting white noise vector field, these different values of the standard deviation generate a field that has more energy in its z-component.However, this intrinsic difference between the three directions propagates throughout the algorithm, and eventually translates in anisotropy between the parallel (zaxis) and perpendicular (x-y plane) direction.The second row of Figure 12 shows the results of the analysis conducted on a field generated using σ = [1, 1, 10].The left panel shows the 1D reduced power spectrum, while middle column shows the 2D power spectrum, both featuring parallel and perpendicular components.The black lines in the 2D power spectrum correspond to isocontours of spectral power levels.The asymmetry of the contours (i.e.slight elongation in the parallel direction) is indeed an indication that the spectral power and the associated power law decay differ between the two directions.The same spectral power level corresponds to larger wavenumbers in the parallel direction than in the perpendicular one.This difference can also be noticed in the 1D power spectra.The right panel shows the current density isocontours for |j| = 0.35.This figure gives visual indications that the resulting field features a preferential direction.In fact, the current density sheets are predominantly extended along the z-axis, contrary to what was shown for the isotropic case.In the preliminary analysis that was conducted for this case, different values of σ z have been analyzed and compared to each other.Here we show the case that produces the maximum level of anisotropy that was observed.
The second parameter considered to construct anisotropic fields is the integration region of the Biot-Savart law from Eq. (3).Although the integration regions for the construction of both the FGF and the magnetic field have been considered in the analysis, only using a vector-valued parameter L R led to more significant results in terms of anisotropic features.Also for this case, as a preliminary survey, we considered different combinations of L R,x , L R,y , and L R,z .Here we show the two most anisotropic fields observed.The results are shown in the third and fourth row of Figure 12 for L R = [0.4,0.4, 0.01] and L R = [0.01,0.01, 0.4], respectively.Recall that we compute the fields on a box of unit length, which makes L R,i = 0.4 an effective maximum value to use in this case.Increasing L R,x and L R,y produces results that are similar to the σ z case: the 1D reduced spectra and the 2D spectrum isocontours elongated in the k ⊥ -direction indicate different spectral power distribution and transfer, j isocontours clearly elongated in the z-direction show the existence of a preferential direction.The results of increasing the L R,z parameters are also similar, but opposite: 1D power spectra shows different powers and spectral slopes, the 2D spectrum isocontours are elongated in the k ∥ -direction and the j isocontours are elongated in the x-y plane.Comparing these results to the ones obtained for σ z , one can notice an increased level of anisotropy in the fields.The 2D power spectrum is more skewed and the current density isocontours follow the preferential direction more closely, assuming an almost bi-dimensional character.However, we also conducted a preliminary analysis concerning the variance anisotropy.Such quantity is computed as: where b ∥ and b ⊥ indicate the fluctuations from the mean value in the parallel and perpendicular direction, respectively, and ⟨⟩ indicates a spatial average.From this analysis we concluded that the most efficient way to introduce variance anisotropy is by varying σ z .Indeed, Σ computed for the case of L R = [0.4,0.4, 0.01] and L R = [0.01,0.01, 0.4] shows slight variation with respect to the isotropic case, while the variance anisotrpy measured for the σ z = 10 case, leads to a difference of 81%.
As it was mentioned at the beginning of this section, this analysis was performed on a rather qualitative level.More in-depth studies on the input parameters should be performed in order to quantify their effect on the resulting anisotropy (both spectral and variance) in the fields.At the same time, this preliminary study effectively shows two techniques through which the BxC generated fields can be rendered anisotropic.

CONCLUSIONS AND OUTLOOK
Synthetic turbulence models represent a promising alternative and/or support to numerical simulations, having the advantage of requiring little computational resources and time.In this work, we have significantly developed our BxC toolkit, improving control over the existing generated turbulent magnetic fields and implementing novel features that are able to render more realistic fields.
The first part of this work focuses on providing the user with precise information on how to control the power spectrum of the generated fields.To do so, an in-depth parameter study has been performed on the parameters that define the modified Biot-Savart laws.After the identification of possible power spectrum features of interest, it has been shown how variations in a restricted set of input parameters can affect and control such features.In particular, the integration region parameter L R and L B control the injection scale in terms of wavenumber and spectral energy respectively.The free power law exponent h B controls the slope of the power spectrum, allowing for a wide range of spectral exponents.The regularization parameter η B controls the dissipation scale.Contrary to what happens at the injection scale, in which the two parameters control separately the two coordinates of the injection point, the regularization parameter affect the dissipation scale both in terms of wavenumber and energy level.This produces a change in the slope of the spectrum, hence an entanglement with the parameter h B .The connection between the two parameters, however, is an advantage that one can exploit in order to produce a change in the dissipation scale only with respect to the wavenumber.This can easily be done appropriately combining changes in both parameters.The parameter study has been conducted in a systematic way, and quantitative relations between parameters and power spectrum features are provided for each parameter considered in the analysis.Table 2 contains all the relations found through the analysis and represents an important tool to generate isotropic magnetic fields with given spectral properties.
The second part of this study aimed at further improving the generated fields in terms of realistic features that can be implemented.In particular, the development focused on two main aspects: generating turbulent fields that are not only purely turbulent, but have a background structure, and introducing anisotropy in the fields.For the first point, it has been shown that it is possible to obtain turbulent fields with a background structure by summing an appropriately constructed purely turbulent field to the desired background magnetic field.Here, "appropriately" means that the turbulent field should be modulated in intensity based on the background field, to avoid the disruption of structures where the background field is weak.It is also suggested that the summing process should be modulated by a constant value, which gives users an additional degree of freedom to regulate the intensity of turbulence with respect to the background field.
Finally, the synthetic model was generalized in order to introduce anisotropy in the fields.With the purpose of keeping the model as simple and efficient as possible, anisotropy has been introduced by changing a set of chosen parameters from single-valued to directiondependent values.Two methods to introduce anisotropy in the fields have been investigated: increasing the values of the white noise standard deviation in one direction and using different values of integration region in the three directions.Both methods successfully reproduced skewed 2D power spectra and current density isocontours extending along a preferential direction.Among these two, different sizes of the integration region produced more significant and effective results.The analysis presented in this work was conducted on a qualitative level and the quantification of the precise effects of these parameters is left for further studies.
The quantification of the effects of the input parameters and the novel features introduced in this work represent an important milestone in the development and application of the BxC code.As of now we can inexpensively generate large data cubes of turbulent magnetic field and turbulent magnetic field structure, with user-controlled power spectrum features.This paves the way for possible astrophysical applications, such as the study of cosmic-ray propagation (e.g.Pucci et al. 2016;Dundovic et al. 2020;Reichherzer et al. 2020;Kuhlen et al. 2022 or see Mertsch (2020) for a review of test particle simulations of cosmic ray in synthetic turbulent fields).Future development options include usercontrolled higher order statistics (i.e.intermittency), as well as generalizing the model to both flow and magnetic field realizations with combined turbulent properties, which would make BxC suitable for many more applications (e.g.study particle acceleration, provide initial conditions for DNSs).

Figure 2 .
Figure 2. Schematic representation of a turbulent power spectrum and its characteristic features.ki, Pi (blue lines) are the wavenumber and spectral power coordinates of the injection scale; k d , P d (red lines) are the wavenumber and spectral power coordinates of the dissipation scale; kmax, Pmax (yellow lines) are the wave number and spectral power coordinates of the maximum point.The slope in the inertial range (ζ) is indicated in green.

Figure 3 .
Figure 3. Effects of varying the LR parameter on ki.Top panel: power spectrum for each value of LR; dashed vertical lines correspond to the (variable) injection scale.Bottom panel: fit performed that shows a linear relation between LR and ki.

Figure 4 .
Figure 4. Effects of varying the LR parameter on kc.Top panel: autocorrelation function for each value of LR.Bottom panel: fit performed that shows a Gaussian relation between LR and kc.

Figure 3
Figure3shows a parameter scan obtained by varying L R between 0.05 and 0.2.The top panel shows the power spectra for each value of L R considered in the analysis.Here, the dashed vertical lines indicate the injection scale k i and the coloring scheme is done according to the parameter values.As expected, we can control the beginning of the inertial range, causing a shift to the left (i.e. to larger scales) as we increase the parameter L R .A fit performed to quantify the influence of this parameter on the injection scale is shown in the bottom panel of Figure3.The relationship between L R and k i is described accurately by a linear function.It can be noted in this figure that the same k i might correspond to multiple consecutive values of L R .This aspect is related to the discretization of the domain, and the fact that the injection scale is generally found at large scales, where fewer wavenumber points are sampled.Exploring this

Figure 5 .
Figure 5. Results for varying the LB parameter.Top panel: power spectrum for each value of LB; dashed horizontal lines correspond to the (variable) injection spectral power.Bottom panel: fits performed to find a quantitative relation between LB and Pi.

Figure 6 .
Figure 6.Results for varying the hB parameter.Top panel: power spectrum for each value of hB; dashed lines indicate the slope given by the minimum and maximum spectral index.Bottom panel: fit performed to find a quantitative relation between hB and ζ.Both plots are colored according to the values of the parameter.

Figure 7 .
Figure 7. Top panel: fit performed to find a quantitative relation between hB and Pi, which is a secondary effect of varying hB.Middle panel: 3D norm of B generated for hB = 1.Bottom panel: comparison between PDF of increments in B for hB = 1 (pink scale) and hB = 2 (gray scale) for different lags.

Figure 8 .
Figure 8. Results for varying the ηB parameter.Top panel: power spectrum for each value of ηB; dashed vertical lines correspond to the (variable) dissipation scale.Bottom panel: fit performed to find a quantitative relation between ηB and k d .Both plots are colored according to the values of the parameter.

Figure 9 .
Figure 9. Fits performed to find quantitative relations for the secondary effects of changing ηB.Top panel: relation between ηB and ζ.Bottom panel: relation between ηB and Pi.

Figure 10 .
Figure 10.Fitting curves describing the relation between hB and ζ, for different values of ηB.

Figure 11 .
Figure 11.Background magnetic field structures implemented in the model and their turbulent realizations for different intensities of turbulence.Top row: uniform magnetic field; middle row: arcade with shear (Terradas et al. 2015); bottom row: cylindrical flux tube with a force-free Gold-Hoyle model (Vandas & Romashets 2017).Left column: background structures; middle column: lightly turbulent realizations; right column: highly turbulent realizations.

Table 1 .
List of input parameters considered in the study together with their reference values and range of variation.LR, LB, ηr, and ηB are all expressed in box units; hr and hB are dimensionless quantities.