Numerical linearized MHD model of flapping oscillations

Kink-like magnetotail flapping oscillations in a Harris-like current sheet with earthward growing normal magnetic field component Bz are studied by means of time-dependent 2D linearized MHD numerical simulations. The dispersion relation and two-dimensional eigenfunctions are obtained. The results are compared with analytical estimates of the double-gradient model, which are found to be reliable for configurations with small Bz up to values ∼0.05 of the lobe magnetic field. Coupled with previous results, present simulations confirm that the earthward/tailward growth direction of the Bz component acts as a switch between stable/unstable regimes of the flapping mode, while the mode dispersion curve is the same in both cases. It is confirmed that flapping oscillations may be triggered by a simple Gaussian initial perturbation of the Vz velocity.


I. INTRODUCTION
We present a numerical study of kink-like flapping oscillations performed by means of time-dependent two-dimensional linearized magnetohydrodynamic (MHD) simulations. Flapping oscillations are long-wavelength perturbations propagating in the magnetotail current sheet in transversal direction, mainly from the central part of the sheet towards the flanks, registered in the Earth, Venus, Jupiter, and Saturn magnetospheres. 1-17 Typically (for Earth magnetotail conditions), flapping waves are observed in the central part of the midtail at distances of 10 À 30 Re, where Re is the Earth radius, propagating with a speed of 20 À 60 km=s ($10 times less than the Alfv en velocity), an amplitude of 1 À 2 Re, and a quasiperiod of several minutes. Flapping motions in the distant magnetotail are reported in Refs. 18 and 19. Thus, these motions are essential in magnetotail physics.
The results of our numerical simulations demonstrate a good agreement with the so-called double-gradient (DG) model, 20,21 which allows us to investigate the model scope.
Even though the DG model is not the only mechanism to explain flapping oscillations, recent case studies have shown that it may suggest the best in-situ data interpretation for some events. 22,23 A detailed discussion of the DG and other MHD [24][25][26] and kinetic [27][28][29][30][31][32][33] models of flapping oscillations is provided in our previous papers. 34,35 To avoid repetitions, here we just outline briefly the essential features of the DG model and its limitations.
The model was developed in Ref. 20 for a current sheet of length L and half-width D, where the magnetic component B x ðzÞ $ z for jzj < D is a linear function within the sheet, and is constant outside: B x ¼ 6B 0 for jzj > D, while B z is a linear function of x. Here and throughout the paper, the xaxis points tailward from the Earth, the y-axis points dawnward, and the z-axis completes the righthanded system, as it is shown in Fig. 1. Considering long-periodic oscillations of the elongated current sheet (D=L ¼ ( 1) with a weak and slowly varying normal magnetic component (maxðB z ÞD= maxðB x ÞL ¼ ( ), the DG model predicts oscillations in the form of a plane wave $ exp½iðky À xtÞ, where the wavenumber k is the y-component of the wave vector k ¼ ð0; k; 0Þ. Two different types of initial perturbation dV z ðzÞ, odd and even, produce the waves of sausage type (lower frequencies) and kink type (higher frequencies), respectively.
In Ref. 21, the analytical solution was obtained also for a Harris-like current sheet (B x $ tanhðz=DÞ) with a small linear normal magnetic component B z ðxÞ. The analytical dispersion curve of the kink mode was derived in the following form: x ¼ x f ffiffiffiffiffiffiffiffiffiffiffiffiffiffi kD kD þ 1 r ; (1) where l 0 is the permeability of free space and x f is the typical frequency called "flapping frequency." The negative product of two magnetic gradients in Eq.
(2) yields imaginary x and, hence, the exponentially growing perturbations, i.e., instability. Vice versa, real x corresponds to the wave propagating in 6y directions. Thus, for the DG mode, the sign of the derivative @B z =@x operates as a switch between stable/unstable regimes, the same as for many other modes, a) Electronic mail: daniil.korovinskiy@oeaw.ac.at e.g., Refs. 26 and 36-43. Noteworthy, the dispersion relations (1) and (2) cover both stable and unstable modes. This resembles the Rayleigh-Taylor instability switching to the propagating wave, when heavy and light liquids swap places. It was also shown that a stable configuration should demonstrate a small gap of the total (gas þ magnetic) pressure, P, in the center of the sheet, while a central maximum of P yields the unstable sheet. 20 The simplicity and clearness of the DG model are achieved by a number of analytical limitations. The linearized MHD equations are solved for (a) incompressible plasma; (b) long-wavelength perturbations as compared to the field line curvature radius, k > R c ; (c) the flapping frequency x f is assumed to be constant; and (d) perturbations are assumed to be independent of x. Thus, the model was developed in the frame of a quasi-one-dimensional approximation; hence, further verification is not redundant.
The model was tested by means of two-dimensional linearized MHD and three-dimensional full MHD simulations of the unstable kink mode, developing in a Harris-like current sheet 44 and at the Pritchett and Coroniti 41 background. 34,35 It is found out that numerical growth rates are close to the analytical prediction with x f ¼ max½x f ðxÞ for 2D simulations and x f ¼ hx f ðxÞi for 3D simulations, where angle brackets mean averaging. It was shown also that an instability does not evolve when curvature radius is too large (R c > L > k). The perturbations were found to be almost uniform on the x-coordinate, in support of the analytical assumptions. 34 Hence, the analytical limitations enumerated above are partly justified. However, several open questions remain and are addressed in the present paper: (a) Does the analytical model match the oscillating regime? (b) How much does compressibility contribute to the development of oscillations? (c) How close are perturbations to the analytical eigenfunctions? (d) How substantial is the contribution of term $B z neglected in the analytical solution? 21 The paper is organized as follows: In Sec. II, we state the problem and analytical approach; the numerical method is described in Sec. III; results are presented in Sec. IV and discussed in Sec. V.

II. PROBLEM AND APPROACH
We study kink oscillations in a Harris-like current sheet with a small additional linearly varying normal magnetic component. Background magnetic field and plasma pressure are fixed as follows: where A ¼ 0.1, B ¼ À0:01, and P 0 ¼ 5. Mass density is assumed to be uniform, q 0 ¼ 1. Here, the magnetic field components are normalized to the field value at "infinity," B 0 ; pressure to B 2 0 =ðl 0 Þ; spatial variables are normalized to the sheet half-width D; and time to the Alfv en time D=V A , where V A is the Alfv en velocity.
Magnetic field lines and the total pressure P for the background configuration are shown in Fig. 2. In the x-range 0 < x < 20; Pðx; zÞ demonstrates a minimum in the center of the sheet (z ¼ 0) for any fixed value of x. However, Pðx; 0Þ is not uniform, the pressure gap deepens from the left to the right, reaching the lowest value (90% of P 0 ) at x ¼ 10 and then flattens out again. In the ranges x < 0 and x > 20, Pðx; zÞ has a bulge in the sheet center; hence, the configuration is unstable there. Thus, the region 0 < x < 20 complies both analytical conditions of the wave mode: positive product of two magnetic gradients and central minimum of the total pressure.
The set of compressional ideal MHD equations 45 is solved numerically using the perturbation method. All variables are represented as a sum of two terms: the equilibrium state U 0 and a small perturbation U 1 ; hence, MHD equations are linearized. The solution U 1 is found in the form of a wave, propagating in y-direction within the current sheet: U 1 ¼ dUðx; z; t; kÞ expðikyÞ, where the perturbation amplitude is dUðx; z; t; kÞ. Therefore, the system of linearized equations for amplitudes takes the conservative form dU ¼ ðdq; fdM i g; fdB i g; dEÞ i¼x;y;z : The variable U denotes the vector of plasma parameters, where M i ¼ qV i is the momentum, E ¼ P=ðj À 1Þ þ 0:5qV 2 þ0:5B 2 is the total energy density, and j is the polytropic index (the value of 5/3 is utilized if not otherwise stated). The expressions for the flux densities F x and F z and for the source function S are given in Appendix. These functions depend on the initial state U 0 , on perturbations dU, on the wave number k, on coordinates x and z, but not on coordinate y. In this way, the initial full MHD system of equations is reduced to a time-dependent two-dimensional problem for complex quantities dUðx; z; t; kÞ, where the physically meaningful quantities are <½dU expðikyÞ. Simulations are seeded with initial perturbation dU 0 . It is particularly convenient to initiate simulations with a perturbation of the normal component of the velocity, dV z0 . For real dV z0 and a background magnetic field component B y ¼ 0, this gives the set of eight non-zero physically meaningful functions: <ðdqÞ cosðkyÞ; <ðdM x Þ cosðkyÞ; =ðdM y Þ sinðkyÞ; <ðdM z Þ cosðkyÞ; <ðdB x Þ cosðkyÞ; =ðdB y Þ sinðkyÞ; <ðdB z Þ cosðkyÞ; <ðdEÞ cosðkyÞ.
Solving the system (6) for different values of the wavenumber k, we aim to obtain and investigate purely oscillating solutions with a frequency xðkÞ. While the exact analytical solutions of (6) cannot be derived, their estimations are provided by the DG model in the frame of a quasi-one-dimensional approximation and eigenmode analysis. These simplified solutions, presented below, form the basis for our approach to the problem.
For perturbations $ exp½iðky À xtÞ, the DG eigenfunctions of the problem are expressed via perturbations of the normal component of velocity, dV z ðz; kÞ Expressions (8)- (14) represent the solution of the system (11)- (15) in Ref. 21 for real dV z and x. Eigenfunctions dV z , in their turn, are expressed via Legendre functions P l k (see Eq. (19), ibid.) l ¼ Àk; For the kink mode, the parity condition requires k þ l ¼ 2m for m ¼ f0; 1; 2; :::g. Here, quantities z, k, and x are normalized, as stated above. Thus, for any fixed value of the wavenumber k, one has an infinite set of decreasing eigenfrequencies x m ðkÞ and corresponding eigenfunctions. For m ¼ 0, Eq. (17) yields the dispersion relation (1) and (2), where the terms D and 1=l 0 disappear due to normalization. Dispersion curves for m ¼ f0; :::; 6g are plotted in Fig. 3, and eigenfunctions dV z ðzÞ for k ¼ 1 and m ¼ f0; 1; 2g are shown in Fig. 4. For arbitrary initial perturbation $ exp½iky, the general solution of the simplified time-dependent system (Eq. (10) in Ref. 21) takes the form dUðx; z; t; kÞ ¼ R 1 m¼0 A m ðkÞdŨ m ðx; z; kÞ exp½Àix m ðkÞt, where dŨ m are the vectors of eigenfunctions and coefficients A m are to be found from an expansion of the initial perturbation. The solution may also be presented in the form of spectrum H½xðkÞ, where the peaking frequency, x peak , and width, 2Dx, depend on the initial perturbation. When the initial perturbation is given by the m-th set of eigenfunctions (8)- (15), the spectrum has the single line at x m ðkÞ.
Generally speaking, these analytical results do not span the true solution of the system (6). However, when the current sheet parameters satisfy the DG model assumptions, we suppose that DG estimates resemble the true solution with a good accuracy. This means that if we initiate the evolution of the sheet parameters with perturbations which are close to the m-th set of eigenfunctions (8)-(15), we may expect that after a transitional process, the solution of Eq. (6) does converge to oscillations with a frequency close to the m-th eigenfrequency of Eq. (17). Hence, the solution spectrum H½xðkÞ is supposed to have a pronounced peaking frequency, x peak ðkÞ, close to x m ðkÞ. The difference between obtained values x peak ðkÞ and x m ðkÞ quantifies the inaccuracy of the DG model approach and predictions. In the current paper, we use the zeroth mode of the solution (8)-(15) to achieve the oscillating solution of the system (6). The obtained dispersion curve is found to be rather close to the highest branch of the analytical dispersion curve, x 0 ðkÞ, given in Eqs.

III. METHOD
Equations (6) are solved numerically by means of the 3rd order central semi-discrete upwind scheme 46 with open boundary conditions @=@n ¼ 0. Semi-discrete schemes allow any appropriate ordinary differential equation (ODE) solver to be used for time integration of the equations for cell averages. In the simulations below, we have used the optimal (in the sense of Courant-Friedrichs-Lewy (CFL) coefficient and the computational cost) strong stability preserving Runge-Kutta method of the 3rd order described in Ref. 47. The integration time step (CFL number ¼ 0:5) is adopted to ensure the convergence of the results with respect to values of the time step. The ðr Á BÞ ¼ 0 constraint is enforced on each time step using the method of projection. 48 To initiate oscillations, two types of initial perturbations are used. The first one is the Gaussian perturbation dV z0 ¼ expðÀz 2 Þ, resembling the analytical eigenfunction V z ðzÞ by Eq. (15) for m ¼ 0, shown in Fig. 4 by the blue curve. Another one represents a full set of analytical eigenfunctions (8)-(15) for m ¼ 0, shown in Fig. 5 for k ¼ f0:5; 1; 2g. The approximate solution (8)- (15) does not match exactly the eigenfunctions of the two-dimensional problem under consideration; however, it is closer than the single Gaussian perturbation of V z . Therefore, the numerical solution reaches the oscillating regime more quickly if analytical eigenfunctions rather than Gaussian perturbation of V z are used to initialize the evolution of (6). In both cases, we observe that after some transitional process, the solution converges to damped harmonic oscillations with a frequency close to the prediction of the DG model. The damping of oscillations is related to numerical dissipation (that can be decreased by reduction of the mesh steps) and to the use of open boundary conditions.
The remnants of transitional processes and the attenuation of oscillations due to numerical damping lead to an inaccuracy of the results. Namely, for any fixed k, instead of a single spectrum, oscillations are characterized by spectrum curves H j ðxÞ, which are calculated individually for each quantity using fast Fourier transform. The more accurate the numerical solution is, the more exact the spectra H j ðxÞ coincide. The spectra widths 2Dx j ðkÞ estimate the "quality" of the DG approximation for a given parameters of the current sheet, though additional spectral broadening results from the numerical dissipation.
These transitional and numerical effects complicate the accurate analysis of the numerically obtained spectra H j ðxÞ. However, when the DG model assumptions are fulfilled, the spectra turn out to be narrow enough, and peaking frequenciesx peak;j coincide to a high accuracy. The symbol "tilde" is used here to separate the true eigenfrequency xðkÞ and its approximate estimate by numerical simulations,x j ðkÞ. Therefore, we have used a simple approach to estimate the frequency of oscillations. Namely, when the numerical Normalization: maxðjdV z jÞ ¼ 1. solution of the system (6) for some fixed value of k is obtained, the frequency of oscillations is computed at one point in the center of the sheet by fitting the quantities <½dqðtÞ; <½dV x ðtÞ; =½dV y ðtÞ; <½dV z ðtÞ; <½dB x ðtÞ; =½dB y ðtÞ; <½dB z ðtÞ; <½dEðtÞ with the functions f j ðtÞ ¼ a j expðc j tÞ sinðx j t þ / 0j Þ; j ¼ f1; 2; :::; 8g: (18) Then, the value of the eigenfrequency x is estimated by averagingx ¼ hx j i6r x . The standard deviation r x quantifies the average width of spectra H j ðxÞ. Numerical attenuation of oscillations is quantified by negative c j .
The numerical results were tested for sensitivity to the box crosswise size and mesh step. In Fig. 6, the wave frequencyx, calculated at the point ðx 0 ¼ 10 À 1=2 5 ; z 0 ¼ À1=2 5 Þ, is plotted for three values of the wavenumber k ¼ f0:1; 1; 10g for different values of the computational box boundary location z max . It is seen that for any value of k, we can place the z-boundary far enough to make the results independent of the boundary location. The lower the wavenumber, the wider should the computational box be, and vice versa. This is due to the fast broadening of the eigenfunctions' profiles across the sheet with decreasing value of k, as it is seen in Fig. 5.
The dependence ofx on the mesh step dx ¼ dz ¼ 1=2 n , where n ¼ f2; 3; 4; 5; 6g, is shown in Fig. 7. It is seen that for a small values of k ($1 or less), we can consider n ¼ 4 as the maximal required grid refinement. For 1 < k < 5, n ¼ 5 seems to be enough, while for k > 5, the furthermost refinement is eligible. However, with the value n ¼ 4, one can expect the error within 10% or less for any k 2 ½0; 10.  18), as a function of wave number (black curve), mesh step (colored asterisks), and polytropic index j (colored crosses). It is seen that the attenuation of oscillations is slower for longer wavelengths and thinner grids. This looks reasonable for pure numerical effect, as a decrease of the wavelength is equivalent to an increase in the mesh step. Conformably, reduction of the mesh step reduces numerical damping.
In the incompressible limit j ! 1, the decrement indicates a trend to negligible values. This means that open boundary conditions do not prevent energy outflow carried away by acoustic mode. Comparison of the numerically obtained frequencies with analytical estimates, presented in Section IV, shows that the effect of these energy losses is minor. Besides, the projection procedure used to enforce the ðr Á BÞ ¼ 0 condition also introduces some amount of dissipation into the numerical scheme.
At the same time, a simple approach adopted in the current paper leaves out all physical mechanisms of the oscillations' damping. In the frame of an ideal MHD, several such mechanisms may be reported. In Ref. 49, it was shown that oscillations may lose their energy in short time (few periods) FIG. 6. Wave frequenciesx ¼ hxi6r x , normalized to the analytically predicted values (1), are calculated in the point ðx 0 ¼ 10 À dx=2; z 0 ¼ Àdz=2Þ and shown as functions of the z boundary location z max for wavenumbers k ¼ 10 (blue), k ¼ 1 (green), and k ¼ 0.1 (red). Simulations are performed in the box ½6; 14 Â ½Àz max ; z max with mesh step dx ¼ dz ¼ 1=2 4 . The values of z max are normalized to z 0 ¼ 64 for k ¼ 0.1, to z 0 ¼ 8 for k ¼ 1, and to z 0 ¼ 4 for k ¼ 10.  (1), are calculated at the point ðx 0 ¼ 10 À dx=2; z 0 ¼ Àdz=2Þ and shown as functions of the mesh step dx ¼ dz ¼ 1=2 n for wavenumbers k ¼ 10 (red), k ¼ 5 (blue), k ¼ 1 (green), and k ¼ 0.5 (black). Simulations are carried out in the x-range ½9; 11.
due to the interaction with the ionosphere. Threedimensional MHD simulations 34 reveal that the flapping frequency is a non-local quantity, which can be affected by changes in remote parts of the current sheet. Furthermore, in Ref. 50, it was demonstrated that in thin current sheets, the Hall effect may suppress the DG mode in the direction opposite to the current.

IV. RESULTS
A series of simulations is carried out in the x-range ½9; 11, covering the deepest total pressure minimum at x ¼ 10, where the background magnetic component jB z j 0:01. In this region, the problem demonstrates the best resemblance to the quasi-one-dimensional analytical model; hence, the best agreement of the results may be expected. Fig. 9 shows the sample solution obtained in the box ½x Â z ¼ ½9; 11 Â ½À64; 64 with mesh step dx ¼ dz ¼ 1=2 4 for k ¼ 1. Normalized quantities <ðdqÞ (blue), <ðdV x Þ (dark green), =ðdV y Þ (red), <ðdV z Þ (cyan), <ðdB x Þ (magenta), =ðdB y Þ (light green), <ðdB z Þ (black), and <ðdEÞ (orange) are plotted in the top panel versus time in the point ð10 À dx=2; Àdz=2Þ. In the bottom panel, the values averaged across the sheet at the cross-section x 0 ¼ 10 À dx=2 are shown. It is seen that the frequency of the sheet oscillations, x aver ¼ 0:0687, is almost the same to the frequency measured in the central point,x point ¼ 0:0685. This means that the entire current sheet oscillates with the same frequency, and oscillations may be registered far enough from the sheet center. The transitional process from initial configuration to the oscillatory solution is seen in the plots of dV y , dV z , and dE at time t < 100.
The dispersion curves of the kink mode, shown in Fig. 10, are calculated in the same box for two types of initial perturbation: eigenfunctions (8)- (15) for all quantities (red), or single dV z0 ¼ expðÀz 2 Þ (blue). Blue points lay almost exactly on the red curve, confirming that the result does not depend on the type of initial perturbation if it bears a resemblance to the zeroth mode eigenfunctions.
One can also see that the standard deviation increases with growing k values, meaning that the spectra broadened, as it is shown in Fig. 11. The Fourier spectra half-width Dx may be estimated as 0.015 for k ¼ 0.4, 0.03 for k ¼ 2, and 0.05 for k ¼ 10. In normalized units Dx=x peak , this gives %f0:3; 0:4; 0:5g, respectively. Fig. 11 illustrates that for small and moderate wavenumbers, the valuex peak is the same for all eight quantities. For high k, valuesx peak;j start to diverge, meaning that the accuracy of the numerical solution diminishes due to the high numerical dissipation.
The numerical dispersion curves of Fig. 10 (red, blue) demonstrate a systematical overshoot of the analytical prediction (green) for k > 4 and undershoot for smaller k values. It turns out that the difference between numerical and analytical solutions displays the contribution of the plasma compressibility, which was neglected in the analytical model, but kept in the numerical one. In Fig. 13, frequenciesx ¼ hx j i 6r x , normalized to the corresponding analytical values, are shown as functions of the polytropic index j for wavenumbers k ¼ 10 (red) and k ¼ 2 (blue). Fitting curves are f 1 ¼ 1=f1 À exp½À2:42ðlgj þ 1Þ 1=3 g and f 2 ¼ tanhf1:44lg½17 þ 4:72jg, respectively. It is seen that numerical values tend to analytical predictions in the incompressible (j ! 1) limit in both cases. The calculated perturbations are compared to the analytical predictions (8)- (15) in Figs. 14 and 15. Profiles of perturbations are computed for the wavenumber k ¼ 1 at x 0 ¼ 10 À 1=2 5 for time t ¼ 200. Perturbations and analytical quantities are normalized to the maxðjdV z jÞ. It is seen that perturbations of velocity components V y and V z demonstrate brilliant agreement, meaning that the divergence-free equation / @dV x =@x þ @dV y =@y þ @dV z =@z ¼ 0 used in Eq. (10) of the work 21 is fulfilled to a high accuracy. Indeed, the perturbation of V x is found to be really weak (the analytical prediction is zero), contrary to the perturbation of the mass density, which turns out to be comparable with dV z .
Perturbations of magnetic components B x and B y demonstrate a good agreement with the analytical predictions in amplitude. The perturbation dB z overshoots the theoretical quantity two times and the energy perturbation occurs two times lower than the analytical estimate, which should be accounted for the contribution of the mass density variation.

062905-8
Korovinskiy et al. seen that with increasing distance from the sheet center (x ¼ 10), where B z ¼ 0, spectra broaden more and more. For x 0 ¼ 4, the normalized spectra half-width Dx=x peak reaches the value of 1. Corresponding profiles of <½dV z ðz; x ¼ x 0 Þ are shown in Fig. 17. With decreasing value of x 0 (growing B z ), the function dV z ðzÞ deviates more and more from the analytical profile. This result makes the major difference with simulations of the DG instability, where perturbation profiles were found almost uniform on x. 34

V. DISCUSSION AND CONCLUSIONS
In this paper, we have investigated long-wavelength kink-like waves propagating within the magnetotail current sheet in transversal direction. Oscillations of that type resemble well-known flapping oscillations. The system of twodimensional time-dependent linearized MHD equations is solved numerically to calculate perturbations of a Harris-like current sheet with a small additional normal magnetic component B z ðxÞ, decreasing in the tailward direction.
Numerical simulations reveal that the analytical DG model 21 provides a reliable solution, when the current sheet parameters satisfy the model assumptions. It was shown previously that this requirement can be fulfilled in the Earth magnetotail. 51 It turns out that the most severe limitation of the DG model is the requirement of a small normal magnetic field component. As can be seen from Figs. 16 and 17, we may estimate the crucial value of B z for our particular configuration as 0:05 B 0 , where B 0 is the lobe magnetic field ($20 nT during the substorm growth phase 52,53 ). At higher values of B z , the analytical profiles of perturbations mismatch the numerical solution, while the half-width of the Fourier spectra of perturbations, Dx, exceeds the peaking frequencyx peak . Nevertheless, the value ofx peak stays close to the DG model prediction.
Thus, while kink-type low-frequency (flapping) oscillations may exist even in current sheets with higher values of B z , the DG model may provide only estimates of their frequencies, but not the profiles of perturbations. This means that the accuracy of the analytical solution for perturbations falls down with growing absolute value of B z , while estimates of the typical frequency of oscillations stay accurate. The latter depends on the basic properties of the current sheet configuration (product of two magnetic gradients and minimum of total pressure). In our particular configuration, the product of magnetic gradients does not change with growing B z , while total pressure minimum is preserved in interval x 2 ð0; 20Þ; hence, the low dependence of the frequency estimates accuracy on jB z j.
At the same time, configurations with B z below 1 À 2 nT are commonly found in THEMIS observations, 53 left), and x 0 ¼ 4 À dx=2 (bottom, right). Simulations are seeded with eigenfunctions (8)- (15) and carried out in the range x 2 ½9; 11 for the top left picture and x 2 ½1; 9 for others. Mesh step dx ¼ dz ¼ 1=2 4 . The analytically predicted frequency is plotted by dashed red lines. Labels of the curves are given in the legend of Fig. 11.
FIG. 17. Normalized profiles of dV z ðzÞ for wavenumber k ¼ 1 are plotted in the cross-sections x 0 ¼ 10 À dx=2 (green), x 0 ¼ 8 À dx=2 (red), x 0 ¼ 6 À dx=2 (black), and x 0 ¼ 4 À dx=2 (blue). Simulations are seeded with eigenfunctions (8)- (15) and carried out in the range x 2 ½9; 11 for the green curve and x 2 ½1; 9 for others. Mesh step dx ¼ dz ¼ 1=2 4 . agreement with the numerical solution. The analytical dispersion curve matches the numerical one with the accuracy up to several percents (see Fig. 10). Simulations with increased value of polytropic index reveal that this minor discrepancy can be explained by plasma compressibility, neglected in the analytical model. Hence, this limitation of the model is found to be not crucial. This supports the conclusion made in Ref. 51, claiming that compressibility may be neglected in studies of low-frequency oscillations x=ðkc s Þ ( 1, where c s is the sound speed (c s % 3 in our simulations with polytropic index j ¼ 5=3).
Then, the condition 1=R c > k > 1=L is satisfied in our simulations for any k 2 ½0:1; 10, as far as R c ¼ jB z = ð@B x =@zÞj z¼0 ¼ jB z j and L $ jmaxðB z Þ=ð@B z =@xÞj $ 10. However, with growing value of k, the difference between analytical and numerical values of x increases, and Fourier spectra of perturbations, shown in Fig. 11, broaden gradually. This is in accordance with our previous findings, 34 stating that the magnetic field line curvature radius is a lower limit for the DG mode wavelength.
A comparison of perturbation profiles with analytical eigenfunctions (Figs. 14 and 15) demonstrates an excellent agreement for quantities dV y and dV z , a good agreement for dB x , and dB y , and a qualitative agreement for dB z and dE. This represents the contribution of the underlined terms of Eq. (10) in Ref. 21, which were neglected in an analytical solution. Noteworthy, the relative magnitude of perturbation dV x amounts to $10 À3 . Taking into consideration the relative phase shift d/ ¼ 111 between dV z and dV x , calculated by Eq. (18), we obtain maxjdV x j % 3 Â 10 À3 , which is much less than maxjdV z j ¼ 1; hence, no plasma motion in the x-direction is found.
Summing up, linearized MHD simulations are found to be a handy tool for studying flapping oscillations in the current sheets with varying B z component. When the latter is small enough, the analytical DG model provides accurate theoretical estimates. Essentially, dispersion curves of oscillations, obtained here, and the corresponding instability, studied in Ref. 44, are close to the same analytical prediction (1) and (2). This supports the wave generation mechanism suggested in Ref. 54. It assumes that the initial perturbation is produced by an instability, developing in the current sheet with tailward growing B z . When @B z =@x changes its sign, the instability switches to the propagating flapping wave. The unique dispersion relation for stable/unstable regimes on the one hand and the frequent change of the B z growth direction in the midtail 55 on the other hand make this scenario probable.
In our simulations, the evolution of system (6) has been initiated either by single Gaussian perturbation of V z or by a set of analytical eigenfunctions (8)-(15) of zero eigenmode. The discrepancy between these starting functions and the solution of (6) grows with increasing value of k. Nevertheless, in both cases and for all k values, the solution converges to the same purely oscillating (up to numerical damping) mode with identical eigenfrequency xðkÞ. This suggests that the flapping oscillations can be triggered by arbitrary initial perturbations bearing some resemblance to the solutions of system (6).
Perturbations similar to a Gaussian perturbation of V z velocity may be caused by many reasons. In particular, one possible mechanism was considered in Ref. 51, where it was shown that bursty bulk flows (BBFs) may produce kink-type oscillations of the current sheet. This mechanism has a substantial advantage, because it may explain the non-zero xcomponent of the wave vector of some flapping wave fronts observed in the magnetotail. 10 On the other hand, the DG model represents a long-wavelength limit of the general ballooning mode. 34 Generally speaking, the wave vector of the latter has k x 6 ¼ 0, while k ¼ ð0; k y ; 0Þ is a limiting case only. In view of that, a generalization of the numerical and analytical models for the case of non-zero k x is the challenging subject for future studies.