Symmetry breaking and universality of decaying MHD Taylor-Green flows
SSymmetry breaking and universality of decaying MHD Taylor-Green flows
V. Dallas ∗ and A. Alexakis Laboratoire de Physique Statistique, ´Ecole Normale Sup´erieure, Universit´e Pierre et Mari´e Curie,Universit´e Paris Diderot, CNRS, 24 rue Lhomond, 75005 Paris, France
We investigate the evolution and stability of a decaying magnetohydrodynamic (MHD) Taylor-Green flow. The chosen flow has been shown to result in a steep total energy spectrum with powerlaw behaviour k − . We investigate the symmetry breaking of this flow by exciting perturbations ofdifferent amplitudes. It is shown that for any finite amplitude perturbation there is a high enoughReynolds number for which the perturbation will grow enough at the peak of dissipation rate result-ing to a non-linear feedback in the flow and subsequently break the Taylor-Green symmetries. Inparticular, we show that symmetry breaking at large scales occurs if the amplitude of the perturba-tion is ρ crit ∼ Re − and at small scales occurs if ρ crit ∼ Re − / . This symmetry breaking modifiesthe scaling laws of the energy spectra at the peak of dissipation rate away from the k − scaling andtowards the classical k − / and k − / power laws. I. INTRODUCTION
In magnetohydrodynamic (MHD) turbulence severalphenomenological theories exist debating for the inter-pretation of the power law of the energy spectrum [1–5].In summary, the power law scaling exponents obtainedin these phenomenologies based on weak and strong tur-bulence arguments both for isotropic and anisotropic en-ergy spectra are − − / − /
2. Numerical simula-tions to date are unable to provide a definitive answer tothis scaling. For example, some direct numerical simu-lations (DNS) obtained energy spectra with k − / whileothers k − / scaling for freely decaying MHD turbulentflows [6, 7]. Astrophysical observations have shown thatthis difference in the power law scaling also exists for themeasured energy spectra of the solar wind [8]. In addi-tion, indications of k − scaling are reported for the mag-netic energy spectrum measured in the magnetosphere ofJupiter [9].Recently, large resolution simulations by Lee et al. [10]demonstrated k − , k − / and k − / total energy spec-trum scalings for different initial conditions of the mag-netic field. Thus, they showed dependence of the energyspectrum at the peak of dissipation on the initial con-ditions. Consequently, this suggests lack of universalityin decaying MHD turbulence. The difference between − / − / − k − spectrum. This scaling of the total energy spectrumwas demonstrated to originate from high shearing regionsthat manifest discontinuities in the magnetic field corre-sponding to strong current sheets [11].All the initial conditions in [10] were satisfying symme- ∗ [email protected] tries of the Taylor-Green (TG) vortex [12]. This propertywas taken into account by numerically enforcing thesesymmetries in order to achieve higher resolutions withless compational cost [10, 13]. The − ∼ O (100)), the TGvortex symmetries did not break within the time intervalof reaching the peak of dissipation. This suggests thatthe TG symmetries are a strong property of the evolutionequations preserved in time. However, Stawarz et al. [15]showed that the TG symmetries can be broken at verylong time scales using runs of low Reynolds numbers dueto round-off error accumulation.Preservation of the TG symmetries hinders the flowfrom exploring all phase space and concequently preventsit from reaching a universal behaviour. Moreover, thebreaking of the TG symmetries can possibly modify thescaling of the energy spectrum by the time of maximumdissipation rate t peak , where the largest inertial range isobtained. Thus, before claiming lack of universality ofspectral exponents for decaying MHD turbulence in pe-riodic boxes, the persistence of the TG symmetries within t peak is an important issue that needs to be resolved.We expect a critical perturbation amplitude to existso that the system transitions from symmetry preserva-tion to symmetry breaking within t peak . The dependenceof this critical amplitude on the Reynolds number andwhether the breaking of the TG symmetries leads to adifferent spectral exponent are the key open questionsthat we address in this work. In order to demonstratelack of universality at the peak of dissipation one needsto show if, at Re (cid:29)
1, there is a finite perturbation ampli-tude below which the power law of the spectrum remainsunchanged. Showing this way that the set of initial con-ditions which lead to a specific behaviour is of non-zeromeasure in the limit of Re → ∞ .In summary, given an infinitesimal perturbation, isthere a Reynolds number such that the symmetries break a r X i v : . [ phy s i c s . f l u - dyn ] O c t within t peak ? Will the breaking of the symmetries lead toa different power law spectrum? Are the discontinuities,which are responsible for the k − spectra, formed due toenforcement of the TG symmetries? Are there universal-ity classes for moderate Reynolds numbers or is there auniversal power law scaling for the high Reynolds num-ber limit? In this work, we investigate these questionsby considering a large set of numerical simulations.The paper is structured as follows. Section II describesthe numerical methodology to solve the governing equa-tions for our decaying MHD turbulent flows and sectionIII provides the necessary details with regards to theTaylor-Green vortex, its symmetries and the measuresof symmetry breaking. In section IV, we define our nu-merical parameters along with our perturbed initial con-ditions. First, we analyse the results from the growthof infinitesimal perturbations (see section V) and thenfrom the finite amplitude perturbations (see section VI)by applying the measures of symmetry breaking. Finally,in section VII we conclude by summarising our findings. II. DNS OF DECAYING MHD TURBULENCE
We consider the three-dimensional, incompressibleMHD equations of fluid velocity u and magnetic induc-tion b to be ∂ t u − ( u × ω ) = − ∇ P + ν ∆ u + ( j × b ) (1) ∂ t b + ( u · ∇ ) b = ( b · ∇ ) u + κ ∆ b (2) ∇ · u = ∇ · b = 0 (3)with ν the kinematic viscosity, κ the magnetic diffusivity, ω ≡ ∇ × u the vorticity, j ≡ ∇ × b the current densityand P = p/ρ + u the fluid pressure, composed by theplasma pressure p divided by ρ the constant mass den-sity plus the hydrodynamic pressure u . Note that themagnetic field has units of Alfv´en velocity, i.e. b / √ ρµ ,where µ = ( κσ ) − is the permeability of free space with σ the electrical conductivity. If ν = κ = 0, the totalenergy E t ≡ (cid:104)| u | + | b | (cid:105) = E u + E b , the magnetic he-licity H b ≡ (cid:104) u · b (cid:105) and the cross helicity H c ≡ (cid:104) a · b (cid:105) areconserved in time (the angle brackets (cid:104) . (cid:105) denote spatialaverages in this study). Here, a is the magnetic poten-tial, which is defined as a ≡ − ∇ − ( ∇ × b ), since one canset b ≡ ∇ × a with ∇ · a = 0.Our numerical method is pseudo-spectral [16], whereeach component of u and b is represented as truncatedGalerkin expansions in terms of the Fourier basis. Thenon-linear terms are initially computed in physical spaceand then transformed to spectral space using fast Fouriertransforms [17]. Aliasing errors are removed using the2/3 dealiasing rule, i.e. wavenumbers k ∈ [1 , N/ N is the number of grid points in each Cartesian coordi-nate of our box of period 2 π . The non-linear terms alongwith the pressure term are computed in such a way that u and b are projected on to a divergence-free space sothat Eqs. (3) are satisfied [18]. The temporal integration of Eqs. (1) and (2) is performed using a second-orderRunge-Kutta method. The code is parallelised using ahybrid parallelisation (MPI-OpenMP) scheme [19]. III. TAYLOR-GREEN VORTEX, SYMMETRIESAND MEASURES OF SYMMETRY BREAKING
The initial conditions that we choose to focus in thisstudy is a magnetic Taylor-Green flow, which results in k − spectra at the peak of dissipation [10, 11, 14, 20]. Inparticular, the initial velocity field is the Taylor-Greenvortex [21] defined as u T G ( x ) = u sin x cos y cos z − cos x sin y cos z , (4)and the initial magnetic field is given by b T G ( x ) = b cos x sin y sin z sin x cos y sin z − x sin y cos z (5)where b and u were chosen so that the norm of the twofields is unity, i.e. (cid:107) u T G (cid:107) = (cid:107) b T G (cid:107) = 1. Here (cid:107) . (cid:107) standsfor the L norm (cid:107) g (cid:107) = V (cid:82) V g · g d x , where g is anarbitrary vector field.Given these initial conditions and in the absence ofany noise the symmetries are preserved by the evolutionequations exactly [13, 14]. In particular, we have reflec-tion (anti)symmetries about the planes x = 0, x = π , y = 0, y = π , z = 0 and z = π as well as ro-tational (anti)symmetries of angle nπ about the axes( x, y, z ) = ( π , y, π ) and ( x, π , π ) and of angle nπ/ π , π , z ) for n ∈ Z . The above mentioned planesthat possess reflection symmetries form the insulatingfaces of the sub-boxes [0 , π ] [12], where the j T G is every-where parallel to these faces. Note that for these partic-ular initial conditions b T G satisfies the same symmetrieswith ω T G and u T G with j T G .It was shown in [11] that the k − spectrum observed inthe numerical simulations originates from the formationof strong current sheets at the reflection symmetry planes x = 0, x = π , y = 0 and y = π . So, we focus on only oneof these symmetries. In particular, we will investigate thereflection symmetry around the plane x = 0. We thendefine the reflection operator R x around the x = 0 planeas R x g x ( x, y, z ) g y ( x, y, z ) g z ( x, y, z ) = − g x ( − x, y, z ) g y ( − x, y, z ) g z ( − x, y, z ) . (6)The TG initial conditions under the action of R x trans-form as follows R x u T G = u T G and R x b T G = − b T G . (7)Given any arbitrary set of fields u , b we define u s and b s as u s = 12 ( u + R x u ) and b s = 12 ( b − R x b ) (8)with u s and b s transforming similar to the TG initialconditions under reflection R x (see Eq. (7)). Similarlywe define u a and b a u a = 12 ( u − R x u ) and b a = 12 ( b + R x b ) (9)as the part of the flow that does not follow the TG sym-metries. Note that u a and b a transform differently underreflection, i.e. R x u a = − u a and R x b a = b a . (10)We will refer to u s , b s as the symmetric part of the flowwhile to u a , b a as the asymmetric part of the flow. Notethat if we start with u = u T G and b = b T G at t = 0,then u a , b a will remain zero throughout the computation.Thus, u a , b a can provide us with a measure of the extentthe symmetries are broken. Here we will focus on twosuch measures. First we consider the ratio of the energiesof asymmetric to the symmetric component of the fields E a /E s , where E a = 12 ( (cid:107) u a (cid:107) + (cid:107) b a (cid:107) ) (11)and E s = 12 ( (cid:107) u s (cid:107) + (cid:107) b s (cid:107) ) . (12)This quantity provides a measure of the degree the TGsymmetries are broken in the large (energy containing)scales. We also focus on the small scales by looking atthe ratio of the dissipation rates (cid:15) a /(cid:15) s , where (cid:15) a = ν (cid:107) ∇ × u a (cid:107) + κ (cid:107) ∇ × b a (cid:107) (13)and (cid:15) s = ν (cid:107) ∇ × u s (cid:107) + κ (cid:107) ∇ × b s (cid:107) . (14) IV. INITIAL CONDITIONS AND SIMULATIONPARAMETERS
To study the stability of the TG symmetries and theirimplications on the energy spectrum a series of numericalsimulations were performed. The simulations were car-ried out on a triple periodic box of size 2 π . The initialconditions were composed by the TG initial conditionsplus small perturbation fields √ ρ u p , √ ρ b p , viz. u = u T G + √ ρ u p and b = b T G + √ ρ b p . (15)The perturbation fields u p , b p were chosen to be a super-position of Fourier modes in spherical shells 2 ≤ | k | ≤ k p . The phases of the Fourier modes were chosen so that R x u p = − u p , R x b p = b p and random otherwise. Inthis way we guarantee that the two perturbation fieldsgive no contribution to u s and b s . The norm of the twofields was set to unity (cid:107) u p (cid:107) = (cid:107) b p (cid:107) = 1 so that theamplitude ρ at t = 0 expresses the ratio of the kineticenergy of the perturbation field to the energy of the TGflow (i.e. ρ ≡ E a | t =0 /E s | t =0 ). Additionally, we define ρ (cid:15) ≡ (cid:15) a | t =0 /(cid:15) s | t =0 as the ratio of the dissipation rate ofthe asymmetric part of the flow to the symmetric part ofthe flow at t = 0.The Reynolds number is defined based on the ve-locity rms value at t = 0 and smallest wavenumber k T G = 1 in the box, i.e. Re ≡ (cid:107) u T G (cid:107) /νk
T G . Withthese scales we can also define the eddy turnover time τ L ≡ ( u T G k T G ) − = 1 at t = 0. The smallest lengthscale in our flows is defined based on Kolmogorov scaling η ≡ ( ν /(cid:15) t ) / , where (cid:15) t = ν (cid:104)| ω | (cid:105) + κ (cid:104)| j | (cid:105) is the totaldissipation rate of energy. In all runs ν = κ and thus thePrandtl number is always unity. The set of parametersfor all the examined runs is given in Table I.TABLE I: Numerical parameters of the DNS. For allruns ν = κ . Note that k max = N/
3, using the 2 / k max η are reported atthe peak of (cid:15) t . ρ = 10 − , k p = 10 ρ = 0 . k p = 4 ρ = 0 . k p = 4 ρ (cid:15) = 1 . · − ρ (cid:15) = 0 . ρ (cid:15) = 0 . Re N k max η Re N k max η Re N k max η
50 128 2.80 50 128 2.80 50 128 2.78100 128 1.69 100 128 1.69 100 128 1.67200 256 2.09 200 256 2.09 200 256 2.04300 256 1.48 300 256 1.47 300 256 1.43500 512 2.25 500 512 2.24 500 512 2.151000 512 1.42 1000 512 1.42 1000 512 1.342000 1024 1.81 2000 1024 1.79 2000 1024 1.665000 2048 1.96
V. GROWTH OF INFINITESIMALPERTURBATIONS
As a first step we look at the temporal evolution offlows with energy ratio ρ = 10 − and dissipation ratio ρ (cid:15) = 1 . · − at t = 0. For this choice and for allReynolds numbers considered here the amplitude of theperturbation (symmetry breaking part of the flow) re-mains much smaller than the symmetric part of the flowat all times of interest. Thus, there is negligible effect ofthe perturbation on the part of the flow that obeys theTG symmetries and u a , b a evolve passively following theMHD equations (1) & (2) linearised around u s , b s .Figure 1 shows the temporal evolution of E s in blue(dark grey) and E a in red (light grey) for the seven dif-ferent Reynolds numbers examined. The lower curves arethe small Re cases while the top curves are the high Re cases. The vertical dashed line indicates t peak , the timethat (cid:15) t is peaked which is the time we are interested in.In this case, very weak variations of t peak were observedwith Re . It is evident that as the Reynolds number is in-creased the growth of the asymmetric part of the energy E a is increased.FIG. 1: (Color online) Evolution of E s and E a as afunction of time for different Reynolds numbers. Thevertical dashed line indicates the time of maximumtotal dissipation rate.Since we are interested in the symmetry breaking atthe peak of dissipation, we plot the ratio E a /E s at t peak as a function of the Reynolds number (see Fig. 2).This energy ratio appears to increase linearly with theFIG. 2: Energy ratio E p /E s at the time of maximumdissipation rate as function of the Reynolds number for ρ = 10 − .Reynolds number. This linear increase of the pertur-bation energy can be understood if we consider that themain source of growth of the perturbation comes from themagnetic shear layer at the x = 0 plane whose strengthincreases with Re. The time scale for the growth of aperturbation in such layers is controlled by the shear rate B /δ where B is the amplitude of the magnetic field inthe layer and δ is the thickness of the shear layer. B has negligible dependence on Re and is determined bythe initial conditions. The thickness of the reconnectionlayer δ is expected to scale like δ ∼ L/ √ S L where L isthe length of the layer that is of the order of the box sizeand S L the Lundquist number S L = B L/κ [22]. For thisproblem S L ∼ Re since (cid:107) u (cid:107) ∼ (cid:107) b (cid:107) and ν = κ . At theshort time scale t peak we expect that transient growthrates will dominate and the growth of the perturbationin time will be linear rather than exponential. Therefore,we expect that the amplitude of the perturbation A p at t peak will increase from the initial value A as: A p ∼ A B δ t peak ∼ A B t peak L Re / (16)from which we conclude that E a E s ∼ A p B ∼ A (cid:18) t peak L (cid:19) Re (17)and hence the linear increase observed in Fig. 2.Figure 3 presents the time evolution of the two dis-sipations (cid:15) s and (cid:15) a divided by the viscosity ν . WhileFIG. 3: (Color online) Evolution of (cid:15) s /ν and (cid:15) a /ν as afunction of time for different Reynolds numbers. Thevertical dashed line indicates the time of maximumtotal dissipation rate.both grow with time, the asymmetric dissipation (cid:15) a /ν increases by roughly four orders of magnitude in twoturnover times for the highest Reynolds number exam-ined.In Fig. 4 we plot the ratio of the two dissipations (cid:15) a /(cid:15) s at t peak . This ratio is increasing faster than linear withReynolds number, indicating that symmetries break evenfaster in the small scales. The scaling observed is closeto (cid:15) a /(cid:15) s ∼ Re / . (18)The fast breaking of the symmetries in the small scalescan also be seen by looking at the energy spectra of theFIG. 4: Dissipation ratio (cid:15) p /(cid:15) s at the time of maximumdissipation rate as function of the Reynolds number for ρ = 10 − .fields u s , u a , b s and b a . Figure 5a and 5b show the sym-metric and asymmetric part of the magnetic ( E b,s , E b,a )and the kinetic ( E u,s , E u,a ) energy spectra, respectively,compensated by k for different times up to t peak . Red(dark grey) curves represent the energy spectra of theasymmetric part of the flow, while blue curves (light grey)the energy spectra of the symmetric part of the flow. Theblack line represents the spectrum of the total field, i.e. E b = E b,s + E b,a and E u = E u,s + E u,a , at the time of themaximum dissipation rate. It is important to note theasymmetric part of the flow is always smaller than thesymmetric part for all scales. Thus, the symmetric partof the flow reproduces the k − energy spectrum, whilethe asymmetric part of the spectrum evolves passively.The results of this section show that no matter howsmall the amplitude ρ of the perturbation added in theTG initial conditions there is a Re for which the per-turbation will grow significantly enough for it to play a(non-linear) dynamical role in the system. This criti-cal amplitude can be estimated from our runs to be ei-ther ρ crit = C Re − if energy estimates are considered or ρ crit = C Re − / if dissipation estimates are considered,where C and C are constants. Consequently, the resultsexplain why symmetry breaking was not observed at t peak in the simulations of [10, 11, 14, 20] due to the presenceof numerical noise. Simulations using single precision ac-curacy introduce perturbations of amplitude ρ ∼ − and thus Re ∼ would be required for the symme-tries to break in the large scales and Re ∼ for thesymmetries to break in the small scales. Simulations atsuch Reynolds numbers cannot be performed on today’slargest supercomputers even at single precision accuracy.Therefore, this would make the observation of symmetrybreaking by numerical noise alone impossible at t peak .We note, however, that symmetry breaking can be ob-served at later times as it was shown in [15]. (a)(b) FIG. 5: (Color online) (a) Magnetic energy and (b)kinetic energy spectra at different times compensatedby k for ρ = 10 − and Re = 2000. The straight linesindicate the proposed spectral slopes k − , k − / , k − / . VI. GROWTH OF FINITE AMPLITUDEPERTURBATIONS
The growth of infinitesimal perturbations gives us alot of information on the growth of symmetry breakingperturbations. However, it does not provide us with anyinformation about a possible change in the spectral expo-nent and the return or not to a universal behaviour. Forthis reason we have performed two series of simulationswith perturbation amplitude ρ = 0 .
01 and ρ = 0 . ρ the perturbation grows suf-ficiently large at t peak to play a dynamical role in the evo-lution of the flow. The wavenumber range of the initialconditions of u p and b p was limited within 2 ≤ | k | ≤ ρ and the dissipation ratio ρ (cid:15) at t = 0are ρ (cid:15) (cid:39) .
54 for the ρ = 0 . ρ (cid:15) (cid:39) .
054 for the ρ = 0 .
01 runs.
A. Temporal behavior
Figure 6a shows the time evolution of Ohmic dissi-pation rate (cid:15) b for the runs with Re = 2000 and threedifferent values of ρ . One can notice that the pertur-bation of amplitude ρ = 0 . (cid:15) b , which has increased in ampli-tude and its peak has been shifted later in time (i.e. t peak /τ L (cid:39) . ρ = 0 .
01 case with a small increase of (cid:15) b ( ∼ t peak /τ L (cid:39) . ρ = 10 − case. Moreover, a new local peak starts to formaround t/τ L (cid:39) . ρ = 10 − run. This new peak, as we will show later, is due to theformation of smaller scales by the breaking of the sym-metries. (a)(b) FIG. 6: (a) Ohmic and (b) viscous dissipation rate as afunction of time for Re = 2000 and three differentvalues of ρ .In Fig. 6b we observe that the viscous dissipation rate (cid:15) u peaks at later times than (cid:15) b at this Re . For the runwith perturbation amplitude ρ = 0 . ρ = 0 .
01 there is a 10% increase with reference to the ρ = 10 − case. It is worth noting that the peak of (cid:15) u forthe ρ = 0 .
01 case coincides with the second local peak of (cid:15) b that takes place at t/τ L (cid:39) . Re and ρ = 0 .
01. For small (a)(b)
FIG. 7: (a) Ohmic and (b) viscous dissipation rate as afunction of time for ρ = 0 .
01 and different values of Re .values of the Reynolds number a single peak appears forthe time evolution of (cid:15) b at t/τ L (cid:39) . t/τ L (cid:39) .
4. The second peak canbe seen clearly only for the Re = 5000 run and it isvery close in amplitude with the first peak at t/τ L (cid:39) . Re , while its effect is muffled atsmaller Re . On the other hand, a single peak is devel-oped in the evolution of the viscous dissipation rate (cid:15) u ,which occurs at t/τ L (cid:39) . Re ≤ t/τ L (cid:39) . Re = 2000 and5000, respectively.The inset in Fig. 7a illustrates the values of (cid:15) b at themoment of the first peak t/τ L = 1 . t/τ L = 2 . (cid:15) b seems to reacheach asymptotic state much faster than the first peak as Re increases with Re = 5000 the transitional point where (cid:15) b | t/τ L =1 . (cid:39) (cid:15) b | t/τ L =2 . . Then, for Re (cid:29) (cid:15) b suggest that the second peak will be-come dominant and Reynolds number independent. Themaximum values of (cid:15) u that are plotted in the inset of Fig.7b as a function of Re also indicate that the viscous dissi-pation rate is far from reaching its asymptotic state evenfor our highest resolution simulations ( Re = 5000), whichare at the limit of the current computational power. It isinteresting that the Ohmic and viscous dissipation obeydifferent high Reynolds number asymptotics. In otherwords, the small scales of the magnetic field seem to reachits universal regime at lower Re than the small scales ofthe velocity field. B. Symmetry breaking and structures
In Figs. 8a and 8b we show the values of the energyratio E a /E s and the dissipation ratio (cid:15) a /(cid:15) s , respectively,at the peak of the energy dissipation rate for all the dif-ferent Re that we consider in this study for runs with ρ = 0 .
01 and 0.1 (see Table I). The effect of non-linearityis evident, since both cases deviate from the scalings E a /E s ∼ Re and (cid:15) a /(cid:15) s ∼ Re / observed in section V.In particular, at high Re the asymmetric part of the en-ergy for the runs with ρ = 0 . E a (cid:39) . E s while in the ρ = 0 .
01 it is significantlysmaller, i.e. E a (cid:39) . E s , with slightly higher value atthe second dissipation peak (see Fig. 8a). This impliesthat at the large, energy containing scales only a mod-est breaking of the symmetries has occurred for ρ = 0 . ρ . For the ρ = 0 . Re = 500 and for ρ = 0 .
01 we have (cid:15) a (cid:39) . (cid:15) s at the first dissipation peak and (cid:15) a (cid:39) . (cid:15) s at the seconddissipation peak for Re = 5000.Symmetry breaking of the resulting structures can bealso realised through visualisations. While three dimen-sional images of the full computational box provide globalinformation, they are sensitive in the choice of iso-contourlevels and very often can be misleading. We have thuschosen to show colour plots of two dimensional slices thatpass through the high current and vorticity density re-gions.Figure 9 illustrates the current (top panels) and vortic-ity (bottom panels) density at the z = π/ ρ = 10 − and Re = 2000 runwhere the perturbation evolved passively. In this case,the strong current and vorticity sheets that appear at the x = 0 , π/ y = 0 , π/ (a)(b) FIG. 8: (a) Energy ratio E a /E s and (b) dissipationratio (cid:15) a /(cid:15) s as functions of Re at the time of maximumdissipation rate for ρ = 0 .
01 and 0.1. The trianglescorrespond to the run with ρ = 0 .
01 for the second peakof dissipation that only appeared at the high Re casesat t/τ L (cid:39) . k − energy spectrum (see [11]). Therefore, the stabilityof these structures is crucial to determine the presenceor absence of universality. The effect of the perturbationon the structures becomes pronounced for the ρ = 10 − and Re = 2000 case (see Fig. 9b). The development ofthe perturbation has lead to the bending and curling ofthe current and vorticity sheets. The basic structures,however, remain unaltered without development of addi-tional features. On the other hand, at higher Reynoldsnumber (i.e. Re = 5000 and ρ = 10 − ) more structuresappear (see Fig. 9c). While the bended current and vor-ticity sheets are still present, the flow at their edge hasbeen fully “randomised” by the development of the in-stability. At this location, structures with no particularorder and reminiscent of flows with random initial condi-tions have formed. Note that these “random” structureshave generated scales smaller than the thickness of the“ordered” current/vorticity sheets. These structures are ( a )( b )( c )( d ) F I G . : ( C o l o r o n li n e ) C u rr e n t(t o pp a n e l s ) a nd v o r t i c i t y ( b o tt o m p a n e l s ) d e n s i t y a tt h e t i m e o f m a x i m u m d i ss i p a t i o n r a t e f o r ( a ) ρ = − , R e = , ( b ) ρ = . , R e = , ( c ) ρ = . , R e = nd ( d ) ρ = . , R e = . possibly responsible for the second peak of the Ohmicdissipation observed in Fig. 7a but also for the peakof the viscous dissipation at later times as Re increases(see Fig. 7b). This observation supports our conjecturethat at Re (cid:29) (cid:15) b will become a globalmaximum. Finally, in Fig. 9d the effect of the pertur-bation is much more pronounced for the ρ = 0 . Re = 2000 run; not only at the edge of the current andvorticity sheets but also away from the symmetry planes,where strong current and vorticity turbulent structureshave emerged. This probably indicates that the ampli-tude of the perturbation ρ = 0 . C. Spectral behaviour
As we stated in the introduction, it is an open questionwhether the breaking of the symmetries will change thepower law exponent of the spectra. Therefore, in thissection we investigate the effect of symmetry breakingon the inertial range scaling of our energy spectra.In Fig. 10 we show the energy spectra for the ρ = 0 . Re = 2000 case. Figure 10a represents the magneticenergy spectrum at the peak of the Ohmic dissipationand Fig. 10b the kinetic energy spectrum at the peakof viscous dissipation. The blue (light-grey) and the red(dark-grey) lines in these figures indicate the symmetricand asymmetric parts of the energy spectra, respectively,whereas the black lines indicate the full spectra, i.e. sym-metric plus asymmetric part. At the peak of Ohmic dissi-pation E b,s and E b,a reach equipartition within the range20 (cid:46) k ≤ k max . The full magnetic energy spectrum com-pensated by k has clearly a positive slope and a linearfit indicates a value close to the k − / scaling. For the ki-netic energy spectrum, equipartition occurs between E u,s and E u,a at all scales. The slope of the full compensatedspectrum k E u is positive and also close to k − / . How-ever, the range of wave numbers that exhibit a powerlaw for the velocity field is shorter than for the magneticfield and a distinction between k − / and k − / is notpossible. Hence, it is clear that for ρ = 0 . k − scaling and it returned to theclassical k − / (or k − / ) turbulence scaling. Note, how-ever, that this case was strongly perturbed at t = 0 witha significant amount of enstrophy and current densityintroduced by the perturbation.The case with ρ = 0 .
01 and Re = 5000 is more in-sightful since both the energy and the enstrophy/currentdensity of the perturbation are significantly smaller thanthe TG initial conditions. The spectra for this case areshown for the two total dissipation peaks in Figs. 11 and12 at times t/τ L = 1 . t/τ L = 2 .
4, respectively. Bothmagnetic and kinetic energy spectra reach equipartition (a)(b)
FIG. 10: (Color online) (a) Magnetic energy and (b)kinetic energy spectra at the peak of Ohmic and viscousdissipation rate, respectively, for ρ = 0 . Re = 2000.between the symmetric and the asymmetric part of theflow at small scales but not at large scales. In particu-lar, at the first total dissipation peak E b,s ∼ E b,a and E u,s ∼ E u,a only for wavenumbers k >
30 (see Fig. 11).The slopes of the compensated full energy spectra k E b and k E u are positive but less than the k − / scaling.Note that for the symmetric part of the magnetic field k E b,s ∼ const . Therefore, k − is still a good scaling for E b,s implying that the change in the slope of E b is dueto the symmetry breaking part of the flow.In the second dissipation peak more scales havereached equipartition between the symmetric and theasymmetric part of the flow (see Fig. 12). The power lawof the full magnetic energy spectrum remains between k − and k − / while the full kinetic energy spectrum iscloser to k − / .The evolution of the scaling exponents for the magneticenergy ( E b ∼ k s b ) and kinetic energy spectra ( E u ∼ k s u )for various cases of Table I are presented in Figs. 13and 14. The scaling exponents were obtained using alinear fit on the energy spectra between wavenumbers0 (a)(b) FIG. 11: (Color online) (a) Magnetic energy and (b)kinetic energy spectra for ρ = 0 .
01 and Re = 5000 atthe first peak of total dissipation rate t/τ L = 1 . < k <
20 for the runs with Re = 1000, 4 < k < Re = 2000 and 4 < k <
40 for the runwith Re = 5000. We note that measured exponents inthis way are sensitive in the choice of the fitting rangeespecially away from t peak . However, the objective hereis to show the time evolution of the exponents and notthe precise value. Our choices were based on Figs. 10to 12 as the most reasonable to our opinion. The effectof the perturbation amplitude on the spectral exponents s b and s u is shown in Fig. 13. The strongly perturbedcase with ρ = 0 . − / s b while no clear saturation can be observedfor s u . This indicates that Reynolds number is not highenough for a clear scaling of the kinetic energy. Evenin the unperturbed case E u does not have a clear powerlaw scaling (see [11]). The weakly perturbed cases forthe magnetic field reach a value close to s b = − s u is close to − (a)(b) FIG. 12: (Color online) (a) Magnetic energy and (b)kinetic energy spectra for ρ = 0 .
01 and Re = 5000 atthe second peak of total dissipation rate t/τ L = 2 . ρ = 0 .
01 case where the high-est Re was obtained. Although s u is close to − t/τ L (cid:39) . − / Re increasesat the peak of the viscous dissipation (see Fig. 14b).The scaling exponent of the magnetic energy spectrumsaturates to a value between − − /
3. It is worthnoting that although the slope at the first peak of Ohmicdissipation ( t/τ L (cid:39) .
7) seems to have saturated as afunction of Re at the second peak of Ohmic dissipation( t/τ L (cid:39) .
4) the slope appears to still increase with Re .We conjecture that this intermediate value of the ex-ponent, i.e. − < s b < − /
3, is a finite Reynolds numbereffect. We expect that as Re increases the small scalesthat break the symmetries will start forming a strongturbulence scaling (i.e. k − / or k − / ). The spectrumin these scales is related to the turbulent fluctuations ob-served in Fig. 9. As Re will increase, more small scaleturbulent fluctuations will be excited and consequentlythe range of validity of this power law scaling will in-crease. Ultimately, this will be the dominant spectrum1 (a)(b) FIG. 13: Scaling exponent of (a) the magnetic energyand (b) the kinetic energy spectrum for different ρ and Re = 2000.as Re → ∞ . This is what we try to depict in Fig. 15 bypresenting schematically the expected energy spectra at Re (cid:29)
1. For the large scales where the symmetries arenot broken, we expect the k − energy spectrum, whichreflects the presence of the strong current/vortex sheets,to persist. The wavenumber k s , that depends on theamplitude of the perturbation, determines the transitionpoint between the k − and the k − / spectrum. If thesymmetries are broken at all scales at Re (cid:29)
1, then the − VII. CONCLUSIONS
In this work we have studied one of the proposed ini-tial conditions in [10] for a freely decaying MHD flow. (a)(b)
FIG. 14: Scaling exponent of (a) the magnetic energyand (b) the kinetic energy spectrum for different Re and ρ = 0 . k − energy spectrum [11] different thanthe more commonly obtained spectra k − / and k − / forrandom initial conditions. Here, we investigated whetherthis behaviour persist when the initial conditions weaklydeviate from the ones proposed in [10] and break the in-volved Taylor-Green symmetries by adding a small per-turbation.We demonstrated that a sufficiently small perturbationevolves passively and it grows at a rate that increaseswith Reynolds number. In particular, it was shown thatthe energy ratio scales as E a /E s ∼ Re the dissipationratio scales as (cid:15) a /(cid:15) s ∼ Re / at the peak of the total dis-sipation rate. Therefore, for any finite amplitude pertur-bation, no matter how small it is, there is a high enoughReynolds number for which the perturbation will growenough at the peak of the total dissipation resulting toa non-linear feedback in the flow and subsequently breakthe TG symmetries.For strong perturbations of amplitude ρ = 0 . , π ] boxes. These new smallscale features change the slope of the energy spectrumfrom k − to the classical turbulence spectrum, i.e k − / or k − / power law scaling.For the smaller amplitude perturbation ρ = 0 .
01 theinitially passive asymmetric part of the flow grows to anamplitude that can play a non-linear role in the MHDequations at large Re . A new dissipation peak appears asa result of the non-linear evolution of the instability. Thestrong shearing regions bend and turbulent structuresappear at the edge of the current sheet causing this newpeak. The scaling exponent of the energy spectrum isclearly larger than − − /
3. Thisintermediate value of the exponent appears because atthe examined Reynolds numbers the small scales havebroken the symmetries and are approaching the strongturbulence scaling, while the large scales still exhibit the k − scaling. Therefore, the measured exponent appearsto take an intermediate value. We argue that the strongturbulence scaling (i.e. k − / or k − / ) will dominate athigher Re .The above results suggest that unless the TG symme-tries are satisfied exactly in periodic boxes they will breakat sufficiently large Re and the strong turbulence scalingof the spectrum will be recovered. Thus, the system willreturn to a universal behaviour up to the distinction be-tween the k − / and the k − / that we could not resolvehere.Another very important issue is why and under whatconditions these large current sheets, which lead to thetransient k − energy spectra form. Even though currentsheets form spontaneously in MHD [22], in this case thespanwise length was the size of the fundamental box andtheir amplitude was strong enough to dominate the en-ergy spectrum something that is not typically observedin random MHD turbulent flows. It is crucial in obser-vations to distinguish the k − spectra that manifest dueto discontinuities in the magnetic field and those due toweak turbulence. Note that these mechanisms are dis-tinctly different even though the energy spectra displaythe same scaling. Some of these questions are going tobe addressed in our future work. ACKNOWLEDGMENTS
V.D. acknowledges the financial support from EU-funded Marie Curie Actions–Intra-European Fellowships(FP7-PEOPLE-2011-IEF, MHDTURB, Project No.299973). The computations were performed using theHPC resources from GENCI-TGCC-CURIE (Project No.x2013056421) and PRACE-FZJ-JUQUEEN (Projectname PRA068). [1] D. Biskamp,
Magnetohydrodynamic turbulence (Cam-bridge University Press, 2003).[2] Y. Zhou, W. H. Matthaeus, and P. Dmitruk, Rev. Mod.Phys. , 1015 (2004).[3] S. Boldyrev, Phys. Rev. Lett. , 115002 (2006).[4] C. S. Ng and A. Bhattacharjee, Phys. Plasmas , 605(1997).[5] S. Galtier, S. V. Nazarenko, A. C. Newell, and A. Pou-quet, J. Plasma Phys. , 447 (2000).[6] W. C. M¨uller and R. Grappin, Phys. Rev. Lett. ,114502 (2005).[7] P. D. Mininni and A. Pouquet, Phys. Rev. Lett. ,254502 (2007).[8] J. J. Podesta, D. A. Roberts, and M. L. Goldstein, TheAstrophysical Journal , 543 (2007).[9] J. Saur, H. Politano, A. Pouquet, and W. H. Matthaeus,Astron. Astrophys. , 699 (2002).[10] E. Lee, M. E. Brachet, A. Pouquet, P. D. Mininni, and D. Rosenberg, Phys. Rev. E , 016318 (2010).[11] V. Dallas and A. Alexakis, Phys. Rev. E (to be pub-lished), e-print arXiv:1306.1380 (2013).[12] M. E. Brachet, D. I. Meiron, S. A. Orszag, B. G. Nickel,R. H. Morf, and U. Frisch, J. Fluid Mech. , 411(1983).[13] E. Lee, M. E. Brachet, A. Pouquet, P. D. Mininni, andD. Rosenberg, Phys. Rev. E , 066401 (2008).[14] V. Dallas and A. Alexakis, Phys. Fluids (accepted), e-print arXiv:1304.0695 (2013).[15] J. E. Stawarz, A. Pouquet, and M.-E. Brachet, Phys.Rev. E , 036307 (2012).[16] D. Gottlieb and S. Orszag, Numerical analysis of spectralmethods: theory and applications , Vol. 26 (SIAM, 1977).[17] M. Frigo and S. Johnson, in