Wavelet-based coherent vorticity sheet and current sheet extraction from three-dimensional homogeneous magnetohydrodynamic turbulence

from three-dimensional homogeneous magnetohydrodynamic turbulence Katsunori Yoshimatsu, Yuji Kondo, Kai Schneider, Naoya Okamoto, Hiroyuki Hagiwara, and Marie Farge Department of Computational Science and Engineering, Nagoya University, Nagoya 464-8603, Japan M2P2–CNRS and CMI, Université de Provence, 39 rue Frédéric Joliot-Curie, 13453 Marseille Cedex 13, France Center for Computational Science, Graduate School of Engineering, Nagoya University, Nagoya 464-8603, Japan LMD-IPSL-CNRS, Ecole Normale Supérieure, 24 rue Lhomond, 75231 Paris Cedex 05, France


I. INTRODUCTION
7][8][9] The extraction and characterization of these organized structures are a long lasting problem.A sound understanding of their dynamics is of fundamental importance and is a prerequisite for accurate and physically based modeling.
Intermittency in MHD turbulence can be attributed to these organized sheetlike structures, noting that flow intermittency is attributed to vortex tubes in HD turbulence as first suggested by Batchelor and Townsend. 10Intermittency can be defined as the lacunarity of such fine-scale activity, in other words, the spatial support decreases in the active regions with scale.A suitable multiscale representation of intermittent fields is beneficial, because it allows us to take into account the lacunarity of the small scale activity.
Wavelet analysis is a powerful tool, as it yields sparse representations of intermittent data due to the fact that it keeps track of both location and scale of the field.2][13][14] Recently, they have been applied to highresolution DNS data of three-dimensional homogeneous isotropic HD turbulence, [15][16][17] obtained through DNS performed on the Earth simulator at resolutions up to 2048 3 grid points. 18he wavelet transform, efficiently computed by the fast wavelet transform, decomposes a given field into scale-space contributions from which it can be perfectly reconstructed.The small-scale contributions have significant values only in well-localized active regions and nonsignificant values in nonactive regions.Thus the number of degrees of freedom of the field can be substantially reduced by keeping the significant contributions of active regions, while discarding the nonsignificant contributions.
The coherent vortex extraction ͑CVE͒ method introduced for two-and three-dimensional hydrodynamic turbulence [19][20][21] is one of the most useful tools for the extraction of coherent vortices from HD turbulence together with a significant reduction in the number of degrees of freedom.It is based on the orthogonal wavelet decomposition of vorticity and allows the wavelet coefficients to be split into two sets, i.e., weak and strong coefficients depending on their magnitude.A subsequent reconstruction in physical space of both coefficient sets yields a separation of the vorticity field into two orthogonal parts, namely, coherent and incoherent vorticity.Coherent vorticity is reconstructed from few wavelet coefficients, whose moduli are above a certain threshold, which is motivated from denoising theory 22 ͑for more details, we refer to Sec.II B͒, it retains the coherent vortex tubes and exhibits statistics similar to those of the total vorticity.The incoherent vorticity reconstructed from the remaining large majority of the wavelet coefficients corresponds to an almost uncorrelated random background flow with quasi-Gaussian one-point statistics.In a previous study, the CVE method was applied to DNS data of threedimensional homogeneous isotropic hydrodynamic turbulence at resolution 256 3 and a Taylor microscale Reynolds number of R = 150. 20It was shown that only 3% of the wavelet coefficients of vorticity represent the coherent vortex tubes, which retain 99% of the energy and 75% of the enstrophy of the flow.Furthermore, the percentage of wavelet coefficients retained by the coherent vorticity decreases with increasing R , as Okamoto et al. 15 have shown for R ranging from 167 to 732.Applications of the CVE method to DNS data of resistive drift-wave turbulence can be found in Ref. 23 and applications to experimental data in turbulent edge plasma in Ref. 24.
The proposed turbulence model called coherent vortex simulation ͑CVS͒ consists of applying the CVE decomposition to the Navier-Stokes equations. 19,25A posteriori tests of this model were performed for two-dimensional flows 26 and three-dimensional turbulent mixing layers. 27CVS is based on the deterministic computation of the coherent flow evolution, using an adaptive wavelet basis and either modeling or neglecting the influence of the incoherent background flow.
In the present paper, we generalize the CVE method to extract coherent vorticity sheets and current sheets out of three-dimensional incompressible homogeneous MHD turbulence.We call this generalized method the coherent vorticity sheet and current sheet extraction ͑CVCE͒ method.To the best of our knowledge, this is the first time that waveletbased extraction methods have been applied to threedimensional MHD turbulence.We study the properties of the extracted coherent vorticity sheets and current sheets using visualization of intense regions of vorticity and current density.The statistical behaviors of the total, coherent, and incoherent fields are analyzed and the dynamics of these fields is studied by considering the energy fluxes.The CVE method has great potential for application to other types of intermittent data, since it is based only on the flow intermittency.It is not clear, however, whether CVCE works as well for MHD turbulence as CVE does for HD turbulence, 15,20,21 because the sheetlike structures in MHD turbulence seem less localized in space than the tubelike structures in HD turbulence.
Homogeneous turbulence was chosen for this paper in order to demonstrate the efficiency of CVCE in the worst possible case where structures are spread all over space in contrast with inhomogeneous turbulence.Indeed, the wavelet representation in CVE is better suited for dealing with inho-mogeneous flows, such as turbulent mixing layers, than for dealing with homogeneous flows.
The outline of the present paper is as follows: in Sec.II, after a short review of orthonormal wavelet analysis, we introduce the CVCE method.In Sec.III, we describe DNS of three-dimensional incompressible homogeneous MHD turbulence without mean magnetic field.In Sec.IV, we apply the CVCE method to the DNS data.We analyze the total, coherent, and incoherent fields in physical space and investigate their statistical behaviors.We also study statistics in spectral space.Finally, conclusions are drawn in Sec.V. We present some perspectives on the coherent vorticity sheet and current sheet simulation method ͑CVCS͒ based on the orthogonal wavelet decomposition.Appendix A discusses the influence of divergence of the vector valued wavelet basis and Appendix B covers the dependence of the number of iterations in the CVCE method.

II. WAVELET-BASED COHERENT VORTICITY AND CURRENT SHEET EXTRACTION
In this section, we fix the notation for the orthogonal wavelet decomposition of a three-dimensional vector valued field and describe main ideas of coherent vorticity sheet and current sheet extraction.For definitions and details on the orthogonal wavelet transform and its extension to higher dimensions, we refer readers to, e.g., Refs.11 and 28.

A. Vector valued orthogonal wavelet decomposition
We consider a three-dimensional 2-periodic vector field v͑x͒ = ͑v 1 ͑x͒ , v 2 ͑x͒ , v 3 ͑x͒͒ with x = ͑x 1 , x 2 , x 3 ͒ ⍀ = ͓0,2͔ 3 ʚ R 3 and v ᐉ L 2 ͑⍀͒ sampled on N =2 3J equidistant grid points, J being the number of octaves in each space direction of the Cartesian coordinate.The three-dimensional orthonormal wavelet transform unfolds v into scale, positions, and seven directions using a function ͑x͒, which is called a three-dimensional mother wavelet and based on a tensor product construction, where ͑x͒ is the one-dimensional mother wavelet and ͑x͒ is the one-dimensional scaling function.A wavelet is well localized in space x R 3 , is oscillating, and is smooth.In other words, a wavelet exhibits a fast decay for ͉x͉ tending to infinity, has at least a vanishing mean value, or better the first m moments of vanish, and its Fourier transform F͓ ͔͑k͒ exhibits fast decay for wavenumbers ͉k͉ tending to infinity.The mother wavelet then generates a family of wavelets , ͑x͒ by dilation and translation, which yields an orthogonal basis of L 2 ͑R 3 ͒ and also of L 2 ͑⍀͒ through the application of a periodization technique. 28he vector field v, having a mean value v, can be decomposed into an orthogonal wavelet series: where the multi-index = ͑j , i 1 , i 2 , i 3 ͒ denotes the scale 2 −j and the position 2 −j i =2 −j ͑i 1 , i 2 , i 3 ͒ of the wavelets for each direction .The index set ⌳ is ⌳ = ͕͑,͉͒ = 1, ... ,7, j = 0, ... ,J − 1, i n = 0, ... ,2 j − 1, and n = 1,2,3͖.
The index set ⌳ can be seen as the octree representation of the orthogonal wavelet coefficients, a generalization of the dyadic tree from one dimension to three dimensions.Consequently, we have 7 ϫ 2 3j wavelet coefficients at scale 2 −j .Due to orthogonality, the wavelet coefficients are given by v ˜, = ͗v , , ͘, where ͗¯, ¯͘ denotes the L 2 -inner product defined by ͗f , g͘ = ͐ ⍀ f͑x͒g͑x͒dx.The coefficients measure fluctuations of v around scale 2 −j and around position 2 −j i in one of the seven directions .The N wavelet coefficients v ˜, are efficiently computed from the N grid point values of v using the fast wavelet transform, which has linear complexity. 28

B. Coherent vorticity sheet and current sheet extraction
The wavelet-based method called CVE was introduced to extract coherent vortices out of two-and threedimensional hydrodynamic turbulent flows. 19,20The CVE method is an efficient tool for extracting organized structures from turbulent flows and reducing the number of degrees of freedom.21]29 Here, we generalize the CVE algorithm and introduce a wavelet-based method for extracting coherent vorticity sheets and current sheets from three-dimensional homogeneous MHD turbulence.
In the CVCE method, we consider the vorticity field ͑=ٌ ϫ u͒ and the current density field j͑=ٌ ϫ b͒ rather than the velocity field u and the magnetic field b, because vorticity preserves the Galilean invariance, whereas velocity does not.Furthermore, as previous studies have shown, 2,5 the regions of intense vorticity overlap well with the regions of intense current density.The magnetic field b is normalized by ͑ 0 0 ͒ 1/2 , where 0 is the permeability of free space and 0 is the fluid density.We choose the orthogonal wavelet representation, where each wavelet coefficient is indexed by its location, scale, and direction, instead of the Fourier representation, where each Fourier coefficient is indexed only by its wavenumber.
As proposed in Ref. 20, we consider the following minimal but hopefully consensual statement that coherent structures are not noise and define coherent structures as what remain after denoising, because there has not been yet a universal definition for them.For the noise, the mathematical definition stating that a noise cannot be compressed in any functional basis is used.Here, as first guess, the simplest possible type of noise, namely, additive, Gaussian, and white ͑uncorrelated͒ noise, is considered.Now, we sketch the procedure of the CVCE algorithm.We apply a nonlinear thresholding method to the wavelet coefficients of the vorticity field and, separately, the current density field j.To achieve this, both fields, and j, are first decomposed into orthogonal wavelet series.Then we split each field into coherent and incoherent contributions in wavelet space by applying nonlinear thresholding.The threshold values for the wavelet coefficients of vorticity ˜, and the wavelet coefficients of current density j ˜, are defined as ˜͑͒ = ͑2 / 3 2 ln N͒ 1/2 and ˜͑j͒ = ͑2 / 3 j 2 ln N͒ 1/2 , respectively, where 2 and j 2 are the variances of the incoherent contributions of ˜, and j ˜, , respectively ͑both of which are a priori unknown͒.Thus the value of ˜͑͒ can be different from that of ˜͑j͒.Both thresholding values depend only on the number of grid points N and on the variances.
The choice of these thresholding values is based on theorems of Donoho and Johnstone 22 proving optimality of the wavelet representation to denoise signals in the presence of Gaussian white noise with given variance.The waveletbased estimator for the denoised signal minimizes the maximal L 2 -error between the original signal and the denoised estimator for functions having a certain inhomogeneous regularity, which thus typically holds for intermittent data.
As the variances of the incoherent contributions are unknown in our case, we overestimate them as a first step by taking the variances of the total vorticity and current density instead.Therewith, for each field, we compute the threshold.Then we can compute the variance of those coefficients being smaller than this threshold, which yields a new improved value for the threshold.The coherent vorticity and coherent current density fields, c and j c , are reconstructed from the wavelet coefficients whose moduli are larger than ˜͑͒ and ˜͑j͒, respectively.The incoherent fields, i and j i , are reconstructed from the wavelet coefficients whose moduli are smaller or equal to ˜͑͒ and ˜͑j͒, respectively.The incoherent fields are also obtained by means of simple subtraction, i = − c and j i = j − j c , which is used in the present paper.The coherent and incoherent contributions for each field thus obtained are orthogonal, which ensures a separation of the total kinetic and magnetic enstrophies, defined as Note that the superscript u stands for the velocity field and b for the magnetic field, corresponding to the vorticity and current density fields, respectively.
The above procedure for determining the thresholding values can be iterated.This yields an improved estimation of the variances of the incoherent fields and thus allows us to compute subsequently improved thresholds.For the present paper, we performed one iteration to privilege good compression rates rather than perfectly denoised contributions.The influence of the number of iterations in the CVCE algorithm is briefly discussed in Appendix B. The mathematical properties for the CVE algorithm were investigated in Ref. 29.The influence of the number of iterations was studied for isotropic hydrodynamic turbulence in Refs.15  ϫ ٌ͑ −2 j͒, are used to reconstruct the coherent velocity field u c , the incoherent velocity field u i , the coherent magnetic field b c , and the incoherent magnetic field b i from the coherent and incoherent vorticities and current densities, respectively.
There is a large collection of possible orthogonal wavelets; the actual choice depends on which properties are preferred, e.g., compact-support, symmetry, smoothness, number of cancellations, computational efficiency. 28Here, we use the Coifman 12 wavelet, which is compactly supported, has four vanishing moments, is quasisymmetric and is defined with a filter of length 12.This filter length leads to a computational cost of the fast wavelet transform in 24N operations, since two filters are used.The dependence of the CVE method on the choice of the wavelet has been tested through the use of different wavelets, e.g., Coifman wavelets of different order, Symmlets, spline wavelets, and Haar wavelets.The results were robust, except in the case of the Haar wavelet, which only has one vanishing moment.A comparison between orthogonal and biorthogonal wavelet decomposition for the CVE applied to three-dimensional isotropic turbulence can be found in Ref. 30.

III. DIRECT NUMERICAL SIMULATION
We performed DNS of incompressible MHD turbulence in a 2 periodic box ⍀ obeying the following equations: where f u and f b are external forces, P is the pressure, is the kinetic viscosity, and is the magnetic diffusivity.The Prandtl number Pr is set to 1, i.e., = .The above equations are solved with a pseudospectral method using a fourth-order Runge-Kutta method for time integration.The aliasing errors are removed by means of a phase shift method.Only modes with wavenumbers satisfying k Ͻ 2 1/2 N 1/3 / 3 are retained, where k = ͉k͉ and N 1/3 is the number of grid points in each direction of the Cartesian coordinate.The wavenumber increment is 1 and the minimum wavenumber is 1.We use the Elsässer variables z Ϯ defined by z Ϯ = u Ϯ b to reduce the number of fast Fourier transforms for computing the nonlinear terms.Random solenoidal external forces F͓f u ͔͑k , t͒ and F͓f b ͔͑k , t͒ introduced by Yoshida and Arimitsu 31 are imposed on the velocity and magnetic fields in a low wavenumber range k Ͻ 2.5 in order to achieve the highest possible Reynolds number under the limitation of available computational resources.Here, F͓ • ͔ expresses the Fourier transform of •.The forcing satisfies the following conditions: where f ᐉ denotes the ᐉth component of f ͑ = u , b͒, k ᐉ is the ᐉth component of k, F is the intensity of f , T c is the correlation time of the external forces, ᐉmn is the third-order alternating tensor, ␦ ᐉm is the Kronecker delta, and ␦ k,k Ј = 1 for k = kЈ.Einstein's summation convention is used for the repeated alphabetical indices and ͗ • ͘ r denotes a random average.The forcing does not generate ensemble averages with cross helicity and magnetic helicity.Readers interested in details of how to generate such random forces are referred to the Appendix of Ref. 31.
The two simulations are performed until the energy dissipation rate per unit mass ͗⑀͘, which is given by ͗⑀͘ = ͗⑀ u ͘ + ͗⑀ b ͘, remains almost constant, which characterizes the behavior of small scales and thus corresponds to a quasistationary state ͑see Fig. 1͒.Here, ͗⑀ u ͘ and ͗⑀ b ͘ are the kinetic and magnetic energy dissipation rates, respectively.The basic parameters of the two DNSs are listed in Table I.
The final state of run 1 is used to initialize run 2. The initial velocity and magnetic fields of run 1 are generated randomly under the following constraints: Here, T i is the initial large eddy turnover time defined by T = L u / u 0 , and L u is the initial integral length scale given by L u = / ͑2u 0 2 ͒͐ 0 k max dkE u ͑k͒ / k.

TABLE I. DNS parameters of runs 1 and 2.
Here, k max is the maximum wavenumber of the retained modes and ⌬t is the time increment.Note that = .
where E u ͑k͒, E b ͑k͒, E C ͑k͒, and E H ͑k͒ are the kinetic energy spectrum, magnetic energy spectrum, cross helicity spectrum, and magnetic helicity spectrum, respectively.Here, we set k p = 2 and the constant C is determined so that the total kinetic and magnetic energy satisfy

IV. CVCE OF ISOTROPIC MHD TURBULENCE AT R M = 154
In this section, we apply the CVCE algorithm, described in Sec.II B, to the DNS data of run 2 ͑DNS 2 ͒ computed with N = 512 3 grid points and a Taylor microscale Reynolds number of R M = 154.

A. Total, coherent, and incoherent vorticities and current density fields
Figure 2 ͑top͒ shows isosurfaces of intense vorticity and current density regions of the total fields ͑green͒ in DNS 2 .We observe vorticity and current sheets as in previous DNS results ͑e.g., Refs. 2 and 4͒.Applying the CVCE algorithm, we then decompose both fields into their coherent and incoherent contributions.Note that the threshold value used in the algorithm for the vorticity field is different from that for the current density field.We show isosurfaces of the coherent vorticity and current density ͑red͒ in Fig. 2 ͑middle͒.The coherent vorticity and current density fields, plotted with the same isosurface values as those of the total fields, are in good agreement with their respective total fields.The coherent vorticity sheets and current sheets are represented by 3.21% of the wavelet coefficients of vorticity and by 3.16% of the wavelet coefficients of current density, respectively.Table III confirms that the coherent parts retain almost all of the kinetic and magnetic energies, 93.2% of the kinetic enstrophy and 93.7% of the magnetic enstrophy.
In contrast, the incoherent vorticity and current density ͑blue͒ is structureless, as shown in Fig. 2 ͑bottom͒; we checked this by zooming in that there is no vorticity and current sheets remain.Note that for the incoherent fields the values of the isosurfaces chosen for visualization are reduced by a factor of 3, as their fluctuations are much smaller than those of the total fields.These incoherent parts, reconstructed The values of the isosurfaces are taken as ͉͉ = ͗ t ͘ +4 for the total and coherent vorticity and ͉j͉ = ͗j t ͘ +4 j for the current density.For the incoherent vorticity and current density fields, the isosurfaces are set to ͉͉ = ͑͗ t ͘ +4 ͒ / 3 and ͉j͉ = ͑͗j t ͘ +4 j ͒ / 3, respectively.Here, ͗ t ͘ and ͗j t ͘ are the mean values of ͉͉ and ͉j͉ for the total vorticity and current density fields, respectively, and and j are standard deviation values of ͉͉ and ͉j͉, respectively.Cubes of size 512 3 are visualized.
TABLE III.Overall compression rates C u and C b , the percentages of the number of the wavelet coefficients retained by the coherent vorticity and current density fields over all scales, respectively, and percentages of kinetic energy, magnetic energy, kinematic enstrophy, and magnetic enstrophy kept by the coherent and incoherent fields.The CVCE method is applied to DNS 2 .from the remaining large majority of the wavelet coefficients of vorticity and current density, retain little of the kinetic and magnetic energies: only 6.8% of the kinetic enstrophy and 6.3% of the magnetic enstrophy ͑see Table III͒.

B. Probability density functions
The probability density functions ͑PDFs͒ of the velocity and magnetic fields of the total, coherent, and incoherent parts are plotted in Fig. 3. Comparison of the total and coherent velocity PDFs ͑two wide PDFs͒ shows that they coincide well.The incoherent velocity PDF is quasi-Gaussian with a strongly reduced variance compared to that of the total field.The same observations hold for the magnetic fields.
In contrast, the vorticity and current density PDFs exhibit different behavior ͑see Fig. 4͒.Although the PDFs of the total and coherent fields of both vorticity and current density almost coincide, they show stretched exponential tails that illustrate the intermittency that is due to the presence of coherent vorticity sheets and current sheets.The PDFs of the incoherent fields have exponential shapes with reduced variances compared to those of the total fields.The skewness and flatness of these total, coherent, and incoherent fields are summarized in Table IV.

C. Energy spectra
Figure 5 shows the kinetic and magnetic energy spectra of total, coherent, and incoherent parts of DNS 2 .The spectra are obtained by integrating energy in threedimensional k-space over spherical shells k = ͉k͉.We observe that the spectra of the total fields exhibit an Iroshnikov and   Kraichnan ͑IK͒ spectrum in the inertial subrange, i.e., E ͑k͒ ϰ k −3/2 ͑ = u , b͒. 32,33 This observation is in accordance with other recent results of forced DNS, 31 which are also confirmed in higher resolution freely decaying DNS. 34or the coherent contributions, we find that the kinetic and magnetic energy spectra are identical with those of the total fields all along the inertial range, respectively.This implies that coherent vorticity sheets and current sheets are responsible for the IK spectrum.For k IK տ 0.5, both spectra of the coherent fields differ from the spectra of the total fields, although the coherent fields still provide significant contributions at scales smaller than k IK Ӎ 0.5.Concerning the incoherent flow, we observe that the scaling of the incoherent kinetic energy spectrum is close to k 2 .This is also the case for the incoherent magnetic field.These k 2 spectra correspond to equipartitions of incoherent kinetic and magnetic energies between all wavenumbers k, respectively.The incoherent velocity and magnetic fields are therefore spatially decorrelated, which is consistent with the observation that incoherent vorticity and current density are structureless ͑Fig.2͒.We also note some energy piling up around the cutoff wavenumber in the total kinetic and magnetic energy spectra.These spikes are retained by the incoherent contributions but not by the coherent contributions.

D. Energy fluxes
Studying the transfer of kinetic and magnetic energy in Fourier space enables us to check the contributions of the coherent and incoherent fields to the nonlinear dynamics.Using the decompositions u = u c + u i and b = b c + b i and Elsässer variables defined by z Ϯ = u Ϯ b, we consider the following eight fluxes of the kinetic and magnetic energy for all possible combinations between coherent and incoherent contributions: where q = ͉q͉, u t = u, and b t = b.Thus, we can define eight fluxes, ⌸ ccc , ⌸ cci , ⌸ cic , ⌸ cii , ⌸ icc , ⌸ ici , ⌸ iic , and ⌸ iii .Figure 6 shows these eight energy fluxes, normalized by the energy dissipation rate ⌸ ␣␤␥ ͑k͒ / ͗⑀͘ versus k IK , together with the total flux denoted by ⌸ ttt .We find that ⌸ ccc coincides with ⌸ ttt all along the inertial range and that the other fluxes are close to zero, which confirms that the nonlinear dynamics is fully captured by the coherent fields.In the dissipative range, the coherent flux ⌸ ccc still dominates, though it begins to depart from the total flux ⌸ ttt , since ⌸ cci and ⌸ icc begin to build up for scales smaller than k IK Ӎ 0.3.The flux ⌸ cci is positive, while ⌸ icc is negative; these fluxes tend to compensate each other.The remaining terms are negligible for all scales.

E. Scale-dependent flatness
Scale-dependent flatness is constructed from the moments of the wavelet coefficients v The scale-dependent flatness of the vector field v at scale 2 −j is a measure of intermittency and quantifies the sparsity of the wavelet coefficients at each scale. 35It is defined by Here, M p,j ͑v͒ is the centered pth order moment of v at scale 2 −j defined by

͑14͒
where M ¯j͑v͒ is the mean value of v ˜, ᐉ at scale 2 −j given by The flatness is equal to 3 at all scales for a Gaussian white noise, which proves that this signal is not intermittent in this case.For more details we refer to Ref. 36.The scale-dependent flatness for the total, coherent, and incoherent velocity fields and magnetic fields, denoted by F j ͑u ␣ ͒ and F j ͑b ␣ ͒, respectively, are shown in Fig. 7 as a function of the dimensionless wavenumber k j IK .Note that ␣ = t, c, i, u t = u and b t = b.The scale 2 −j can be related to the wavenumber k j via the relation k j =2 j / 1.3, where 1/1.3 is the centroid wavenumber of the Coifman 12 wavelet.At low wavenumbers the statistical quantities of the incoherent contributions yield erroneous results, since too few wavelet coefficients of the incoherent vorticity and current density fields represent their contribution at large scales.Therefore, we only consider the five smallest scales and omit the large ones.
Figure 7 shows that the flatness factors of the total velocity and magnetic fields increase with k j IK .The increase in F j ͑b͒ ͑up to 220͒ for smaller scales is greater than that of F j ͑u͒ ͑up to 150͒.This stronger increase proves that the magnetic field is more intermittent than the velocity field, which is consistent with the conclusions drawn in previous works 37,38 showing that the pth-order structure function scaling exponent of the velocity field is larger than that of the magnetic field for integer p͑Ն4͒.
The flatness of the coherent velocity and magnetic fields is comparable to that of the total velocity and magnetic fields, respectively, which illustrates that the coherent vorticity sheets and current sheets are responsible for the intermittency of the velocity and magnetic fields.The situation for the incoherent contributions of the velocity and magnetic fields differs significantly: almost no dependence of flatness on scale is observed and the values are much smaller for k IK տ 0.2, compared to the values of the total and coherent fields.Nevertheless, these values are not equal to 3, the value obtained for a Gaussian white noise.

F. Scale-dependent compression rate
Figure 8 shows plots of the scale-dependent compression rates for the vorticity and current density fields, defined by versus the normalized wavenumber k j IK , where N c,j is the number of the wavelet coefficients at scale 2 −j retained by the coherent part of , and N t,j is the total number of wavelet coefficients at scale 2 −j , given by N t,j =7ϫ 2 j−1 .We observe that both compression rates decrease with increasing k j IK .For k j IK Շ 0.1, almost all coefficients are retained by the coherent fields, while the percentages of the retained coefficients decreases for k j IK Ͼ 0.1; thus the compression rates are improved.Note that due to the octree representation of the wavelet coefficients, the number N j is very large for small scales, i.e., for large j.This implies that the overall compression rate C is dominated by C j at small scales.The scale-dependent compression rates C j b and C j u are comparable to each other for a given scale.Now, let us discuss how small wavelet coefficients are retained by either coherent fields c or j c for a given scale.The percentage C j ഫ is defined by where N c,j ഫ is the number of directions and locations at each scale where the coefficients are retained by either the coherent part of or that of j.As shown in Fig. 8, C j ഫ exhibits a behavior similar to the scale-dependent compression rate C j of both fields.Note that the overall union rate, i.e., the ratio of the number of the wavelet coefficients retained by either the coherent vorticity field or the coherent current density field to the total number of the wavelet coefficients N, denoted by C ഫ , is 4.13%.

G. Scale-dependent overlap rate
Next, in Fig. 9, we examine scale-dependent overlap rate defined by which is a measure for characterizing the degree of spacescale overlap between coherent vorticity sheets and current sheets.Here, N c,j പ is the number of directions and locations at each scale where the wavelet coefficients, ˜, and j ˜, , belong to both coherent parts of and j.For k j IK Շ 0.1, the coherent vorticity sheets and current sheets overlap well.The scale-dependent overlap rate O j decreases with increasing k j IK for k j IK Ͼ 0.1, showing loglike decay for k j IK տ 0.2.
The observation of O j may be consistent with intuitive impression given from the visualization of intense vorticity and current regions shown in Fig. 2 ͑top͒.The overall overlap rate, i.e., the ratio of the number of the wavelet coefficients retained by both of the coherent fields to the number of the wavelet coefficients retained by either the coherent vorticity field or the coherent current density field for all scales, is 54.2%.

V. CONCLUSION AND PERSPECTIVES
We introduced a method for extracting coherent vorticity sheets and current sheets out of three-dimensional homogeneous MHD turbulence.This technique generalizes the CVE method, previously introduced for HD turbulence. 19,20It is based on the orthogonal wavelet decomposition of the vorticity and current density fields.Thresholding the wavelet coefficients allows us to split each of the fields into two orthogonal contributions: a coherent and organized component and an incoherent and random component.Coherent vorticity sheets and current sheets have thus been extracted from DNS data of forced incompressible homogeneous MHD turbulence without mean magnetic field, computed with a fully dealiased Fourier pseudospectral method at resolution 512 3 and a Taylor microscale Reynolds number R M of 154.We have shown that few wavelet coefficients are sufficient to represent the coherent vorticity sheets and current sheets, while the large majority of the coefficients correspond to incoherent background fields, which are structureless and contain no vorticity and current sheets.These findings contradict those of Politano et al. 2 who argued that intense vorticity and current sheets are strongly dissipative and thus cannot be viewed as coherent structures.According to our definition, the extracted sheets are indeed coherent.
The statistics of the coherent velocity and magnetic fields are similar to those of the total velocity and magnetic fields, respectively.The coherent fields contain most of the kinetic and magnetic energies and enstrophies of the total fields.The coherent kinetic and magnetic energy spectra coincide with the spectra of the total fields all along the inertial range and they differ only in the dissipative range.In contrast, the spectra of the incoherent fields are close to the spectra corresponding to energy equipartitions.
Concerning higher order statistics, we found that the PDFs of the total and coherent vorticity and current density have stretched exponential tails and coincide almost perfectly, while the incoherent contributions have reduced variances and exhibit exponential shapes.The PDFs of the incoherent velocity and magnetic fields yield a quasi-Gaussian distribution with strongly reduced variances.
Studying the flux of the kinetic and magnetic energy confirms that the nonlinear dynamics is fully captured by the coherent fields.We also introduced a scale-dependent overlap rate, a measure characterizing the degree of space-scale overlap between coherent vorticity sheets and current sheets.We found that the coherent vorticity sheets and current sheets overlap well for k j IK Շ 0.1, while the degree of overlap monotonically decreases with increasing k j IK for k j IK Ͼ 0.1.In the wavenumber range k j IK տ 0. Rosenberg et al. 39,40 for two-dimensional MHD flow; the former is a generalization of the CVS method.CVS is based on the deterministic computation of the coherent flow evolution using an adaptive wavelet basis and modeling the influence of the incoherent background flow, which was proposed for HD turbulence. 25Applications of CVS to twodimensional flows and to three-dimensional turbulent mixing layers can be found in Refs.26 and 27, respectively.Thus, CVCS is promising and should be pursued in future studies.

ACKNOWLEDGMENTS
The computations were carried out on the HPC2500 system at the Information Technology Center of Nagoya University.
This work was supported by a Grant-in-Aid for Young Scientists ͑B͒ No. 20760055 from the Ministry of Education, Culture, Sports, Science and Technology of Japan and by the 21st Century COE program "Frontiers of Computational Science."M.F. and K.S. acknowledge financial support from the Agence Nationale de la Recherche, project "M2TFPЉ and from the Association CEA-Euratom.The authors are grateful to K. Yoshida for fruitful discussion about DNS of MHD turbulence.

APPENDIX A: DIVERGENCE PROBLEM
The coherent and incoherent vorticity fields and current density fields obtained by the CVCE are not perfectly divergence-free, because the orthogonal wavelet transform does not commute with the divergence operator and the vector-valued wavelet basis used is not divergence-free.In practice, however, this is not crucial.For the coherent vorticity and current density fields, the divergent components of the coherent fields are 1.18% of the total kinetic enstrophy and 1.11% of the total magnetic enstrophy, respectively, in DNS 2 .They appear mostly in the dissipative range, as shown below.This observation is qualitatively similar to what is found for the CVE in hydrodynamic isotropic turbulence. 15,21n Fig. 10 we plot the enstrophy and current density spectra of the coherent and incoherent parts together with the divergence-free counterparts, denoted by a c Ј, a i Ј ͑a = , j͒, and their corresponding divergent parts, ٌ a .We observe that the divergent parts contribute little to the kinetic and magnetic enstrophies.Therefore, it is suggested that the lack of  perfect incompressibility of the coherent vorticity and current density fields is not a key issue in CVCS.

APPENDIX B: INFLUENCE OF THE NUMBER OF ITERATIONS IN CVCE
Here, we briefly discuss how the overall compression rates C u and C b depend on the number of iterations I in the CVCE algorithm.Figure 11 shows that C u and C b monotonically increase with I, tending to 14.5% and 15.4%, respectively.The coherent and incoherent fields of the vorticity and current density obtained after convergence ͑I = 100͒ retain the vorticity sheets and the current density sheets better than the coherent fields obtained with one iteration, I =1, as expected.The PDFs of the incoherent velocity and magnetic fields are quasi-Gaussian and the incoherent energy spectra exhibit k 2 slopes in a certain wavenumber range around k IK = 0.5 ͑figure omitted͒.These properties for the incoherent fields are the same as those obtained with I =1.
Figure 12 shows that scale-dependent overlap rates O j logarithmically decrease as k j IK increases in the range k j IK տ 0.4 for the coherent fields obtained for both I = 1 and I = 100.We find that both decay rates are the same.

FIG. 1 .
FIG.1.The evolution of the energy dissipation rate per unit mass, ͗⑀͘, in runs 1 and 2 performed up to t =8͑ϳ3.3Ti ͒ and t = 4.5͑ϳ2.7Ti ͒, respectively.Here, T i is the initial large eddy turnover time defined by T = L u / u 0 , and L u is the initial integral length scale given by L u = / ͑2u 0 2 ͒͐ 0 k max dkE u ͑k͒ / k.

FIG. 3 .
FIG.3.͑Color online͒ PDFs of the ᐉth components of ͑a͒ the velocity and ͑b͒ the magnetic fields for the total, coherent, and incoherent fields in DNS 2 .Dotted lines: The Gaussian distributions with zero mean and the same standard deviation as each corresponding incoherent component.

FIG. 4 .
FIG.4.͑Color online͒ PDFs of the ᐉth components of ͑a͒ the vorticity and ͑b͒ the current density fields for the total, coherent, and incoherent fields in DNS 2 .

FIG. 11 .FIG. 12 .
FIG. 11.Compression rates C u and C b vs number of iterations I in DNS 2 for = u and b.

TABLE II .
Characteristics of the two DNS runs at t = t f .

TABLE IV .
Statistical properties of the total, coherent, and incoherent fields obtained by CVCE applied to DNS 2 .
Author complimentary copy.Redistribution subject to AIP license or copyright, see http://php.aip.org/php/copyright.jsp 2, the decay is logarithmic.The above findings motivate the further development of adaptive numerical methods for MHD turbulence, such as the CVCS and adaptive grid refinement techniques with spectral elements.The latter has been developed by