# Dynamical Model for Nonlinear Mirror Modes near Threshold

E.A. Kuznetsov, T. Passot, P.L. Sulem

<figure><img src="image_1.png" /></figure>

## To cite this version:

E.A. Kuznetsov, T. Passot, P.L. Sulem. Dynamical Model for Nonlinear Mirror Modes near Threshold. Physical Review Letters, 2007, 98, pp.235003. 10.1103/PhysRevLett.98.235003. hal-00388125

HAL Id: hal-00388125

https://hal.science/hal-00388125v1

Submitted on 13 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Dynamical Model for Nonlinear Mirror Modes near Threshold

E. A. Kuznetsov,¹,² T. Passot,³ and P. L. Sulem³

¹L. D. Landau Institute of Theoretical Physics, 2 Kosygin Street, 119334 Moscow, Russia

²P. N. Lebedev Physical Institute, 53 Leninsky Ave., 119991 Moscow, Russia

³CNRS, Observatoire de la Côte d'Azur, BP 4229, 06304 Nice Cedex 4, France

(Received 12 January 2007; published 8 June 2007)

Using a reductive perturbative expansion of the Vlasov-Maxwell (VM) equations for magnetized plasmas, a pseudodifferential equation of gradient type is derived for the nonlinear dynamics of mirror modes near the instability threshold. This model, where kinetic effects arise at a linear level only, develops a finite-time singularity, indicating the existence of a subcritical bifurcation. A saturation mechanism based on the local variations of the ion Larmor radius, is then phenomenologically supplemented. In contrast with previous models where saturation is due to the cooling of a population of trapped particles, the resulting equation correctly reproduces results of numerical simulations of VM equations, such as the development of magnetic humps from an initial noise, and the existence of stable large-amplitude magnetic holes both below and slightly above threshold.

DOI: 10.1103/PhysRevLett.98.235003

PACS numbers: 52.35.Py, 52.25.Xz, 94.05.-a, 94.30.cj

Magnetic holes associated with maxima of density and pressure are often encountered in planetary magnetosheaths close to the magnetopause, and in the solar wind (see, e.g., [1,2]). These structures are often viewed as ultralow frequency (ULF) waves resulting from the mirror instability [3] that develops in a collisionless plasma with ion temperature anisotropy and a beta exceeding a few units. Such conditions can be met under the effect of the plasma compression in front of the magnetopause [4], where the effect of ULF waves on the energy and mass transfer remains an important open problem in the context of the interactions between the solar wind and the magnetosphere. Mirror structures are also observed when the plasma is linearly stable [5,6], which may be viewed as the signature of a bistability regime resulting from a subcritical bifurcation, whose existence was interpreted on the basis of a simple energetic argument within the simplified description of anisotropic magnetohydrodynamics [7].

The linear mirror instability has been extensively studied analytically, at least near threshold ([8] and references therein), and also by means of particle-in-cell (PIC) simulations [9]. It turns out that wave-particle resonance plays a central role in driving the instability, while finite Larmor radius (FLR) effects are at the origin of the quenching of the instability at small scales. In contrast, a theoretical understanding of the nonlinear phase remains limited, based on a phenomenological modeling of particle trapping [10,11] that hardly reproduces simulations of Vlasov-Maxwell equations [12]. The aim of this Letter is to provide a nonlinear model for the mirror mode dynamics, based on a reductive perturbative expansion of the Vlasov-Maxwell (VM) equations. For the sake of simplicity, we consider an electron-proton plasma with cold electrons.

From VM equations, one easily derives the equation governing the mean proton velocity $u$ that, when neglecting electron inertia, reads

$$ \rho \frac{du_i}{dt} = -\partial_j \left[ \left( p_{\perp} + \frac{|B|^2}{8\pi} \right) \delta_{ij} - \left( p_{\perp} - p_{\|} + \frac{|B|^2}{4\pi} \right) \hat{b}_i \hat{b}_j + \Pi_{ij} \right]. \quad (1) $$

Here $\Pi$ denotes the nongyrotropic contribution to the pressure tensor and $\hat{b} = B/|B|$. The other notations are standard. When the control parameter (temperature anisotropy or beta of the plasma) is close to the instability threshold, the left hand side of Eq. (1) becomes subdominant. The equation governing the mirror dynamics is then obtained perturbatively by expanding the projection of the resulting pressure balance on the plane perpendicular to the local field. In this approach, the ion pressure tensor elements are computed from the VM system, near a bi-Maxwellian equilibrium state characterized by temperatures $T_{\perp}^{(0)}$ and $T_{\|}^{(0)}$ and a constant ambient magnetic field of amplitude $B_0$ taken along the $z$ direction. At a linear level, the instability growth rate takes the form

$$ \gamma = |k_z| v_{\text{th}\parallel} \frac{\beta_{\|}}{\sqrt{\pi} \beta_{\perp}} \left[ \frac{\beta_{\perp}}{\beta_{\|}} - 1 - \frac{1}{\beta_{\perp}} - \frac{k_z^2}{k_{\perp}^2 \beta_{\perp}} \left( 1 + \frac{\beta_{\perp} - \beta_{\|}}{2} \right) - \frac{3}{4 \beta_{\perp}} k_{\perp}^2 \rho_L^2 \right], \quad (2) $$

which involves the ion parallel thermal velocity $v_{\text{th}\parallel} = (2T_{\|}^{(0)}/m)^{1/2}$, the ion Larmor radius $\rho_L = (2T_{\perp}^{(0)}/m)^{1/2}/\Omega$ [where $\Omega = eB_0/(mc)$ is the ion gyrofrequency], and the parameters $\beta_{\perp} = 8\pi p_{\perp}/B_0^2$ and $\beta_{\|} = 8\pi p_{\|}/B_0^2$. This growth rate can be recovered by expanding the general expression given in [8], in the limit of small transverse wave numbers. The instability takes place when $\beta_{\perp}/\beta_{\|} > 1 + \beta_{\perp}^{-1}$ and, near threshold, develops in quasiperpendicular directions, making the longitudinal magnetic perturbation dominant.
Because of the lack of space, the reductive perturbative expansion performed on the VM system to estimate the components of the pressure tensor [to be substituted in Eq. (1)], cannot be detailed here and will be published elsewhere. A few remarks are, however, in order. The linear part of the resulting Eq. (5) written below corresponds to the instability growth rate given by Eq. (2). Note that the time derivative is due to the resonance between the nonpropagating mirror wave and particles with almost zero parallel velocity. FLR corrections originate not only from the gyroviscous tensor Π in Eq. (1) but also include contributions from the gyrotropic pressures. Furthermore, the nonlinear term is not affected by FLR corrections (nor by resonances) and will be, in the framework of this Letter, more conveniently computed within the drift kinetic approximation [13] that assumes scales large enough for nongyrotropic effects to be negligible. For this purpose, noting the absence of parallel electric field in the case of cold electrons, and the subdominant character of the electric drift, we write the drift kinetic equation for the ions in the form

$$\partial_t f + v_{\parallel} \hat{b} \cdot \nabla f - \mu \hat{b} \cdot \nabla |B| \partial_{v_{\parallel}} f = 0. \quad (3)$$

In this equation, the adiabatic invariant $\mu = v_{\parallel}^2/(2|B|)$ plays the role of a parameter. The pressures $p_{\parallel}$ and $p_{\perp}$ are given by the integrals $p_{\parallel} = 2\pi m B \int v_{\parallel}^2 f d\mu dv_{\parallel}$ and $p_{\perp} = 2\pi m B^2 \int \mu f d\mu dv_{\parallel}$. Measuring the distance to threshold by $\epsilon \equiv \beta_{\perp} \chi^{-1} |\beta_{\perp}/\beta_{\parallel} - 1 - \beta_{\perp}^{-1}| \ll 1$ with $\chi = 1 + (\beta_{\perp} - \beta_{\parallel})/2$, we rescale parallel (perpendicular) space derivatives by $\epsilon$ (respectively, $\epsilon^{1/2}$) and time derivative by $\epsilon^2$, consistently with the linear growth rate. This scaling assumes that the instability is not early arrested by quasilinear effects and particle trapping, a point discussed below. We also expand $f$ and $B_z$ up to second order in powers of $\epsilon$, denoting by a superscript (p) the corresponding coefficient of $\epsilon^p$. We note that the parallel current is subdominant, leading to $B_{\perp}^{(3/2)} = (-\Delta_{\perp})^{-1} \nabla_{\perp} \partial_z B_z^{(1)}$.

Defining $\tilde{p}_{\perp}^{(n)} = \pi m \int v_{\perp}^2 f^{(n)} v_{\perp} dv_{\perp} dv_{\parallel}$, we have

$$p_{\perp}^{(2)} = (B_z^{(1)}/B_0)^2 p_{\perp}^{(0)} + 2(B_z^{(1)}/B_0) \tilde{p}_{\perp}^{(1)} + \tilde{p}_{\perp}^{(2)},$$

up to an additional contribution proportional to $B_z^{(2)}$ that cancels out in the final equation due to the threshold condition.

On the considered time scale, the effect of nonlinear Landau resonance is negligible in the nonlinear contributions that can thus be estimated from the equations

$$v_{\parallel} \partial_z f^{(1)} - \mu \partial_z B_z^{(1)} \partial_{v_{\parallel}} f^{(0)} = 0, \\ v_{\parallel} \partial_z f^{(2)} + (2\mu^2/v_{\text{th}\parallel}^2) B_z^{(1)} \partial_z B_z^{(1)} \partial_{v_{\parallel}} f^{(0)} = 0.$$

For an equilibrium bi-Maxwellian distribution $f^{(0)} \propto \exp(-v_{\parallel}^2/v_{\text{th}\parallel}^2 - 2\mu B_0/v_{\text{th}\perp}^2)$, we have $f^{(1)} = -(2\mu/v_{\text{th}\parallel}^2)B_z^{(1)}f^{(0)}$ and $f^{(2)} = (2\mu^2/v_{\text{th}\parallel}^4)(B_z^{(1)})^2 f^{(0)}$. As a consequence, because of the vicinity to threshold we obtain

$$p_{\perp}^{(2)} + \frac{(B_z^{(1)})^2}{8\pi} = \frac{3(B_z^{(1)})^2}{8\pi} \left(1 + \frac{1}{\beta_{\perp}}\right) > 0. \quad (4)$$

We finally get the closed equation for the longitudinal magnetic perturbation that, after rescaling, reads

$$\partial_T U = -\mathcal{H} \partial_Z [(\sigma - \Delta_{\perp}^{-1} \partial_{ZZ} + \Delta_{\perp})U - 3U^2], \quad (5)$$

with $\sigma = \pm 1$, depending of the positive or negative sign of $\beta_{\perp}/\beta_{\parallel} - 1 - \beta_{\perp}^{-1}$. Here $\mathcal{H}$ is the Hilbert transform along the direction of the ambient magnetic field. The rescaled variables are defined as $Z = (2/\sqrt{3})\epsilon z \rho_L^{-1} \chi^{1/2}$, $\mathbf{R}_{\perp} = (2/\sqrt{3})\sqrt{\epsilon} \mathbf{r}_{\perp} \rho_L^{-1} \chi^{1/2}$, and $T = (2/\sqrt{3})\epsilon^2 \Omega t (\sqrt{\pi} \beta_{\perp})^{-1} \times (\chi \beta_{\parallel}/\beta_{\perp})^{3/2}$, together with $B_z^{(1)}/B_0 = \epsilon 2 \chi \beta_{\perp} (1 + \beta_{\perp})^{-1} U$.

Equation (5) possesses the remarkable property of being of the form $\partial_T U = -\hat{K}_Z \delta F / \delta U$, where $\hat{K}_Z = -\mathcal{H} \partial_Z$ is a positive definite operator (whose Fourier transform is $|K_Z|$), and

$$F = \int \left[ -\frac{\sigma}{2} U^2 + \frac{U}{2} \Delta_{\perp}^{-1} \partial_{ZZ} U + \frac{1}{2} (\nabla_{\perp} U)^2 + U^3 \right] d^3 R$$

has the meaning of a free energy or a Lyapunov functional. This quantity can only decrease in time, since

$$\frac{dF}{dT} = \int \frac{\delta F}{\delta U} \partial_T U d^3 R = - \int \frac{\delta F}{\delta U} \hat{K}_Z \frac{\delta F}{\delta U} d^3 R \le 0. \quad (6)$$

The time derivative of $F$ identically vanishes at the stationary localized solutions defined by the equation

$$-\frac{\delta F}{\delta U} = (\sigma - \Delta_{\perp}^{-1} \partial_{ZZ} + \Delta_{\perp})U - 3U^2 = 0,$$

for which one easily checks that nonzero solutions do not exist when $\sigma = +1$. As a consequence, $dF/dT$ is strictly negative above threshold. In the nonlinear regime, negativity of this derivative implies $\int U^3 d^3 R < 0$, which corresponds to the formation of magnetic holes.

Equation (5) simplifies when the spatial variations are limited to a direction making a fixed angle with the ambient magnetic field. After a simple rescaling, one gets

$$\partial_T U = \hat{K}_{\Xi} [(\sigma + \partial_{\Xi\Xi})U - 3U^2], \quad (7)$$

where $\Xi$ is the coordinate along the direction of variation. This equation can be referred to as a “dissipative Kortveg–de Vries (KdV) equation”, since its stationary solutions coincide with those of the usual KdV equation. The presence of the Hilbert transform in Eq. (7) nevertheless leads to a dynamics significantly different from that described by soliton equations.

Equation (7) was integrated using a spectral method above threshold with an initial condition in the form of a sine function involving several periods within the computational domain. After an initial phase associated with the linear instability, we observe the formation of several magnetic holes whose number is progressively reduced to
one. After a while, the solution blows up. Its profile at a time close to collapse is displayed in Fig. 1.

In order to characterize the nature of the singularity of Eq. (7), it is convenient to introduce the similarity variables $\xi = (T_0 - T)^{-1/3}\Xi$, $\tau = -\log(T_0 - T)$, and to look for a solution in the form $U = (T_0 - T)^{-2/3}g(\xi, \tau)$, where $g(\xi, \tau)$ satisfies the equation

$$\partial_{\tau}g + (2/3)g + (\xi/3)\partial_{\xi}g = \hat{K}_{\xi}[\partial_{\xi\xi}g - 3g^2] + e^{-\tau}\hat{K}_{\xi}g.$$

At large $|\xi|$, that corresponds to the limit $T \to T_0$, the asymptotic solution $\tilde{g}$ obeys $2\tilde{g} + \xi\partial\tilde{g}/\partial\xi = C\xi^{-2}$, where $C = \frac{2}{\pi} \int_{-\infty}^{\infty} g^2(\xi')d\xi' > 0$, and has the form $\tilde{g} = C\xi^{-2}\log|\xi/\xi_0|$. For $U$, it gives the asymptotic solution $U_{asymp} = (C/\Xi^2)\log|\Xi/\Xi_0(t)|$ with $\Xi_0(t) = (T_0 - T)^{1/3}\xi_0$, that, as $T \to T_0$, has an almost time independent tail. For $|\Xi| < (T_0 - T)^{1/3}|\xi_0|$ the solution is negative and becomes singular as $\Xi$ approaches the origin. The inset in Fig. 1 that displays $U_{asymp}$ for $\Xi_0 = 5$ and $C = 10$, shows a good agreement with the simulation of the initial value problem.

Asymptotically self-similar solutions can also be constructed in several dimensions, when rescaling the longitudinal coordinate by $(T_0 - T)^{1/2}$, the transverse ones by $(T_0 - T)^{1/4}$ and the amplitude of the solution by $(T_0 - T)^{-1/2}$. Existence of a finite-time singularity for the initial value problem can be established for initial conditions for which the functional $F$ is negative, when the term involving $\sigma$ can be neglected, an approximation consistent with the dynamics.

At the level of the primitive VM equations, the singularity can be viewed as the signature of a subcritical bifurcation that cannot be described in the framework of the present asymptotics where no saturation originating from nonlinear kinetic effects is retained.

Among these effects, we concentrate on particle trapping and nonlinear FLR corrections. As argued in [14],

<figure><img src="image_3.png" /><figcaption>FIG. 1 (color online). Profile of the solution of Eq. (7) in the linearly unstable regime, close to the collapse time (here and in the other figures, the horizontal axis is labeled by grid point numbers); insert: Asymptotic profile $U_{asymp}$ for $\Xi_0 = 5$ and $C = 10$.</figcaption></figure>

saturation results from both particle trapping (efficient at large $v_\perp$) and quasilinear effects [15] (for small $v_\perp$), the latter mechanism leading, in the absence of coherent structures, to amplitudes of order $\epsilon^2$. In fact, another branch of solutions can be reached, where the hydrodynamic nonlinearity that leads to the finite time collapse of Eq. (5) can rapidly become important and dominate over the former effects. The singularity indeed induces a dynamics much faster than in the linear regime, thus permitting the amplitude to reach a level large enough for all previously subdominant nonlinearities to become relevant. No systematic approach being presently possible, it is reasonable to address these effects separately.

Models involving saturation processes based on the cooling of a population of trapped particles, but neglecting FLR corrections, were previously discussed [10,11]. These models mainly predict the formation of deep magnetic holes, while simulations of the VM equations, using either PIC codes [12,16] or an Eulerian description [17] display the formation of magnetic humps above threshold together with a phenomenon of bistability, associated with the existence of stable large-amplitude magnetic holes both below and above threshold.

It turns out that a saturation mechanism based on the simple effect of the variation of the local ion Larmor radius $\rho_L$ can be phenomenologically supplemented to our asymptotic model, making it consistent with VM simulations. For this purpose, writing that $\rho_L^2 \propto T_\perp/|B|^2 \propto 1/|B| \approx 1/B_\xi$ (because of the conservation of $\mu$), this variation can be modeled in Eq. (5) by replacing the term $\Delta_\perp U$ by $[1/(1 + \alpha U)]\Delta_\perp U$. Furthermore, in addition to the Laplacian which results from the leading order expansion of a nonlocal operator associated with FLR corrections [8], we retain the next order contribution in the form $(4/9)[\nu/(1 + \alpha U)^2]\Delta_\perp^2 U$. Equation (7) is modified accordingly. The coefficient $\nu$ is related to the size of the computational domain and is here taken equal to $10^{-2}$ in all the simulations, while $\alpha$ is a combination of the distance to threshold and of the value of $\beta_\perp$. Such terms that become important as small scales are formed, lead to a saturation of the singularity with a minimum of $U$ that

<figure><img src="image_4.png" /><figcaption>FIG. 2. Quasistatic solution of the saturated equation for $\sigma = 1$, $\nu = 0.01$, and $\alpha = 0.2$, obtained from a small-amplitude noise, by the development of the mirror instability.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 3. Same as Fig. 2, but with $\alpha = 1.2$.</figcaption></figure>

scales like $1/\alpha$, making the minimum of $B_z$ independent of $\alpha$. Retaining higher order terms does not drastically affect the value of this minimum, while keeping only the Laplacian would make it almost zero.

Figure 2 displays the profile of the solution obtained by time integration of the modified version of Eq. (7) with $\sigma = +1$, assuming $\alpha = 0.2$ and a random initial noise of maximal amplitude 0.01. We observe the suppression of the singularity and the formation of flat minima. For this value of $\alpha$, associated with small $\beta_\perp$ and/or proximity to threshold, we observe both holes and peaks, with a slight preference for holes, a tendency being more pronounced at smaller values of $\alpha$.

Increasing $\alpha$ significantly affects the geometry of the structure with, for $\alpha = 1.2$, a clear domination of magnetic humps (Fig. 3), as observed in PIC and Vlasov simulations. In the present simulation, after the arrest of the linear instability, six peaks were formed, which gradually disappear, through a coarsening phenomenon previously observed in hybrid PIC simulations [12].

Below threshold, when integrating the modified equation with $\alpha = 0.32$ and random initial conditions with a relatively large maximal amplitude equal to 1.5, the formation of spiky magnetic holes are observed, as exemplified in Fig. 4. For larger values of $\alpha$, such solutions no longer exist.

<figure><img src="image_6.png" /><figcaption>FIG. 4. Quasistatic solution of the saturated equation for $\sigma = -1$, $\nu = 0.01$, and $\alpha = 0.32$, obtained with large initial perturbations.</figcaption></figure>

In conclusion, we demonstrate that the hydrodynamic nonlinearity leads to a subcritical bifurcation for mirror modes near threshold, which provides a dynamic interpretation of the observation of magnetic holes below the instability threshold. Above onset, in addition to the classical quasilinear solutions, another branch of solutions with much larger amplitudes is obtained. Within the proposed model, the obtained structures (holes or humps) are consistent with the statistics of mirror structures in the terrestrial magnetosheath obtained from CLUSTER satellite data [6].

We thank F. Califano and P. Hellinger for communication of their results before publication, and to G. Belmont, V. Génot, and O. A. Pokhotelov for fruitful discussions. This work was performed in the framework of ISSI team “The effect of ULF turbulence and flow chaotization on plasma energy and mass transfers at the magnetopause.” The work of E. K. was supported by RFBR (Grant No. 06-01-00665) and by Landau-CNRS agreement during his visits to Nice; that of T. P. and P. L. S. by “Programme National Soleil Terre” of CNRS.

[1] E. A. Lucek, M. W. Dunlop, T. S. Horbury, A. Balogh, P. Brown, P. Cargill, C. Carr, K. H. Fornacon, E. Georgescu, and T. Oddy, Ann. Geophys. **19**, 1421 (2001).

[2] K. Sperveslage, F. M. Neubauer, K. Baumgärtel, and N. F. Ness, Nonlin. Processes Geophys. **7**, 191 (2000).

[3] A. A. Vedenov and R. Z. Sagdeev, in *Plasma Physics and Problem of Controlled Thermonuclear Reactions*, edited by M. A. Leontovich (Pergamon, NY, 1958), Vol. III, p. 332.

[4] P. Hellinger and P. Travnicek, J. Geophys. Res. **110**, A04210 (2005).

[5] G. Erdös and A. Balogh, J. Geophys. Res. **101**, 1 (1996).

[6] V. Génot, E. Budnik, C. Jacquey, J. Sauvaud, I. Dandouras, and E. Lucek, *AGU Fall Meeting* (AGU, Washington, DC, 2006), Abstract No. C1412+.

[7] T. Passot, V. Ruban, and P. L. Sulem, Phys. Plasmas **13**, 102310 (2006).

[8] O. A. Pokhotelov, R. Z. Sagdeev, M. A. Balikhin, and R. A. Treumann, J. Geophys. Res. **109**, A09213 (2004).

[9] S. P. Gary, J. Geophys. Res. **97**, 8519 (1992).

[10] M. G. Kivelson and D. S. Southwood, J. Geophys. Res. **101**, 17365 (1996).

[11] P. G. E. Pantellini, J. Geophys. Res. **103**, 4789 (1998).

[12] K. Baumgärtel, K. Sauer, and E. Dubinin, Geophys. Res. Lett. **30**, 1761 (2003).

[13] R. M. Kulsrud, in *Handbook of Plasma Physics*, edited by M. N. Rosenbluth and R. Z. Sagdeev, Basic Plasma Physics Vol. 1, edited by A. A. Galeev and R. N. Sudan (1983), p. 115.

[14] F. G. E. Pantellini, D. Burgers, and S. J. Schwartz, Adv. Space Res. **15**, 341 (1995).

[15] V. D. Shapiro and V. I. Shevchenko, Sov. Phys. JETP **18**, 1109 (1964).

[16] P. Hellinger (private communication).

[17] F. Califano (private communication).