## 5 Numerical Simulations

Numerical simulations currently represent one of the main source of information about non-linear evolution of fluid flows. The actual super-computers are now powerful enough to simulate equations (NS or MHD) that describe turbulent flows with Reynolds numbers of the order of 10^{4}in two-dimensional configurations, or 10

^{3}in three-dimensional one. Of course, we are far from achieving realistic values, but now we are able to investigate turbulence with an inertial range extended for more than one decade. Rather the main source of difficulties to get results from numerical simulations is the fact that they are made under some obvious constraints (say boundary conditions, equations to be simulated, etc.), mainly dictated by the limited physical description that we are able to use when numerical simulations are made, compared with the extreme richness of the phenomena involved: numerical simulations, even in standard conditions, are used tout court as models for the solar wind behavior. Perhaps the only exception, to our knowledge, is the attempt to describe the effects of the solar wind expansion on turbulence evolution like, for example, in the papers by Velli et al. (1989, 1990); Hellinger and Trávníček (2008). Even with this far too pessimistic point of view, used here solely as a few words of caution, simulations in some cases were able to reproduce some phenomena observed in the solar wind.

Nevertheless, numerical simulations have been playing a key role, and will continue to do so in our seeking an understanding of turbulent flows. Numerical simulations allows us to get information that cannot be obtained in laboratory. For example, high resolution numerical simulations provide information at every point on a grid and, for some times, about basic vector quantities and their derivatives. The number of degree of freedom required to resolve the smaller scales is proportional to a power of the Reynolds number, say to , although the dynamically relevant number of modes may be much less. Then one of the main challenge remaining is how to handle and analyze the huge data files produced by large simulations (of the order of Terabytes). Actually a lot of papers appeared in literature on computer simulations related to MHD turbulence. The interested reader can look at the book by Biskamp (1993) and the reviews by Pouquet (1993, 1996).

### 5.1 Local production of Alfvénic turbulence in the ecliptic

The discovery of the strong correlation between velocity and magnetic field fluctuations has represented the motivation for some MHD numerical simulations, aimed to confirm the conjecture by Dobrowolny et al. (1980b). The high level of correlation seems to be due to a kind of self-organization (dynamical alignment) of MHD turbulence, generated by the natural evolution of MHD towards the strongest attractive fixed point of equations (Ting et al., 1986; Carbone and Veltri, 1987, 1992). Numerical simulations (Carbone and Veltri, 1992; Ting et al., 1986) confirmed this conjecture, say MHD turbulence spontaneously can tends towards a state were correlation increases, that is, the quantity , where is the cross-helicity and the total energy of the flow (see Appendix B.1), tends to be maximal.

The picture of the evolution of incompressible MHD turbulence, which comes out is rather nice but solar wind turbulence displays a more complicated behavior. In particular, as we have reported above, observations seems to point out that solar wind evolves in the opposite way. The correlation is high near the Sun, at larger radial distances, from 1 to 10 AU the correlation is progressively lower, while the level in fluctuations of mass density and magnetic field intensity increases. What is more difficult to understand is why correlation is progressively destroyed in the solar wind, while the natural evolution of MHD is towards a state of maximal normalized cross-helicity. A possible solution can be found in the fact that solar wind is neither incompressible nor statistically homogeneous, and some efforts to tentatively take into account more sophisticated effects have been made.

A mechanism, responsible for the radial evolution of turbulence, was suggested by Roberts and Goldstein (1988); Goldstein et al. (1989); and Roberts et al. (1991, 1992) and was based on velocity shear generation. The suggestion to adopt such a mechanism came from a detailed analysis made by Roberts et al. (1987a,b) of Helios and Voyager interplanetary observations of the radial evolution of the normalized cross-helicity at different time scales. Moreover, Voyager’s observations showed that plasma regions, which had not experienced dynamical interactions with neighboring plasma, kept the Alfvénic character of the fluctuations at distances as far as 8 AU (Roberts et al., 1987b). In particular, the vicinity of Helios trajectory to the interplanetary current sheet, characterized by low velocity flow, suggested Roberts et al. (1991) to include in his simulations a narrow low speed flow surrounded by two high speed flows. The idea was to mimic the slow, equatorial solar wind between north and south fast polar wind. Magnetic field profile and velocity shear were reconstructed using the six lowest Fourier modes as shown in Figure 67. An initial population of purely outward propagating Alfvénic fluctuations () was added at large and was characterized by a spectral slope of . No inward modes were present in the same range. Results of Figure 67 show that the time evolution of spectrum is quite rapid at the beginning, towards a steeper spectrum, and slows down successively. At the same time, modes are created by the generation mechanism at higher and higher but, along a Kolmogorov-type slope .

These results, although obtained from simulations performed using 2D incompressible spectral and pseudo-spectral codes, with fairly small Reynolds number of , were similar to the spectral evolution observed in the solar wind (Marsch and Tu, 1990a). Moreover, spatial averages across the simulation box revealed a strong cross-helicity depletion right across the slow wind, representing the heliospheric current sheet. However, magnetic field inversions and even relatively small velocity shears would largely affect an initially high Alfvénic flow (Roberts et al., 1992). However, Bavassano and Bruno (1992) studied an interaction region, repeatedly observed between 0.3 and 0.9 AU, characterized by a large velocity shear and previously thought to be a good candidate for shear generation (Bavassano and Bruno, 1989). They concluded that, even in the hypothesis of a very fast growth of the instability, inward modes would not have had enough time to fill up the whole region as observed by Helios 2.

The above simulations by Roberts et al. (1991) were successively implemented with a compressive
pseudo-spectral code (Ghosh and Matthaeus, 1990) which provided evidence that, during this turbulence
evolution, clear correlations between magnetic field magnitude and density fluctuations, and between
and density fluctuations should arise. However, such a clear correlation, by-product of the non-linear
evolution, was not found in solar wind data (Marsch and Tu, 1993b; Bruno et al., 1996). Moreover,
their results did not show the flattening of e^{–} spectrum at higher frequency, as observed by
Helios (Tu et al., 1989b). As a consequence, velocity shear alone cannot explain the whole
phenomenon, other mechanisms must also play a relevant role in the evolution of interplanetary
turbulence.

Compressible numerical simulations have been performed by Veltri et al. (1992) and Malara et al. (1996, 2000) which invoked the interactions between small scale waves and large scale magnetic field gradients and the parametric instability, as characteristic effects to reduce correlations. In a compressible, statistically inhomogeneous medium such as the heliosphere, there are many processes which tend to destroy the natural evolution toward a maximal correlation, typical of standard MHD. In such a medium an Alfvén wave is subject to parametric decay instability (Viñas and Goldstein, 1991; Del Zanna et al., 2001; Del Zanna, 2001), which means that the mother wave decays in two modes: i) a compressive mode that dissipates energy because of the steepening effect, and ii) a backscattered Alfvénic mode with lower amplitude and frequency. Malara et al. (1996) showed that in a compressible medium, the correlation between the velocity and the magnetic field fluctuations is reduced because of the generation of the backward propagating Alfvénic fluctuations, and of a compressive component of turbulence, characterized by density fluctuations and magnetic intensity fluctuations .

From a technical point of view it is worthwhile to remark that, when a large scale field which varies on a narrow region is introduced (typically a -like field), periodic boundaries conditions should be used with some care. Roberts et al. (1991, 1992) used a double shear layer, while Malara et al. (1992) introduced an interesting numerical technique based on both the glue between two simulation boxes and a Chebyshev expansion, to maintain a single shear layer, say non periodic boundary conditions, and an increased resolution where the shear layer exists.

Grappin et al. (1992) observed that the solar wind expansion increases the lengths normal to the radial
direction, thus producing an effect similar to a kind of inverse energy cascade. This effect perhaps should be
able to compete with the turbulent cascade which transfers energy to small scales, thus stopping the
non-linear interactions. In absence of non-linear interactions, the natural tendency towards an increase of
is stopped. These inferences have been corroborated by further studies like those by Grappin and
Velli (1996) and Goldstein and Roberts (1999). A numerical model treating the evolution of
e^{+} and e^{–}, including parametric decay of e^{+}, was presented by Marsch and Tu (1993a). The
parametric decay source term was added in order to reproduce the decreasing cross-helicity
observed during the wind expansion. As a matter of fact, the cascade process, when spectral
equations for both e^{+} and e^{–} are included and solved self-consistently, can only steepen the
spectra at high frequency. Results from this model, shown in Figure 68, partially reproduce the
observed evolution of the normalized cross-helicity. While the radial evolution of e^{+} is correctly
reproduced, the behavior of e^{–} shows an over-production of inward modes between 0.6 and 0.8 AU
probably due to an overestimation of the strength of the pump-wave. However, the model is
applied to the situation observed by Helios at 0.3 AU where a rather flat e^{–} spectrum already
exists.

### 5.2 Local production of Alfvénic turbulence at high latitude

An interesting solution to the radial behavior of the minority modes might be represented by local generation mechanisms, like parametric decay (Malara et al., 2001a; Del Zanna et al., 2001), which might saturate and be inhibited beyond 2.5 AU.

Parametric instability has been studied in a variety of situations depending on the value of the plasma (among others Sagdeev and Galeev, 1969; Goldstein, 1978; Hoshino and Goldstein, 1989; Malara and Velli, 1996). Malara et al. (2000) and Del Zanna et al. (2001) recently studied the non-linear growth of parametric decay of a broadband Alfvén wave, and showed that the final state strongly depends on the value of the plasma (thermal to magnetic pressure ratio). For the instability completely destroys the initial Alfvénic correlation. For (a value close to solar wind conditions) the instability is not able to go beyond some limit in the disruption of the initial correlation between velocity and magnetic field fluctuations, and the final state is as observed in the solar wind (see Section 4.2).

These authors solved numerically the fully compressible, non-linear MHD equations in a one-dimensional
configuration using a pseudo-spectral numerical code. The simulation starts with a non-monochromatic,
large amplitude Alfvén wave polarized on the plane, propagating in a uniform background magnetic
field. Successively, the instability was triggered by adding some noise of the order 10^{–6} to the initial density
level.

During the first part of the evolution of the instability the amplitude of unstable modes is small and,
consequently, non-linear couplings are negligible. A subsequent exponential growth, predicted by the linear
theory, increases the level of both e^{–} and density compressive fluctuations. During the second part of the
development of the instability, non-linear couplings are not longer disregardable and their effect is firstly to
slow down the exponential growth of unstable modes and then to saturate the instability to a level that
depends on the value of the plasma .

Spectra of are shown in Figure 69 for different times during the development of the instability. At
the beginning the spectrum of the mother-wave is peaked at , and before the instability saturation
the back-scattered e^{–} and the density fluctuations are peaked at and ,
respectively. After saturation, as the run goes on, the spectrum of e^{–} approaches that of e^{+} towards a
common final state characterized by a Kolmogorov-like spectrum and e^{+} slightly larger than
e^{–}.

The behavior of outward and inward modes, density and magnetic magnitude variances and the normalized cross-helicity is summarized in the left column of Figure 70. The evolution of , when the instability reaches saturation, can be qualitatively compared with Ulysses observations (courtesy of B. Bavassano) in the right panel of the same figure, which shows a similar trend.

Obviously, making this comparison, one has to take into account that this model has strong limitations
like the presence of a peak in e^{+} not observed in real polar turbulence. Another limitation, partly due to
dissipation that has to be included in the model, is that the spectra obtained at the end of the
instability growth are steeper than those observed in the solar wind. Finally, a further limitation
is represented by the fact that this code is 1D. However, although for an incompressible 1-D
simulation we do not expect to have turbulence development, in this case, since parametric
decay is based on compressive phenomena, an energy transfer along the spectrum might be at
work.

In addition, Umeki and Terasawa (1992) studying the non-linear evolution of a large-amplitude incoherent Alfvén wave via 1D magnetohydrodynamic simulations, reported that while in a low beta plasma () the growth of backscattered Alfvén waves, which are opposite in helicity and propagation direction from the original Alfvén waves, could be clearly detected, in a high beta plasma () there was no production of backscattered Alfvén waves. Consequently, although numerical results obtained by Malara et al. (2001b) are very encouraging, the high beta plasma (), characteristic of fast polar wind at solar minimum, plays against a relevant role of parametric instability in developing solar wind turbulence as observed by Ulysses. However, these simulations do remain an important step forward towards the understanding of turbulent evolution in the polar wind until other mechanisms will be found to be active enough to justify the observations shown in Figure 63.