ViDA: a Vlasov-DArwin solver for plasma physics at electron scales

TL;DR

ViDA is a specialized Vlasov-Darwin solver designed for high-accuracy, small-scale plasma physics simulations, enabling detailed studies of electron-scale phenomena like magnetic reconnection and turbulence.

Contribution

The paper introduces ViDA, a novel numerical code that efficiently simulates electron-scale plasma dynamics using a Vlasov-Darwin approach with low noise and high accuracy.

Findings

Successfully reproduces linear and nonlinear wave propagation

Captures magnetic reconnection physics accurately

Demonstrates efficient parallelization on high-performance computing clusters

Abstract

We present a Vlasov-DArwin numerical code (ViDA) specifically designed to address plasma physics problems, where small-scale high accuracy is requested even during the non linear regime to guarantee a clean description of the plasma dynamics at fine spatial scales. The algorithm provides a low-noise description of proton and electron kinetic dynamics, by splitting in time the multi-advection Vlasov equation in phase space. Maxwell equations for the electric and magnetic fields are reorganized according to Darwin approximation to remove light waves. Several numerical tests show that ViDA successfully reproduces the propagation of linear and nonlinear waves and captures the physics of magnetic reconnection. We also discuss preliminary tests of the parallelization algorithm efficiency, performed at CINECA on the Marconi-KNL cluster. ViDA will allow to run Eulerian simulations of a non-relativistic fully-kinetic collisionless plasma and it is expected to provide relevant insights on important problems of plasma astrophysics such as, for instance, the development of the turbulent cascade at electron scales and the structure and dynamics of electron-scale magnetic reconnection, such as the electron diffusion region.

Equations

$${\partial }_t{f}_{\alpha }+\mathbf{v}\cdot \nabla f_{\alpha }+\frac{Z_{\alpha }e}{m_{\alpha }}\left(\mathbf{E}+\frac{\mathbf{v}}{c}\times \mathbf{B}\right)\cdot {\nabla }_{\mathbf{v}}{f}_{\alpha }=0$$

$$\nabla \cdot {\mathbf{E}}_L=4\pi {\rho }_c$$

$$\nabla \cdot \mathbf{B}=0$$

$$\nabla \times {\mathbf{E}}_T=-\frac{1}{c}{\partial }_t\mathbf{B}$$

$$\nabla \times \mathbf{B}=\frac{1}{c}{\partial }_t{\mathbf{E}}_L+\frac{4\pi }{c}\mathbf{j}$$

$${\rho }_c=e\sum _{\alpha }{Z}_{\alpha }{n}_{\alpha }=e\sum _{\alpha }{Z}_{\alpha }\int d\mathbf{v}{f}_{\alpha }$$

$$\mathbf{j}=\sum _{\alpha }{\mathbf{j}}_{\alpha }=e\sum _{\alpha }{Z}_{\alpha }{n}_{\alpha }{\mathbf{V}}_{\alpha }=e\sum _{\alpha }{Z}_{\alpha }\int d\mathbf{v}\mathbf{v}{f}_{\alpha }$$

$${\partial }_t{f}_{\alpha }+\left(\mathbf{v}\cdot \nabla \right)f_{\alpha }+\frac{Z_{\alpha }}{{\mu }_{\alpha }}\left(\mathbf{E}+\mathbf{v}\times \mathbf{B}\right)\cdot {\nabla }_{\mathbf{v}}{f}_{\alpha }=0$$

$${\nabla }^2\phi =-{\zeta }^2\sum Z_{\alpha }{n}_{\alpha }{\mathbf{E}}_L=-\nabla \varphi$$

$${\nabla }^2\mathbf{B}=-{\bar{u}}^2{\zeta }^2\nabla \times \mathbf{j}$$

$${\nabla }^2{\hat{\mathbf{E}}}_T-{\bar{u}}^2{\zeta }^2\sum _{\alpha }\frac{Z_{\alpha }^2{n}_{\alpha ,0}}{{\mu }_{\alpha }}{\hat{\mathbf{E}}}_T={\bar{u}}^2{\zeta }^2[-\nabla \cdot \sum _{\alpha }{Z}_{\alpha }⟨\mathbf{v}\mathbf{v}{⟩}_{\alpha }+$$

$$+\sum _{\alpha }\frac{Z_{\alpha }^2}{{\mu }_{\alpha }}\left(n_{\alpha }{\mathbf{E}}_L+⟨\mathbf{v}{⟩}_{\alpha }\times \mathbf{B}\right)]$$

$${\nabla }^2\Theta =\nabla \cdot {\hat{\mathbf{E}}}_T{\mathbf{E}}_T={\hat{\mathbf{E}}}_T-\nabla \Theta$$

$$\nabla \cdot \mathbf{B}=0$$

$${\mathcal{E}}_{kin}+{\mathcal{E}}_{th}+{\mathcal{E}}_{mag}+{\mathcal{E}}_{el}=const$$

$${\partial }_t{f}_{\alpha }+\left(\mathbf{v}\cdot \nabla \right)f_{\alpha }=0$$

$${\partial }_t{f}_{\alpha }+\frac{Z_{\alpha }}{{\mu }_{\alpha }}\left(\mathbf{E}+\mathbf{v}\times \mathbf{B}\right)\cdot {\nabla }_{\mathbf{v}}{f}_{\alpha }=0$$

$$f_{\alpha }(\mathbf{x},\mathbf{v},t_N)={\left[{\Lambda }_{\mathbf{x}}\left(\frac{\Delta t}{2}\right){\Lambda }_{\mathbf{v}}\left(\Delta t\right){\Lambda }_{\mathbf{x}}\left(\frac{\Delta t}{2}\right)\right]}^N{f}_{\alpha }(\mathbf{x},\mathbf{v},0)$$

$$\delta u_{e,y}=A\sin (kx)$$

$$\delta u_{e,z}=A\cos (kx)$$

$$\delta u_{p,y}=A\sin (kx)$$

$$\delta u_{p,z}=A\cos (kx)$$

$$\delta u_{e,y}=Asin(kx)$$

$$\delta u_{e,z}=Acos(kx)$$

$$B_{x}(y)=B_{0}\biggl{[}\tanh\biggl{(}\frac{y-L_{y}/2}{L_{1}}\biggr{)}-\tanh\biggl{(}\frac{y}{L_{2}}\biggr{)}-\tanh\biggl{(}\frac{y-L_{y}}{L_{2}}\biggr{)}\biggr{]}.$$

$$n_0(T_e+T_p)=\frac{B_0^2}{2}$$

$$\frac{u_{e,0}}{T_e}=-\frac{u_{p,0}}{T_p}$$

$$\frac{j_z(y)}{n(y)}\equiv u_0=u_{p,0}-u_{e,0}$$

$$u_{e,0}=-\frac{u_0}{1+T_p/T_e}$$

$$u_{p,0}=\frac{u_0}{1+T_e/T_p}$$

$$\delta \psi (x,y)={\psi }_0\cos (2\pi x/L_x)\cos (2\pi y/L_y)$$

Full text

1 Introduction

Despite being studied with great efforts for about a century, natural and laboratory plasmas exhibit several complex phenomena that still need to be understood, mainly because of the strongly non-linear 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 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 (MHD, 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, 2005; Califano & Mangeney, 2008; Bruno & Carbone, 2016) as well as fusion devices dynamics, where collisions can become effective at scales smaller than the electron kinetic scales (Falchetto et al., 2008).

As a result, the plasma is allowed to freely access the entire phase space and to manifest dynamical states far from thermal equilibrium (Valentini et al., 2005; Galeotti & Califano, 2005; Marsch, 2006; Franci et al., 2015; Servidio et al., 2015, 2017; Franci et al., 2018; Sorriso-Valvo et al., 2018b; Pezzi et al., 2018; Cerri et al., 2018; Sorriso-Valvo et al., 2019). 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 to be relevant (Califano et al., 2007). Within this context, the Vlasov equation for each particle species, self-consistently coupled to Maxwell equations for 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., 2016; Pezzi et al., 2016, 2019). 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 reduced phase-space geometry. A numerical approach is therefore mandatory to describe the dynamics of collisionless magnetized plasmas in fully nonlinear regime.

As of today, numerical simulations have provided significant insights on 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., 2007). 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., 2009; Valentini et al., 2010; Servidio et al., 2012; Franci et al., 2015; Servidio et al., 2015; Franci et al., 2016; Valentini et al., 2016; Cerri et al., 2017; Cerri & Califano, 2017; Valentini, F. et al., 2017; Pezzi et al., 2017b, a; Franci et al., 2018; Perrone et al., 2018; Cerri et al., 2018; Sorriso-Valvo et al., 2018a), the dynamo effect in turbulent plasmas (Rincon et al., 2016), the interaction of solar wind and Earth’s magnetosphere at global scales (Palmroth et al., 2013; Kempf et al., 2013; Pokhotelov et al., 2013; Von Alfthan et al., 2014; Hoilijoki et al., 2016; Pfau-Kempf et al., 2018; Palmroth et al., 2018) and the dynamics of magnetic reconnection (Birn et al., 2001; Shay et al., 2001; Pritchett, 2008; Califano et al., 2018). To reduce the computational cost of the simulation reduced models –such as the gyro-kinetics (Howes et al., 2006, 2008a, 2008b; Schekochihin et al., 2008; Tatsuno et al., 2009; Howes et al., 2011; TenBarge et al., 2013; Told et al., 2015; Howes, 2016) or the finite Larmor radius Landau fluid ones (Passot & Sulem, 2007; Sulem & Passot, 2015; Sulem et al., 2016; Kobayashi et al., 2017)– 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., 2016; Fuselier et al., 2016) allowed, for the first time, to investigate the plasma dynamics at electron scale. The MMS mission focuses primarily on kinetic processes occurring in the electron diffusion region of magnetic reconnection  (Burch et al., 2016; Torbert et al., 2016, 2018) 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., 2016; Breuillard et al., 2018; Chasapis et al., 2018). Magnetic reconnection often takes place within a turbulent environment where coherent structures –such as current sheets and X-points– naturally develop (Retinò et al., 2007; Servidio et al., 2009, 2010; Haggerty et al., 2017; Phan et al., 2018). At the same time, plasma jets generated by magnetic reconnection can provide energy for sustaining the turbulence itself (Pucci et al., 2017; Cerri et al., 2017; Pucci et al., 2018). 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, 2003; Uzdensky, 2011) as well as for tokamaks, being a major cause of loss of plasma confinement and plasma heating (Helander et al., 2002; Tanabe et al., 2015).

In order to properly combine and compare the experimental evidences at electron scales with theoretical investigations (Hesse et al., 2016), a huge numerical effort needs to be done yet. To this end, only few numerical algorithms which retain both proton and electron kinetic physics are nowadays available. Most of them are PIC codes (Markidis et al., 2010; Daughton et al., 2011; Zeiler et al., 2002; Camporeale & Burgess, 2011; Karimabadi et al., 2013; Leonardis et al., 2013; Divin et al., 2015; Wan et al., 2015; Lapenta et al., 2015; Grošelj et al., 2017; Yang et al., 2017; Parashar et al., 2018; Shay et al., 2018; Lapenta et al., 2019; González et al., 2019), 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 behavior 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 plasma dynamics via Eulerian fully-kinetic codes have became affordable, thanks to the improved supercomputer capabilities (Schmitz & Grauer, 2006a; Umeda et al., 2009, 2010; Tronci & Camporeale, 2015; Delzanno, 2015; Umeda & Wada, 2016, 2017; Ghizzo et al., 2017; Juno et al., 2018; Roytershteyn et al., 2019; Skoutnev et al., 2019). 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 in applying the so-called Darwin approximation (Kaufman & Rostler, 1971; Birdsall & Langdon, 2004; Schmitz & Grauer, 2006a, b) to the Maxwell equations based on the expansion of the Maxwell system in the small parameter $v^{2}/c^{2}$ (Mangeney et al., 2002) ($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_{\phi}\sim c$, $v_{\phi}$ being the wave phase speed); by doing so, the numerical stability condition for the timestep 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., 2007) (hereafter referred as HVM code) and has been extended specifically to include electron kinetic dynamics. The paper is organized as follows. In Sect. 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 Sect. 3 we present first results obtained through this algorithm, concerning the propagation of i) electrostatic Langmuir waves, ii) whistler waves and iii) Alfvén waves. In Sect. 4, we describe the onset of the electron Weibel instability which is a plasma instability driven by the presence of a electron temperature anisotropy (Weibel, 1959). In Sect. 5 we present preliminary results concerning one of the main potential applications of ViDA: the magnetic reconnection process at electron scales. Then, in Sect. 6, we discuss the performances of the algorithm. Finally, we conclude and summarize in Sect. 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 timesteps for numerical integration (Kaufman & Rostler, 1971; Mangeney et al., 2002; Birdsall & Langdon, 2004; Schmitz & Grauer, 2006a). Indeed, since Maxwell equations allow for the propagation of waves at the light speed $c$, the timestep $\Delta t$ of any explicit numerical scheme solving these equations would be limited by the Courant-Friedrichs-Lewy (CFL) condition, $\Delta t\lesssim\Delta x/c$ (Peyret & Taylor, 1986). 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 by electrons and protons. The Vlasov-Darwin system of equations reads (in CGS units):

[TABLE]

where $f_{\alpha}({\bm{x}},{\bm{v}},t)$ is the distribution function (DF) of the $\alpha=p,e$ species, $m_{\alpha}$ and $Z_{\alpha}$ are respectively the mass and charge number of the $\alpha$ species and $c$ is the light speed. $\partial_{t}$, $\nabla$ and $\nabla_{{\bm{v}}}$ indicate the derivatives with respect to the time $t$, the spatial coordinates ${\bm{x}}$ and the velocity coordinates ${\bm{v}}$, respectively. ${\bm{E}}({\bm{x}},t)={\bm{E}}_{L}({\bm{x}},t)+{\bm{E}}_{T}({\bm{x}},t)$ and ${\bm{B}}({\bm{x}},t)$ are the electric and magnetic field, respectively. The electric field has been decomposed into a longitudinal (irrotational, $\nabla\times{\bm{E}}_{L}=0$) and a transverse (solenoidal, $\nabla\cdot{\bm{E}}_{T}=0$) component (Griffiths, 1962). According to the Darwin approximation, the transverse component of the displacement current has been neglected in the Ampere’s law [Eq. (5)] (Birdsall & Langdon, 2004; Schmitz & Grauer, 2006a). 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., 2007; Tronci & Camporeale, 2015).

The plasma charge density $\rho_{c}({\bm{x}},t)$ and the current density ${\bm{j}}({\bm{x}},t)$ are defined through the first two velocity moments of the particle DFs:

[TABLE]

Equations (1–5) can be further simplified to obtain a set of Helmholtz-like equations of state without explicit time derivatives (see Birdsall & Langdon (2004); Schmitz & Grauer (2006a) 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_{\alpha,0}=\bar{n}/\bar{U}^{3}$ (being $\bar{n}=\bar{L}^{-3}$ the equilibrium density), it is straightforward to get the dimensionless Vlasov-Darwin system of equations:

[TABLE]

where $\langle h\rangle_{\alpha}=\int d\bm{v}f_{\alpha}h$. In Eqs. (8–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 $\mu_{\alpha}=m_{\alpha}/\bar{m}$, $\bar{u}=\bar{U}/c$ and $\zeta=\bar{\omega}_{p}\bar{t}$, being $\bar{\omega}_{p}=\sqrt{4\pi e^{2}\bar{n}/\bar{m}}$. Note that in Eq. (11) we have omitted a term $~{}\bar{u}^{2}\nabla\partial_{tt}\phi$ which could generate, in principle, an irrotational component, and we have introduced Eqs. (12) to preserve the solenoidality of ${\bm{E}}_{T}$ (Schmitz & Grauer, 2006a). The spatial dependence of $n_{\alpha}$ on the left-hand side of Eq. (11) has been neglected ($n_{\alpha}\simeq n_{\alpha,0}$) to let coefficients be constant (Valentini et al., 2007).

2.1 Conservation properties

It is straightforward to verify that Eqs. (8–13) satisfy the mass $\int d{\bm{x}}d\bm{v}f_{\alpha}$ and entropy $S_{\alpha}=\int d\bm{x}d\bm{v}f_{\alpha}\log f_{\alpha}$ conservation. The energy conservation equation, obtained by multiplying Eq. (8) by $m_{\alpha}{\bm{v}}^{2}/2$, integrating over the phase-space volume $\int d\bm{x}d\bm{v}$ and summing over the species reads:

[TABLE]

where the kinetic energy is $\mathcal{E}_{kin}=\sum_{\alpha}(m_{\alpha}/2)\int d\bm{x}\,n_{\alpha}{\bm{u}}_{\alpha}^{2}$, the thermal energy is $\mathcal{E}_{th}=\sum_{\alpha}(3/2)\int d\bm{x}\,n_{\alpha}T_{\alpha}$, the magnetic energy is $\mathcal{E}_{mag}=\sum_{\alpha}(\bar{m}/2\bar{u}^{2}\zeta^{2})\int d\bm{x}\,{\bm{B}}^{2}$ and the electrostatic energy is $\mathcal{E}_{el}=\sum_{\alpha}(\bar{m}/2\zeta^{2})\int d\bm{x}\,{\bm{E}}_{L}^{2}$. Note that the temperature of the $\alpha$-species is defined as $3n_{\alpha}T_{\alpha}/m_{\alpha}=\int d\bm{v}({\bm{v}}-{\bm{u}}_{\alpha})^{2}f_{\alpha}$ and, to get Eq. (14), we have used $\int d\bm{x}\,{\bm{w}}_{T}\cdot{\bm{w}}_{L}=0$, ${\bm{w}}_{T}$ and ${\bm{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 (1976) in the electrostatic limit and later extended to the full electromagnetic case (Mangeney et al., 2002). Darwin equations are solved through standard Fast Fourier Transform (FFT) algorithms.

In our case, the splitting algorithm for Eq. (8) reads as follows:

[TABLE]

In the first equation ${\bm{v}}$ is a parameter, while in the second equation ${\bm{x}}$, ${\bm{E}}$ and ${\bm{B}}$ are parameters. At time $t$, the solutions of Eqs. (15–16) can be written as $\Lambda_{\bm{x}}(t)F_{\alpha}({\bm{x}},{\bm{v}})$ and $\Lambda_{\bm{v}}F_{\alpha}({\bm{x}},{\bm{v}})$, respectively. In last expressions, $\Lambda_{\bm{x}}$ and $\Lambda_{\bm{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. (2002); while $F_{\alpha}({\bm{x}},{\bm{v}})=f_{\alpha}({\bm{x}},{\bm{v}},t=0)$ is the initial condition. Note that $\Lambda_{\bm{v}}$ depends also on the particle species $\alpha$ and, for the sake of simplicity, we avoid to explicitly report such dependence.

The splitting scheme is a symplettic, second order accurate in time (see Mangeney et al. (2002) where the stability condition of the advection operator is also discussed) and the numerical solution at $t=t_{N}=N\Delta t$ is given by:

[TABLE]

At $t=0$, the distribution function $F_{\alpha}({\bm{x}},{\bm{v}})=f_{\alpha}({\bm{x}},{\bm{v}},t=0)$ is first advected in physical space by a half time-step, obtaining ${\tilde{f}}_{\alpha}({\bm{x}},{\bm{v}},\Delta t/2)$. Then, the following structure is executed:

Computing the moments of ${\tilde{f}}_{\alpha}$ and evaluating the electromagnetic fields ${\bm{E}}_{L}$, ${\bm{E}}_{T}$ and ${\bm{B}}$, at $t=\Delta t/2$, through Eqs. (9–13); 2. 2.

Performing a time-step advection in velocity space: ${\hat{f}}_{\alpha}=\Lambda_{\bm{v}}(\Delta t){\tilde{f}}_{\alpha}$. 3. 3.

Performing a time-step advection in physical space: ${\tilde{f}}_{\alpha}=\Lambda_{\bm{x}}(\Delta t){\hat{f}}_{\alpha}$.

This last structure is repeated in the algorithm, according to Eq. (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_{\bm{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_{\bm{v},\alpha}=[-v^{max}_{\alpha,x},v^{max}_{\alpha,x}]\times[v^{max}_{\alpha,y},v^{max}_{\alpha,y}]\times[-v^{max}_{\alpha,z},v^{max}_{\alpha,z}]$ is discretized by $(2N_{\alpha,v_{x}}+1)\times(2N_{\alpha,v_{y}}+1)\times(2N_{\alpha,v_{z}}+1)$ grid points. Velocity-space boundary conditions impose $f_{\alpha}(|v_{i}|>v^{max}_{\alpha,i})=0$ ($i=x,y,z$). In order to ensure mass conservation, $v^{max}_{\alpha,i}$ is typically set to be a large multiple of the thermal speed $v_{th,\alpha}=\sqrt{T_{\alpha}/m_{\alpha}}$.

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 ($1D$, $2D$, or $3D$), 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 ${\bm{B}}$ and ${\bm{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 VDF by Mangeney et al. (2002), 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 (2017) and it is essential to perform high-resolution Eulerian Vlasov simulations, in particular in $3D$. Preliminary tests on performance and scalability are reported in Sect. 6.

2.3 Normalizations of the Vlasov-Darwin equations

In order to normalize Eqs. (8–13), three possible choices have been implemented in ViDA:

Electrostatic normalization. Characteristic quantities are: length $\bar{L}=\lambda_{D,e}$, time $\bar{t}=\omega_{p,e}^{-1}$, velocity $\bar{U}=v_{th,e}$ and mass $\bar{m}=m_{e}$. Here $\lambda_{D,e}=\sqrt{T_{e}/4\pi n_{e}e^{2}}$, $\omega_{p,e}=(\sqrt{4\pi n_{e}e^{2}/m_{e}})^{-1}$, $v_{th,e}=\sqrt{T_{e}/m_{e}}=\lambda_{D,e}\omega_{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. 2. 2.

Electromagnetic normalization. Characteristic quantities are: length $\bar{L}=d_{e}$, time $\bar{t}=\omega_{p,e}^{-1}$, velocity $\bar{U}=c$ and mass $\bar{m}=m_{e}$, where $d_{e}=c/\omega_{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. 3. 3.

Hybrid normalization. Characteristic quantities are: length $\bar{L}=d_{p}$, time $\bar{t}=\Omega_{c,p}^{-1}$, velocity $\bar{U}=v_{A}$ and mass $\bar{m}=m_{p}$. In previous expressions $\Omega_{cp}=eB_{0}/m_{p}c$, $v_{A}=B_{0}/\sqrt{4\pi n_{p}m_{p}}$, $d_{p}=v_{A}/\Omega_{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 behavior and reliability of the code has been tested against the propagation of Langmuir waves, in both linear and nonlinear regime, 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 $\delta n_{e}/n_{e,0}=A\sin(kx)$, $A=10^{-4}$ and $k=k_{1}=2\pi/L$ being the amplitude and the wave-number, respectively. The box length is $L=18\lambda_{D,e}$ ($k_{1}=0.35$) and $v^{max}_{e,i}=5v_{th,e}$ ($i=x,y,z$). The system evolution is reproduced up to a maximum time $t_{max}=100\omega_{p,e}^{-1}$, while the numerical recurrence time is $t_{rec}=2\pi/k\Delta v\simeq 180\omega_{p,e}^{-1}$(Cheng & Knorr, 1976; Galeotti et al., 2006; Pezzi et al., 2016).

The left panel of Fig. 1 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, 1946); the observed damping rate shows a very good agreement with the theoretical prediction $\gamma_{L}=-3.37\times 10^{-2}\omega_{p,e}$ (red-dashed line), evaluated through a numerical linear solver for the roots of the electrostatic Vlasov dielectric function. In the right panel of Fig. 1 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},\omega_{R})$ as a function of the pulsation $\omega_{R}$. As expected, the resonant curve displays a well-defined frequency peak in correspondence of a value of the pulsation $\omega_{R}=1.22\omega_{p,e}$. In the right panel of the figure, the vertical red-dashed line represents the value of the theoretical resonant pulsation $\omega_{R,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 $\Delta\omega_{R}=2\pi/t_{max}\simeq 0.063\omega_{p,e}$. Again, numerical results are in very good agreement with theoretical predictions.

In order to show the dependence of real $\omega_{R}$ and imaginary $\omega_{I}$ part 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\pi/L$ ($L=2\pi 10\lambda_{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 $\omega_{R}$ (left) and $\omega_{I}$ (right), in units of $\omega_{pe}$, as a function of the wave number $k\lambda_{D,e}$. A very good agreement with theoretical expectations (red-dashed curves) is recovered for both real and imaginary parts of the complex frequency.

We conclude this section by focusing on the nonlinear regime of the Langmuir wave dynamics (see, for example, Brunetti et al. (2000) and refs. 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_{max}=400\omega_{pe}^{-1}$. The box length is $L_{x}=18\lambda_{D,e}$, while $v^{max}_{e,i}=5v_{th,e}$ ($i=x,y,z$). As reported in the left panel of Fig. 3, the time evolution of the electric-field Fourier component shows an early linear damping phase (Landau, 1946), until particle trapping arrests the damping and produces oscillations of the signal envelope (O’Neil, 1965). At large times, a phase-space vortex is observed in the electron DF in the vicinity of the wave phase speed, as reported in the center panel of Fig. 3(center). As it has been shown in O’Neil (1965), the nonlinear trapping time $\tau_{nl}$ depends on the saturation amplitude $A_{1}$ of the electric oscillations as $\tau_{nl}\sim A_{1}^{-1/2}$. For each of the ten simulations, we evaluated $A_{1}$ and $\tau_{nl}$ at the time of the first peak of the electric envelope oscillations. The right panel of Fig. 3 shows in log-log plot $\tau_{nl}$ as a function of $A_{1}$ (stars), compared to the theoretical expectation (red line), showing a very nice agreement.

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 ${\bm{B}_{0}}=B_{0}{\bm{e}}_{x}$ ($B_{0}=1$) and protons are not fixed. The box length is $L_{x}=2\pi 10d_{e}$, while $v^{max}_{e,i}=10v_{th,e}$ and $v^{max}_{p,i}=7v_{th,p}$ ($i=x,y,z$). Note that increasing the value of $v^{max}_{p(e)}$ 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_{th,e}/c=10^{-3}$ and $\zeta=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:

[TABLE]

where $A=10^{-3}$ and $k=k_{1}=2\pi/L_{x}$.

By solving Darwin equations, these perturbations generate a current density and then magnetic fluctuations. Figure 4 reports the time evolution of $|B_{y}|(k_{1},t)$ (left) and the frequency peak of the spectral magnetic energy $|B_{y}|^{2}(k_{1},\omega_{R})$ as a function of the pulsation $\omega_{R}$ (right). The magnetic field clearly oscillates at the correct frequency $\omega_{R,th}=0.91\omega_{p,e}$, that can be evaluated from the linear dispersion relation for whistler waves (obtained by assuming motionless protons and cold electrons): $\omega_{R,th}(k)=B_{0}k^{2}/(1+k^{2})$ (Krall & Trivelpiece, 1973). Note that, for the considered wave-number, a negligible damping of whistler waves is expected. In the left panel of Fig. 4, red-dashed line represents the value of the resonant pulsation from above expression for $\omega_{R,th}$, while the two vertical red dot-dashed lines indicate the interval of numerical uncertainty $\Delta\omega_{R}=2\pi/t_{max}\simeq 0.063\omega_{p,e}$.

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 ${\bm{B}_{0}}=B_{0}{\bm{e}}_{x}$ and $B_{0}=1$. The box length is $L_{x}=2\pi 50d_{p}$, while $v^{max}_{e,i}=5v_{th,e}$ and $v^{max}_{p,i}=5v_{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 $m_{p}/m_{e}=25$, thus avoiding extremely small timesteps, while $T_{e}/T_{p}=1$, $\bar{u}=v_{A}/c=10^{-3}$ and $\zeta=c/v_{A}=10^{3}$. Initial proton $\beta_{p}$ is $\beta_{p}=2v_{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:

[TABLE]

where $A=10^{-4}$ and $k=k_{1}=2\pi/L_{x}$. Figure 5 shows the time evolution of $|B_{z}|(k_{1},t)$ (left) and the magnetic spectral energy $|B_{z}|^{2}(k_{1},\omega_{R})$ as a function of $\omega_{R}$ (right). The recovered resonant peak is in agreement with the theoretical pulsation, evaluated through a fully-kinetic linear solver of the dispersion relation (Camporeale & Burgess, 2017). Moreover, no Landau damping is observed, since it occurs at much smaller scales (Barnes, 1966; Vàsconez et al., 2014; Camporeale & Burgess, 2017). 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.

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, 1959; Gary & Karimabadi, 2006; Califano et al., 2008; Palodhi et al., 2009, 2010; Chen & Chacón, 2014, 2015; Camporeale & Zimbardo, 2015).

Here we focus on the development of the electron Weibel instability, that produces electromagnetic fluctuations transverse to the wavevector ${\bm{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 ${\bm{k}}=k\hat{{\bm{e}}}_{x}$, although we have verified that instability is triggered in the same way also in the different phase space configuration ($2D$ and $3D$). 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_{th,e,x}=2.5\times 10^{-2}c$ and $v_{th,e,y}=v_{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_{th,p}=v_{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^{max}_{e(p)}=5v_{th,e(p)}$ in each velocity directions.

The initial equilibrium has been perturbed through a sinusoidal, transverse perturbation, imposed on the electron bulk speed:

[TABLE]

where $A=2\times 10^{-5}$ and $k=k_{1}=2\pi/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 $\omega_{I}^{th}\simeq 4\times 10^{-3}\omega_{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, 2014).

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 GEM challenge (Birn et al., 2001), in order to allow for a direct comparison to previous studies (Birn et al., 2001; Schmitz & Grauer, 2006b). For this reason, we have also chosen the hybrid normalization (see Sect. 2.3).

The equilibrium is set by adapting the Harris equilibrium (Harris, 1962) 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:

[TABLE]

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 (1962) and $L_{1}$ is the corresponding current sheet thickness. The second and third hyperbolic tangent in Eq.(24) 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 uniform at the initial time and the density $n(y)$ is defined in order to satisfy pressure balance. Then, from Eq. (5) and considering $\partial_{t}{\bf E}_{L}=0$ at the initial time, we get the initial current density ${\bf j}=(0,\ 0,\ j_{z}(y))$.

Following the prescriptions of the Harris equilibrium we get, in normalized units,

[TABLE]

Eq. (26) 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 Eqs.(26)–(27) we have:

[TABLE]

It is worth to point 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., 2001), fluctuations are superposed to the initial magnetic field in order to obtain a single magnetic island at the center of the space domain at the initial time. In particular, $\delta\mathbf{B}=\nabla\delta\psi\times\hat{z}$ and

[TABLE]

where, as already stated, $L_{x}$ and $L_{y}$ are the lengths of the spatial domain in $x$ and $y$ direction respectively. According to GEM challenge, in scaled units, $\psi_{0}$ is set to $0.1$.

By using the relation $\delta\mathbf{B}(x,y)=\nabla\delta\psi(x,y)\times\hat{z}$ and Eq.(10), we derive the expression for the current density fluctuations $\delta\mathbf{j}(x,y)$ consistent with $\delta\psi(x,y)$. In particular, it is possible to define $\delta\mathbf{j}(x,y)=(0,0,\delta 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}^{max}=5\ v_{th,e}$ and $v_{p}^{max}=5\ v_{th,p}$, where the normalized $v_{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. (2001).

In Figure 7 we show the evolution of the reconnected flux given by the difference $\Delta\psi$ between the magnetic flux $\psi$ evaluated at the X point and at the O point. Accordingly to the initial perturbation, the X and the O point are initially located at $(L_{x}/2,L_{y}/2)$ and $(0,L_{y}/2)$ and their location does not significantly change throughout the simulation run. The behavior of $\Delta\psi$ is very similar to the evolution of the reconnected flux in Ref. (Birn et al., 2001). Reconnection evolves with a reconnected flux that remains close to zero until $t\sim 15\ \Omega_{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\ \Omega_{c,p}^{-1}$.

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 $\psi$ 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 et al., 1994; Uzdensky & Kulsrud, 2006). The quadrupolar structure that we find is analogous to the one obtained with other kinetic codes, both Eulerian (Schmitz & Grauer, 2006b, see Fig. 2) and Lagrangian (Pritchett, 2001, 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. (2001) (see Fig.6(d)). The electron current at the X line has a thickness comparable to $d_{p}$ which corresponds to $5\ d_{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 Ref. (Schmitz & Grauer, 2006b).

In Figure 9 we show the reconnection outflow of protons and electrons at $t^{*}=18.13\ \Omega_{c,p}$. In particular, we note that at $x=3\ d_{p}$ (panel (a)), corresponding to a distance of $9.8\ d_{p}$ from the X-point located at $L_{x}/2=12.8\ d_{p}$, the electron velocity is characterized by two peaks corresponding to the separatrices, while the proton velocity is concentrated in the center 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 Fig. 8(b). Fig. 9(b) shows the same quantities of Fig.9(a) at a distance of $2.3\ d_{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}$.

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 center (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 by 1KNL processor (68 cores, 1.40GHz), 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.

ViDA numerically integrates VD equations in a six-dimensional phase-space ($3D$–$3V$: $x,y,z,v_{x},v_{y},v_{z}$). Only the $3D$ physical space is parallelized using cubic cells (squared in the $2D$ configuration). For implementing these tests, we have chosen $51^{3}$ velocity gridpoints for each particle VDF (protons and electrons), which represent a typical value adopted in production runs, and we have performed about $100$ timesteps per test (note that changing the step number does not act on the code scalability). Note that two VDFs 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 analyzed the parallel performance in the $2D$–$3V$ configuration, adopting a physical-space grid with $1024\times 1024$ points. This setup requires about 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 Fig. 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 Fig. 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 $3D$ configuration with $128\times 128\times 64$ grid points in the spatial domain. The strong scaling from $128$ to $1024$ cores is shown in Fig. 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 Fig 10(d), it can be appreciated that the parallel efficiency is high only up to several hundreds cores.

These preliminary tests show a reasonable parallel efficiency on KNL architecture, at least up to some hundreds 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 to finally highlight 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 about twice as computationally expensive as the HVM code (Valentini et al., 2007), which has been recently used for 3D simulations of plasma turbulence (see for instance Cerri et al. (2018)). 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 Marconi supercomputer using 16 nodes and 512 MPI processes. On the other hand, being ViDA a code for a new piece of physics, it is difficult to foresee for the exact cost of a 3D 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 3D 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.

for a new piece of physics, it is difficult to foresee for the exact cost of a 3D 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.

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 timestep 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 $\lambda_{D}$. As typically the case for space plasmas, but often also in the laboratory, inter-particle collisions are not described, since collisional scales are assumed to be smaller than other characteristics dynamical scales.

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 at today a strong computational challenge and a frontier problem for the understanding of the 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 the Ohm’s law and the mechanisms of electron heating, which are among the main targets of satellite MMS data analysis (Torbert et al., 2016; Genestreti et al., 2018; Cozzani et al., 2019). 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., 2016).

Acknowledgements.

This project (FC, AR, FV) has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 776262 (AIDA, www.aida-space.eu). OP 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”. EC 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).