1 Introduction
Despite being studied with great effort for approximately a century, natural and laboratory plasmas exhibit several complex phenomena that still need to be understood, mainly because of the strongly nonlinear interactions and the presence of kinetic effects. In this context, investigating plasma dynamics is decisive for understanding fundamental processes occurring in different systems, ranging from very-far astrophysical objects to the near-Earth environment and laboratory fusion devices. These systems routinely present a strongly nonlinear dynamics, which develops on a large range of spatial and time scales, including the ones associated with kinetic processes. In such systems, energy is typically injected at large fluid scales and cascades towards smaller scales, driving the system to cross three different physical regimes, ranging from fluid (magnetohydrodynamics, Hall-MHD) to ion kinetic and eventually to electron kinetic scales. This multi-scale physics is the direct consequence of the weak plasma collisionality, that characterizes solar-wind and astrophysical plasmas (Kulsrud Reference Kulsrud2005; Califano & Mangeney Reference Califano and Mangeney2008; Bruno & Carbone Reference Bruno and Carbone2016) as well as fusion device dynamics, where collisions can become effective at scales smaller than the electron kinetic scales (Falchetto et al. Reference Falchetto, Scott, Angelino, Bottino, Dannert, Grandgirard, Janhunen, Jenko, Jolliet and Kendl2008).
As a result, the plasma is allowed to freely access the entire phase space and to manifest dynamical states far from thermal equilibrium (Galeotti & Califano Reference Galeotti and Califano2005; Valentini et al. Reference Valentini, Carbone, Veltri and Mangeney2005; Marsch Reference Marsch2006; Franci et al. Reference Franci, Verdini, Matteini, Landi and Hellinger2015; Servidio et al. Reference Servidio, Valentini, Perrone, Greco, Califano, Matthaeus and Veltri2015, Reference Servidio, Chasapis, Matthaeus, Perrone, Valentini, Parashar, Veltri, Gershman, Russell and Giles2017; Cerri, Kunz & Califano Reference Cerri, Kunz and Califano2018; Franci et al. Reference Franci, Landi, Verdini, Matteini and Hellinger2018; Pezzi et al. Reference Pezzi, Servidio, Perrone, Valentini, Sorriso-Valvo, Greco, Matthaeus and Veltri2018; Sorriso-Valvo et al. Reference Sorriso-Valvo, Perrone, Pezzi, Valentini, Servidio, Zouganelis and Veltri2018b , Reference Sorriso-Valvo, Catapano, Retinò, Le Contel, Perrone, Roberts, Coburn, Panebianco, Valentini and Perri2019). As an example, we highlight here the fundamental role of the collisionless magnetic reconnection, that – even within a fluid theory framework – drives a strongly nonlinear dynamics (at both ion and electron scales), without collisions being relevant (Califano, Faganello & Pegoraro Reference Califano, Faganello and Pegoraro2007). Within this context, the Vlasov equation for each particle species, self-consistently coupled to Maxwell equations for the fields, provides a complete description of the system dynamics, although in some cases the role of weak collisions should be also considered (Navarro et al. Reference Navarro, Teaca, Told, Groselj, Crandall and Jenko2016; Pezzi, Valentini & Veltri Reference Pezzi, Valentini and Veltri2016; Pezzi et al. Reference Pezzi, Perrone, Servidio, Valentini, Sorriso-Valvo and Veltri2019). The Vlasov–Maxwell model is an nonlinear integro-differential set of equations in multi-dimensional phase space, whose analytic solutions are only available in a few simplified cases and in a reduced phase-space geometry. A numerical approach is therefore mandatory to describe the dynamics of collisionless magnetized plasmas in the fully nonlinear regime.
As of today, numerical simulations have provided significant insights into the plasma dynamics at proton and sub-proton spatial scales, where proton kinetic effects are dominant, while electrons can be approximated as an isothermal fluid (hybrid framework) (Valentini et al. Reference Valentini, Trávníček, Califano, Hellinger and Mangeney2007). In this range of scales, both particle-in-cell (PIC) and Eulerian hybrid codes have been extensively employed to investigate in detail a variety of physical phenomena such as, for instance, the development of the intermittent cascade of turbulent fluctuations (Parashar et al. Reference Parashar, Shay, Cassak and Matthaeus2009; Valentini, Califano & Veltri Reference Valentini, Califano and Veltri2010; Servidio et al. Reference Servidio, Valentini, Califano and Veltri2012; Franci et al. Reference Franci, Verdini, Matteini, Landi and Hellinger2015; Servidio et al. Reference Servidio, Valentini, Perrone, Greco, Califano, Matthaeus and Veltri2015; Franci et al. Reference Franci, Landi, Matteini, Verdini and Hellinger2016; Valentini et al. Reference Valentini, Perrone, Stabile, Pezzi, Servidio, Marco, Marcucci, Bruno, Lavraud and Keyser2016; Cerri & Califano Reference Cerri and Califano2017; Cerri, Servidio & Califano Reference Cerri, Servidio and Califano2017; Valentini et al. Reference Valentini, Vásconez, Pezzi, Servidio, Malara and Pucci2017; Pezzi et al. Reference Pezzi, Parashar, Servidio, Valentini, Vásconez, Yang, Malara, Matthaeus and Veltri2017a ,Reference Pezzi, Parashar, Servidio, Valentini, Vásconez, Yang, Malara, Matthaeus and Veltri b ; Cerri et al. Reference Cerri, Kunz and Califano2018; Franci et al. Reference Franci, Landi, Verdini, Matteini and Hellinger2018; Perrone et al. Reference Perrone, Passot, Laveder, Valentini, Sulem, Zouganelis, Veltri and Servidio2018; Sorriso-Valvo et al. Reference Sorriso-Valvo, Carbone, Perri, Greco, Marino and Bruno2018a ), the dynamo effect in turbulent plasmas (Rincon et al. Reference Rincon, Califano, Schekochihin and Valentini2016), the interaction of solar wind and the Earth’s magnetosphere at global scales (Kempf et al. Reference Kempf, Pokhotelov, von Alfthan, Vaivads, Palmroth and Koskinen2013; Palmroth et al. Reference Palmroth, Honkonen, Sandroos, Kempf, von Alfthan and Pokhotelov2013; Pokhotelov et al. Reference Pokhotelov, von Alfthan, Kempf, Vainio, Koskinen and Palmroth2013; Von Alfthan et al. Reference Von Alfthan, Pokhotelov, Kempf, Hoilijoki, Honkonen, Sandroos and Palmroth2014; Hoilijoki et al. Reference Hoilijoki, Palmroth, Walsh, Pfau-Kempf, von Alfthan, Ganse, Hannuksela and Vainio2016; Palmroth et al. Reference Palmroth, Ganse, Pfau-Kempf, Battarbee, Turc, Brito, Grand in, Hoilijoki, Sandroos and von Alfthan2018; Pfau-Kempf et al. Reference Pfau-Kempf, Battarbee, Ganse, Hoilijoki, Turc, von Alfthan, Vainio and Palmroth2018) and the dynamics of magnetic reconnection (Birn et al. Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001; Shay et al. Reference Shay, Drake, Rogers and Denton2001; Pritchett Reference Pritchett2008; Califano et al. Reference Califano, Cerri, Faganello, Laveder and Kunz2018). To reduce the computational cost of the simulation, reduced models – such as the gyro-kinetic (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006, Reference Howes, Dorland, Cowley, Hammett, Quataert, Schekochihin and Tatsuno2008a ,Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin b ; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Plunk, Quataert and Tatsuno2008; Tatsuno et al. Reference Tatsuno, Dorland, Schekochihin, Plunk, Barnes, Cowley and Howes2009; Howes et al. Reference Howes, Tenbarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011; TenBarge, Howes & Dorland Reference TenBarge, Howes and Dorland2013; Told et al. Reference Told, Jenko, TenBarge, Howes and Hammett2015; Howes Reference Howes2016) or the finite Larmor radius Landau fluid ones (Passot & Sulem Reference Passot and Sulem2007; Sulem & Passot Reference Sulem and Passot2015; Sulem et al. Reference Sulem, Passot, Laveder and Borgogno2016; Kobayashi et al. Reference Kobayashi, Sahraoui, Passot, Laveder, Sulem, Huang, Henri and Smets2017) – have been also widely adopted to describe plasma dynamics at kinetic scales.
Within the context of space plasmas, recent high-resolution observations conducted by the magnetospheric multi-scale (MMS) mission (Burch et al. Reference Burch, Moore, Torbert and Giles2016a ; Fuselier et al. Reference Fuselier, Lewis, Schiff, Ergun, Burch, Petrinec and Trattner2016) allowed us, for the first time, to investigate the plasma dynamics at the electron scale. The MMS mission focuses primarily on kinetic processes occurring in the electron diffusion region of magnetic reconnection (Burch et al. Reference Burch, Torbert, Phan, Chen, Moore, Ergun, Eastwood, Gershman, Cassak and Argall2016b ; Torbert et al. Reference Torbert, Burch, Giles, Gershman, Pollock, Dorelli, Avanov, Argall, Shuster and Strangeway2016, Reference Torbert, Burch, Phan, Hesse, Argall, Shuster, Ergun, Alm, Nakamura and Genestreti2018) and its unprecedented high-resolution observations confirm a very complex picture where several mechanisms can be at work in producing small-scale fluctuations (Le Contel et al. Reference Le Contel, Retinò, Breuillard, Mirioni, Robert, Chasapis, Lavraud, Chust, Rezeau and Wilder2016; Breuillard et al. Reference Breuillard, Matteini, Argall, Sahraoui, Andriopoulou, Contel, Retinò, Mirioni, Huang and Gershman2018; Chasapis et al. Reference Chasapis, Matthaeus, Parashar, Wan, Haggerty, Pollock, Giles, Paterson, Dorelli and Gershman2018). Magnetic reconnection often takes place within a turbulent environment where coherent structures – such as current sheets and X-points – naturally develop (Retinò et al. Reference Retinò, Sundkvist, Vaivads, Mozer, André and Owen2007; Servidio et al. Reference Servidio, Matthaeus, Shay, Cassak and Dmitruk2009, Reference Servidio, Matthaeus, Shay, Dmitruk, Cassak and Wan2010; Haggerty et al. Reference Haggerty, Parashar, Matthaeus, Shay, Yang, Wan, Wu and Servidio2017; Phan et al. Reference Phan, Eastwood, Shay, Drake, Sonnerup, Fujimoto, Cassak, Øieroset, Burch and Torbert2018). At the same time, plasma jets generated by magnetic reconnection can provide energy for sustaining the turbulence itself (Cerri et al. Reference Cerri, Servidio and Califano2017; Pucci et al. Reference Pucci, Servidio, Sorriso-Valvo, Olshevsky, Matthaeus, Malara, Goldman, Newman and Lapenta2017, Reference Pucci, Matthaeus, Chasapis, Servidio, Sorriso-Valvo, Olshevsky, Newman, Goldman and Lapenta2018). Reconnection is important for space and astrophysical plasmas as it is responsible for major plasma heating and particle acceleration in solar and stellar coronae, magnetars, accretion disks and astrophysical jets (Lyutikov Reference Lyutikov2003; Uzdensky Reference Uzdensky2011) as well as for tokamaks, being a major cause of loss of plasma confinement and plasma heating (Helander, Eriksson & Andersson Reference Helander, Eriksson and Andersson2002; Tanabe et al. Reference Tanabe, Yamada, Watanabe, Gi, Kadowaki, Inomoto, Imazawa, Gryaznevich, Michael and Crowley2015).
In order to properly combine and compare the experimental evidence at electron scales with theoretical investigations (Hesse et al.
Reference Hesse, Aunai, Birn, Cassak, Denton, Drake, Gombosi, Hoshino, Matthaeus and Sibeck2016), a huge numerical effort needs still to be made. To this end, only few numerical algorithms which retain both proton and electron kinetic physics are nowadays available. Most of them are PIC codes (Zeiler et al.
Reference Zeiler, Biskamp, Drake, Rogers, Shay and Scholer2002; Markidis, Lapenta & Rizwan-uddin Reference Markidis and Lapenta2010; Camporeale & Burgess Reference Camporeale and Burgess2011; Daughton et al.
Reference Daughton, Roytershteyn, Karimabadi, Yin, Albright, Bergen and Bowers2011; Karimabadi et al.
Reference Karimabadi, Roytershteyn, Wan, Matthaeus, Daughton, Wu, Shay, Loring, Borovsky and Leonardis2013; Leonardis et al.
Reference Leonardis, Chapman, Daughton, Roytershteyn and Karimabadi2013; Divin et al.
Reference Divin, Khotyaintsev, Vaivads, Andr, Markidis and Lapenta2015; Lapenta et al.
Reference Lapenta, Markidis, Goldman and Newman2015; Wan et al.
Reference Wan, Matthaeus, Roytershteyn, Karimabadi, Parashar, Wu and Shay2015; Grošelj et al.
Reference Grošelj, Cerri, Navarro, Willmott, Told, Loureiro, Califano and Jenko2017; Yang et al.
Reference Yang, Matthaeus, Parashar, Haggerty, Roytershteyn, Daughton, Wan, Shi and Chen2017; Parashar, Matthaeus & Shay Reference Parashar, Matthaeus and Shay2018; Shay et al.
Reference Shay, Haggerty, Matthaeus, Parashar, Wan and Wu2018; González et al.
Reference González, Parashar, Gomez, Matthaeus and Dmitruk2019; Lapenta et al.
Reference Lapenta, Pucci, Goldman and Newman2019), which capture the full dynamics (including electron scales) since their computational cost is smaller with respect to low-noise Eulerian (Vlasov) codes. However, at variance with noise-free Eulerian algorithms, PIC codes fail in providing a clean description of small-scale fluctuations (e.g. the electric-field behaviour around the X-point) and particle distribution functions in phase space, since they suffer from intrinsic statistical noise. Only very recently the first attempts to describe the plasma dynamics via Eulerian fully kinetic codes have become affordable, thanks to the improved supercomputer capabilities (Schmitz & Grauer Reference Schmitz and Grauer2006a
; Umeda, Togano & Ogino Reference Umeda, Togano and Ogino2009; Umeda et al.
Reference Umeda, Miwa, Matsumoto, Nakamura, Togano, Fukazawa and Shinohara2010; Delzanno Reference Delzanno2015; Tronci & Camporeale Reference Tronci and Camporeale2015; Umeda & Wada Reference Umeda and Wada2016; Ghizzo, Sarrat & Del Sarto Reference Ghizzo, Sarrat and Del Sarto2017; Umeda & Wada Reference Umeda and Wada2017; Juno et al.
Reference Juno, Hakim, TenBarge, Shi and Dorland2018; Roytershteyn et al.
Reference Roytershteyn, Boldyrev, Delzanno, Chen, Grošelj and Loureiro2019; Skoutnev et al.
Reference Skoutnev, Hakim, Juno and TenBarge2019). As stated above, Eulerian algorithms generally require a computational cost significantly large as compared to PIC codes. A way to reduce the computational cost of a fully kinetic Eulerian simulation consists of applying the so-called Darwin approximation (Kaufman & Rostler Reference Kaufman and Rostler1971; Birdsall & Langdon Reference Birdsall and Langdon2004; Schmitz & Grauer Reference Schmitz and Grauer2006a
,Reference Schmitz and Grauer
b
) to the Maxwell equations based on the expansion of the Maxwell system in the small parameter
$v^{2}/c^{2}$
(Mangeney et al.
Reference Mangeney, Califano, Cavazzoni and Travnicek2002) (
$v$
being the typical plasma bulk speed). Within this approximation, all wave modes (including those triggered by charge separation) are retained except for light waves (
$v_{\unicode[STIX]{x1D719}}\sim c$
,
$v_{\unicode[STIX]{x1D719}}$
being the wave phase speed); by doing so, the numerical stability condition for the time step can be significantly relaxed.
In the present work we present a newly developed fully kinetic Eulerian Vlasov–DArwin algorithm (ViDA) which integrates numerically the kinetic equations for a non-relativistic globally neutral plasma composed of protons and electrons. Equations are discretized on a fixed-in-time grid in phase space with periodic boundary conditions in the physical domain. ViDA originates from the hybrid Vlasov–Maxwell code (Valentini et al. Reference Valentini, Trávníček, Califano, Hellinger and Mangeney2007) (hereafter referred to as the HVM code) and has been extended specifically to include electron kinetic dynamics. The paper is organized as follows. In § 2 we revisit the Darwin approximation and describe the system of equations that is numerically integrated through ViDA. We discuss in detail the strategy of the numerical integration of the Vlasov equation for each species and we show that the Darwin version of the Maxwell equations can be written as a set of Helmholtz and Poisson-like equations, solvable through a spectral method. In the same section, we also provide a description of the algorithm design. Then, in § 3 we present the first results obtained through this algorithm, concerning the propagation of (i) electrostatic Langmuir waves, (ii) whistler waves and (iii) Alfvén waves. In § 4, we describe the onset of the electron Weibel instability which is a plasma instability driven by the presence of an electron temperature anisotropy (Weibel Reference Weibel1959). In § 5 we present preliminary results concerning one of the main potential applications of ViDA: the magnetic reconnection process at electron scales. Then, in § 6, we discuss the performance of the algorithm. Finally, we conclude and summarize in § 7.
2 The Vlasov–Darwin (VD) model
The Darwin approximation, that we briefly revisit in the current section, has been adopted to reduce the limitations on the time steps for numerical integration (Kaufman & Rostler Reference Kaufman and Rostler1971; Mangeney et al.
Reference Mangeney, Califano, Cavazzoni and Travnicek2002; Birdsall & Langdon Reference Birdsall and Langdon2004; Schmitz & Grauer Reference Schmitz and Grauer2006a
). Indeed, since Maxwell equations allow for the propagation of waves at the light speed
$c$
, the time step
$\unicode[STIX]{x0394}t$
of any explicit numerical scheme solving these equations would be limited by the Courant–Friedrichs–Lewy (CFL) condition,
$\unicode[STIX]{x0394}t\lesssim \unicode[STIX]{x0394}x/c$
(Peyret & Taylor Reference Peyret and Taylor1986). The Darwin approximation, by dropping the transverse displacement current term, rules out the transverse light waves (i.e. the fastest waves in the system that propagate at phase speed
$c$
) and significantly relaxes the CFL condition.
We consider a non-relativistic, collisionless, fully kinetic plasma composed of electrons and protons. The Vlasov–Darwin system of equations reads (in CGS units),
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\,f_{\unicode[STIX]{x1D6FC}}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}f_{\unicode[STIX]{x1D6FC}}+\frac{Z_{\unicode[STIX]{x1D6FC}}e}{m_{\unicode[STIX]{x1D6FC}}}\left(\boldsymbol{E}+\frac{\boldsymbol{v}}{c}\times \boldsymbol{B}\right)\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{v}}f_{\unicode[STIX]{x1D6FC}}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}_{L}=4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{c}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\times \boldsymbol{E}_{T}=-\frac{1}{c}\unicode[STIX]{x2202}_{t}\boldsymbol{B}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\times \boldsymbol{B}=\frac{1}{c}\unicode[STIX]{x2202}_{t}\boldsymbol{E}_{L}+\frac{4\unicode[STIX]{x03C0}}{c}\boldsymbol{j}, & \displaystyle\end{eqnarray}$$
where
$f_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v},t)$
is the distribution function (DF) of the
$\unicode[STIX]{x1D6FC}=p,e$
species,
$m_{\unicode[STIX]{x1D6FC}}$
and
$Z_{\unicode[STIX]{x1D6FC}}$
are respectively the mass and charge number of the
$\unicode[STIX]{x1D6FC}$
species and
$c$
is the light speed;
$\unicode[STIX]{x2202}_{t}$
,
$\unicode[STIX]{x1D735}$
and
$\unicode[STIX]{x1D735}_{\boldsymbol{v}}$
indicate the derivatives with respect to the time
$t$
, the spatial coordinates
$\boldsymbol{x}$
and the velocity coordinates
$\boldsymbol{v}$
, respectively;
$\boldsymbol{E}(\boldsymbol{x},t)=\boldsymbol{E}_{L}(\boldsymbol{x},t)+\boldsymbol{E}_{T}(\boldsymbol{x},t)$
and
$\boldsymbol{B}(\boldsymbol{x},t)$
are the electric and magnetic field, respectively. The electric field has been decomposed into a longitudinal (irrotational,
$\unicode[STIX]{x1D735}\times \boldsymbol{E}_{L}=0$
) and a transverse (solenoidal,
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{E}_{T}=0$
) component (Griffiths Reference Griffiths1962). According to the Darwin approximation, the transverse component of the displacement current has been neglected in Ampere’s law (2.5) (Birdsall & Langdon Reference Birdsall and Langdon2004; Schmitz & Grauer Reference Schmitz and Grauer2006a
). The Darwin model, that retains at least the longitudinal component of the displacement current, is generally closer to the full Maxwell system with respect to models where the displacement current is completely neglected (Valentini et al.
Reference Valentini, Trávníček, Califano, Hellinger and Mangeney2007; Tronci & Camporeale Reference Tronci and Camporeale2015).
The plasma charge density
$\unicode[STIX]{x1D70C}_{c}(\boldsymbol{x},t)$
and the current density
$\boldsymbol{j}(\boldsymbol{x},t)$
are defined through the first two velocity moments of the particle DFs,
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{c}=e\mathop{\sum }_{\unicode[STIX]{x1D6FC}}Z_{\unicode[STIX]{x1D6FC}}n_{\unicode[STIX]{x1D6FC}}=e\mathop{\sum }_{\unicode[STIX]{x1D6FC}}Z_{\unicode[STIX]{x1D6FC}}\int \text{d}\boldsymbol{v}\,f_{\unicode[STIX]{x1D6FC}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}=\mathop{\sum }_{\unicode[STIX]{x1D6FC}}\boldsymbol{j}_{\unicode[STIX]{x1D6FC}}=e\mathop{\sum }_{\unicode[STIX]{x1D6FC}}Z_{\unicode[STIX]{x1D6FC}}n_{\unicode[STIX]{x1D6FC}}\boldsymbol{V}_{\unicode[STIX]{x1D6FC}}=e\mathop{\sum }_{\unicode[STIX]{x1D6FC}}Z_{\unicode[STIX]{x1D6FC}}\int \text{d}\boldsymbol{v}\boldsymbol{v}\,f_{\unicode[STIX]{x1D6FC}}. & \displaystyle\end{eqnarray}$$
Equations (2.1)–(2.5) can be further simplified to obtain a set of Helmholtz-like equations of state without explicit time derivatives (see Birdsall & Langdon (Reference Birdsall and Langdon2004) and Schmitz & Grauer (Reference Schmitz and Grauer2006a
) for details). By normalizing equations using a characteristic length
$\bar{L}$
, time
$\bar{t}$
, velocity
$\bar{U}=\bar{L}/\bar{t}$
, mass
$\bar{m}$
and distribution function
$f_{\unicode[STIX]{x1D6FC},0}=\bar{n}/\bar{U}^{3}$
(with
$\bar{n}=\bar{L}^{-3}$
the equilibrium density), it is straightforward to get the dimensionless Vlasov–Darwin system of equations,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\,f_{\unicode[STIX]{x1D6FC}}+(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\,f_{\unicode[STIX]{x1D6FC}}+\frac{Z_{\unicode[STIX]{x1D6FC}}}{\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}}(\boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{v}}\,f_{\unicode[STIX]{x1D6FC}}=0,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D711}=-\unicode[STIX]{x1D701}^{2}\sum Z_{\unicode[STIX]{x1D6FC}}n_{\unicode[STIX]{x1D6FC}},\quad \boldsymbol{E}_{L}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D719},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\boldsymbol{B}=-\bar{u}^{2}\unicode[STIX]{x1D701}^{2}\unicode[STIX]{x1D735}\times \boldsymbol{j}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\hat{\boldsymbol{E}}_{T}-\bar{u}^{2}\unicode[STIX]{x1D701}^{2}\mathop{\sum }_{\unicode[STIX]{x1D6FC}}\frac{Z_{\unicode[STIX]{x1D6FC}}^{2}n_{\unicode[STIX]{x1D6FC},0}}{\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}}\hat{\boldsymbol{E}}_{T}=\bar{u}^{2}\unicode[STIX]{x1D701}^{2}\left[-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\mathop{\sum }_{\unicode[STIX]{x1D6FC}}Z_{\unicode[STIX]{x1D6FC}}\langle \boldsymbol{v}\boldsymbol{v}\rangle _{\unicode[STIX]{x1D6FC}}+\mathop{\sum }_{\unicode[STIX]{x1D6FC}}\frac{Z_{\unicode[STIX]{x1D6FC}}^{2}}{\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}}(n_{\unicode[STIX]{x1D6FC}}\boldsymbol{E}_{L}+\langle \boldsymbol{v}\rangle _{\unicode[STIX]{x1D6FC}}\times \boldsymbol{B})\right], & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D6E9}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{E}}_{T},\quad \boldsymbol{E}_{T}=\hat{\boldsymbol{E}}_{T}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6E9},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0, & & \displaystyle\end{eqnarray}$$
where
$\langle h\rangle _{\unicode[STIX]{x1D6FC}}=\int \text{d}\boldsymbol{v}f_{\unicode[STIX]{x1D6FC}}h$
. In (2.8)–(2.13), the electric and magnetic fields are normalized to
$\bar{E}=\bar{m}\bar{U}/e\bar{t}$
and
$\bar{B}=\bar{m}c/e\bar{t}$
, respectively. Note also that we set
$k_{B}=1$
. Non-dimensional parameters are
$\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}=m_{\unicode[STIX]{x1D6FC}}/\bar{m}$
,
$\bar{u}=\bar{U}/c$
and
$\unicode[STIX]{x1D701}=\bar{\unicode[STIX]{x1D714}}_{p}\bar{t}$
, with
$\bar{\unicode[STIX]{x1D714}}_{p}=\sqrt{4\unicode[STIX]{x03C0}e^{2}\bar{n}/\bar{m}}$
. Note that in (2.11) we have omitted a term
$\bar{u}^{2}\unicode[STIX]{x1D735}\unicode[STIX]{x2202}_{tt}\unicode[STIX]{x1D719}$
which could generate, in principle, an irrotational component, and we have introduced (2.12a,b
) to preserve the solenoidality of
$\boldsymbol{E}_{T}$
(Schmitz & Grauer Reference Schmitz and Grauer2006a
). The spatial dependence of
$n_{\unicode[STIX]{x1D6FC}}$
on the left-hand side of (2.11) has been neglected (
$n_{\unicode[STIX]{x1D6FC}}\simeq n_{\unicode[STIX]{x1D6FC},0}$
) to let the coefficients be constant (Valentini et al.
Reference Valentini, Trávníček, Califano, Hellinger and Mangeney2007).
2.1 Conservation properties
It is straightforward to verify that (2.8)–(2.13) satisfy mass
$\int \text{d}\boldsymbol{x}\,\text{d}\boldsymbol{v}\,f_{\unicode[STIX]{x1D6FC}}$
and entropy
$S_{\unicode[STIX]{x1D6FC}}=\int \text{d}\boldsymbol{x}\,\text{d}\boldsymbol{v}\,f_{\unicode[STIX]{x1D6FC}}\log f_{\unicode[STIX]{x1D6FC}}$
conservation. The energy conservation equation, obtained by multiplying (2.8) by
$m_{\unicode[STIX]{x1D6FC}}\boldsymbol{v}^{2}/2$
, integrating over the phase-space volume
$\int \text{d}\boldsymbol{x}\,\text{d}\boldsymbol{v}$
and summing over the species reads,
$$\begin{eqnarray}{\mathcal{E}}_{\text{kin}}+{\mathcal{E}}_{\text{th}}+{\mathcal{E}}_{\text{mag}}+{\mathcal{E}}_{\text{el}}=\text{const.},\end{eqnarray}$$
where the kinetic energy is
${\mathcal{E}}_{\text{kin}}=\sum _{\unicode[STIX]{x1D6FC}}(m_{\unicode[STIX]{x1D6FC}}/2)\int \text{d}\boldsymbol{x}n_{\unicode[STIX]{x1D6FC}}\boldsymbol{u}_{\unicode[STIX]{x1D6FC}}^{2}$
, the thermal energy is
${\mathcal{E}}_{\text{th}}=\sum _{\unicode[STIX]{x1D6FC}}(3/2)\int \text{d}\boldsymbol{x}n_{\unicode[STIX]{x1D6FC}}T_{\unicode[STIX]{x1D6FC}}$
, the magnetic energy is
${\mathcal{E}}_{\text{mag}}=\sum _{\unicode[STIX]{x1D6FC}}(\bar{m}/2\bar{u}^{2}\unicode[STIX]{x1D701}^{2})\int \text{d}\boldsymbol{x}\boldsymbol{B}^{2}$
and the electrostatic energy is
${\mathcal{E}}_{\text{el}}=\sum _{\unicode[STIX]{x1D6FC}}(\bar{m}/2\unicode[STIX]{x1D701}^{2})\int \text{d}\boldsymbol{x}\boldsymbol{E}_{L}^{2}$
. Note that the temperature of the
$\unicode[STIX]{x1D6FC}$
-species is defined as
$3n_{\unicode[STIX]{x1D6FC}}T_{\unicode[STIX]{x1D6FC}}/m_{\unicode[STIX]{x1D6FC}}=\int \text{d}\boldsymbol{v}(\boldsymbol{v}-\boldsymbol{u}_{\unicode[STIX]{x1D6FC}})^{2}f_{\unicode[STIX]{x1D6FC}}$
and, to get (2.14), we have used
$\int \text{d}\boldsymbol{x}\boldsymbol{w}_{T}\boldsymbol{\cdot }\boldsymbol{w}_{L}=0$
,
$\boldsymbol{w}_{T}$
and
$\boldsymbol{w}_{L}$
being a generic transverse and longitudinal vector, respectively. In each of the tests described in the present work we have checked the conservation of these quantities: their variation with respect to initial values is always smaller than the 1 %.
2.2 ViDA algorithm and code design
The Vlasov equation for each species is integrated numerically by employing the time splitting method first proposed by Cheng & Knorr (Reference Cheng and Knorr1976) in the electrostatic limit and later extended to the full electromagnetic case (Mangeney et al. Reference Mangeney, Califano, Cavazzoni and Travnicek2002). Darwin equations are solved through standard fast Fourier transform (FFT) algorithms.
In our case, the splitting algorithm for (2.8) reads as follows:
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{t}\,f_{\unicode[STIX]{x1D6FC}}+(\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\,f_{\unicode[STIX]{x1D6FC}}=0 & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\,f_{\unicode[STIX]{x1D6FC}}+\frac{Z_{\unicode[STIX]{x1D6FC}}}{\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}}(\boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{v}}\,f_{\unicode[STIX]{x1D6FC}}=0. & \displaystyle\end{eqnarray}$$
In the first equation
$\boldsymbol{v}$
is a parameter, while in the second equation
$\boldsymbol{x}$
,
$\boldsymbol{E}$
and
$\boldsymbol{B}$
are parameters. At time
$t$
, the solutions of (2.15)–(2.16) can be written as
$\unicode[STIX]{x1D6EC}_{\boldsymbol{x}}(t)F_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v})$
and
$\unicode[STIX]{x1D6EC}_{\boldsymbol{v}}F_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v})$
, respectively. In last expressions,
$\unicode[STIX]{x1D6EC}_{\boldsymbol{x}}$
and
$\unicode[STIX]{x1D6EC}_{\boldsymbol{v}}$
are the advection operators in physical and velocity space, whose explicit definition, based on the third-order van Leer scheme, can be found in Mangeney et al. (Reference Mangeney, Califano, Cavazzoni and Travnicek2002); while
$F_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v})=f_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v},t=0)$
is the initial condition. Note that
$\unicode[STIX]{x1D6EC}_{\boldsymbol{v}}$
depends also on the particle species
$\unicode[STIX]{x1D6FC}$
and, for the sake of simplicity, we avoid to explicitly report such dependence.
The splitting scheme is a symplectic, second order accurate in time (see Mangeney et al. (Reference Mangeney, Califano, Cavazzoni and Travnicek2002) where the stability condition of the advection operator is also discussed) and the numerical solution at
$t=t_{N}=N\unicode[STIX]{x0394}t$
is given by,
$$\begin{eqnarray}f_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v},t_{N})=\left[\unicode[STIX]{x1D6EC}_{\boldsymbol{x}}\left(\frac{\unicode[STIX]{x0394}t}{2}\right)\unicode[STIX]{x1D6EC}_{\boldsymbol{v}}(\unicode[STIX]{x0394}t)\unicode[STIX]{x1D6EC}_{\boldsymbol{x}}\left(\frac{\unicode[STIX]{x0394}t}{2}\right)\right]^{N}f_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v},0).\end{eqnarray}$$
At
$t=0$
, the distribution function
$F_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v})=f_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v},t=0)$
is first advected in physical space by half of a time step, obtaining
$\tilde{f}_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v},\unicode[STIX]{x0394}t/2)$
. Then, the following structure is executed:
(i) computing the moments of
$\tilde{f}_{\unicode[STIX]{x1D6FC}}$
and evaluating the electromagnetic fields
$\boldsymbol{E}_{L}$
,
$\boldsymbol{E}_{T}$
and
$\boldsymbol{B}$
, at
$t=\unicode[STIX]{x0394}t/2$
, through (2.9a,b
)–(2.13);(ii) performing a time-step advection in velocity space:
$\hat{f}_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x1D6EC}_{\boldsymbol{v}}(\unicode[STIX]{x0394}t)\tilde{f}_{\unicode[STIX]{x1D6FC}}$
;(iii) performing a time-step advection in physical space:
$\tilde{f}_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x1D6EC}_{\boldsymbol{x}}(\unicode[STIX]{x0394}t)\hat{f}_{\unicode[STIX]{x1D6FC}}$
.
This last structure is repeated in the algorithm, according to (2.17), in order to get the evolved distribution function at any time instant.
The phase-space domain is discretized as follows. The physical space
$D_{\boldsymbol{x}}=[0,L_{x}]$
$\times [0,L_{y}]\times [0,L_{z}]$
is discretized with
$N_{x}\times N_{y}\times N_{z}$
gridpoints and periodic boundary conditions are used. The velocity space
$D_{\boldsymbol{v},\unicode[STIX]{x1D6FC}}=[-v_{\unicode[STIX]{x1D6FC},x}^{\text{max}},v_{\unicode[STIX]{x1D6FC},x}^{\text{max}}]\times [v_{\unicode[STIX]{x1D6FC},y}^{\text{max}},v_{\unicode[STIX]{x1D6FC},y}^{\text{max}}]\times [-v_{\unicode[STIX]{x1D6FC},z}^{\text{max}},v_{\unicode[STIX]{x1D6FC},z}^{\text{max}}]$
is discretized by
$(2N_{\unicode[STIX]{x1D6FC},v_{x}}+1)\times (2N_{\unicode[STIX]{x1D6FC},v_{y}}+1)\times (2N_{\unicode[STIX]{x1D6FC},v_{z}}+1)$
grid points. Velocity-space boundary conditions impose
$f_{\unicode[STIX]{x1D6FC}}(|v_{i}|>v_{\unicode[STIX]{x1D6FC},i}^{\text{max}})=0$
(
$i=x,y,z$
). In order to ensure mass conservation,
$v_{\unicode[STIX]{x1D6FC},i}^{\text{max}}$
is typically set to be a large multiple of the thermal speed
$v_{\text{th},\unicode[STIX]{x1D6FC}}=\sqrt{T_{\unicode[STIX]{x1D6FC}}/m_{\unicode[STIX]{x1D6FC}}}$
.
The ViDA algorithm has been designed in such a way that the user can select (i) different normalizations of the model equations, (ii) the possibility of setting motionless protons and (iii) different dimensionalities of the physical-space domain (one, two or three dimensions), the velocity-space domain being always three-dimensional (
$3V$
). Within ViDA spatial vectors always have three components and can be function of one, two or three spatial variables, depending on the physical-space dimensionality. Since Darwin equations are a set of Helmholtz-like equations, initial perturbations have to be introduced through the particle DFs (and their moments): this represents a difference with respect to standard codes where also magnetic perturbations can be introduced. A check on the solenoidality of
$\boldsymbol{B}$
and
$\boldsymbol{E}_{T}$
is also implemented at each time step.
The computational effort necessary to solve VD equations is significant and a massive parallelization, based on both MPI and OpenMP paradigms is implemented. The MPI paradigm, first introduced for the DF by Mangeney et al. (Reference Mangeney, Califano, Cavazzoni and Travnicek2002), is adopted to parallelize the physical-space computational domain for both particle DFs (and their moments) and electromagnetic field. Hence each MPI thread accesses a finite portion of phase space, composed by a sub-portion of physical space and by the whole velocity space. Within each MPI thread, the OpenMP directives are adopted to parallelize the velocity-space cycles. The parallelization of the electromagnetic field is a new feature recently introduced in the HVM code in Cerri & Califano (Reference Cerri and Califano2017) and it is essential to perform high-resolution Eulerian Vlasov simulations, in particular in three dimensions. Preliminary tests on performance and scalability are reported in § 6.
2.3 Normalizations of the Vlasov–Darwin equations
In order to normalize (2.8)–(2.13), three possible choices have been implemented in ViDA:
(i) Electrostatic normalization. Characteristic quantities are: length
$\bar{L}=\unicode[STIX]{x1D706}_{D,e}$
, time
$\bar{t}=\unicode[STIX]{x1D714}_{p,e}^{-1}$
, velocity
$\bar{U}=v_{\text{th},e}$
and mass
$\bar{m}=m_{e}$
. Here
$\unicode[STIX]{x1D706}_{D,e}=\sqrt{T_{e}/4\unicode[STIX]{x03C0}n_{e}e^{2}}$
,
$\unicode[STIX]{x1D714}_{p,e}=(\sqrt{4\unicode[STIX]{x03C0}n_{e}e^{2}/m_{e}})^{-1}$
,
$v_{\text{th},e}=\sqrt{T_{e}/m_{e}}=\unicode[STIX]{x1D706}_{D,e}\unicode[STIX]{x1D714}_{p,e}$
and
$m_{e}$
are the electron Debye length, the electron plasma frequency, the electron thermal speed and the electron mass, respectively. This normalization is appropriate for describing phenomena occurring at electron scales, such as the propagation of electrostatic plasma waves.(ii) Electromagnetic normalization. Characteristic quantities are: length
$\bar{L}=d_{e}$
, time
$\bar{t}=\unicode[STIX]{x1D714}_{p,e}^{-1}$
, velocity
$\bar{U}=c$
and mass
$\bar{m}=m_{e}$
, where
$d_{e}=c/\unicode[STIX]{x1D714}_{p,e}$
is the electron skin depth. This normalization can be adopted for describing electromagnetic phenomena, where both protons and electrons are involved, such as magnetic reconnection and plasma turbulence at kinetic scales.(iii) Hybrid normalization. Characteristic quantities are: length
$\bar{L}=d_{p}$
, time
$\bar{t}=\unicode[STIX]{x1D6FA}_{c,p}^{-1}$
, velocity
$\bar{U}=v_{A}$
and mass
$\bar{m}=m_{p}$
. In previous expressions
$\unicode[STIX]{x1D6FA}_{cp}=eB_{0}/m_{p}c$
,
$v_{A}=B_{0}/\sqrt{4\unicode[STIX]{x03C0}n_{p}m_{p}}$
,
$d_{p}=v_{A}/\unicode[STIX]{x1D6FA}_{cp}$
and
$m_{p}$
are the proton cyclotron frequency, the proton Alfvén speed, the proton skin depth and the proton mass, respectively. This normalization is useful for investigating the turbulent cascade in the sub-proton range, where electron physics starts to play a role.
These three normalizations can be adopted to describe, in a more natural way (i.e. characteristic scales close to unity), phenomena where electrostatic, electromagnetic or proton inertial effects dominate, respectively.
3 Numerical tests of ViDA
In this section we report the results of several tests performed to evaluate the capabilities of ViDA in describing basic collisionless plasma physics dynamics. The proper behaviour and reliability of the code have been tested against the propagation of Langmuir waves, in both linear and nonlinear regimes, whistler waves and Alfvén waves.
3.1 Propagation and damping of Langmuir waves
For these tests we adopted the electrostatic normalization. We discuss results of simulations performed with motionless protons in 1D–3V phase-space configuration, where Langmuir waves propagate along the
$x$
direction. Physical and velocity space have been discretized with
$N_{x}=128$
and
$N_{e,v_{x}}=50$
,
$N_{e,v_{y}}=N_{e,v_{z}}=15$
grid points, respectively. In the case of mobile protons (
$m_{p}/m_{e}=1836$
and
$T_{p}/T_{e}=1$
), the propagation of Langmuir waves has been reproduced with lower phase-space resolution, but the results are quantitatively similar to those discussed in the following. We have also separately tested the propagation of Langmuir waves along
$y$
and
$z$
directions by carrying out 2D–3V and 3D–3V runs.
The initial equilibrium is given by an electron Maxwellian distribution spatially homogeneous. The plasma is unmagnetized, the initial electron temperature is
$T_{e}=1$
(in scaled units). At
$t=0$
, the electron number density is perturbed through a sinusoidal perturbation
$\unicode[STIX]{x1D6FF}n_{e}/n_{e,0}=A\sin (kx)$
,
$A=10^{-4}$
and
$k=k_{1}=2\unicode[STIX]{x03C0}/L$
being the amplitude and the wavenumber, respectively. The box length is
$L=18\unicode[STIX]{x1D706}_{D,e}$
(
$k_{1}=0.35$
) and
$v_{e,i}^{\text{max}}=5v_{\text{th},e}$
(
$i=x,y,z$
). The system evolution is reproduced up to a maximum time
$t_{\text{max}}=100\unicode[STIX]{x1D714}_{p,e}^{-1}$
, while the numerical recurrence time is
$t_{\text{rec}}=2\unicode[STIX]{x03C0}/k\unicode[STIX]{x0394}v\simeq 180\unicode[STIX]{x1D714}_{p,e}^{-1}$
(Cheng & Knorr Reference Cheng and Knorr1976; Galeotti, Califano & Pegoraro Reference Galeotti, Califano and Pegoraro2006; Pezzi et al.
Reference Pezzi, Valentini and Veltri2016).
Figure 1(a) shows the time evolution of the amplitude of the
$k_{x}=k_{1}$
Fourier component of the electric field
$|E_{x}|(k_{1},t)$
, in a semi-logarithmic plot. The electric field undergoes Landau damping (Landau Reference Landau1946); the observed damping rate shows a very good agreement with the theoretical prediction
$\unicode[STIX]{x1D6FE}_{L}=-3.37\times 10^{-2}\unicode[STIX]{x1D714}_{p,e}$
(red-dashed line), evaluated through a numerical linear solver for the roots of the electrostatic Vlasov dielectric function. In figure 1(b) we report the resonant curve, obtained by Fourier transforming the electric signal in space and time; we plot the spectral electric energy
$|E_{x}|^{2}(k_{1},\unicode[STIX]{x1D714}_{R})$
as a function of the pulsation
$\unicode[STIX]{x1D714}_{R}$
. As expected, the resonant curve displays a well-defined frequency peak in correspondence of a value of the pulsation
$\unicode[STIX]{x1D714}_{R}=1.22\unicode[STIX]{x1D714}_{p,e}$
. In figure 1(b), the vertical red-dashed line represents the value of the theoretical resonant pulsation
$\unicode[STIX]{x1D714}_{R,\text{th}}$
obtained through the linear solver, while the two vertical red-dot-dashed lines indicate the interval of uncertainty of the numerical code, due to the time discretization
$\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R}=2\unicode[STIX]{x03C0}/t_{\text{max}}\simeq 0.063\unicode[STIX]{x1D714}_{p,e}$
. Again, numerical results are in very good agreement with theoretical predictions.
(a) Time evolution of
$|E_{x}|(k_{1},t)$
. The red dashed line indicates the predicted Landau damping rate
$\unicode[STIX]{x1D6FE}_{L}$
. (b) Energy peak
$|E_{x}|^{2}(k_{1},\unicode[STIX]{x1D714}_{R})$
as a function of the pulsation
$\unicode[STIX]{x1D714}_{R}$
. The red dashed line indicates the theoretical wave frequency
$\unicode[STIX]{x1D714}_{R,\text{th}}$
, while red dot-dashed lines show the
$\unicode[STIX]{x1D714}_{R}$
– resolution, i.e.
$\unicode[STIX]{x1D714}_{R,\text{th}}\pm \unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R}/2$
. The theoretical expectations for the Langmuir wave damping and the pulsation have been obtained with a numerical solver of the linear dispersion relation.
In order to show the dependence of the real
$\unicode[STIX]{x1D714}_{R}$
and imaginary
$\unicode[STIX]{x1D714}_{I}$
parts of the frequency as a function of the wavenumber, we have performed an additional
$1D$
–
$3V$
run, in which the initial perturbation is a superposition of the first six wavenumbers
$k_{x}=[k_{1},6k_{1}]$
, where
$k_{1}=2\unicode[STIX]{x03C0}/L$
(
$L=2\unicode[STIX]{x03C0}10\unicode[STIX]{x1D706}_{D,e}$
); the other parameters are the same as in the previous run. To avoid numerical recurrence, phase space has been discretized with
$N_{x}=128$
,
$N_{e,v_{x}}=100$
and
$N_{e,v_{y}}=N_{e,v_{z}}=15$
. Figure 2 reports by stars the dependence of
$\unicode[STIX]{x1D714}_{R}$
(a) and
$\unicode[STIX]{x1D714}_{I}$
(b), in units of
$\unicode[STIX]{x1D714}_{pe}$
, as a function of the wavenumber
$k\unicode[STIX]{x1D706}_{D,e}$
. A very good agreement with theoretical expectations (red-dashed curves) is recovered for both real and imaginary parts of the complex frequency.
Pulsation
$\unicode[STIX]{x1D714}_{R}$
(a) and damping rate
$\unicode[STIX]{x1D6FE}_{I}$
(b), in units of
$\unicode[STIX]{x1D714}_{p,e}$
, from the simulation (black dots) and from the linear numerical solver (red dashed line) as a function of the wavenumber
$k\unicode[STIX]{x1D706}_{D,e}$
.
We conclude this section by focusing on the nonlinear regime of the Langmuir wave dynamics (see, for example, Brunetti, Califano & Pegoraro (Reference Brunetti, Califano and Pegoraro2000) and references therein). We have performed ten different runs, varying the amplitude of the initial density perturbation in the range
$A=[8\times 10^{-3},8\times 10^{-2}]$
. In this case, phase space has been discretized with
$N_{x}=128$
,
$N_{e,v_{x}}=150$
and
$N_{e,v_{y}}=N_{e,v_{z}}=15$
, while
$t_{\text{max}}=400\unicode[STIX]{x1D714}_{pe}^{-1}$
. The box length is
$L_{x}=18\unicode[STIX]{x1D706}_{D,e}$
, while
$v_{e,i}^{\text{max}}=5v_{\text{th},e}$
(
$i=x,y,z$
). As reported in figure 3(a), the time evolution of the electric-field Fourier component shows an early linear damping phase (Landau Reference Landau1946), until particle trapping arrests the damping and produces oscillations of the signal envelope (O’Neil Reference O’Neil1965). At large times, a phase-space vortex is observed in the electron DF in the vicinity of the wave phase speed, as reported in figure 3(b). As it has been shown in O’Neil (Reference O’Neil1965), the nonlinear trapping time
$\unicode[STIX]{x1D70F}_{\text{nl}}$
depends on the saturation amplitude
$A_{1}$
of the electric oscillations as
$\unicode[STIX]{x1D70F}_{\text{nl}}\sim A_{1}^{-1/2}$
. For each of the ten simulations, we evaluated
$A_{1}$
and
$\unicode[STIX]{x1D70F}_{\text{nl}}$
at the time of the first peak of the electric envelope oscillations. Figure 3(c) shows in log–log plot
$\unicode[STIX]{x1D70F}_{\text{nl}}$
as a function of
$A_{1}$
(stars), compared to the theoretical expectation (red line), showing a very nice agreement.
(a) Time evolution of
$|E_{x}(k_{1},t)|$
for the simulation with
$A=8\times 10^{-2}$
. (b) Contour plot of
$\tilde{f}_{e}(x,v_{x})=\int \text{d}v_{y}\,\text{d}v_{z}\,f_{e}(x,\boldsymbol{v})$
in the plan
$x$
–
$v_{x}$
. (c) Nonlinear time
$\unicode[STIX]{x1D70F}_{\text{nl}}$
as a function of the first peak amplitude
$A_{1}$
. The red line reports the predicted scaling
${\sim}-0.5$
(the result of the linear fit is
${\sim}-0.48$
).
3.2 Propagation of whistler waves
To reproduce the propagation of whistler waves at electron scales, the electromagnetic normalization has been employed. Protons are assumed just as a fixed neutralizing background. Again, we have verified that the ViDA code behaves exactly in the same manner in the three spatial directions. Hence, we discuss here the result of a
$1D$
–
$3V$
run, where
$\boldsymbol{B}_{0}=B_{0}\boldsymbol{e}_{x}$
(
$B_{0}=1$
) and protons are not fixed. The box length is
$L_{x}=2\unicode[STIX]{x03C0}10d_{e}$
, while
$v_{e,i}^{\text{max}}=10v_{\text{th},e}$
and
$v_{p,i}^{\text{max}}=7v_{\text{th},p}$
(
$i=x,y,z$
). Note that increasing the value of
$v_{p(e)}^{\text{max}}$
has been necessary to ensure mass conservation. The phase space has been discretized with
$N_{x}=128$
,
$N_{e,v_{x}}=N_{e,v_{y}}=N_{e,v_{z}}=50$
and
$N_{p,v_{x}}=N_{p,v_{y}}=N_{p,v_{z}}=35$
. We also set
$m_{p}/m_{e}=1836$
,
$T_{e}/T_{p}=1$
,
$\bar{u}=v_{\text{th},e}/c=10^{-3}$
and
$\unicode[STIX]{x1D701}=1$
.
The equilibrium is composed of Maxwellian velocity distributions for both protons and electrons and homogeneous density. The initial equilibrium is then perturbed with the following electron bulk-speed perturbations:
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}u_{e,y}=A\sin (kx), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}u_{e,z}=A\cos (kx), & \displaystyle\end{eqnarray}$$
where
$A=10^{-3}$
and
$k=k_{1}=2\unicode[STIX]{x03C0}/L_{x}$
.
By solving the Darwin equations, these perturbations generate a current density and then magnetic fluctuations. Figure 4 reports the time evolution of
$|B_{y}|(k_{1},t)$
(a) and the frequency peak of the spectral magnetic energy
$|B_{y}|^{2}(k_{1},\unicode[STIX]{x1D714}_{R})$
as a function of the pulsation
$\unicode[STIX]{x1D714}_{R}$
(b). The magnetic field clearly oscillates at the correct frequency
$\unicode[STIX]{x1D714}_{R,\text{th}}=0.91\unicode[STIX]{x1D714}_{p,e}$
, that can be evaluated from the linear dispersion relation for whistler waves (obtained by assuming motionless protons and cold electrons):
$\unicode[STIX]{x1D714}_{R,\text{th}}(k)=B_{0}k^{2}/(1+k^{2})$
(Krall & Trivelpiece Reference Krall and Trivelpiece1973). Note that, for the considered wavenumber, a negligible damping of whistler waves is expected. In figure 4(a), the red-dashed line represents the value of the resonant pulsation from the above expression for
$\unicode[STIX]{x1D714}_{R,\text{th}}$
, while the two vertical red dot-dashed lines indicate the interval of numerical uncertainty
$\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R}=2\unicode[STIX]{x03C0}/t_{\text{max}}\simeq 0.063\unicode[STIX]{x1D714}_{p,e}$
.
(a) Time evolution of
$|B_{y}|(k_{1},t)$
. (b) Spectral magnetic energy
$|B_{y}|^{2}(k_{1},\unicode[STIX]{x1D714}_{R})$
as a function of the pulsation
$\unicode[STIX]{x1D714}_{R}$
. The red dashed line indicates the theoretical wave frequency
$=\unicode[STIX]{x1D714}_{R,\text{th}}$
, while the red dot-dashed lines show the
$\unicode[STIX]{x1D714}_{R}$
-resolution
$\unicode[STIX]{x1D714}_{R,\text{th}}\pm \unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R}/2$
.
3.3 Propagation of Alfvén waves
Here we show numerical results concerning the propagation of Alfvén waves along a background magnetic field. The adopted normalization for these tests is the hybrid one. We perform a
$1D$
–
$3V$
run, where
$\boldsymbol{B}_{0}=B_{0}\boldsymbol{e}_{x}$
and
$B_{0}=1$
. The box length is
$L_{x}=2\unicode[STIX]{x03C0}50d_{p}$
, while
$v_{e,i}^{\text{max}}=5v_{\text{th},e}$
and
$v_{p,i}^{\text{max}}=5v_{\text{th},p}$
(
$i=x,y,z$
). The phase space has been discretized with
$N_{x}=32$
,
$N_{e,v_{x}}=N_{e,v_{y}}=N_{e,v_{z}}=25$
,
$N_{p,v_{x}}=40$
and
$N_{p,v_{y}}=N_{p,v_{z}}=35$
gridpoints. The mass ratio has been artificially set to
$m_{p}/m_{e}=25$
, thus avoiding extremely small time steps, while
$T_{e}/T_{p}=1$
,
$\bar{u}=v_{A}/c=10^{-3}$
,
$\unicode[STIX]{x1D701}=c/v_{A}=10^{3}$
and
$\unicode[STIX]{x1D6FD}_{p}=2v_{\text{th},p}^{2}/v_{A}^{2}=1$
.
The initial equilibrium, composed of spatially homogeneous Maxwellian protons and electrons, has been perturbed with the following proton bulk-speed perturbations:
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}u_{p,y}=A\sin (kx), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}u_{p,z}=A\cos (kx), & \displaystyle\end{eqnarray}$$
where
$A=10^{-4}$
and
$k=k_{1}=2\unicode[STIX]{x03C0}/L_{x}$
. Figure 5 shows the time evolution of
$|B_{z}|(k_{1},t)$
(a) and the magnetic spectral energy
$|B_{z}|^{2}(k_{1},\unicode[STIX]{x1D714}_{R})$
as a function of
$\unicode[STIX]{x1D714}_{R}$
(b). The recovered resonant peak is in agreement with the theoretical pulsation, evaluated through a fully kinetic linear solver of the dispersion relation (Camporeale & Burgess Reference Camporeale and Burgess2017). Moreover, no Landau damping is observed, since it occurs at much smaller scales (Barnes Reference Barnes1966; Vàsconez et al.
Reference Vàsconez, Valentini, Camporeale and Veltri2014; Camporeale & Burgess Reference Camporeale and Burgess2017). This test represents the first attempt towards a general description of Alfvén waves, where electron physics is also taken into account. Since including electron physics is currently too computational demanding, we plan to continue the investigation in a separate, future work.
(a) Time evolution of
$|B_{z}|(k_{1},t)$
. (b) Magnetic spectral energy
$|B_{z}|^{2}(k_{1},\unicode[STIX]{x1D714}_{R})$
as a function of
$\unicode[STIX]{x1D714}_{R}$
. The red solid line indicates the theoretical wave frequency
$\unicode[STIX]{x1D714}_{R,\text{th}}$
, while red dashed lines show the
$\unicode[STIX]{x1D714}_{R}$
-resolution
$\unicode[STIX]{x1D714}_{R,\text{th}}\pm \unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{R}/2$
.
4 Temperature anisotropy driven instability: electron Weibel instability
Another class of interesting numerical tests, which can be performed to point out the reliability of the ViDA code, concerns the onset of micro-instabilities, such as whistler, mirror and Weibel instabilities driven by a temperature anisotropy (Weibel Reference Weibel1959; Gary & Karimabadi Reference Gary and Karimabadi2006; Califano et al. Reference Califano, Hellinger, Kuznetsov, Passot, Sulem and Trvnek2008; Palodhi, Califano & Pegoraro Reference Palodhi, Califano and Pegoraro2009, Reference Palodhi, Califano and Pegoraro2010; Chen & Chacón Reference Chen and Chacón2014; Camporeale & Zimbardo Reference Camporeale and Zimbardo2015; Chen & Chacón Reference Chen and Chacón2015).
Here we focus on the development of the electron Weibel instability, that produces electromagnetic fluctuations transverse to the wavevector
$\boldsymbol{k}$
. The most suitable normalization to perform this analysis is the electromagnetic one. In particular we discuss the results of a 1D–3V run with
$\boldsymbol{k}=k\hat{\boldsymbol{e}}_{x}$
, although we have verified that instability is triggered in the same way also in the different phase-space configuration (two and three dimensions). The mass ratio is
$m_{p}/m_{e}=100$
, while
$T_{e}/T_{p}=0.01$
. Electrons are initialized with a bi-Maxwellian distribution function, with thermal speeds
$v_{\text{th},e,x}=2.5\times 10^{-2}c$
and
$v_{\text{th},e,y}=v_{\text{th},e,z}=4\times 10^{-2}c$
, this giving a temperature anisotropy
$A=T_{y(z)}/T_{x}=2.56$
. Protons have a Maxwellian velocity distribution at
$t=0$
, with a thermal speed
$v_{\text{th},p}=v_{\text{th},e,x}$
and homogeneous density. However, in this case, protons mainly act just as a neutralizing background, not being involved in the dynamics during the linear stage (i.e. during the instability development). No background magnetic field has been introduced. Physical space, whose length is
$L_{x}=32d_{e}$
, has been discretized with
$N_{x}=64$
gridpoints. Velocity space has been discretized with
$51^{3}$
gridpoints for both protons and electrons and
$v_{e(p)}^{\text{max}}=5v_{\text{th},e(p)}$
in each velocity directions.
The initial equilibrium has been perturbed through a sinusoidal, transverse perturbation, imposed on the electron bulk speed,
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}u_{e,y}=A\sin (kx), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}u_{e,z}=A\cos (kx), & \displaystyle\end{eqnarray}$$
where
$A=2\times 10^{-5}$
and
$k=k_{1}=2\unicode[STIX]{x03C0}/L_{x}$
. Such bulk-speed perturbations produce a current density, which in turn generates magnetic fluctuations. Figure 6 reports the time evolution of the magnetic spectral energy density
$W_{B}(k_{1},t)=(|B_{y}|^{2}(k_{1},t)+|B_{z}|^{2}(k_{1},t))/2$
. The red-dashed line indicates the expected linear instability growth rate
$\unicode[STIX]{x1D714}_{I}^{\text{th}}\simeq 4\times 10^{-3}\unicode[STIX]{x1D714}_{p,e}$
, evaluated through a linear solver for the roots of the kinetic electromagnetic dielectric function. In the early stage of the simulation,
$W_{B}$
increases exponentially with a growth rate in very good agreement with the expected one. Then, oscillations saturate at a nearly constant value in the nonlinear regime of evolution (Chen & Chacón Reference Chen and Chacón2014).
Time evolution of
$W_{B}(k_{1})=(|B_{y}|^{2}(k_{1})+|B_{z}|^{2}(k_{1}))/2$
. The red dashed line indicates the linear instability growth
$\unicode[STIX]{x1D714}_{I}^{\text{th}}\simeq 4\times 10^{-3}\unicode[STIX]{x1D714}_{p,e}$
, calculated with a linear solver of the fully kinetic dispersion relation (Camporeale & Burgess Reference Camporeale and Burgess2017).
5 Dynamics of magnetic reconnection
In this section we present the results of a magnetic reconnection simulation. Generally speaking, Vlasov simulations of magnetic reconnection represent a strong numerical challenge because of the huge memory and CPU time required by Eulerian algorithms. This approach, if successful, would certainly provide a crucial contribution to the understanding of the magnetic reconnection process especially at electron scales, thanks to the fact that Eulerian algorithms allow for an almost noise-free description of fields and particle DFs. A noise-free description is crucial to properly understand e.g. which electromagnetic fluctuations contribute to the reconnection electric field in the form of anomalous resistivity and how distribution functions are modified leading to electron heating.
We have performed a
$2D$
–
$3V$
symmetric magnetic reconnection simulation. Reconnection is symmetric when the values of magnetic field and density are equal on the two opposite sides of the current sheet. The initial condition of our simulation is the one adopted in the geospace environmental modeling (GEM) challenge (Birn et al.
Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001), in order to allow for a direct comparison to previous studies (Birn et al.
Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001; Schmitz & Grauer Reference Schmitz and Grauer2006b
). For this reason, we have also chosen the hybrid normalization (see § 2.3).
The equilibrium is set by adapting the Harris equilibrium (Harris Reference Harris1962) to the periodic boundary conditions in the spatial domain. In particular, the component of the magnetic field
$B_{x}(y)$
corresponding to the double current sheet profile reads,
$$\begin{eqnarray}B_{x}(y)=B_{0}\left[\tanh \left(\frac{y-L_{y}/2}{L_{1}}\right)-\tanh \left(\frac{y}{L_{2}}\right)-\tanh \left(\frac{y-L_{y}}{L_{2}}\right)\right].\end{eqnarray}$$
This profile is characterized by the presence of two gradients (the current sheets) varying as an hyperbolic tangent and located at
$y=L_{y}/2$
and
$y=0$
(and so at
$y=L_{y}$
) where
$L_{y}$
is the length of the spatial domain in the
$y$
direction. The first hyperbolic tangent is the one defined in Harris (Reference Harris1962) and
$L_{1}$
is the corresponding current sheet thickness. The second and third hyperbolic tangents in (5.1) have been included to satisfy the spatial periodicity; the value of
$L_{2}$
is taken sufficiently large compared to
$L_{1}$
to slow down the development of reconnection in the second current sheet with respect to the main one. The electron and ion temperature are set as uniform at the initial time and the density
$n(y)$
is defined in order to satisfy pressure balance. Then, from (2.5) and considering
$\unicode[STIX]{x2202}_{t}\boldsymbol{E}_{L}=0$
at the initial time, we get the initial current density
$\boldsymbol{j}=(0,0,j_{z}(y))$
.
Following the prescriptions of the Harris equilibrium we get, in normalized units,
$$\begin{eqnarray}\displaystyle & \displaystyle n_{0}(T_{e}+T_{p})=\frac{B_{0}^{2}}{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{u_{e,0}}{T_{e}}=-\frac{u_{p,0}}{T_{p}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{j_{z}(y)}{n(y)}\equiv u_{0}=u_{p,0}-u_{e,0}. & \displaystyle\end{eqnarray}$$
Equation (5.3) corresponds to the no charge separation condition of the Harris equilibrium so that quasi-neutrality is imposed,
$n_{e}(y)=n_{p}(y)=n(y)$
. In other words, the electric field is zero at the initial time. Moreover, from (5.3)–(5.4) we have,
$$\begin{eqnarray}\displaystyle & \displaystyle u_{e,0}=-\frac{u_{0}}{1+T_{p}/T_{e}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{p,0}=\frac{u_{0}}{1+T_{e}/T_{p}}. & \displaystyle\end{eqnarray}$$
It is worth pointing out that this is not an exact Vlasov kinetic equilibrium. In particular, it differs from the equilibrium presented by Harris since in this simulation the spatial domain is periodic in the varying
$y$
-direction. On the other hand, the initial configuration is in force balance and we have checked that the initial equilibrium is not significantly affected by, for example, ballistic effects within the time scale of reconnection considered here. As for the GEM challenge (Birn et al.
Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001), fluctuations are superposed on the initial magnetic field in order to obtain a single magnetic island at the centre of the space domain at the initial time. In particular,
$\unicode[STIX]{x1D6FF}\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}\times \hat{z}$
and
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}(x,y)=\unicode[STIX]{x1D713}_{0}\cos (2\unicode[STIX]{x03C0}x/L_{x})\cos (2\unicode[STIX]{x03C0}y/L_{y}),\end{eqnarray}$$
where, as already stated,
$L_{x}$
and
$L_{y}$
are the lengths of the spatial domain in
$x$
and
$y$
directions, respectively. According to GEM challenge, in scaled units,
$\unicode[STIX]{x1D713}_{0}$
is set to 0.1.
By using the relation
$\unicode[STIX]{x1D6FF}\boldsymbol{B}(x,y)=\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}(x,y)\times \hat{z}$
and (2.10), we derive the expression for the current density fluctuations
$\unicode[STIX]{x1D6FF}\boldsymbol{j}(x,y)$
consistent with
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D713}(x,y)$
. In particular, it is possible to define
$\unicode[STIX]{x1D6FF}\boldsymbol{j}(x,y)=(0,0,\unicode[STIX]{x1D6FF}j_{z}(x,y))$
. Finally, the initial electron and proton distribution functions are shifted Maxwellian distributions with drift velocities along the
$z$
direction and temperature
$T_{e}$
and
$T_{p}$
.
The phase space has been discretized with
$N_{x}\times N_{y}=512\times 512$
gridpoints in the spatial domain,
$N_{e,v_{x}}\times N_{e,v_{y}}\times N_{e,v_{z}}=41\times 41\times 81$
gridpoints in the velocity domain for electrons and
$N_{p,v_{x}}\times N_{p,v_{y}}\times N_{p,v_{z}}=31\times 31\times 31$
gridpoints in the velocity domain for protons. We also set
$v_{e}^{\text{max}}=5v_{\text{th},e}$
and
$v_{p}^{\text{max}}=5v_{\text{th},p}$
, where the normalized
$v_{\text{th},p}$
is set to
$1$
. Other simulation parameters are
$L_{1}=0.5d_{p}$
,
$L_{2}=2.5d_{p}$
,
$m_{p}/m_{e}=25$
,
$n_{\infty }=0.2$
,
$T_{e}/T_{p}=0.2$
,
$L_{x}=L_{y}=25.6d_{p}$
. Also, we set
$B_{0}=1$
and
$n_{0}=1$
. All parameters are chosen to be as close as possible to the simulation parameters listed in Birn et al. (Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001).
In figure 7 we show the evolution of the reconnected flux given by the difference
$\unicode[STIX]{x0394}\unicode[STIX]{x1D713}$
between the magnetic flux
$\unicode[STIX]{x1D713}$
evaluated at the X point and at the O point. Accordingly to the initial perturbation, the X-point and the O point are initially located at
$(L_{x}/2,L_{y}/2)$
and
$(0,L_{y}/2)$
and their locations do not significantly change throughout the simulation run. The behaviour of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D713}$
is very similar to the evolution of the reconnected flux in Birn et al. (Reference Birn, Drake, Shay, Rogers, Denton, Hesse, Kuznetsova, Ma, Bhattacharjee and Otto2001). Reconnection evolves with a reconnected flux that remains close to zero until
$t\sim 15\unicode[STIX]{x1D6FA}_{c,p}^{-1}$
, when a sharp increase is observed. Then, the reconnection rate stays relatively constant until the reconnected flux begins to saturate at
$t\sim 30\unicode[STIX]{x1D6FA}_{c,p}^{-1}$
.
Time evolution of the reconnected magnetic flux
$\unicode[STIX]{x0394}\unicode[STIX]{x1D713}$
.
Contour plots of
$B_{z}$
(a); out-of-plane electron current density
$j_{e,z}$
(b); out-of-plane proton current density
$j_{p,z}$
(c); and electron number density
$n_{e}$
(d). The quantities are shown at the time
$t^{\ast }=15.27\unicode[STIX]{x1D6FA}_{c,p}^{-1}$
. At that time
$\unicode[STIX]{x0394}\unicode[STIX]{x1D713}=1.18$
. All the panels are zoomed in
$y$
in the interval
$[6d_{p},19d_{p}]$
.
In figure 8 we show the contour plots of the out-of-plane magnetic field
$B_{z}$
(a), of the electron current density in the
$z$
-direction
$j_{e,z}$
(b), of the proton current density in the
$z$
-direction
$j_{p,z}$
(c) and of the electron number density
$n_{e}$
(d). In each panel, the contour lines of the magnetic flux
$\unicode[STIX]{x1D713}$
are superposed;
$B_{z}$
exhibits the typical Hall quadrupolar pattern usually observed during symmetric magnetic reconnection. This magnetic signature indicates that the ions are demagnetized while the electrons are still frozen to the magnetic field. The difference in their dynamics produces the out-of-plane
$B_{z}$
(Mandt, Denton & Drake Reference Mandt, Denton and Drake1994; Uzdensky & Kulsrud Reference Uzdensky and Kulsrud2006). The quadrupolar structure that we find is analogous to the one obtained with other kinetic codes, both Eulerian (Schmitz & Grauer Reference Schmitz and Grauer2006b
, see figure 2) and Lagrangian (Pritchett Reference Pritchett2001, see Plate 1(b)). We note that the
$j_{p,z}$
pattern closely follows the density pattern (
$n_{e}\simeq n_{p}$
) so that
$j_{p,z}$
is depleted at the X point while it reaches its maximum value within the magnetic island. On the other hand,
$j_{e,z}$
is enhanced at the X point and the region of strong
$j_{e,z}$
is elongated along
$x$
. Away from the X point,
$j_{e,z}$
splits into two branches that identify the separatrices, as it was also observed by Shay et al. (Reference Shay, Drake, Rogers and Denton2001) (see their figure 6d). The electron current at the X line has a thickness comparable to
$d_{p}$
which corresponds to
$5d_{e}$
. The maximum value of the normalized
$B_{z}$
is 0.09 while the maximum values of
$j_{p,z}$
and
$j_{e,z}$
are 0.39 and 1.49, respectively. These values are overall slightly smaller than the values found in a similar Vlasov–Darwin simulation described in Schmitz & Grauer (Reference Schmitz and Grauer2006b
).
In figure 9 we show the reconnection outflow of protons and electrons at
$t^{\ast }=18.13\unicode[STIX]{x1D6FA}_{c,p}$
. In particular, we note that at
$x=3d_{p}$
(panel (a)), corresponding to a distance of
$9.8d_{p}$
from the X-point located at
$L_{x}/2=12.8d_{p}$
, the electron velocity is characterized by two peaks corresponding to the separatrices, while the proton velocity is concentrated in the centre of the outflow region and it reaches lower values, as expected. The presence of the two peaks is consistent with the
$j_{e,z}$
pattern shown in figure 8(b). Figure 9(b) shows the same quantities of figure 9(a) at a distance of
$2.3d_{p}$
from the X-point where the outflow is still developing and we note that
$u_{e,x}$
is rather similar in shape and value to
$u_{p,x}$
.
(a) The
$x$
component of the electron velocity
$u_{e}$
(black line) and of the proton velocity
$u_{p}$
(red line) at
$x^{\ast }=3.00d_{p}$
and (b) at
$x^{\ast }=10.52d_{p}$
. The quantities are shown at the time
$t^{\ast }=18.13\unicode[STIX]{x1D6FA}_{c,p}^{-1}$
.
6 Performance test on the ViDA code
In this section, we present preliminary performance tests of ViDA implemented on the Marconi-KNL cluster at the CINECA supercomputing centre (Casalecchio di Reno (BO), Italy). The Marconi cluster is equipped with 3600 Lenovo Adam Pass nodes, interconnected through the Intel OmniPath network and each one composed of a 1KNL processor (68 cores, 1.40 GHz), formally 96 GB of RAM (effective 83 GB) and 16 GB of MCDRAM. The tests have been performed on a simple equilibrium configuration (Maxwellian DFs with no perturbations). We remark however that this choice does not affect the code performance.
(a) Strong scaling from 128 to 1024 cores on Marconi KNL using 2 OpenMP threads per Task MPI on the 2.5-dimensional set-up. (b) Strong scaling using 64 MPI task and from 2 to 16 OpenMP Threads per task on the 2.5-dimensional set-up. (c) Strong scaling from 128 to 1024 cores on Marconi KNL using 2 OpenMP threads per Task MPI on the three-dimensional set-up. (d) Weak scaling from 64 to 1024 cores on Marconi KNL using 2 OpenMP threads per Task MPI on the three-dimensional set-up.
ViDA numerically integrates the VD equations in a six-dimensional phase space (
$3D$
–
$3V$
:
$x,y,z,v_{x},v_{y},v_{z}$
). Only the three-dimensional physical space is parallelized using cubic cells (squared in the two-dimensional configuration). For implementing these tests, we have chosen
$51^{3}$
velocity gridpoints for each particle DF (protons and electrons), which represent a typical value adopted in production runs, and we have performed approximately 100 time steps per test (note that changing the step number does not act on the code scalability). Note that two DFs are advanced in time through the ViDA algorithm, this limiting the number of spatial grid points per single processor and hence increasing the number of communications.
As a first step, we have analysed the parallel performance in the
$2D$
–
$3V$
configuration, adopting a physical-space grid with
$1024\times 1024$
points. This set-up requires approximately 6 TB of RAM, corresponding to, at least, 64 Marconi-KNL nodes. We have performed a strong scaling test by reducing the number of MPI tasks per node from 8 to 1 and maintaining the same number of two OpenMP threads per task. Results are presented in figure 10(a): the parallel efficiency scales efficiently up to 512 cores. As the number of cores increases, the efficiency is degraded owing to the more significant weight of MPI communications. This is mainly due to the huge memory request of the code combined with the Marconi KNL architecture. The code performance would strongly benefit from using a computer architecture with a larger RAM and a lower number of cores per node. We have also verified that the performance degradation cannot be handled by using an OpenMP strategy, as shown in figure 10(b), since the code performance is not affected by increasing the number of threads per node. In summary, within the current parallelization, the best performance is achieved with 32 MPI threads and 2 OpenMP tasks per node on a KNL system.
A slightly better performance is achieved using a full three-dimensional configuration with
$128\times 128\times 64$
grid points in the spatial domain. The strong scaling from 128 to 1024 cores is shown in figure 10(c). A weak scaling test has been also performed by multiplying the number of spatial points and the number of cores (nodes) by the same factor. From the results, presented in figure 10(d), it can be appreciated that the parallel efficiency is high only up to several hundred cores.
These preliminary tests show a reasonable parallel efficiency on KNL architecture, at least up to some hundreds of cores. We are presently working to increase the code efficiency, in particular optimizing the communication pattern of the ViDA algorithm. We note that these results in part depend on the employed architecture. It is worth finally highlighting that, for instance, by using a Skylake machine with 192 GB and 48 cores per node, we would be able to increase by a factor of 3 the number of spatial gridpoints per node, thus increasing the parallel efficiency of the code, as the number of communications would strongly decrease.
Concerning the computational costs, the ViDA code is approximately twice as computationally expensive as the HVM code (Valentini et al. Reference Valentini, Trávníček, Califano, Hellinger and Mangeney2007), which has been recently used for three-dimensional simulations of plasma turbulence (see for instance Cerri et al. (Reference Cerri, Kunz and Califano2018)). More specifically, the reconnection run presented here – which is the most expensive test in this paper in terms of required computational resources – has a cost of slightly less than 0.1 Mh on a Marconi supercomputer using 16 nodes and 512 MPI processes. On the other hand, with ViDA being a code for a new piece of physics, it is difficult to foresee the exact cost of a three-dimensional reconnection (or turbulence) run because the numerical and physical parameters, as well as the duration of the run, can vary significantly with respect to the standard ones used with the HVM code. Based on the experience with the HVM code, we may suggest that a high-resolution three-dimensional run of magnetic reconnection focusing on the electron physics would take from a few to a few tens of Mh. Such a significant allocation of computing time can be obtained, for example, in the framework of a PRACE project.
7 Conclusion
In this paper we have presented a fully kinetic code (ViDA) based on a Vlasov–Darwin algorithm, where only light waves are excluded in order to relax the constraint on the time step advancement. This approach is particularly suited for the investigation of the kinetic dynamics from sub-ion scales down to the electron kinetic scales
$d_{e}$
and to the Debye length
$\unicode[STIX]{x1D706}_{D}$
. As typically the case for space plasmas, but often also in the laboratory, inter-particle collisions are not described, since collisional effects are assumed to be insignificant.
ViDA has been tested against several waves modes, in particular Alfvén, whistlers and plasma waves. The development of the Weibel instability and reconnection, both in a regime where the main dynamics is driven by the electrons, has been also reproduced. These tests represent typical regimes of interest for studying the electron-scale kinetic dynamics representing today a strong computational challenge and a frontier problem for the understanding of electron plasma physics.
One of the main future objectives of ViDA will be the study of the structure and dynamics of the electron diffusion region, including the role of anomalous resistivity in Ohm’s law and the mechanisms of electron heating, which are among the main targets of satellite MMS data analysis (Torbert et al. Reference Torbert, Burch, Giles, Gershman, Pollock, Dorelli, Avanov, Argall, Shuster and Strangeway2016; Genestreti et al. Reference Genestreti, Varsani, Burch, Cassak, Torbert, Nakamura, Ergun, Phan, Toledo-Redondo and Hesse2018; Cozzani et al. Reference Cozzani, Retinò, Califano, Alexandrova, Le Contel, Khotyaintsev, Vaivads, Fu, Catapano and Breuillard2019). Last but not least, we will make use of the ViDA code for the study of the plasma turbulent dynamics focusing on the problem of the ‘dissipative’ scale, of primary interest in the context of the solar-wind turbulent heating at kinetic scales (Vaivads et al. Reference Vaivads, Retinò, Soucek, Khotyaintsev, Valentini, Escoubet, Alexandrova, André, Bale and Balikhin2016).
Acknowledgements
This project (F.C., A.R., F.V.) has received funding from the European Unions Horizon 2020 research and innovation programme under grant agreement no. 776262 (AIDA, www.aida-space.eu). O.P. is partly supported by the International Space Science Institute (ISSI) in the framework of the International Team 405 entitled ‘Current Sheets, Turbulence, Structures and Particle Acceleration in the Heliosphere’. E.C. is partially supported by NWO Vidi grant 639.072.716. Numerical simulations discussed here have been performed on the Marconi cluster at CINECA (Italy), under the ISCRA initiative (IsC53_MRVDS and IsB16_VDMMS) and under the INAF-CINECA initiative (INA17_C2A16 and INA17_C2A16b).
