Chaotic behavior of Eulerian MHD turbulence

# Chaotic behavior of Eulerian MHD turbulence

Richard D.J.G. Ho    Arjun Berera    Daniel Clark SUPA, School of Physics and Astronomy, University of Edinburgh,
JCMB, King’s Buildings, Peter Guthrie Tait Road, Edinburgh EH9 3FD, United Kingdom
July 19, 2019
###### Abstract

We study the chaos of a turbulent conducting fluid using direct numerical simulation in the Eulerian frame. We predict that the Lyapunov exponent, which measures the exponential separation of initially close solutions of the magnetohydrodynamic equations, is proportional to the inverse of the Kolmogorov microscale time and also obtain new results for this relation in hydrodynamic turbulence, specifically deriving a previously unknown co-efficient. These predictions agree with simulation results. The simulations also show a diminution of chaos from the introduction of magnetic helicity, which is expected to be eliminated at maximum helicity. Linear growth of the difference between fields was recently found in hydrodynamics and we find here that it extends to the magnetic and velocity fields, with growth rates dependent on the dissipation rate of the relevant field. We infer that the chaos in the system is totally dominated by the velocity field and connect this work to real magnetic systems such as solar weather and confined plasmas.

## I Introduction

Turbulence displays chaotic dynamics. A small change in initial conditions will result in a large difference in the state at later times, with this difference growing exponentially. This exponential growth of error puts limits on the predictability of the system, and understanding these limits is important for forecasting the evolution of fluids governed by turbulence. Turbulence, or its underlying equations, is interesting from a dynamical systems point of view. In this context, it is just another dynamical system of which we wish to know the chaotic properties.

In homogeneous isotropic turbulence (HIT), recent results have shown a relationship between the Reynolds number, Re, and the maximal Lyapunov exponent, , for chaos in a Eulerian (considering trajectories in phase space) description Berera and Ho (2018); Boffetta and Musacchio (2017) consistent with theoretical predictions by Ruelle Ruelle (1979); Crisanti et al. (1993). In these simulations, two initially close fields were evolved concurrently and their difference quantified. Surprisingly, a limit was found on the growth of this difference which is proportional to the dissipation rate. In this paper the analysis is extended to magneto-hydrodynamic (MHD) turbulence, also known as hydromagnetic turbulence. First there is a derivation of the predicted behavior for Navier-Stokes turbulence which is somewhat different to that done previously Ruelle (1979); Crisanti et al. (1993). Here a new result is obtained, namely we explicitly find the prefactor for the relation of Lyapunov exponent and Kolmogorov microscale time. Afterwards, this analysis is applied to MHD.

The evolution of an uncharged fluid is described by the Navier-Stokes equation, whilst the evolution of an electrically conducting fluid is described by the MHD equations. The incompressible MHD equations are

 ∂tu =−1ρ∇P−(u⋅∇)u+νΔu+1ρ(∇×b)×b , (1) ∂tb =(b⋅∇)u−(u⋅∇)b+ηΔb , ∇⋅u=∇⋅b=0 , (2)

with velocity field , magnetic field , density , pressure , viscosity , and magnetic diffusivity . Like the Navier-Stokes equations which they modify, they also display turbulence Biskamp (2003). In modifying the Navier-Stokes equation, new ideal invariants are introduced. One of the striking changes from this modification is that the conservation of one of these invariants, magnetic helicity, introduces an inverse transfer of energy from the small to large scale Brandenburg and Subramanian (2005).

The literature investigating the intersection of MHD and chaos is sparse. The relative dispersion of charged particles is predicted to grow exponentially for intermediate times when they are initially close Misguich and Balescu (1982); Misguich et al. (1987). The technique of extracting a Lyapunov-exponent spectrum from a time series Wolf et al. (1985) has been applied to experimental data from an undriven plasma system Huang et al. (1994). Use of this technique reveals a transition from quasiperiodicity to chaos and allows the evaluation of the Lyapunov spectrum. In solar physics, a low dimensional attractor is suggested to be responsible for the data seen in pulsation events of solar radio emissions Kurths and Herzel (1986). Scalar models of turbulence applied to MHD with Pr have found that the maximal Lyapunov exponent obeys Grappin et al. (1986). Later results for turbulent models of Navier-Stokes found a similar scaling Yamada and Ohkitani (1987), which is roughly the scaling found in DNS Berera and Ho (2018). This is the first analysis of chaos in MHD turbulence in an Eulerian context.

The way helicity affects chaotic behaviour has also been investigated. Numerical simulations of ABC flow and the non-linear dynamo problem have indicated that an increase in helicity results in a diminution of Lagrangian (considering particles in real space) chaos in MHD Zienicke et al. (1998). Simulations have also shown that an increase in the strength of the magnetic field corresponds to a diminution of chaos in helically forced ABC flows Rempel et al. (2012, 2013). In ABC flows, an increase of magnetic strength may be due to an increase of helicity in the system, so these simulations are consistent with an increase in helicity causing a diminution of chaos. Physical experiment has shown a similar diminution of chaos with the increase of helicity in reverse field pinch plasmas Escande et al. (2000).

The relationship between MHD and chaos is interesting from the viewpoint of predictability. Just as it is important to quantify the predictability of weather forecasts to understand the time horizons over which a prediction can be considered accurate, it is important to understand the predictability of forecasts where the governing equations are the MHD equations, such as in space weather Mirmomeni and Lucas (2009), solar physics Karak and Nandy (2012), and high latitude ground magnetic fields Weigel et al. (2003). For instance, much effort is put into understanding solar flares, a phenomena governed by the MHD equations, because of the damage they can cause to artificial satellites and thus global communications Amari et al. (2014). Whilst HIT is an ideal description of turbulence, it nonetheless should represent the behaviour of a turbulent system far from boundaries or at small scales. The work here may have practical use in these situations but also in tokamak reactors and fusion research where boundaries are present.

## Ii Theory for Navier-Stokes turbulence

In a chaotic system, two states which initially differ by (these can either be particle positions in which case we have a Lagrangian description, or a measure of the difference between two fields in which case we have an Eulerian description), the separation will grow as . In this case is the maximal Lyapunov exponent. According to the theory of Ruelle Ruelle (1979), the maximal Lyapunov exponent for Navier-Stokes turbulence is given by

 λ∼1τ , (3)

where is the Kolmogorov microscale time and the dissipation. The argument used by Ruelle requires the existence of an inertial range in order to justify a characteristic exponent, an inverse timescale, which is dependent only on the dissipation. In hydrodynamic turbulence this relation becomes

 λ=1T0Reα , (4)

where is the large eddy turnover time and Re is the Reynolds number. The Kolmogorov theory predicts that Crisanti et al. (1993); Kolmogorov (1941). Intermittency corrections predict that . However, DNS results have shown that, in a Eulerian description, Berera and Ho (2018); Boffetta and Musacchio (2017); Mohan et al. (2017) whilst in a Lagrangian description Biferale et al. (2005). These DNS results also show a corresponding behaviour for which either rises (for ) or falls (for ) with Re.

The dissipation anomaly implies that, if the relation is accurate, then it is not possible to have . Considering the kinetic energy equation

 ε∼ν∫VSijSij dV (5)

where is the strain rate tensor, and does not become zero as Re . This is because whilst , the vortex lines become more and more stretched, generating ever smaller scales. For fixed large scale parameters and , whilst varying the viscosity, Re and so if we define such that

 ν∫VSijSij dV∼Reβ , (6)

then otherwise the dissipation will vanish asymptotically, which is not the case due to the dissipation anomaly. This means that the Komogorov time obeys

 1τ2∼εν∼∫VSijSij dV∼Re2α , (7)

implying that and . The fact that this is not the case in a Lagrangian sense means that the Lyapunov exponent cannot be proportional to .

The argument used by Ruelle to predict the Lyapunov exponent is that it should be proportional to the largest inverse time, and is the shortest time scale in the system. However, in MHD there are many other time scales which may be smaller than . Before turning to MHD, here we re-derive the relation of Ruelle for hydrodynamics to decide which of the turbulent time scales is relevant.

### ii.1 The maximal Lyapunov exponent of tracer particles in turbulence

Consider two passive tracer particles at points and separated by . After a short time they will be separated by . If the flow is chaotic then

 |δxδt|=|r|exp(λδt)=|xb−xa|exp(λδt) , (8)

where is the maximal Lyapunov exponent. This is usually defined in a large time limit, but is valid provided our initial perturbation is small. If the velocity is at point and at point then

 δxδt=r+(ub−ua)δt . (9)

By squaring both sides

 (r⋅r)exp(2λδt)=(r⋅r) +2r⋅(ub−ua)δt+(ub−ua)⋅(ub−ua)δt2 . (10)

Expanding at small and neglecting terms of gives

 (r⋅r)λ =r⋅(ub−ua) (11) λ =r⋅(ub−ua)(r⋅r) . (12)

The velocities and can be split into their longitudinal, , and normal, , components with respect to their displacement , which has magnitude . In this case and . Thus . This analysis is done for a fixed distance . Considering the Lyapunov exponent again

 λ=uLb−uLar , (13)

which in principle is dependent on position and time, and so a global averaging must be taken. This means

 λ2 =uLbuLb−2uLauLb+uLauLar2 (14) ⟨λ2⟩ =⟨uLbuLb⟩−2⟨uLauLb⟩+⟨uLauLa⟩r2 . (15)

Now recall the definition of the longitudinal correlation function Batchelor (1953), which obeys

 CLL=f(r)u2 , (16)

where is the rms velocity and is the dimensionless longitudinal correlation function, which is only a function of in HIT. Under the Taylor hypothesis

 f(r)≃1−r22l2 , (17)

where is the Taylor microscale. Now

 ⟨λ2⟩ =2u2−2CLL(r)r2=u2l2 , (18)

which is a prediction of the mean squared Lyapunov exponent. A very similar derivation to this is done by de Devitiis de Divitiis (2018), where it is shown that . This gives a prediction for the Lyapunov exponent for pairs of tracer particles, dropping the angled bracket notation, of

 λ=u2l=12√15τ , (19)

where the last equality comes from the relation , which was not noted by de Devitiis. Whilst this is the same functional dependence as found Ruelle, and whilst has been noted before, this is the first time that this explicit value for the prefactor has been pointed out in the literature, to the knowledge of the authors. The argument by Ruelle Ruelle (1979) rests on different assumptions than our own and never specifies homogeneous isotropic turbulence. The arguments here give an equivalent result.

This prefactor is an explicit prediction. We note that and this compares well with the relation found in DNS for tracer particles in Biferale et al. (2005) of . Similarly, comparing it to the theory of two-particle dispersion in the dissipation subrange, , with and being the Batchelor constant. Our prediction has , which is similar to simulations and experiment Salazar and Collins (2009); Yeung and Pope (1989); Guala et al. (2007). However, this method should be applicable to other turbulent configurations which have different functional dependencies for the longitudinal correlation functions. As well, this still puts a bound on the prediction for the Lyapunov exponent given that variance is positive definite, i.e. .

### ii.2 Eulerian and Lagrangian maximal Lyapunov exponents

The above derivation works for in a Lagrangian sense, which is labelled for the current discussion. It can be related to the Eulerian Lyapunov exponent .

The velocity of a particle , , is equal to the rate of change of its position . A small perturbation to the velocity has

 ui+δui=∂tri+δui≡∂tri+∂tδri . (20)

The statistics of the small displacement should be similar to those of two particles initially separated by a small distance under the conditions of homogeneity and isotropy. This is because, after letting time advance slightly, the two realisations of the same particle are now separated by a small distance and should themselves have negligible effect on the rest of the flow.

Using the language of momentum space, we have introduced a perturbation to the momentum space of one of the particles and said this is equivalent to a perturbation in physical space of two particles.

The difference between the two velocities of the one particle in momentum space can be labelled as and the separation between the two particles as . Thus

 δv ∼∂tδr (21) ∂t(δv)2 ∼2˙δr⋅¨δr=13∂3t(δr)2−23...δr⋅δr . (22)

The last equality is just an identity. In this way we relate the growth of a perturbation between one particle in two different systems to the growth of a perturbation between two particles in one system.

What was done is the same as using the Eulerian to Lagrangian transformation. The Eulerian description of fluid flow assigns a velocity at each point in space to the flow. This can be denoted , which is the velocity at position at time . The Lagrangian description of fluid flow assigns a velocity to each parcel of the fluid which is then tracked. This is denoted , which is the position of a parcel which was at position at time 0, i.e . These have the following relation

 u(r(x0,t),t)=∂tr(x0,t) , (23)

where is the partial derivative with respect to time.

Note that

 ∂tEd =12∫∂t(δv)2dx (24) ∼16∂3t∫(δr)2dx−13∫...δr⋅δrdx . (25)

So , whilst . The second integral should thus vanish unless there is some process which provides jerk preferentially in a fixed direction. Even in a system with a mean flow, this should not be the case. It may be related to a closure approximation of the system. As such,

 ∂tEd∼∂3t⟨δr2⟩ , (26)

which relates the rate of change of the error and the change in time of the mean squared displacement, , between particles. For instance, if grows exponentially, then so will , and with the same exponent.

The Lagrangian Lyapunov exponent is defined by and so . Thus, if , then and so . This assumption also requires homogeneity, isotropy, and that the statistics of the field are not significantly changed by introducing the perturbation, that is to say there is a statistically steady flow. Indeed, looking at the supplemental material of Berera and Ho (2018), for hydrodynamic turbulence it is found that . At very large , this is very close to , which is our prediction of . As well, in other DNS and large-eddy simulations a relation of has also been found Nastac et al. (2017).

Previous predictions for have not distinguished between the Lagrangian frame and Eulerian frame for turbulence, such as the argument by Ruelle which makes no such distinction Ruelle (1979). We have just argued that in homogeneous isotropic turbulence. However, if Re, there are different values of in a Lagrangian and Eulerian sense. Specifically, for Lagrangian tracers Biferale et al. (2005) whilst for Eulerian trajectories Berera and Ho (2018); Boffetta and Musacchio (2017). We believe this may be due to the averaging procedure. The Lagrangian Lyapunov exponent is not uniform through space despite homogeneity and its distribution will be heavily affected by the strength of the turbulence. In this way, a smaller Reynolds number would allow larger patches of low Lagrangian chaos, which could reduce the average Lyapunov exponent. The difference may also come from the dependence of the relation on Reynolds number.

## Iii Theory for λ in MHD

The previous theory has been applied to Navier-Stokes turbulence and finds new relations. We now look at how it applies to MHD. Since , the subscript notation will be dropped. In MHD, the charged particles should follow the velocity evolution equation Eq. (1) passively. In MHD the functional dependence of is expected to change, but if it does not, then for MHD should be the same as for hydrodynamics and so it would obey Eq. (19). In any case, it should only depend on the longitudinal correlation function of the velocity field. The longitudinal correlation function for the velocity field in MHD is defined , specifically

 ⟨λ2⟩ =2u2−2CLL,MHD(r)r2 . (27)

If the previous assumptions are true, then this prediction follows. Using the above arguments, we predict that the Lyapunov exponent for the particles should also hold for the velocity field itself.

Measuring the Lyapunov exponent is also an implicit measure of , and thus the structure functions, in MHD. The fact that the data generated is more noisy would also imply that the longitudinal correlation functions in MHD are also more noisy than in hydrodynamics. Because there is little theory on the form of , we use the assumptions from hydrodynamics. If we accept the assumptions which lead Taylor to his hypothesis about the correlation functions, namely a prediction that the transverse correlation functions behave as a parabola, then the exact same relation should also hold for MHD turbulence. This assumption should be most tenable in the case of zero magnetic and cross helicity. In effect we assume that . Thus we say that the Ruelle prediction that , where is the kinetic dissipation, should also hold for MHD in a Eulerian sense. This prediction that is also backed up by the findings of scalar turbulence models mentioned previously Grappin et al. (1986); Yamada and Ohkitani (1987).

This extension only holds for the velocity part of the MHD equations, from this argument no prediction can be made about the chaotic nature of the magnetic field. However, we note that there is only an indirect non-linearity in the equation for the evolution of the magnetic field (via the velocity field) rather than the direct non-linearity of the velocity field with itself. As such, we also predict that the chaos due to the velocity field is dominant over that for the magnetic field.

In extending findings of Eulerian chaos in hydrodynamics to MHD, there are other complications that must be tested, such as whether has any dependence on Re or Pr = , the magnetic Prandtl number, which has important effects on dissipation rates the presence of inverse spectral transfer McKay et al. (2018). This paper relies on this previous foundational work.

Care must be taken in the consideration of the two fields and , as either separate or to some union of the two fields. As well, in thinking about the saturation of error growth found in hydrodynamics, there may be a dependence on individual dissipation rates.

Results from DNS simulations suggest that the ratio of and depends on the Prandtl number according to the relationship Pr Brandenburg (2014); McKay et al. (2018) where depends on the presence of helicity in the system with and so should become totally dominant over . The limit for the growth of the difference between the two fields in hydrodynamics Berera and Ho (2018) may carry over to MHD, and the specific dependence on either or needs to be tested. These are examined in Section V.

## Iv Direct numerical simulation

We performed direct numerical simulations of forced HIT on the incompressible MHD equations using a fully de-aliased pseudo-spectral code in a periodic cube of length with unit density. An external forcing function was applied to maintain energy in the system. The code and forcing are fully described in Yoffe (2012); Linkmann (2016); McKay et al. (2017). The field is initialised with a set of random variables following a Gaussian distribution with zero mean. The initial kinetic and magnetic energy spectra were where is the peak wave number. The primary forcing used was an adjustable helicity forcing, forced for wave numbers . This allowed maintenance of a desired magnetic helicity, which was zero except where noted, and cross helicity, which was always zero. The necessary benchmarking for the simulations was done in McKay et al. (2018). This included making sure that all simulations were fully resolved in both fields. This was key in allowing greater confidence in the accuracy of the final results presented in this paper.

After a statistically steady state was reached, a copy of evolved fields and were made. At one time point, these fields were perturbed by adding a white noise to them of size , creating fields and . This introduced a difference between the fields, with both sets of fields then evolved independently.

The Lyapunov exponents was obtained, by tracking the two independently evolved fields. The difference spectra and were measured as

 Eud(k,t) =12∫|k|=kdk|^u1(k,t)−^u2(k,t)|2 , (28) Ebd(k,t) =12∫|k|=kdk|^b1(k,t)−^b2(k,t)|2 , (29)

, and . We model the growth of and by an exponential with and . In our simulations we found that the exponential growth of both fields was the same and occurred with approximately the same rate, and so for all runs. For reference this is termed the direct method.

As a way to cross-check our main method, we also use another common approach, finite time Lyapunov exponents (FTLEs) Ott (2002). This method does the following: At a time after the initial perturbation, and at every subsequent time interval of the field is rescaled according to the rule

 u2→u1+ (δ0δu)(u2−u1) , (30) b2→b1+ (δ0δb)(b2−b1) . (31)

Here . For each time interval an FTLE, , was defined

 γ=1δtln(δδ0) . (32)

The Lyapunov exponent is then the average of these .

If the fields were treated separately, with and , the FTLEs measured were different for both fields. The FTLE for the kinetic field was similar to that found using the direct method, whilst that for the magnetic field averaged to zero. Using the direct method, which is the more physically natural description, the magnetic fields diverged exponentially and so we kept as coming from the union of the two fields. The fact that splitting the fields showed that the field had a FTLE of almost zero is in agreement with our theoretical argument in Section III that the magnetic field is only weakly chaotic and that the most important contribution to chaos comes from the velocity field. Indeed, this finding provides extra support for that conclusion.

Using the FTLEs allowed a check on the direct method. The time over which the FTLEs were averaged was not very long due to simulation restrictions, however, they did agree broadly with the other method and allowed us to be more sure of the validity of the results. None of the results were qualitatively changed if the FTLE method was used as opposed to the direct method, and the quantitative changes were small. The averaging procedure of the FTLE tends to produce more stable results, but at the same time this averaging may miss features such as the long term linear behavior. However, it is the most effective to use to test against the direct method we are using.

## V Results

### v.1 Re and Pr Dependence

We use the adjustable helical forcing (with the magnetic helicity set to zero) for a set of parameters of and , which vary independently from 0.01 to 0.0003125. After a statistically steady state was reached, the Lyapunov exponent was measured using the direct method. Although the data is noisy, this noise appeared in both methods used to measure the Lyapunov exponent, as such we expect that this noise is an inherent issue with MHD turbulence as opposed to a measurement issue. The noise may arise from the more complicated dynamics of MHD. Specifically, although total magnetic helicity sums to zero, there are likely regions of opposite magnetic helicity which can reduce the chaos and cannot totally be eliminated.

A grid of viscosities and diffusivities is shown in Table 1. The corresponding grid point is the maximal Lyapunov exponent. What is seen very easily in the table is that those simulations with very low Pr (in the bottom left corner of the table) have much higher than the high Pr simulations, which is also seen using the FTLE method. Since low Pr means that the kinetic turbulence is dominant, this suggests that more kinetic dominance results in more chaotic dynamics.

In this data, Re varies from 50 to 2000, and Pr varies from 1/32 to 32. The values of forcing resulted in a total dissipation of . A full table of simulation parameters is shown in Table 2. The wave number dependence of the field, and , were similar to those found in pure hydrodynamics and similar to each other Berera and Ho (2018).

Using this data is plotted against the prediction of as is shown in Fig. 1. The inset shows against Re, and has no clear dependent trend. Meanwhile, the main data shows broad agreement that , although there is an offset which may be related to the onset of chaos in the system only after a certain level of turbulence is reached. In the figure, the solid (green) line shows the best fit to the data, whilst the dashed black line shows the relationship between and for hydrodynamics found in Berera and Ho (2018). By eye, the relationship from pure hydrodynamic turbulence is not noticeably worse than that of the fit to actual MHD simulations. Specifically, the proportionality constant from the data is (whilst for the FTLE method it was ), which is not far from our theoretical prediction.

The data for the larger in MHD is noisier (in general at higher Re) than for hydrodynamics, whilst for low Re the data is comparatively clear. However, it seems that the chaos is entirely dependent on the chaos due to the velocity field as also shown by the results from removing the Lorentz force later in the paper. The higher results are especially noisy, whilst up to this scale, the data has less noise. This is probably due to the fact that these higher simulations are on the edge of our resolution range, but the fit does not weigh these heavily.

The data in Table 1 can be split by Prandtl number, and then calculate the dependence on Reynolds number separately. Each was approximated as following Eq. (4), giving a different Pr. This is equivalent to saying that in Eq. (4), which is 0.53 in Navier-Stokes turbulence (or 0.64 in other simulations Boffetta and Musacchio (2017); Mohan et al. (2017)), depends on the Prandtl number as

 λ=1T0Ref(Pr) , (33)

where is some function of the Prandtl number, attaining roughly 1/2 at Pr = 1. A plot of Pr is shown in Fig. 2. The solid black line shows a power law behaviour Pr, whilst the dashed (green) line shows a constant . For Pr = 1, we find (the FTLE method produced an equivalent result within one standard deviation).

Given that there is no theory to suggest either option, and that has the greatest predictive quality, the data is consistent with the idea that any dependence of on Pr comes indirectly through its dependence on . The data for Pr is consistent with there being no fundamental dependence at all. If, instead, it were a function of Rm, we note that Rm = PrRe, and so the dependence would be the same as if Pr were kept constant whilst calculating the exponent, . That the results become independent of Pr is similar to results found in McKay et al. (2018) of dissipation rates.

Indeed, if we want to be guided by the principle of universality, that at sufficiently high Re or Pr we should have only a function of Pr or Re respectively, then such a power law dependence for would be precluded. For example, this universality was observed in Brandenburg (2014) and put on a firm footing in fully resolved simulations of MHD in McKay et al. (2018). These simulations show that there is complete dominance of the dissipation by the kinetic dissipation, ie at sufficient Prandtl number and Rm. At increasing Rm, the ratio of , and depending on the degree of accuracy required, the Lyapunov exponent can be approximated as solely dependent on and , which would be consistent with an extension to MHD of the first Kolmogorov hypothesis of similarity Kolmogorov (1941).

Given that there are also many other time scales that could be used, such as the Kolmogorov microscale time for the magnetic field or the large eddy turnover time for the magnetic field, we did check the relationship between these and and found none that were better related than . In any case, none of the others had theoretical justification. As such, we suggest that our extension of the Ruelle theory to MHD has the greatest predictive power consistent with our data. This may be due to the nature of the MHD equations themselves. We note that there is a direct non-linearity in the evolution of , whilst is only non-linear indirectly through the evolution of . There is also an indirect non-linearity of through the evolution of . The direct non-linearity seems to be the most important for the chaos.

Since the exponential growth of and had the same exponent, we surmise that if one field becomes sufficiently different, it drags the other one along with it. The data suggests that the difference in drives the difference in . This is backed by the fact that the relationship between and for hydrodynamics fits the data. This would be the case if the difference in the magnetic field merely reacts to that in the velocity field, which is the true driver of chaos in the system.

Given the relationship to and not any MHD quantities, we suggest that this direct non-linearity is the most important feature and driver of the MHD chaos. However, in some MHD dominated by magnetic reconnection it cannot be precluded that the most important time scale would be the reconnection rate, which becomes independent of , the Lundquist number, for sufficiently high values of Loureiro et al. (2012).

### v.2 Helicity

Given past findings that suggest that helicity results in a diminution of chaos Zienicke et al. (1998); Escande et al. (2000), and the fact that we can directly control and measure the relative helicity of our system, we test this directly in the Eulerian sense. The relative helicity is defined

 ρ2b=⟨b⋅a⟩⟨b⋅a⟩⟨b⋅b⟩⟨a⋅a⟩ , (34)

with vector potential . We run a set of simulations where we vary the quantity by injecting helicity into the system with the adjustible helicity forcing described in McKay et al. (2017). All simulations were done with hypoviscosity, which took energy out at the large scale (low wave number) to ensure that the simulations did not blow up due to the inverse transfer present in helical simulations. After a statistically steady state was reached, including having roughly constant , the FTLE method was used to measure the Lyapunov exponent, except for the inset of Fig. 4 which was done using the direct method. The results of this are shown in Fig. 4; note in the main plot the Reynolds number was 120. What is seen is that increasing does indeed create a diminution of chaos. The data suggest that as , . Thus, we have assumed a parameter dependence of , where is the Lyapunov exponent at and is some fit parameter, here it is equal to . This is shown in grey (blue) in the plot. However, the uncertainty on means that many functional dependencies on are consistent, but these would need a general tendency for at maximal magnetic helicity.

The inset of Fig. 4 shows a run with Re 320. In the run without helicity, there is a regular exponential increase in whilst for the run with , after the initial perturbation, the two realisations remain close for many large eddy turnovers, . This effect is rather dramatic. Helicity has some organising effect on the fluid Brandenburg and Subramanian (2005). As helicity increases, the organising effect becomes greater and the chaos diminishes. We believe that the diminution of chaos cannot only be a result of the inverse transfer present in helical MHD. For instance, in two-dimensional hydrodynamics there is an inverse transfer of energy, but still a robust amount of chaos is seen Boffetta and Musacchio (2001).

Given our assumption that the Lyapunov exponent for Eulerian and Lagrangian chaos are equal, and that the exponent for Lagrangian chaos is intimately related to the longitudinal correlation function, we would expect that depends on the magnetic helicity in the system. However, would require from Eq. (18), such that the whole field is correlated. This would violate the assumption of homogeneity and isotropy on which Eq. (18) relies. As such, there will be a breakdown in this relationship, and the organisational property of the magnetic helicity suppresses the chaos coming from the velocity field. As well, if and are equal, we also expect that helicity diminishes chaos even in a Lagrangian sense, which was seen in ABC flows of MHD Zienicke et al. (1998).

It has been seen before that helicity plays an important role in MHD turbulent dynamics. Specifically, even without initial helicity, the system is attracted to states that are helical and laminar even with a large Re Dallas and Alexakis (2015). These Beltrami states have greatly reduced non-linearity Servidio et al. (2008). This attractor behaviour is also seen in pure hydrodynamics Linkmann and Morozov (2015), where these laminar states also exist. The laminar states in hydrodynamics are purely attractive, such that once the system becomes laminar, it does not delaminarize. However, in MHD, the helical states can be exited spontaneously, but with a tendency to stay near these states longer as helicity is increased. This may explain the difference in behaviour of perturbations made in the laminar states in hydrodynamics, where decays exponentially, and those made here, where remains roughly constant. In this sense, the MHD simulations act like the hydrodynamic simulations in Berera and Ho (2018) which have very low Re, but which have still not relaminarized.

### v.3 Linear Growth

In hydrodynamic turbulence, there is a limit on the growth rate of , which eventually seems to grows linearly in time with rate equal to the dissipation rate as has been found recently Berera and Ho (2018). We find an analogous effect in MHD with important consequences for the predictability of MHD systems with high Pr, such as galactic scale magnetic fields. The total dissipation is defined, , as the sum of the dissipations due to and field respectively. A negative damping forcing could be used to maintain a constant dissipation rate whilst forcing the velocity to maintain a roughly zero relative helicity McKay et al. (2017).

Eventually, both and do enter a linear growth phase. However, the linear rate for both fields is different. Specifically, for the rate is determined by . In Fig. 4, the growth of is shown with time for three different runs, one with Re = 120 and = 0.11, one with Re = 80 and = 0.036, and one with Re = 45 and = 0.0012. In the plot is divided by and after this normalisation, the linear growth rate for each is roughly the same. Eventually . The values for used here differ by a factor of 100 so we feel relatively confident that this linear relation is based upon the magnetic dissipation, which resolves the question raised in Section III.

For the growth of it is difficult to distinguish whether the behaviour is based on or because the two values tended to be close. We suspect it is the former, and this would be more appealing symmetrically, but the difference is less pronounced in the values of the data. This is because it is difficult in practice to get .

The fact that the exponential growth of both fields is the same, whilst the linear growth at late times is different suggest that they are caused by fundamentally different processes. Or at least they come form different aspects of dynamical systems theory, and that these processes are controlled by different variables.

The differing linear growth rates is a new and very interesting result in our opinion, with important implications for predictability. Specifically, if , as happens for high Pr Brandenburg (2014); McKay et al. (2018), whilst the magnetic field still has a significant amount of energy in it, the predictability time would diverge. It also helps to justify the kinetic approximation (where there is no feedback of the magnetic field to the velocity field) in some instances. This is because the kinetic field, diverging rapidly, must have many different fields which correspond to the same evolution of the magnetic field.

### v.4 FSLEs

As an alternative to the direct method, the level of chaos can be quantified using finite size Lyapunov exponents (FSLEs). These are defined by the time that it takes for an error of size to grow by a factor . A set of five simulations was performed for three different Reynolds numbers, each with Pr = 1. After an initial small perturbation, the time taken for the perturbation to grow by a factor was measured and an average of these times was taken across the set of five simulations. Using this, an FSLE Lyapunov exponent can be defined . At very low . The velocity and magnetic fields were treated separately and the results for each is shown in Fig. 5 and Fig. 6 respectively. The values are normalised by the rms quantity of the relevant field and a relevant timescale , where is the energy of that field and the dissipation due to that field.

Lorenz has suggested that for fluid turbulence a disturbance in the inertial range will grow with rate equal to the local eddy turnover time Lorenz (1969). This eventually predicts that Boffetta and Musacchio (2017). However, because the interplay of different fields and the potential presence of non-local interactions between the and fields, there is no great justification for this in MHD. As such, the data is presented here to facilitate comparison with hydrodynamic results. The inset of the figures shows a linear fit to data which is closer over a larger range of values than the fit.

The data is nosier than for hydrodynamics, this is again like all the previous data we have found. In both linear and logarithmic plots, the FSLE collapse onto a slope which is roughly independent of Re. This collapse is more consistent for the velocity field. The parameters used did not ensure a constant ratio , which varies from 2 for Re = 80 to 1/2 for Re = 670. As such, we are fairly confident that the relevant dissipation for determining the dynamics is that of the respective field and not the total dissipation. In this way, the fields, despite being non-linearly connected, have linear behaviour dictated by their own evolution.

In the region where is a linear function of , where and are positive constants. This last equality comes from a small expansion of about 1 in the definition of . This differential equation is solved by

 δ=cm+(cδ0−m)e−ct , (35)

where is the separation at . This gives a sigmoid shape, with maximum rate of growth when is at half of its maximum. The maximum of is . At . Remembering the normalisation of and we redefine and , where for the velocity field and for the magnetic field from Fig. 5 and Fig. 6. In both cases, . The rms values are defined and . Thus, the maximum rate of growth for the relevant difference has

 ∂tEid(t)≤9a48εi . (36)

For , this predicts a maximum of , where , which is very close to the 1.12 previous found in hydrodynamic turbulence Berera and Ho (2018). This also implies that the magnetic field has a maximum growth of difference of .

A study of dynamics of systems with different timescales, of which MHD turbulence is definitely an example, looked at the study of both scales through the use of FSLE Boffetta et al. (1998). It concludes with a note that parametrization of the fast scales is not crucial and that the slow mode dynamics are dominant. In MHD turbulence, the dynamics of the velocity field are equivalent to the fast modes and the magnetic field to the slow modes. As such, we also predict that the velocity field can be successfully parametrized whilst still capturing the most important magnetic field dynamics, as is done in the static field approximation.

### v.5 Influence of the Lorentz force

The MHD equations Eq. (1-2) can be changed to omit the action of the Lorentz force on the velocity field. Although this means that the fields no longer conserve energy, we use this as a diagnostic tool to disentangle the amount of chaos that comes from perturbations to each field. This has been performed in prior simulations of MHD turbulence investigating the inverse transfer of energy, finding that it is not necessary for inverse transfer to persist Berera and Linkmann (2014). In our simulations, if the energy spectra when decreasing with wave number were approximated as following , the for both fields was somewhat smaller than for an equivalent simulation with a Lorentz force.

We perform a simulation using the helical forcing as above for Pr = 1 with and Re = 74, where we have removed the Lorentz force . This makes the magnetic field equation linear. Even with the Lorentz force, the non-linearity for in the MHD equations comes entirely from . In this case, the magnetic field alone cannot produce chaos or turbulence, but this gets input from .

Either the magnetic or kinetic field can be perturbed individually. Both perturbations were made from identical instances of the fields. A magnetic field perturbation would never induce a kinetic perturbation but a kinetic perturbation can induce a magnetic perturbation.

Whilst perturbing the kinetic field causes an exponential divergence, a magnetic field perturbation can remain relatively stable for many large eddy turnover times. The ad hoc removal of the Lorentz force is sometimes done in seed dynamo simulations, where it is argued that the magnetic field is too small to affect the velocity field. Because of this tendency, the Reynolds number used was a modest 74, otherwise the magnetic field has a tendency to accumulate energy indefinitely. However, in simulations which did have an unphysical exponential growth of magnetic energy Berera and Linkmann (2014), the perturbation in the magnetic field grew no faster than the magnetic field itself did. As such, we find the presence or absence of dynamo action should be irrelevant for chaos.

For the velocity perturbation, both magnetic and velocity fields diverged at the same rate, which agrees with our interpretation that it is the velocity field which drives the chaos. The Lyapunov exponent found was roughly in agreement with our theoretical prediction from Eq. (19), being dependent on the kinetic dissipation rate.

Although the results in Fig. 1 are noisier than for hydrodynamic turbulence, this diagnostic tool of removing the Lorentz force is another result which supports our theoretical argument in Section III. Namely, it also shows that the chaotic properties are dominated by the velocity field properties, here being the dependence on as opposed to any other timescales. This diagnostic tool also shows that the magnetic field itself is not chaotic. This is important for simulations where the Lorentz force is removed by affecting the statistics, such as in seed dynamo simulations.

Looking at the prediction from Eq. (27), if there is no Lorentz term, then exactly. In the absence of a Lorentz force, the velocity field behaves the same as in hydrodynamic turbulence. Since, by removing the Lorentz term, our results are not significantly changed, we can state that the chaos in the system with a magnetic field is not significantly different to that without it. As such, our prior assumptions that the magnetic field should not affect the chaos greatly and that the longitudinal correlation functions are also not greatly affected are in agreement with our results, and are strengthened by them.

## Vi Discussion and Conclusion

The data presented here include Pr which cover three orders of magnitude. Although many of the results here are independent of magnetic Prandtl number, there is direct comparison with many physical systems. Here we look at some of them in turn and see how our measurements can be useful for their further study.

In toroidal plasmas, such as those found in tokamak reactors, Pr is expected to be on the order of 100 Itoh et al. (1993); Mendonca et al. (2018). However, the turbulence is strongly confined and should be greatly affected by the boundary conditions, as such, an assumption of homogeneity and isotropy will not capture the full dynamics. There may be other practical reasons not to apply MHD as a model for the plasmas in tokamak reactors. Even so, the level of chaos in the system should be related to the generation of instabilities in the flow. These instabilities can affect the confinement of the plasma. As well, the level of instability should be affected by the level of inverse cascade, which is greater at increased Rm.

Accretion disks around black holes and neutron stars are predicted to have regions where Pr is close to unity Balbus and Henri (2008). Though these should be affected by relativistic effects, if the results in this paper extend beyond the Prandtl number range of our simulations they may also apply to more typical accretion disks. By understanding the ratio of the dissipations, we can understand the relative importance of ion and electron heating in these systems which has an effect on the luminosity and thus their observational characteristics. In regions with greater chaos of the particles, those regions with lower Kolmogorov microscale time and lower magnetic helicity, the process of accretion should be reduced. If particles are, on average, moving together and not drifting apart exponentially, there should be a longer time for them to accrete. As such, in any accretion disk, we expect that the accretion rate should be increased by magnetic helicity and . Indeed, simulations have shown already that accretion disks are associated with turbulent dynamo action, which is also associated with magnetic helicity Brandenburg and Subramanian (2005). Thus, an increased accretion rate should be associated with diminution of chaos, here consistent with increased magnetic helicity.

Liquid sodium experiments are the natural venue for testing the prediction that in a Lagrangian sense. As well, if and are related, then an increase in magnetic helicity should be reflected by a lessening of chaos in Lagrangian tracers.

Understanding the chaos seen in MHD turbulence is important for understanding the scope for prediction of turbulent conducting fluids. The findings here may also be applicable to coupled dynamical system where there is a direct non-linearity for only one of the fields and multiple relevant timescales.

Our extension of the Ruelle prediction to MHD that is most consistent with our own data. We conclude that this is because the MHD equations are directly non-linear only in and indirectly through the coupling for . We have also confirmed previous findings that show magnetic helicity results in a diminution of chaos and further predict that a fully helical system should have zero Lyapunov exponent, although this was not done for DNS Zienicke et al. (1998); Escande et al. (2000).

One of the most interesting results is the new finding that the growth of both and becomes linear with rates that depend on the dissipation rate of the relevant field. This may apply more generally to non-linearly coupled fields. Specifically in the case of high Pr, where becomes dominant and very small, then this has important implications for the long term predictability of galactic plasmas and magnetic fields. For high Re, we should expect to be very large and that any small scale error should grow in size very quickly and so the the separation between the two fields will enter the linear regime very quickly. Thus, any predictability time will be dominated by for the relevant field. For magnetic fields, this can become extremely large, and if there is a large amount of helicity in the system, then the predictability time of the magnetic fields can become very long.

###### Acknowledgements.
We thank Mairi E. McKay for helpful discussion. This work has used resources from ARCHER ARCHER () via the Director’s Time budget. This work used the Cirrus UK National Tier-2 HPC Service at EPCC Cirrus () funded by the University of Edinburgh and EPSRC (EP/P020267/1). R.D.J.G.H is supported by the U.K. Engineering and Physical Sciences Research Council (EP/M506515/1), D.C. is supported by the University of Edinburgh. A.B. acknowledges funding from the U.K. Science and Technology Facilities Council.

## References

• Berera and Ho (2018) A. Berera and R. D. Ho, Phys. Rev. Lett. 120, 024101 (2018).
• Boffetta and Musacchio (2017) G. Boffetta and S. Musacchio, Phys. Rev. Lett. 119, 054102 (2017).
• Ruelle (1979) D. Ruelle, Phys. Lett. 72A, 81 (1979).
• Crisanti et al. (1993) A. Crisanti, M. H. Jensen, G. Paladin,  and A. Vulpiani, J. Phys. A: Math. Gen. 26, 6943 (1993).
• Biskamp (2003) D. Biskamp, Magnetohydrodynamic turbulence (Cambridge Univ. Press, 2003).
• Brandenburg and Subramanian (2005) A. Brandenburg and K. Subramanian, Phys. Reports 417, 1 (2005).
• Misguich and Balescu (1982) J. H. Misguich and R. Balescu, Plasma Phys. 24, 289 (1982).
• Misguich et al. (1987) J. H. Misguich, R. B. H. L., Pecsell, T. Mikkelsen, S. E. Larsen,  and Q. Xiaoming, Plasma Phys. Controlled Fusion 29, 825 (1987).
• Wolf et al. (1985) A. Wolf, J. B. Swift, H. L. Swinney,  and J. A. Vastano, Physica D 16, 285 (1985).
• Huang et al. (1994) W. Huang, W. X. Ding, D. L. Feng,  and C. X. Yu, Phys. Rev. E 50, 1062 (1994).
• Kurths and Herzel (1986) J. Kurths and H. Herzel, Solar Phys. 107, 39 (1986).
• Grappin et al. (1986) R. Grappin, J. Leorat,  and A. Pouquet, J. de Physique 47, 1127 (1986).
• Yamada and Ohkitani (1987) M. Yamada and K. Ohkitani, J. Phys. Soc. Jpn. 56, 4210 (1987).
• Zienicke et al. (1998) E. Zienicke, H. Politano,  and A. Pouquet, Phys. Rev. Lett. 81, 4640 (1998).
• Rempel et al. (2012) E. L. Rempel, A. C. L. Chian,  and A. Brandenburg, Physica Scripta 86, 018405 (2012).
• Rempel et al. (2013) E. L. Rempel, A. C. L. Chian, A. Brandenburg, P. R. Muñoz,  and S. C. Shadden, J. Fluid Mech. 729, 309 (2013).
• Escande et al. (2000) D. F. Escande, P. Martin, S. Ortolani, A. Buffa, P. Franz, L. Marrelli, E. Martines, G. Spizzo, S. Cappello, A. Murari, et al., Phys. Rev. Lett. 85, 1662 (2000).
• Mirmomeni and Lucas (2009) M. Mirmomeni and C. Lucas, Space Weather 7, S07002 (2009).
• Karak and Nandy (2012) B. B. Karak and D. Nandy, ApJ Lett. 761, L13 (2012).
• Weigel et al. (2003) R. S. Weigel, A. J. Klimas,  and D. Vassiliadis, J. Geophys. Res. 108, 1298 (2003).
• Amari et al. (2014) T. Amari, A. Canou,  and J. Aly, Nature 514, 465 (2014).
• Kolmogorov (1941) A. N. Kolmogorov, Dokl. Akad. Nauk SSSR,  30, 301 (1941).
• Mohan et al. (2017) P. Mohan, N. Fitzsimmons,  and R. D. Moser, Phys. Rev. Fluids 2, 114606 (2017).
• Biferale et al. (2005) L. Biferale, G. Boffeta, A. Celani, B. J. Devenish, A. Lanotte,  and F. Toschi, Phys. Fluids 17, 115101 (2005).
• Batchelor (1953) G. K. Batchelor, The Theory of Homogeneous Turbulence (Cambridge Univ. Press, 1953).
• de Divitiis (2018) N. de Divitiis, arXiv preprint arXiv:1802.08081  (2018).
• Salazar and Collins (2009) J. P. L. C. Salazar and L. R. Collins, Annu. Rev. Fluid Mech. 41, 405 (2009).
• Yeung and Pope (1989) P. K. Yeung and S. B. Pope, J. Fluid Mech. 207, 531 (1989).
• Guala et al. (2007) M. Guala, A. Liberzon, A. Tsinober,  and W. Kinzelbach, J. Fluid Mech. 574, 405 (2007).
• Nastac et al. (2017) G. Nastac, J. W. Labahn, L. Magri,  and M. Ihme, Phys. Rev. Fluids 2, 094606 (2017).
• McKay et al. (2018) M. McKay, A. Berera,  and R. Ho, “Fully-resolved array of simulations investigating the influence of the magnetic prandtl number,” arXiv preprint arXiv:1806.10202 (2018).
• Brandenburg (2014) A. Brandenburg, ApJ 791, 12 (2014).
• Yoffe (2012) S. R. Yoffe, Ph.D. thesis, University of Edinburgh (2012), arXiv:1306.3408.
• Linkmann (2016) M. F. Linkmann, Ph.D. thesis, University of Edinburgh (2016), http://hdl.handle.net/1842/19572.
• McKay et al. (2017) M. E. McKay, M. Linkmann, D. Clark, A. A. Chalupa,  and A. Berera, Phys. Rev. Fluids 2, 114604 (2017).
• Ott (2002) E. Ott, Chaos in Dynamical Systems (Cambridge Univ. Press, 2002).
• (37) The data is publically available, see http://dx.doi.org/10.7488/ds/2410.
• Loureiro et al. (2012) N. F. Loureiro, R. Samtaney, A. A. Schekochihin,  and D. A. Uzdensky, Phys. Plasmas 19, 042303 (2012).
• Boffetta and Musacchio (2001) G. Boffetta and S. Musacchio, Phys. Fluids 13, 1060 (2001).
• Dallas and Alexakis (2015) V. Dallas and A. Alexakis, Phys. Fluids 27, 045105 (2015).
• Servidio et al. (2008) S. Servidio, W. H. Matthaeus,  and P. Dmitruk, Phys. Rev. Lett. 100, 095005 (2008).
• Linkmann and Morozov (2015) M. F. Linkmann and A. Morozov, Phys. Rev. Lett. 115, 134502 (2015).
• Lorenz (1969) E. N. Lorenz, Tellus 21, 289 (1969).
• Boffetta et al. (1998) G. Boffetta, A. Crisanti, F. Paparella, A. Provenzale,  and A. Vulpiani, Physica D 116, 301 (1998).
• Berera and Linkmann (2014) A. Berera and M. Linkmann, Phys. Rev. E 90, 041003 (2014).
• Itoh et al. (1993) K. Itoh, S. I. Itoh, A. Fukuyama, M. Yagi,  and M. Azumi, J. Phys. Soc. Jpn. 62, 4269 (1993).
• Mendonca et al. (2018) J. Mendonca, D. Chandra, A. Sen,  and A. Thyagaraja, Phys. Plasmas 25, 022504 (2018).
• Balbus and Henri (2008) S. A. Balbus and P. Henri, ApJ 674, 408 (2008).
• (49) ARCHER, http://www.archer.ac.uk.
• (50) Cirrus, http://www.cirrus.ac.uk.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters