Modeling the Solar Wind Turbulent Cascade Including Cross Helicity: With and Without Expansion [original]

Modeling the Solar Wind Turbulent Cascade Including Cross Helicity: With and

Without Expansion

Roland Grappin

1,2

, Andrea Verdini

3,4

, and W.-C. Müller

5,6

Laboratoire de Physique des Plasmas (LPP), École Polytechnique, Rte de Saclay, 91120 Palaiseau, France; [email protected]

CNRS, Observatoire de Paris, Sorbonne Université, Université Paris Saclay, Ecole Polytechnique, Institut Polytechnique de Paris, F-91120 Palaiseau, France

Università di Firenze, Dipartimento di Fisica e Astronomia, Firenze, Italy

INAF, OAA, Firenze, Italy

Technische Universität Berlin, ER 3-2, Hardenbergstr. 36a, D-10623 Berlin, Germany

Max-Planck/Princeton Center for Plasma Physics, Hofgartenstr., 8 80539 Munich, Germany/Princeton University P.O. Box 430 Princeton, NJ 08544-1019, USA

Received 2022 February 7; revised 2022 April 7; accepted 2022 April 27; published 2022 July 19

Abstract

Simulations of the turbulent cascade forming in the solar wind, including cross helicity, commonly adopt a

homogeneous setup, not taking into account wind expansion. Here we want to assess the predictions of decaying

3D compressible (low Mach number)MHD simulations, respectively homogeneous and with expansion, in order

to examine which is the most fruitful approach to understanding the turbulent cascade in the solar wind. We follow

turbulent evolution during 10 nonlinear turnover times, considering several initial values of the initial spectral

slope and cross helicity. In the expanding case, the transverse sizes of the plasma volume are stretched by a factor

of 5 during the simulation, corresponding to traveling from 0.2 up to 1 au. In homogeneous simulations, the

relative cross helicity rises, and the Elsässer spectra E

show “pinning,”with a steep dominant spectrum and ﬂat

subdominant spectrum, the ﬁnal spectral indices depending on cross helicity but not initial indices. With

expansion, the relative cross helicity decreases, and dominant and subdominant spectra share the same index, with

the index relaxing to an asymptotic value that generally depends on the initial index. The absence of pinning, as

well as the decrease of relative cross helicity, probably both rely on the permanent injection by expansion of an

excess of magnetic energy at the largest scales, equivalent to injecting subdominant energy. Also, spectra generally

steepen when initially starting ﬂatter than k

−5/3

but stop evolving at a ﬁnite time/distance.

Uniﬁed Astronomy Thesaurus concepts: Interplanetary turbulence (830);Solar wind (1534);Space plasmas (1544);

Alfven waves (23);Interplanetary medium (825);Magnetohydrodynamics (1964)

1. Introduction

Turbulent dissipation and heating play an important role in

the evolution of the interplanetary plasma and the acceleration

of the solar wind, as demonstrated in solar wind models

(Cranmer et al. 2007; Verdini et al. 2010; Chandran et al. 2011;

van der Holst et al. 2014; Shoda et al. 2018,2019; Réville et al.

2020; Matsumoto 2021). At the base of this, there is the

question of the nature of the turbulent cascade at work in the

solar wind.

To proceed, one can compare direct predictions of cascade

theories to observations, possibly with the help of direct

simulations. A good example is that of local spectral anisotropy,

computed in a frame attached to the local mean magnetic ﬁeld and

using structure functions. Assuming axial symmetry around the

direction of the local mean ﬁeld, work by Horbury et al. (2008)

ﬁrst recovered the so-called critical balance results, in which the

perpendicular cascade is strong and the spectrum along the mean

ﬁeld is generated by linear transport along it. Later, work by Chen

et al. (2012)showed important deviations from axisymmetry

around the local mean ﬁeld. This was reproduced by 3D MHD

simulations showing that expansion breaks the local axisymmetry

in the whole inertial range (Verdini & Grappin 2015).

A second example is that of the origin of the global spectral

anisotropy, that is, the origin of the so-called Maltese cross

(Matthaeus et al. 1990)obtained in a frame attached to the

mean magnetic ﬁeld computed during a ﬁxed interval of time.

We identiﬁed the two branches of the Maltese cross as

corresponding to the signature of, respectively, weak or strong

expansion. We could reproduce the two conﬁgurations via

simulations with expansion and in both cases recover the

observed slow plasma cooling with distance (Montagud-Camps

et al. 2018,2020). A last interesting example is that of the

origin of the switchback formation in the early solar wind,

which seems to be possibly attributed to the expansion of the

plasma (Squire et al. 2020).

Despite these successes, the large majority of work in which

a cascade theory is proposed for the solar wind neglects

expansion (e.g., Boldyrev 2005; Lithwick et al. 2007; Perez &

Boldyrev 2010; Chandran et al. 2015). These works predict

speciﬁc spectral scalings for z

(see deﬁnition below)or

velocity and magnetic ﬂuctuations covering the observed range

[−5/3, −3/2]. However, in situ observations show that the

spectral indices, at scales above the sub-ion scale, are

distributed in a large interval of values at a given distance

(Grappin et al. 1991; Boldyrev et al. 2011), and that the

spectral index varies with distance (Bavassano et al. 1982;

Chen et al. 2020)and the average properties of the ﬂow, such

as proton temperature (Grappin et al. 1991), solar wind speed,

and cross helicity (Chen et al. 2013). In this paper, we

investigate such properties via numerical simulations, compar-

ing the results obtained with and without expansion and

insisting on Alfvénic (i.e., high cross helicity)streams, which

are ubiquitous in the solar wind.

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 https://doi.org/10.3847/1538-4357/ac6ba4

Original content from this work may be used under the terms

of the Creative Commons Attribution 4.0 licence. Any further

distribution of this work must maintain attribution to the author(s)and the title

of the work, journal citation and DOI.

Note that a whole family of works dealing with the evolution

of turbulence in the solar wind (e.g., Adhikari et al. 2020)does

not describe the details of the cascade, such as the evolution of

turbulent energy spectra, but instead directly proposes

approximate evolution equations for moments, i.e., energies

and correlations, including some known properties associated

with expansion.

In the present work, we instead follow in detail the formation

and decay of the turbulent cascade by solving the primitive

MHD equations in two versions: the standard compressive

MHD equations (without expansion)on one side and MHD

equations with ﬁnite expansion on the other side, also called

EBM equations (see below). The plasma evolution is followed

during 10 nonlinear times, which in the expanding case

corresponds to the heliocentric distance of the plasma volume

increasing by a factor of 5, e.g., varying from 0.2 to 1 au. We

consider in this work cases dominated by outward-propagating

Alfvén waves, as well as cases with equal populations

propagating in both directions, as both cases are important in

the solar wind.

We run weakly compressive simulations, i.e., with a turbulent

Mach number equal to 0.12, characteristic of the quasi-

incompressible ﬂuctuations at hour scales (Matthaeus et al.

1991; Bavassano & Bruno 1995;Ofman2010). We adopt in our

simulations with expansion an expansion rate ò=0.4 (see

deﬁnition below; Equation (7)), a value for which expansion

terms are comparable to nonlinear terms at the large scales. Such

an expansion rate is characteristic of scales of the order of 1 day

(Montagud-Camps et al. 2020); a more realistic study would

require using a smaller expansion rate and/or a very high

resolution in order to include scales smaller than 1 hr, but this

would require important numerical resources, and we think that

the present choice of parameters is a reasonable compromise for a

ﬁrst study. These choices allow one to study a quasi-

incompressible turbulent state and clearly reveal the differences

between homogeneous turbulence and turbulence with expansion.

The plan of the present paper is as follows. Section 2

describes the equations, basic parameters, and initial conditions

of the simulations. Section 3gives the results. Section 4is a

discussion, and Section 5is a conclusion.

2. Method: Equations

2.1. Equations and Deﬁnitions

Here we brieﬂy derive the basic equations and method (for a

more detailed derivation, see, e.g., Grappin et al. 1993; Dong

et al. 2014; Montagud-Camps et al. 2018).

We start with the MHD equations for the density ρ,(isotropic)

pressure P, velocity ﬂuctuation ˆ

uUUe

=- (where Uis the

total velocity and U

is the mean radial ﬂow amplitude),and

magnetic ﬁeld B. Consider a Cartesian frame with XYZ

coordinates, the X-axis parallel to the radial passing through

the middle of the box, and change to a Galilean frame moving

with the mean wind along the radial coordinate. In this frame,

the plasma volume is uniformly stretched in the transverse

directions (see Figure 1(b)), thus neglecting curvature terms.

All ﬁelds are assumed periodic in the comobile coordinates

x,y, and z:

()t,1t=

() ()xXU a,2

x0t=-

() ( )yYat,3=

() ( )zZat.4=

The parameter

()aLL LL 5

yxz

is the initial aspect ratio of the domain (L

, and L

being the

size of the domain in the three directions and the sufﬁx“0”

denoting the initial value). Note that the radial size L

is a

constant, while the other sizes of the domain increase linearly

with time. The parameter a(t)measures the transverse

expansion. It is deﬁned as the heliospheric distance R(t)of

the barycenter of the plasma domain, normalized by the initial

distance R



() ( )aRtR tLL LL1.6

yyzz

000

==+==

In this equation, ò=da/dt is the expansion parameter,deﬁned

as the initial ratio between the characteristic expansion and

turnover times in the transverse directions (perpendicular to the

radial),



()

ku ,7

exp

rms

where

LL22

yyz

00 0 0

pp==is the initial minimum wave-

number in the transverse directions, and u

rms

is the initial rms

value of the velocity ﬂuctuations. The EBM equations ﬁnally

read, with dissipation terms omitted,



() ( ) ()ua2, 8

trr r¶+ =-



() () ()uuPPP Pa..2,9

tgg¶+ + =-



() ( ) ()

()

uuu BBPB a.2.,

t2rr¶+ + + - =-





() ()BuBBuBu a.. . ,11

¶+-+=-



()PT,12r=

where (

)

uu0, ,

=and ()



BBB2, ,

xyz

=. The nabla

operator that appears in the previous equations is written in

terms of comobile coordinates as ∇=(1/(a

)∂

,(1/a(t))∂

(1/a(t))∂

). The plasma is transported radially, which implies

that, in a local Cartesian coordinate system (x,y,z), where x

Figure 1. Sketch of the initial and ﬁnal plasma volumes after 10 nonlinear times in (a)standard homogeneous MHD simulations and (b)simulations with expansion

(see text).

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller

represents the local radial direction, it expands in directions y

and zperpendicular to the local radial direction (Figure 1(b)).

Expansion modiﬁes the evolution in two ways: (i)by

damping the different ﬁelds’amplitudes (see the right-hand

side terms in Equations (8)–(12)) and (ii)by damping the

gradients perpendicular to the radial due to the 1/a(t)factor in

the nabla expression. Explicit viscoresistive terms (not shown)

are added to the previous equations in order to dissipate the

energy that is transported along the spectrum down to the

smallest available scales. This energy is given to the internal

energy. Note that in the expanding case, the viscoresistive

terms are modiﬁed in order to minimize the decrease of the

Reynolds number with time (Montagud-Camps et al. 2018).

Note that if we choose ò=0 in Equations (8)–(12),we

recover the standard compressible MHD equations.

Finally, we deﬁne the so-called Elsässer variables,

() ()zu BBsign , 13

0dr=



where δB=B−B

.Wedeﬁne the normalized cross helicity

(hereafter simply called “cross helicity”)as the relative energy

excess in the z

mode:

()

.14

s=-

+=áñ-áñ

áñ+áñ

2.2. Numerics, Initial Conditions

Periodic boundary conditions are assumed in all directions of

a cube of side 2πin the comobile variables x,y, and z. We use a

standard pseudospectral code to numerically solve the above

3D EBM equations with a guide ﬁeld. All results are from runs

with a resolution of N

=512

grid points. A third-order

Runge–Kutta time integration scheme is used.

The initial state consists of nonzero amplitudes for the z

(k)

and z

−

(k)ﬁelds, excited within shells 1 kΔk, with k=|k|

(in units of 2πL

/l, with lthe wavelength), and weighted so that

the 1D spectral index is approximately m

. Initial incompressi-

bility is obtained by imposing k.z

(k)=0. Density and pressure

are uniform. The velocity ﬁelds and magnetic ﬂuctuations

are then deduced from the z

ﬁelds using the deﬁnition

(Equation (13)). Random phases are chosen for both z

and

−

ﬁelds, which implies that z

(k)and z

−

(k)are uncorrelated,

thus leading (since 〈z

−

〉=u

−b

/ρ)to ub

rms rms r.

Amplitudes are chosen so as to lead to the chosen initial cross

helicity c

Figure 1shows in a schematic way the evolution of the

plasma volume from start to end (after 10 nonlinear times)in

the two kinds of simulations: (a)homogeneous and (b)

expanding. In the latter case, the plasma volume is transported

from the initial heliocentric distance R

to 5R

, with

=0.2 au.

In the homogeneous case (a), the mean magnetic ﬁeld is

parallel to the xaxis and remains so. The aspect ratio of the

domain is adapted to the amplitude of the mean ﬁeld so that

rms

In the expanding case (b),wedonotstartwithameanﬁeld

aligned with the radial, as this is not the standard situation.

Instead, we adopt a slightly rotated conﬁguration (

B0.2

=).

In this way, due to the expansion of the plasma volume in the

transverse directions during transport, the conservation of

magnetic ﬂux through the expanding faces of the volume leads

=;i.e.,themeanﬁeld progressively rotates to an

angle of 45

with respect to the radial, a standard value at 1 au.

The initial domain (as well as the initial eddies)is chosen in most

of our simulations to be elongated in the radial direction (i.e., with

an aspect ratio a

=5), so that at the end, the domain aspect ratio

is unity. The interest of this choice is rationalized in Section 4.2.

Table 1presents the list of runs considered here. Parameters

are the expansion parameter ò(Equation (7)); the initial spectral

index m

, which is imposed in a wavenumber interval Δk

starting from the minimum wavenumber; the initial value of the

normalized cross helicity σ

,deﬁned as the relative excess of

dominant energy; the initial aspect ratio of the domain

; the mean magnetic ﬁeld Bt0

=; and the

dynamic viscosity μ

As already mentioned, the initial turbulent Mach number is

M=u

rms

/c=0.12, with c=((5/3)(P/ρ))

1/2

, where Pis the

gas pressure and ρis the density. Such a low Mach number

implies quasi-incompressibility, which allows one to consider

the two Elsässer energies as quasi-inviscid invariants, i.e.,

leading to separate cascades of the two energies. Most runs

Table 1

Initial Conditions for the Simulations

Run òm

Δkc

[10

−4

]

hA0, 5, 8 0 −5/3 16 0, 0.5, 0.8 [1, 0, 0]12

hB0, 5, 8 0 −1 16 0, 0.5, 0.8 [1, 0, 0]12

eA0, 5, 8 0.4 −5/3 64 0, 0.5, 0.83 [2, 0.4, 0]5 1.3

eB0, 5, 8 0.4 −1 64 0, 0.5, 0.83 [2, 0.4, 0]5 1.3

eC5, 8 0.4 −3/2 64 0.5, 0.83 [2, 0.4, 0]5 1.3

eE8 0.4 −1.25 64 0.83 [2, 0.4, 0]5 1.3

hA8k40−5/3 4 0.8 [1, 0, 0]11

eA8k16 0.4 −5/3 16 0.83 [2, 0.4, 0]5 1.3

eC8a1 0.4 −3/2 64 0.83 [2, 0.4, 0]1 1.3

eC8a3 0.4 −3/2 64 0.83 [2, 0.4, 0]3 1.3

hA8b50−5/3 16 0.8 [5, 0, 0]52

eC8b6 0.4 −3/2 64 0.83 [6, 1.2, 0]5 1.3

e2A0, 5, 8 0.2 −5/3 64 0, 0.5, 0.83 [2, 0.4, 0]5 1.3

e2B0, 5, 8 0.2 −1 64 0, 0.5, 0.83 [2, 0.4, 0]5 1.3

Note. From left to right: name of the run, expansion parameter ò, initial spectral index m

, initial power-law extent Δk, normalized cross helicity c

, initial mean

magnetic ﬁeld vector Bt

=, aspect ratio of the numerical domain a

, and dynamic viscosity μ

. Note that the ﬁrst four rows group together runs with three different

values of cross helicities. Resolution is 512

for all runs. All runs are integrated during 10 nonlinear times, except run eB8, which ends at t=20.

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller

have ò=0.4, which corresponds to about a 1 day timescale in

Helios data at 1 au, at solar minimum, while the value ò=0.2

corresponds to an 8 hr scale (Montagud-Camps et al. 2020).

With one exception, the value of the initial spectral slope m

−1, −5/3, or −3/2; as we will see, this leads to nonnegligible

differences in the evolution of the spectra. The effect of varying

the spectral extent Δk, aspect ratio a

(Equation (5)), and mean

magnetic ﬁeld will be discussed in detail in Section 4.

2.3. Diagnostic Tools, Initial Spectra

Simulations provide 3D energy spectra E

, from which we

deduce 1D energy spectra in the three xyz directions,

() ( ) ( )E k E k k k dk dk,, , 15

xDxyz yz

and similarly for the other two spectra, E(k

)and E(k

). The 1D

spectra will show either total energy (E

), that is, kinetic +

magnetic energy per unit mass, or, most of the time, the E

and

−

spectra separately. In the expanding case, the 0xdirection is

that of the radial, with the plane x0ycontaining the mean

magnetic ﬁeld direction.

In the following, we use either physical or normalized

wavenumbers. For physical wavenumbers (k

,ork

), we use

as a unit wavenumber the initial smallest transverse wave-

number,

()kL2, 16

00 0

where L

is the transverse size of the initial numerical domain.

Due to differences in both initial conditions and anisotropy

evolution, it will also prove useful to use normalized

wavenumbers when comparing spectra with and without

expansion. This simply consists of normalizing the abscissa

(either k

,ork

)by its minimum value when plotting

spectra.

Spectral indices will be deﬁned as

(( ) ( ) ( ) ( )mEkkEkkkklog log , 17

ii2121

== =



where k

can be either k

, or an equivalent of k

, for which

one takes the average of the two spectra E(k

)and E(k

). Note

that the interval [k

]will be ﬁxed in normalized coordinates.

We present in Figure 2the initial total energy spectra versus

physical wavenumbers k

(solid line),k

(dotted), and k

(dashed)for the four series of runs hAs,hBs,eAs, and eBs. All

spectra are compensated by k

−5/3

. One sees that for

homogeneous runs hAs and hBs, the 1D spectra occupy the

same wavenumber range in all three directions; this corre-

sponds to the domain aspect ratio a

=1. For runs with

expansion, the radial (E(k

)) spectrum is shifted with respect to

the two spectra E(k

)and E(k

), which corresponds to the

domain aspect ratio a

=5.

As seen in the ﬁgure, another difference between the

homogeneous and expanding simulations is the initial extent of

the excited spectrum, which is Δk=16 in the homogeneous

runs and Δk=64 in the expanding runs (see Table 1). It is seen

that in the ﬁgure, the scaling is constant only in the ﬁrst half of

the interval Δk, while the spectrum is becoming steeper in the

second half of the interval. Also, even in the ﬁrst half of the

excited interval, the spectra are actually somewhat steeper than

the imposed scaling

m0, as is clearly seen from the two cases

with m

=−5/3. The value m

is thus to be taken as a plain

indication, while meaningful evaluations of the spectral indices

(see Equation (17)) will be given when analyzing the time

evolution of the spectra in the next section.

3. Results

3.1. Varying c

and m

We ﬁrst consider the effect of varying the initial cross

helicity c

and spectral index m

3.1.1. Global Quantities

Here we consider global quantities versus time: the two

Elsässer energies

=á ñ

, the ratio b

rms

/|B

|, and the

normalized cross helicity σ

(t).

Figure 3shows the decay of the dominant and subdominant

energies E

and E

−

versus time. Results for the homogeneous

(hA

0,5,8

)and expanding (eA

0,5,8

)runs are given in the left and

right panels, respectively. In each panel, solid, dotted, and

dashed lines indicate, respectively, c

=0, 0.5, and 0.8 (for

homogeneous runs)or 0.83 (in the expanding case).As

expected, when c

is nonzero, cross helicity is seen to grow

with time for homogeneous runs and decay for expanding runs.

We consider in Figure 4the decay of energies E

versus radial

distance Rin the expanding case; in the left panel, we show the

three runs eA

0,5,8

with m

=−5/3, and in the right panel, we

show the three runs eB

0,5,8

with m

=−1. The 1/RWKB

energy decay appears as a short thin solid line. The two long

thin solid lines correspond to the mean E

energy decay found

by Bavassano et al. (2000,2002)in Helios data. Run eB8is

closer to these observations than run eA8; clearly, the energy

decay curves for both E

and E

−

depend on the initial

conditions.

Figure 5shows b

rms

/|B

|and σ

versus time stacked on one

another for runs hAs,hBs,eAs, and eBs, i.e., varying c

and m

Figure 2. Initial 1D total energy spectra for runs hAs,hBs,eAs, and eBs. Each

panel shows the 1D spectra in the three directions: E

)(solid),E

)

(dotted), and E

)(dashed). All spectra are compensated by k

−5/3

. As a rule,

)and E

)are almost indistinguishable.

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller

In each case, c

=0, 0.5, and 0.8 (or 0.83 in the expanding

case)are shown, respectively, with solid, dotted, and dashed

lines. For the homogeneous runs (top panels), the mean

magnetic ﬁeld is constant; thus, as expected, b

rms

/|B

|is

decaying with time. Also, increasing c

leads to slower

magnetic energy decay. This is not so surprising, since an

increase in c

leads to a smaller z

−

and thus a longer turnover

time for z

(and thus for b

rms

:()kz1

+-. On the contrary,

runs with expansion (bottom panels)show increasing b

rms

/|B

because expansion more strongly damps the radial component

of the mean ﬁeld than ﬂuctuations due to the linear Alfvén

coupling.

Note that, for the expanding case, the three curves do not

depend much on cross helicity. This means that the main source

of damping is expansion, not turbulent decay, as could be

suspected by remarking that the decay of the dominant energy

curve in Figure 4is not very different from the WKB 1/Rlaw.

Note also that the time evolution of b

rms

depends strongly

on the initial angle of the mean magnetic ﬁeld with the radial; it

is clear that starting with a large angle will lead to a much

slower decrease of B

, while the contrary will occur if the initial

angle is close to zero. The angle we considered here has been

chosen to match the average observed angle at 1 au, i.e., 45

Consider now cross helicity versus time. When starting with

a nonzero c

, homogeneous ﬂows (top panels)show a clear

increase of σ

due to the larger relative decay rate of the minor

species z

−

. This is explained by the fact that the minor species

cascade is driven by the dominant amplitude, thus leading to a

rapid decay of the minor species, called “dynamic alignment”

(Dobrowolny et al. 1980; Matthaeus & Montgomery 1980;

Grappin et al. 1982). On the contrary, ﬂows with expansion

(bottom panels)show a decrease of σ

with time (and thus

distance). Such a cross-helicity decrease was ﬁrst observed in

Helios and Voyager observations by Roberts et al. (1987).At

last, we note (and discuss later below)that the cross-helicity

decrease is slower when starting with ﬂatter spectra (compare

the two bottom panels, runs eA5,8 and eB5,8).

3.1.2. Spectral Quantities

We present in Figure 6a cut at k

=0 of the 3D E

spectra at

times zero and 10 for the two runs with initial large cross

helicity, with the homogeneous in the top panels (hA8)and

expansion in the bottom panels (eA8). The direction of the

mean magnetic ﬁeld is indicated by a straight line, which in the

expanding case rotates progressively so that it is 45

off the

radial at the end of run eA8. We see that the direction of the

cascade is perpendicular to B

in both runs, and that it rotates

with B

in the expanding case. Note that other anisotropies are

possible in the expanding case; starting with a quasi-isotropic

spectrum and a large enough expanding rate generally leads to

spectra aligned with the radial direction, a source of the second

branch in the so-called Maltese cross (Matthaeus et al. 1990;

Verdini & Grappin 2016; Montagud-Camps et al. 2020).

We next consider in Figure 7the 1D reduced spectra E

(k)

averaged during the interval 8 t10 for the four families of

runs hA0, 5, 8; hB0, 5, 8; eA0, 5, 8; and eB0, 5, 8. The spectra

(k)and E

−

(k)are compensated by k

−3/2

. Ordinates are

arbitrary, chosen to optimize readability. Recall that different

directions are used for the homogeneous and expanding runs.

For the homogeneous runs, we plot the average of the spectra

E(k

)and E(k

), i.e, spectra in directions perpendicular to the

mean magnetic ﬁeld, while for the runs with expansion, we plot

E(k

), that is, in the radial direction. In both cases, we

normalize the wavenumbers by their minimum values.

We added vertical dotted lines at k=4 and 10, which we

will consider from now on as markers of the inertial range

Figure 3. Varying cross helicity: total energies E

vs. time for homogeneous

runs hAs and runs with expansion eAs. Here c

=0, 0.5, and 0.8 (0.83 for the

run with expansion)are shown by solid, dotted, and dashed lines, respectively.

Figure 4. Varying cross helicity: total energies E

vs. heliocentric distance R

for runs with expansion eAs and eBs. Styles are the same as in Figure 3. The

two long thin solid lines are the mean energies E

and E

−

as found by

Bavassano et al. (2000,2002)from Helios data, and the short thin solid line is

the 1/R(WKB)decay.

Figure 5. Varying cross helicity and initial spectral index b

rms

(top)and

cross helicity σ

vs. time. Styles are the same as in Figure 3.

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller

boundaries, since spectra can be approximated by power laws

in this range.

Examining the different spectra shows interesting differences

in the distribution of cross helicity along the spectral range and

the dependence (or not)on initial conditions. First, in runs with

nonzero σ

, cross helicity is concentrated on the largest scales

in the homogeneous runs, while in the expanding runs, it is

either uniformly distributed or showing a minimum at the

largest scale. Second, in the homogeneous runs, the ﬁnal slope

does not depend on the initial slope, while in the runs with

expansion, it does.

These properties are well visible in Figure 8, where we show

the evolution of the two spectral indices m

and m

−

with time,

computed in the interval 4 k10. Cross helicity grows from

left to right, with the homogeneous runs appearing in the top

panels and runs with expansion in the bottom panels. Each

panel allows one to compare the evolution of indices m

(solid

lines)and m

−

(dotted lines)with two different starting indices:

=−5/3 in double thick green lines and −1 in thick black

lines.

Homogeneous and expanding runs show completely differ-

ent evolutions. Expanding runs are almost independent of cross

helicity and show simple properties: (i)the two indices

converge toward values that depend strongly on their initial

values, and (ii)the two indices m

and m

−

remain about equal.

For homogeneous runs, we see that (i)m

and m

−

converge

toward values that do not depend on their initial value, and (ii)

when cross helicity is substantial, indices become unequal, with

|>|m

−

|; i.e., E

becomes steeper than E

−

. This behavior

was ﬁrst found in closure calculations by Grappin et al. (1983)

and commonly denoted as “spectral pinning”(Lithwick &

Goldreich 2003), as discussed below.

Figure 6. Cuts at k

=0of3DE

=0)spectra at time t=0 and 10.

Straight lines indicate the direction of the mean magnetic ﬁeld. Top panels:

homogeneous run hA8. Bottom panels: run with expansion eA8. The unit

wavenumber is L

p, where L

is the initial transverse size of the domain.

Figure 7. Varying c

and m

. Shown are the 1D reduced energy spectra E

(k)

and E

−

(k)in normalized coordinates, compensated by k

−3/2

and averaged in

the time interval 8 t10. The ordinates are arbitrary to maximize

readability. The runs are hAs,hBs,eAs, and eBs. Styles are the same as in

Figure 3. Vertical dotted lines (k=4 and 10 in normalized coordinates)

indicate the wavenumber interval chosen to compute spectral indices in the

following.

Figure 8. Varying c

and m

. Shown are indices m

and m

−

vs. time (solid and

dotted lines, respectively)computed in the normalized wavenumber interval

4k10. Each panel shows two runs with a given c

(from left to right, 0,

0.5, 0.8, 0, 0.5, and 0.83). Runs with m

=−5/3 appear as double thick green

lines, and runs with m

=−1 appear as thick black lines. Solid horizontal lines

mark spectral indices −3/2 and −5/3.

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller

3.2. Spectral Extent Δk, Mean Magnetic Field, and Expansion

Parameter ò

In series hAs,hBs,eAs, and eBs analyzed above, we varied

and m

. A look at these runs in Table 1also shows that a

Δk, and B

are systematically different in the homogeneous

and expanding runs. We thus consider here the successive

effect of varying Δk,Bt0

=, as well as an intermediate value of

the expansion parameter òbetween 0 and 0.4. We postpone the

case of the initial aspect ratio of the domain a

to the Section 4.

In Figure 9, we consider two runs, one homogeneous (hA8)

and one with expansion (eA8)with high cross helicity and

=−5/3. Each one is compared with a similar run with Δk

four times smaller (hA8k4 and eA8k16). The left panels show

the spectra at time t=10, and the right panels show the spectral

indices m

versus time. Homogeneous runs are shown in the

top panels, expanding runs in the bottom panels. Runs with

large Δkappear in green lines, and runs with reduced Δk

appear in black lines.

Decreasing Δkproduces initially steeper spectra, as seen in

the right panels (black lines), for both the homogeneous and

expanding runs. This is due to the initial lack of excitation in

the reference interval [4, 10]. However, at t=10, the ﬁnal

slopes become comparable to those of runs with the standard

values of Δk, as seen in the right panels. We conclude that the

choice of different values of Δk(16 and 64)for the

homogeneous and expanding runs are of minor importance,

as relaxation toward the same turbulent state is ensured

whatever the value of Δk, albeit more or less rapidly.

In Figure 10, we show the effect of increasing the initial

mean magnetic ﬁeld from B

=(1, 0, 0)to (5, 0, 0)in the

homogeneous case (runs hA8 and hA8b5)and from B

=(2,

0.4, 0)to (6, 1.2, 0)in the expanding case (runs eC8 and

eC8b6; see Table 1).

In both the homogeneous and expanding cases, it is seen that

increasing the mean magnetic ﬁeld leads to steeper spectra. The

effect is stronger in the homogeneous case. Such an effect has

been observed in several works investigating the transition

from strong to weak anisotropic cascade, albeit with zero cross

helicity (Dmitruk et al. 2003; Rappazzo et al. 2007; Perez &

Boldyrev 2008). Note also that the spectral pinning (i.e.,

|>|m

−

|)we have observed previously in homogeneous

runs is seen to appear here in expanding runs, as well when the

mean ﬁeld is large enough.

We ﬁnally show in Figure 11 the indices m

(t)with half our

standard value of the expansion rate, i.e., ò=0.2. We consider

the usual three values of cross helicity, 0, 0.5, and 0.83, and the

two starting indices m

=−5/3(green lines)and −1(black

lines). The resulting curves are globally similar to those shown

in the bottom panels of Figure 8; however, one sees for runs

with nonzero c

and m

=−5/3 the appearance of some

pinning (|m

|>|m

−

|), albeit more moderate than previously

seen in homogeneous runs (top panels in Figure 8).

4. Discussion

We have compared in this paper the results obtained on the

turbulent evolution using two models: (i)homogeneous MHD

and (ii)MHD with expansion. We ﬁrst summarize the

properties predicted by the two models, and then we discuss

the validity of the parameters chosen.

Figure 9. Varying the initial spectral extent Δk:E

spectra compensated by

3/2

at t=10 (left panels)and spectral indices (right panels). Top panels:

homogeneous runs hA4(green lines)and hA8k4(with reduced Δk; black lines).

Bottom panels: expanding runs eA8(green lines)and eA8k16 (with reduced

Δk; black lines). Spectra are compensated by k

−3/2

, with vertical dotted lines

marking the inertial range. Solid and dotted lines show the m

and m

−

indices,

respectively.

Figure 10. Increasing the initial mean magnetic ﬁeld B

: spectral indices m

(solid lines)and m

−

(dotted lines)vs. time. Starting at the top left panel, we

show (i)homogeneous runs hA8 with B

=(1, 0, 0),(ii)hA8b5 with B

=(5,

0, 0),(iii)runs with expansion eC8 with B

=(2, 0.4, 0), and (iv)eC8b6 with

=(6, 1.2, 0). Homogeneous runs have 0.8

=, and runs with expansion

have 0.8

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller

4.1. Contrasted Properties of the Two Models

In homogeneous MHD simulations, spectral pinning is

found, that is, a steeper dominant component E

, with the

difference between the two indices m

and m

−

increasing when

cross helicity is larger. In other words, the spectral indices

depend on total cross helicity. This is not compatible with the

usual picture proposed by the strong perpendicular cascade

(that is, perpendicular to the mean magnetic ﬁeld), for which

both spectra E

and E

−

are parallel, with m

−

=−5/3

(Lithwick et al. 2007); instead, it is compatible with the weak

isotropic cascade prediction of the generalized Iroshnikov–

Kraichnan theory (Grappin et al. 1983), for which m

−

−3. This value is found in our homogeneous runs when the

mean magnetic ﬁeld is not too much larger than unity. The

pinning phenomenon has been found in previous MHD

simulations. First, in 2D incompressible MHD simulations

(with no mean ﬁeld), the weak isotropic version is found

−

=−3)in Pouquet et al. (1988). More recently, Perez

& Boldyrev (2010)found it but, to explain why solar wind

observations show no spectral pinning, proposed that the solar

wind turbulence are forced to exhibit parallel E

spectra due to

a quasi-inﬁnite inertial range resulting from an inﬁnite

Reynolds number.

We have seen that spectral pinning is absent from our

simulations with expansion as soon as the expansion parameter

is large enough and the mean magnetic ﬁeld is not much larger

than the b

rms

amplitude. In that case, the two spectra E

are

indeed more or less parallel, i.e., m

−

, as observed in solar

wind spectra (Grappin et al. 1990; Chen et al. 2020). As a rule,

initially steep spectra (e.g., m

;−5/3)remain steep at larger

distances, while initially ﬂat spectra (e.g., m;−1)steepen and

relax to intermediate index values, e.g., −1.5. Again, a similar

behavior is observed in solar wind data (Chen et al. 2020).In

the homogeneous case, on the contrary, the relaxed spectral

index is independent of the initial spectral index, as expected in

a turbulent state.

In the rest of the discussion below, we consider in turn the

following issues:

(i)the choice of the initial domain aspect ratio for the

expanding cases;

(ii)the difference between radial and perpendicular directions

in the expanding cases;

(iii)the relation between spectral index and cross correlation;

(iv)the origin of the absence of pinning in the solar wind; and

(v)the reason that the spectral index does not always

converge to −5/3 in expanding runs.

4.2. Domain Aspect Ratio

In the simulations with expansion examined up to now, we

started with a numerical domain extended in the radial direction

by a factor of 5 (a

=5), thus arriving at the end at a cubic

domain. A cubic domain is a priori supposed to favor the local

nonlinear interactions, which are necessary to ensure a

“normal”turbulent cascade. However, the opposite choice

=1)is also possible; in that case, local interactions are

maximum at the beginning of the evolution, not at the end. We

compare in Figure 12 the two cases at the end of the simulation,

t=10, further adding a run with the intermediate value a

=3.

The ﬁgure shows E

spectra along k

(solid lines)and k

(dotted lines)for the runs with expansion (from left to right):

eC8a1, eC8a3, and eC8 with a

=1, 3, and 5, respectively. The

inertial range is marked by dotted lines as previously. Note that

we use physical wavenumbers, so that the inertial range is

given by 4/5k10/5(the factor 1/5 is due to the fact that

the ﬁnal transverse size is ﬁve times the initial one in all

three runs).

We see that, for both a

=5 and 3, turbulence at t=10 is

quasi-isotropic, which means that the coupling between the

radial direction and the perpendicular plane is efﬁcient, while it

is much less so when a

=1, which is probably due to the too-

large importance of nonlocal interactions in this case, at least at

the end of the computation. From this, we conclude that

choosing a

between 3 and 5 is a reasonable choice. See

Montagud-Camps et al. (2020)for a similar discussion in the

context of turbulent dissipation with a large turbulent Mach

number.

4.3. Radial versus zDirection

In simulations with expansion, we have chosen to measure

spectral scaling in the radial direction to allow direct

comparison with observations. It is of interest to also examine

the scaling in directions perpendicular to the mean ﬁeld for the

sake of comparison with the scaling obtained in homogeneous

runs and because the preferred cascade develops perpendicu-

larly to the mean magnetic ﬁeld (Figure 6).

In Figure 13, we compare the spectral evolution along

directions k

(left panels)and k

(right panels)at times t=0,

2, K, 10 in the two expanding runs eA8 and eB8. The E

spectra are compensated by k

−3/2

We see that the interval [4, 10](in normalized coordinates)is

a reasonable candidate for the inertial range in both the k

and

Figure 11. Intermediate expansion rate ò=0.2. Shown is the time evolution of

indices m

and m

−

, computed in 4 k/k

10 for six runs, with growing

initial cross helicity from left to right. See Figure 8for styles. Green lines

denote runs with m

=−5/3, and black lines are runs with m

=−1.

Figure 12. Varying the domain aspect ratio a

in expanding runs eC8a1,

eC8a3, and eC8. Shown are the 1D spectra at t=10 in different directions

using physical wavenumbers E

)(solid lines)and E

)(dotted lines).

From left to right, a

=1, 3, and 5. The two vertical dotted lines mark the

boundaries of the inertial range adopted above; see text. The unit length

is ()

0p.

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller

directions. Actually, the spectral indices do not seem very

different. Figure 14 shows the spectral indices for the two runs

separately: left, computed in the k

direction, and right, in the k

direction. One sees that in the k

case, the indices m

show

larger ﬂuctuations, which is clearly related to the spectra being

more irregular. We thus conclude that the differences we have

observed up to now in spectral behavior between the

homogeneous and expanding runs are not due to our choice

of selecting the k

direction, i.e., not a direction perpendicular

to the mean ﬁeld, in the expanding case.

4.4. Comparison with Helios Data

We attempt here to compare the evolution of some properties

of the solar wind turbulence as observed during the Helios 1

mission with those shown by MHD with expansion. A ﬁrst

property is the combined decay of cross helicity with distance

and that of the spectral slope. The second is the relation

between cross helicity and the spectral index.

The observed turbulent evolution is illustrated in Figure 15

by four averaged spectra E

and E

−

, compensated by f

−5/3

The left panel shows fast stream averages (>500 km s

−1

), and

the right panel shows slow stream averages (<500 km s

−1

Each panel shows two couples of spectra E

, one close to the

Sun (0.3 au <R<0.5 au; solid lines)and one far from the Sun

(0.8 au <R<1 au; dotted lines).

Visual inspection shows that (i)cross helicity decreases with

distance in either fast or slow streams; (ii)in either fast or slow

streams, spectra steepen with distance; (iii)either close to the

Sun or (see Chen et al. 2013)close to 1 au, the spectra are

ﬂatter when cross helicity is higher; and (iv)the previous

property is not built during transport from 0.3 to 1 au but is

already present at 0.3 au. Hence, this relation between spectral

index and cross helicity is a by-product of processes acting

closer to the Sun, either the low corona or the accelerating

region.

We now turn to the simulations with expansion. We show in

Figure 16 the index m

versus time for two series of runs. In

the left panel, we show runs with high cross helicity σ

=0.83,

and we show lower cross helicity 0.5

=in the right panel. In

each case, we vary the initial spectral index m

from −1to−5/

3. The same style is adopted in the two panels for each pair of

runs with the same value of m

. More precisely, the runs

considered for high cross helicity in the left panel are eB8, eE8,

eC8, and cA8, and those for moderate σ

in the right panel are

eB5, eC5, and eA5.

The insert in each panel shows the evolution of cross helicity

with time for the speciﬁc runs eE8 and eC5, denoted

respectively in the ﬁgures by the letters F and S for “fast”

and “slow,”to be explained below. Note that the cross-helicity

curves shown are representative of the other runs in each of the

two panels. So, as we already knew, property (i)is shared with

Helios data.

We can also see that m

(t)is systematically decreasing with

time, except when starting from −5/3 or close to it. This

conclusion might have been guessed already from Figure 8;

however, here the larger ensemble of initial indices allows a

more secure conclusion. So, property (ii)is also shared by

Helios data.

Figure 13. The E

(thick lines)and E

−

(thin lines)spectra along the radial (k

;

left panels)and perpendicular (k

; right panels)directions for two runs with

expansion: eA8(top)and eB8(bottom). Times t=0, 2, K, 10, spectra are

compensated by k

−3/2

. Dotted lines are the boundaries of the interval used to

compute the spectral indices.

Figure 14. Spectral indices m

(solid lines)and m

−

(dotted lines)vs. time

computed in the radial (4k

10; left panel)and z(4k

10; right panel)

directions. Shown are runs eA8(thick green lines)and eB8(thin black lines).

Figure 15. The E

spectra of the Helios 1 mission during the ﬁrst 4 months:

relation between spectral index, distance, and cross helicity/wind speed. Fast

streams are shown in the left panel, and slow streams are shown in the right

panel. Spectra are averaged separately close to the Sun (solid lines)and close to

1au(dotted lines); see text for details. Vertical dotted lines mark the frequency

interval [f

, 2.5f

]with f

=3×10

−4

Hz (1/f

;1hr).

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller

Note that assuming a wind that would be made with equal

weights of all runs shown in the ﬁgure would produce no

correlation at all between the spectral index and the cross

helicity at either t=0 or 10; property (iii)should thus not be

recovered. So, let us introduce the correlation at t=0, for

instance, by considering the distribution made of just the two

runs F and S, deﬁned previously. Run F has a σ

larger than

that of run S, as well as a spectrum ﬂatter than that of run S, at

both t=0 and 10 (and in between). As we see in the ﬁgure,

runs F and S at time t=10 still satisfy the required property.

So, property (iii)is probably satisﬁed, assuming that the solar

wind turbulent stream is made of a family of streams with

initial properties close to those of runs F and S.

4.5. Why Does Pinning Not Appear in Solar Wind Turbulence?

Figure 17 suggests the solution. It shows in the two left

panels the familiar E

spectra at t=10 for the homogeneous

run hA8(top)and the expanding run eA8(bottom), both

compensated by k

−5/3

. In the two right panels, we plot the

magnetic (thick line)and kinetic spectra (thin line). In the

expanding run (and not in the homogeneous run), there is a

clear magnetic excess at the largest scale. Correspondingly, in

the expanding run, one sees that the E

−

spectrum also has an

“excess”at the largest scale. The interpretation is simple: at the

largest scale, the magnetic excess algebraically leads to an

excess of E

−

. There is thus a permanent large-scale injection of

zero cross helicity at the largest scale.

Our ﬁnal interpretation is that this large-scale injection of a

magnetic excess is the source of (i)the decrease with time of

the relative cross helicity in the expanding plasma and (ii)the

absence of pinning.

Now, how is the magnetic energy excess injected at large

scales? A critical scale is the WKB scale, below which

magnetic and velocity ﬂuctuations decrease at the same rate

with distance (δu=δb/ρ

1/2

;1/R

1/2

). Above this critical

WKB scale, expansion dominates; i.e., the Alfvén coupling is

inefﬁcient, and the velocity ﬂuctuations damp globally faster

than the magnetic ﬂuctuations due to the conservation of

magnetic ﬂux, angular momentum, and radial momentum

(Zhou & Matthaeus 1990; Grappin et al. 1993; Oughton &

Matthaeus 1995; Dong et al. 2014). This is sufﬁcient to

generate the magnetic excess as soon as the expansion is large

enough for the largest scale of the simulation to indeed be non-

WKB. Figures 10 and 11 conﬁrm this, as they show that partial

pinning reappears when either the expansion parameter is

decreased or the mean magnetic ﬁeld is increased.

On the contrary, whenever the large non-WKB scales gain

importance, either because the initial mean magnetic ﬁeld is

decreased or as the expansion rate is increased, one ﬁnds (not

shown)a large-scale increase of the magnetic excess and an

associated acceleration of the normalized cross-helicity

decrease.

Note that other large-scale sources of low or zero cross

helicity are present in the solar wind, such as shear ﬂows

associated with stream structures (Roberts et al. 1991). Stream

structures certainly contribute to the decrease of cross helicity

with distance, but locally, while the expansion contribution is

ubiquitous in the solar wind.

4.6. Why Does the Spectral Index Not Always Converge to

−5/3?

When expansion is present, we have seen that the spectral

index converges to a value that depends on the initial index m

except when the initial slope is close enough to the value −5/3.

A possible explanation is that, since expansion damps the

turbulent amplitude in addition to turbulent dissipation, the

nonlinear time is increased; thus, the evolution is slowed down

compared to the homogeneous case.

The evolution of the spectral indices m

is shown in

Figure 18 up to time t=20 for run eB8. One sees that the two

indices show no clear drift toward −5/3 in the interval

6t20, remaining at values close to −3/2. To examine

whether the above explanation is acceptable, we deﬁne the

“age”of the turbulence as the number Nof nonlinear times t

accumulated since the beginning of the calculation. It can be

expressed as

() () () ( )Nt k t u t dt,18

0rms

=¢¢¢

Figure 16. Spectral index m

evolution for high and low cross helicity:

0.8

=in the left panel and 0.5 in the right panel. In each panel, a series of

initial indices m

is considered. Inserts show cross-helicity evolution for runs F

(left)and S (right).

Figure 17. Origin of the suppression of pinning in the expanding case:

homogeneous run hA8(top row)and expanding run eA8(bottom). The left

panels show E

spectra at t=10, compensated by k

−5/3

. The right panels

show Ev (thin lines)and Eb (thick lines)spectra.

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller

where k

and u

rms

are, respectively, the smallest wavenumber

and the rms velocity ﬂuctuation.

The age versus time is shown in Figure 19 for four runs with

large cross helicity without and with expansion: hA8, hB8, eA8,

and eB8 with black, blue, red, and green solid lines,

respectively. The difference between the two homogeneous

runs and the two runs with expansion is drastic. The two

homogeneous runs have their age growing almost linearly with

time, from t=0 up to 10, due to (i)the constancy of the largest

scale and (ii)the limited decay of the largest eddies with time.

The two runs with expansion, on the contrary, show a quasi-

saturation of their age after time t;8. This is due to both a

higher-amplitude decay rate and the linear expansion of the

plasma with time in the two directions transverse to radial

(



()kt a RR t11

µ= µ+; Equation (6)).

To evaluate the relative importance of these two effects on

the freezing, we plotted with dotted lines the age obtained when

taking into account only the turbulence amplitude (i.e.,

replacing ()

0with ()

0=)for runs eA8 and eB8. The

resulting curves show a decrease compared to the two curves

for homogeneous runs but clearly do not correspond to a

saturation of the age. One may thus conclude that the main

origin of the quasi-freezing of the spectral evolution is due

principally to the expansion of the largest eddies and

secondarily to the indirect effect of expansion, which is the

increased decay of the turbulent amplitude.

Note that in the simple, nonturbulent case of the shock

formation of an acoustic wave with a wavevector lying within

the perpendicular plane, one analytically ﬁnds that the increase

of the effective expansion, when initially large enough, stops

the wave from steepening before the shock is formed (Grappin

et al. 1993).

5. Conclusion

In this paper, we considered decaying turbulence, starting

with random ﬂuctuations of velocity and magnetic ﬁeld, and

varying initial cross helicity and spectral index. Our key results

can be summarized as follows (we denote results obtained

using homogeneous runs with (H)and results obtained using

runs with expansion with (E)).

1. The m

indices converge to values independent of the

initial conditions (H); on the contrary, the indices keep a

memory of the initial conditions (E).

2. Pinning occurs with nonzero c

; as a result, |m

|/|m

−

relaxes to a value that grows with c

(H). On the contrary,

as a rule, m

−

, probably due to a permanent

injection of the subdominant mode E

−

at the largest

scale by expansion (E).

3. The memory of initial conditions allows one to under-

stand the conservation of the correlation between high

cross helicity and ﬂat spectra during the transport

between 0.2 and 1 au, but the generation of the

correlation probably occurs close to the Sun, requiring

physical effects that are absent from the present

simulations (E).

4. When increasing B

, pinning persists, but both E

spectra

steepen, as already known in the 0

=case (H). Some

pinning appears, and the spectra also steepen (E).

It would be interesting in future work to consider smaller

values of the expansion parameter in order to be directly

applicable to the inertial range. Also, increasing the resolution

would allow one to describe the break separating the large-

scale 1/frange from the true inertial range.

This work was granted access to the HPC resources of

IDRIS under allocation 2020-A0090407683 made by GENCI.

R.G. acknowledges fruitful discussions on the manuscript with

Thierry Passot, Romain Meyrand, Victor Montagud-Camps,

Nicolas Aunai, and Davide Manzini.

ORCID iDs

Roland Grappin https://orcid.org/0000-0001-7847-3586

References

Adhikari, L., Zank, G. P., Zhao, L. L., et al. 2020, ApJS,246, 38

Bavassano, B., & Bruno, R. 1995, JGR,100, 9475

Bavassano, B., Dobrowolny, M., Mariani, F., & Ness, N. F. 1982, JGR,

87, 3617

Bavassano, B., Pietropaolo, E., & Bruno, R. 2000, JGR,105, 15959

Bavassano, B., Pietropaolo, E., & Bruno, R. 2002, JGR,107, 1452

Boldyrev, S. 2005, ApJL,626, L37

Boldyrev, S., Perez, J. C., Borovsky, J. E., & Podesta, J. J. 2011, ApJL,

741, L19

Figure 18. Evolution of spectral indices up to time t=20, run eB8. Figure 19. Age of turbulence for four runs with cross helicity 0.8

=. Run

hA8 is black, hB8 blue, eA8 red, and eB8 green. Solid lines correspond to age

deﬁned as in Equation (18), and dotted lines correspond to replacing

wavenumber (

0with its initial value in Equation (18), i.e., not taking

expansion into account in the age evaluation.

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller

Chandran, B. D. G., Dennis, T. J., Quataert, E., & Bale, S. D. 2011, ApJ,

743, 197

Chandran, B. D. G., Schekochihin, A. A., & Mallet, A. 2015, ApJ,807, 1

Chen, C., Bale, S., Bonnell, J. W., et al. 2020, ApJS,246, 53

Chen, C., Bale, S., Salem, C. S., & Maruca, B. A. 2013, ApJ,770, 125

Chen, C., Mallet, A., Schekochihin, A. A., et al. 2012, ApJ,758, 120

Cranmer, S. R., van Ballegooijen, A. A., & Edgar, R. J. 2007, ApJS,171, 520

Dmitruk, P., Gomez, D. O., & Matthaeus, W. H. 2003, PhPl,10, 3584

Dobrowolny, M., Mangeney, A., & Veltri, P. 1980, PhRvL,45, 144

Dong, Y., Verdini, A., & Grappin, R. 2014, ApJ,793, 118

Grappin, R., Frisch, U., Pouquet, A. G., & Léorat, J. 1982, A&A, 105, 6

Grappin, R., Léorat, J., & Pouquet, A. G. 1983, A&A, 126, 51

Grappin, R., Mangeney, A., & Marsch, E. 1990, JGR,95, 8197

Grappin, R., Velli, M., & Mangeney, A. 1991, AnGeo, 9, 416

Grappin, R., Velli, M., & Mangeney, A. 1993, PhRvL,70, 2190

Horbury, T. S., Forman, M., & Oughton, S. 2008, PhRvL,101, 175005

Lithwick, Y., & Goldreich, P. 2003, ApJ,582, 1220

Lithwick, Y., Goldreich, P., & Sridhar, S. 2007, ApJ,655, 269

Matsumoto, T. 2021, MNRAS,500, 4779

Matthaeus, W. H., Goldstein, M. L., & Roberts, D. A. 1990, JGR,95, 20673

Matthaeus, W. H., Klein, L. W., Ghosh, S., & Brown, M. R. 1991, JGR,

96, 5421

Matthaeus, W. H., & Montgomery, D. 1980, in Int. Conf. on Nonlinear

Dynamics 357 (New York: New York Academy of Sciences),203

Montagud-Camps, V., Grappin, R., & Verdini, A. 2018, ApJ,853, 153

Montagud-Camps, V., Grappin, R., & Verdini, A. 2020, ApJ,902, 34

Ofman, L. 2010, LRSP,7, 4

Oughton, S., & Matthaeus, W. H. 1995, JGR,100, 14783

Perez, J. C., & Boldyrev, S. 2008, ApJ,672, L61

Perez, J. C., & Boldyrev, S. 2010, ApJL,710, L63

Pouquet, A. G., Sulem, P., & Meneguzzi, M. 1988, PhFl,31, 2635

Rappazzo, A. F., Velli, M., Einaudi, G., & Dahlburg, R. B. 2007, ApJ,

657, L47

Réville, V., Velli, M., Panasenco, O., et al. 2020, ApJS,246, 24

Roberts, D. A., Ghosh, S., Goldstein, M. L., & Matthaeus, W. H. 1991, PhRvL,

67, 3741

Roberts, D. A., Goldstein, M. L., Klein, L. W., & Matthaeus, W. H. 1987, JGR,

92, 12023

Shoda, M., Suzuki, T. K., Asgari-Targhi, M., & Yokoyama, T. 2019, ApJL,

880, L2

Shoda, M., Yokoyama, T., & Suzuki, T. K. 2018, ApJ,853, 190

Squire, J., Chandran, B. D. G., & Meyrand, R. 2020, ApJL,891, L2

van der Holst, B., Sokolov, I. V., Meng, X., et al. 2014, ApJ,782, 81

Verdini, A., & Grappin, R. 2015, ApJL,808, L34

Verdini, A., & Grappin, R. 2016, ApJ,831, 179

Verdini, A., Velli, M., Matthaeus, W. H., Oughton, S., & Dmitruk, P. 2010,

ApJL,708, L116

Zhou, Y., & Matthaeus, W. H. 1990, JGR,95, 10291

The Astrophysical Journal, 933:246 (12pp), 2022 July 10 Grappin, Verdini, & Müller