Menura: a code for simulating the interaction between a turbulent solar wind and solar system bodies

Despite the close relationship between planetary science and plasma physics, few advanced numerical tools allow to bridge the two topics. The code Menura proposes a breakthrough towards the self-consistent modelling of these overlapping field, in a novel 2-step approach allowing for the global simulation of the interaction between a fully turbulent solar wind and various bodies of the solar system. This article introduces the new code and its 2-step global algorithm, illustrated by a first example: the interaction between a turbulent solar wind and a comet.


Introduction
For about a century, three main research fields have taken an interest in the various space plasma environments found around the Sun. On the one hand, two of them, namely planetary science and solar physics, have been exploring the solar system, to understand the functioning and history of its central star, and of its myriad of orbiting bodies. On the other hand, the third one, namely fundamental plasma physics, has been using the solar wind as a handy wind tunnel which allows researchers to study fundamental plasma phenomena not easily reproducible on the ground in laboratories. During the last two decades, bridges between these communities have been developing, as the growing knowledge of each community was bringing the fields ever closer, to a point where overlapping topics made the communities work together. For instance, if initially planetary scientists were studying the interaction between solar system bodies and a steady, ideally laminar solar wind, they soon had to consider the event-full and turbulent nature of the solar wind to go further in the in situ space data analysis, further in their understanding of the interactions at various obstacles. If plasma physicists were originally interested in a pristine solar wind unaffected by the presence of obstacles, they realised that the environment close to these obstacles could provide combinations of plasma parameters otherwise not accessible to their measurements in the unaffected solar wind. For a while now, we have seen planetary studies focusing on the effects of solar wind transient effects (such as Coronal mass Ejection CME or Co-rotational Interaction Region CIR) on planetary plasma environments, at Mars [1], Mercury [2], Venus [3] and comet 67P/C-G [4,5] to only cite a few, the effect of large scale fluctuations in the upstream flow on Earth's magnetosphere [6], and more generally the effect of solar wind turbulence on Earth magnetosphere and ionosphere [7,8]. Similarly, plasma physicists have developed comprehensive knowledge of plasma waves and plasma turbulence in the Earth magnetosheath, presenting relatively high particle densities and electromagnetic field strengths favourable for space instrumentation, in a region more easily accessible to space probes than regions of unaffected solar wind [9,10]. More recently, the same community took an interest in various planet magnetospheres, depicting plasma turbulence in various locations and of various parameters [11] and all references therein.
Various numerical codes have been used for the global simulation of the interaction between a laminar solar wind and solar system bodies, using MHD [12], hybrid [13], or fully kinetic [14] solvers. Similarly, solar wind turbulence in the absence of an obstacle has also been simulated using similar MHD [15], hybrid [16], and fully kinetic [17] solvers. In this context, we identify the lack of a numerical approach for the study of the interaction between a turbulent plasma flow (such as the solar wind) and an obstacle (such as a magnetosphere, either intrinsic or induced). Such a tool would provide the first global picture of these complex interactions. By shading new lights on the long-lasting dilemma between intrinsic phenomena and phenomena originating from the upstream flow, it would allow invaluable comparisons between self-consistent, global, numerical results, and the worth of observational results provided by the various past, current and future exploratory space missions in our solar system.
The main points of interest and main questions motivating such a model can be organised as such: • Macroscopic effects of turbulence on the obstacle shape and position of the plasma boundaries (e.g. bow shock, magnetopause), large scale magnetic reconnection, atmospheric escape, dynamical evolution of the magnestosphere.
• Microscopic physics and instabilities within the interaction region, induced by upstream turbulence energy transport by plasma waves, energy conversion by wave-particle interactions, energy transfers by instabilities.
• The way incoming turbulence is processed by planetary plasma boundaries sudden change of spatial and temporal scales, change of spectral properties, existence of a memory of turbulence downstream magnetospheric boundaries.
Indirectly, because of the high numerical resolution required to properly simulate plasma turbulence, this numerical experiment will provide an exploration of the various obstacles with the same high resolution in both turbulent and laminar runs, resolutions that have rarely been reached for planetary simulations, except for Earth's magnetosphere.
Menura, the new code presented in this publication, splits the numerical modelling of the interaction into two steps.
Step 1 is a decaying turbulence simulation, in which electromagnetic energies initially injected at the large spatial scales of the simulation box cascades towards smaller scales.
Step 2 uses the output of Step 1 to introduce an obstacle moving through this turbulent solar wind.
The code is written in c++ and uses CUDA APIs for running its solver exclusively on Graphics Processing Units (GPUs). Section 2 introduces the solver, which is tested against classical plasma phenomena in Section 3. Sections 4 and 5 tackle the first and second step of the new numerical modelling approach, illustrating the decaying turbulence phase, and introducing the algorithm for combining the output of Step 1 together with the modelling of an obstacle (Step 2). Section 6 presents the first global result of Menura, providing a glimpse of the potential of this numerical approach, and introducing the forthcoming studies.
Menura source code is open source, available under the GNU General Public License license.

The solver
In order to (i) achieve global simulations of the interactions while (ii) modelling the plasma kinetic behaviour, with regard to the computation capabilities currently available, a hybrid Particle-In-Cell (PIC) solver has been chosen for Menura. This well-established type of models resolves the Vlasov equation for the ions by discretising the ion distribution function as macro-particles characterized by discrete positions in phase space, and electrons as a fluid, with characteristics evaluated at the nodes 1 of a grid, together with ion moments and electromagnetic fields. The fundamental computational steps of a hybrid PIC solver are: • Particles' position advancement, or "push".
• Particles' moments mapping, or "gathering": density, current, eventually higher order, as required by the chosen Ohm's law.
• Electromagnetic field advancement, using either an ideal, resistive or generalised Ohm's law and Faraday's law.
These steps are summarised in Figure 6. Details about these classical principles can be found in [18] and references therein. The bottleneck of PIC solvers is the particles' treatment, especially the velocity advancement and the moments computation (namely density and current). The simulation of plasma turbulence especially requires large amounts of macro-particles per grid nodes. We therefore want to minimise both the amount of operations done on the particles and the number of particles itself. A popular method which minimises the amount of these computational passes through all particles is the Current Advance Method (CAM) [19], for instance used for the hybrid modelling of turbulence by [16]. Figure 1 presents Menura's solver algorithm, built around the CAM, similar to the implementation of [13]. In this scheme, only four passes through all particles are performed, one position and one velocity pushes and two particle moments mappings. The second moment mapping in Figure 1, i.e. step 2, also produces the two pseudo-moments Λ and Γ used to advance the current as: with E * the estimated electric field after the magnetic field advancement of step 4. W (r n+1 ) is the shape function, which attributes different weights for each node surrounding the macro-particle [18].
Central finite differences for evaluating derivatives and second order interpolations are used throughout the solver. The grid covering the physical simulation domain has an additional 2-node wide band, the guard or ghost nodes, allowing to solve derivatives using (central) finite differences at the very edge of the physical domain. For periodic boundary conditions, as used along all directions during Step 1 of the simulation, the value at the opposite edge of the physical domain are copied in the guard nodes. Other boundary conditions will be discussed later when introduced. r n Step 1 Step 0 --Position push --Moments, Free-streaming J - Step 2 --Moments, Free-streaming J + Step 3 --Moments averaging Step 4 --Faraday/Ohm's laws Step 5 --Generalised Ohm's law, E estimate Step 6 --Current advance Step 7 Step 8 --Velocity push, Boris v n+3/2 Step 8 Step Step 5 (B n+1 , J + , ρ n+1 ) Step 7 Step 4 Step 6 Step 3 Step 0 Step 2 --Generalised Ohm's law, E corrected The mapping of the particle moments are done using an order-two, triangular shape function: one macro-particle contributes to 9 grid nodes in 2D space ( respectively 27 in 3D space), using 9 (respectively 27) different weights. The interpolation of the fields values from the nodes to the macroparticles' positions uses the exact same weights, with 9 (respectively 27) neighbouring nodes contributing to the fields values at a particle position.
As illustrated in Figure 1, the position and velocity advancements are done at interleaved times, similarly as a classical second order leap-frog scheme. However, since the positions of the particles are needed to evaluate their acceleration, the CAM scheme is not strictly speaking a leap-frog integration scheme. Another difference in this implementation is that velocities are advanced using the Boris method [20].
The Ohm's law is at the heart of the hybrid modelling of plasmas. Menura uses the following form of the law, here given in SI units. In this formulation, the electron inertia is neglected, and the quasi-neutral approximation n ∼ n i ∼ n e is used [17]. Additionally, neglecting the time derivative of the electric field in the Ampere-Maxwell's law (Darwin's hypothesis), one get the total current through the curl of the magnetic field. This formulation highlights the need for only three variables to be followed through time, namely the magnetic field, and the particles position and velocity, while all other variables can be reconstructed from these three.
The Faraday's law is used for advancing the magnetic field in time: The electron pressure is obtained assuming it results from a polytropic process, with an arbitrary index κ, to be carefully chosen by the user.
Using much less memory than the particles' variables, the fields can be advanced in time using a smaller time step and another leap-frog-like approach, as illustrated in Figure 1, step 4 [19]. Additional spurious high-frequency oscillations are the default behaviour of finite differences schemes. Two main families of methods are used to filter out these features, the first being an additional step of field smoothing, the second using the direct inclusion of a diffusive term in the differential equation of the system, acting as a filter [21]. For Menura, we have retained the second approach, implementing a term of hyper-resistivity in the Ohm's law.
The stability of hybrid solvers is sensitive to low ion densities. We use a threshold value equal to a few percent of the background density, 5% in the following examples, threshold below which a node is considered as a vacuum node, and only the resistive terms of the generalised Ohm's law of Equation 4 are solved using a higher value of resistivity η h vacuum [22]. This way, terms proportional to 1/n do not exhibit nonphysical values where the density may get locally very low, due to the thermal noise of the PIC macro-particle discretisation.
All variables in the code are normalised using the background magnetic field amplitude B 0 and the background plasma density n 0 . All variables are then expressed in terms of either these two background values, or equivalently in terms of the proton gyrofrequency ω ci0 and the Alfven velocity v A0 . We define normalised variablesã as obtained by dividing its physical value by its "background" value:ã = a a 0 All background values are given in Table 1, and the normalised equations of the solver are given in Appendix A.

Physical tests
In this section, the code is tested against well-known, collisionless plasma processes, and their solutions given by the linear full kinetic solver WHAMP [23]. We first explore MHD scales, simulating Alfvénic and magnetosonic modes. We use a 2-dimensional spatial domain with one preferential dimension chosen as x. A sum of six cosine modes in the component of the magnetic field along the x-direction direction are initialised, corresponding to the first six harmonics of this periodic box. The amplitude of these modes is 0.05 times the background magnetic field B 0 , which is taken either along (Alfvén mode) or across (magnetosonic mode) the propagation direction x. Data are recorded along time and along the main spatial dimension x, resulting in the 2D field B(x, t). The 2-dimensional Fourier transform of this field is given in Figure 2. In this (ω, k)-plane, each mode can be identified as a point of higher power, six points for six initial modes. The solutions given by WHAMP for the same plasma parameters are shown by the solid lines, and a perfect match is found between the two models. Close to the ion scale k.d i = 1, WHAMP and Menura display two different branches that originate from the Alfvén mode, splitting for higher frequencies into the whistler and the ion cyclotron branches. The magnetosonic modes were also tested using different polytropic indices, resulting in a shift of the dispersion relation along the ω-axis. Changing polytropic indices for both models resulted in the same agreement.
With the MHD scales down to ion inertial scales now validated, we explore the ability of the solver to account for further ion kinetic phenomena, first with the classical case of the two-stream instability (also known as the ionbeam instability, given the following configuration). Two Maxwellian ion beams are initialised propagating with opposite velocities along the main dimension x. A velocity separation of 15v th is used in order to excite only one unstable mode. The linear kinetic solver WHAMP is used to identify the expected growth rate associated to the linear phase of the instability, before both beams get strongly distorted and mixed in phase space during the nonlinear phase of the instability (not capture by WHAMP). During this linear phase, Menura results in a growing circularly polarised wave, and the amplitude's growth of the wave is shown in Figure 3. Both growths match perfectly.
Finally, we push the capacities of the model to the case of the damping of an ion acoustic wave through Landau resonance. A very high number of macro-particles per grid node is required to resolve this phenomenon, so enough resonant particles take part in the interaction with the wave. The amplitude of the initial, single acoustic mode is taken as 0.01 times the background density, taken along the main spatial dimension of the box. The decrease in the density fluctuation through time, spatially averaged, is shown in Figure 3, with again the corresponding solution from WHAMP. A satisfying agreement is found during the first 6 oscillations, before the noise in the hybrid solver output (likely associated to the macroparticle thermal noise) takes over. Admittedly, the amount of particle per node necessary to well resolve this phenomenon is not practical for the global simulations which Menura (together with all global PIC simulations) aims for. For the classical tests presented above, spanning over MHD and ion kinetic scales tests, Menura agrees with theoretical and linear results. In the next section, the simulation of a decaying turbulent cascade provides one final physical validation of the solver, through all these scales at once.

Step 1: Decaying turbulence
We use Menura to simulate plasma turbulence using a decaying turbulent cascade approach: at initial time t = 0, a sum of sine modes with various wave vectors k, spanning over the largest spatial scales of the simulation domain, are added to both the homogeneous background magnetic field B 0 and the particles velocities u i . Particle velocities are initialised according to a Maxwellian distribution, with a thermal speed equal to one Alfvén speed. Without any other forcing later-on, this initial energy cascades, as time advances, towards lower spatial and temporal scales, while forming vortices and reconnecting current sheets [16]. Using such Alfvénic perturbation is motivated by the predominantly Alfvénic nature of the solar wind turbulence measured at 1 au [24].
In this 2-dimensional set-up, B 0 is taken along the z-direction, perpendicular to the simulated spatial domain (x, y), whereas all initial perturbations are defined within the simulation plane. Their amplitude is 0.5 B 0 , while their wave vectors are taken with amplitude between k inj min = 0.01 d i0 and k inj max = 0.1 d i0 , so energy is only injected in MHD scales, in the inertial range. Because we need these perturbation fields to be periodic along both directions, the k x and k y of each mode corresponds to harmonics of the simulation box dimensions. Therefore, a finite number of wave vector directions is initialised. Along these constrained directions, each mode in both fields has two different, random phases. The magnetic field is initialised such that is it divergence-free.
For this example, the box is chosen to be 500 d i0 wide on both dimensions, subdivided by a grid 1000 2 nodes wide. The corresponding ∆x is 0.5 d i0 , and spatial frequencies are resolved over the range [0.0062, 6.2] d i0 . The time step is 0.05 ω −1 ci0 . 2000 particles per grid node are initialised with a thermal speed of 1 v A . The temperature is isotropic and a plasma beta of 1 is chosen for both the ion macro-particles and the electronic massless fluid.
At time t = 500 ω −1 ci0 , the perpendicular (in-plane) fluctuations of the magnetic field have reached the state displayed in Figure 4, left-hand panel. Vortices and current sheets give a maximum B ⊥ /B 0 of about 1, a result consistent with solar wind turbulence observed at 1 au [24]. The omnidirectional power spectra of both the in-plane magnetic field fluctuations and the in-plane ion bulk velocity fluctuations are shown in the right-hand panel of the same figure. Omni-directional spectra are computed as follows, withf the (2D) Fourier transform of f : These spectra are not further normalised and are given in arbitrary units. We then compute a binned statistics over this 2-dimensional array to sum up its values within the chosen bins of k ⊥ , which correspond to rings in the (k ⊥ , k )-plane. The width of the rings is arbitrarily chosen so the resulting 1-dimensional spectrum is well resolved (not too few bins), and not too noisy (not too many bins).
For a vector field such as B ⊥ = (B x , B y ), the spectrum is computed as the sum of the spectra of each field component: The perpendicular magnetic and kinetic energy spectra exhibit power laws over the inertial (MHD) range consistent with spectral indexes -5/3 and -3/2, respectively, between k inj max = 0.1 d i0 and break points around 0.5 d i0 . We remind that a spectral index -5/3 is consistent with the Goldreich-Sridhar strong turbulence phenomenology [25] that leads to a Kolmogorov-like scaling in the plane perpendicular to the background magnetic field, while a B 0 Table 2: Initial parameters of the decaying turbulence run spectral index -3/2 is consistent with the Iroshnikov-Kraichnan scaling [26]. These spectral sloped are themselves consistent with observations of magnetic and kinetic energy spectra associated to solar wind turbulence [27,28]. For higher wavenumbers, both spectral slopes get much steeper, and after a transition region within [0.5, 1.] d i0 get to a value of about -3.2 and -4.5 when reaching the proton kinetic scales for respectively the perpendicular magnetic and kinetic energies, consistent with spectral index found at subion scales by previous authors [16, 29, e.g.]. Additionally, the initial spectra of the magnetic field and bulk velocity perturbations are over-plotted, to show respectively where the energy is injected in spatial frequencies and the level of noise introduced by the finite number of particles per node used.

Step 2: Obstacle
Menura has shown satisfactory results on plasma turbulence, over three orders of magnitude in wavenumbers. We now start the second phase of the simulation, resuming it at t = 500 ω −1 ci0 , corresponding to the snapshot studied in the previous section. We keep all parameters unchanged, but add an obstacle with a relative velocity with regard to the frame used in the first phase, evolving through this developed turbulence. Particles and field are advanced with the exact same time and spatial resolutions as previously, so the interaction between this obstacle and the already-present turbulence is solved with the same self-consistency as in the first phase, with only one ingredient added: the obstacle.

A comet
This obstacle is chosen here to be an intermediate activity comet, meaning that its neutral outgassing rate is typical of an icy nucleus at a distance of about 2 au from the Sun. Comprehensive knowledge on this particular orbital phase of comets has recently been generated by the European Rosetta mission, which orbited its host comet for two years [30]. The first and foremost interest of such an object for this study is its size, which can be evaluated using the gyroradius of water ions in the solar wind at 2 au. The expected size of the interaction region is about 4 times this gyro-radius [31], and with the characteristic physical parameters of Table 2, the estimated size of the interaction region is 480 d i0 . In other words, the interaction region spans exactly over the range of spatial scales probed during the first phase of the simulation, including MHD and ion kinetic scales. The second interest of a comet is its relatively simple numerical implementation. Without a solid body, without gravity and without an intrinsic magnetic field, the obstacle is only made of cometary neutral particles being photo-ionised within the solar wind. Over the scales of interest for this study, the neutral atmosphere can be modelled by a 1/r 2 radial density profile, and considering the coma to be optically thin, ions are injected in the system with a rate following the same profile. This is the Haser Model [32], and simulating a comet over scales of hundreds of d i0 only requires to inject cold cometary ions at each time step with the rate with r the distance from the comet nucleus of negligible size, ν i the ionisation rate of cometary neutral molecules, n 0 the neutral cometary density, Q the neutral outgassing rate, u 0 the radial expansion speed of the neutral atmosphere.

Reference frame
The first phase of the simulation, the decaying turbulence phase, was done in the plasma frame, in which the average ion bulk velocity is 0. Classically, planetary plasma simulation are done in the planet reference frame: the obstacle is static and the wind flows through the simulation domain. In this case, a global plasma reference frame is most of the time not defined. In Menura, we have implemented the second phase of the simulation -the interaction phase -in the exact same frame as the first phase, which then corresponds to the plasma frame of the upstream solar wind, before interaction. In other words, the turbulent solar wind plasma is kept "static", and the obstacle is moving through this plasma. The reason motivating this choice is to keep the turbulent solar wind "pristine", by continuing its resolution over the exact same grid as in phase one. Another motivation for working in the solar wind reference frame is illustrated in Figure 5, in which we compare the exact same simulation done in each frame, using a laminar upstream flow. If the macroscopic result remains unchanged between the two frames, we find strong small scale numerical artifacts propagating upstream of the interaction in the comet reference frame, absent in the solar wind reference frame. Small scale oscillations are current in hybrid PIC simulations, and are usually filtered with either resistivity and/or hyper-resistivity, or with an ad-hoc smoothing method. Note that none of these methods are used in the present example. We demonstrate here the role of the reference frame in the production of one type of small scale oscillations, and insure that their influence over the spectral content of upstream turbulence is minimised, already without the implemented hyper-resistivity.
To summarise, by keeping the same reference frame during Step 1 and 2, the only effective difference between the two phases is the addition of sunward moving cometary macro-particle.  Figure 5: The interaction between a comet and a laminar flow, in the object rest frame (right) and the upstream solar wind reference frame (left). The magnetic field amplitude is shown.

Algorithm
By working in the solar wind reference frame, the obstacle is moving within the simulation domain. Eventually, the obstacle would reach the boundaries of the box, before steady-state is reached. We therefore need to somehow keep the obstacle close to the centre of the simulation domain. This is done by shifting all particles and fields of n∆x every m iterations, n, m ∈ N, as illustrated in Figure 6. Using integers, the shift of the field is simply a side-way copy of themselves without the need of any interpolation, and the shift of the particles is simply the subtraction of n∆x to their xcoordinate. Field values as well are particles ending up downstream of the simulation domain are discarded. This leaves only the injection boundary to be dealt with. There, we simply inject a slice of fields and particles picked from the output of Step 1, using the right slice index in order to inject the continuous turbulent solution, as shown in Figure 6. These slices are n∆x wide.
With idx it the index of the iteration, the algorithm illustrated in Figure  6 is then: • Inject cometary ions according to q i (r) (cf. Eq. 11) • Advance particles and fields (cf. Figure 1) • If idx it%m=0 -Shift particles and fields of -n∆x -Discard downstream values -Inject upstream slice idx slice from Step 1 output This approach has one constraint, we cannot fine-tune the relative speed v 0 between the wind and the obstacle, which has to be v 0 = n m ∆x ∆t (12) in order for the obstacle to come back to its position every m iterations, and therefore not drift up-or downstream of the simulation domain.

CUDA and MPI implementation, performances
The computation done by Menura's solver (Figure 1) is entirely executed on multiple GPUs (Graphics Processing Units), written in c++ in conjunction with the CUDA programming model and the MPI standard, which allows to split the problem and distribute it over multiple cards (i.e. processors). GPUs can run simultaneously thousands of threads, and can therefore tremendously accelerate such applications. The first implementation of a hybrid-PIC model on such devices was done by [33]. However their still limited memory (up to 80GB at the time of writing) is a clear constraint for large problems, especially for turbulence simulation which requires a large range of spatial scales and a very large number of particles per grid node. The use of multiple cards becomes then unavoidable, and the communications between them is implemented using a CUDA-aware version of MPI, the Message Parsing Interface. The division of the simulation domain in the current version of Menura is kept very simple, with equal size rectangular sub-domains distributed along the direction perpendicular to the motion of the obstacle: one sub-domain spans the entire domain along the x-axis with its major dimension, as shown in Figure 6. MPI communications are done for particles after each position advancement, and for fields after each solution  of the Ohm's law and the Faraday's law. But since the shift of fields and particles described in the previous section is done purely along the obstacle motion direction, no MPI communication is needed after the shifts, thanks to the orientation of the sub-domains.
Another limitation in using GPUs is the data transfer time between the CPU and the cards. In Menura, all variables are initialised on the CPU, and are saved from the CPU. Data transfers are then unavoidable, before starting the main loop, and every time we want to save the current state of the variables. During Step 2 of the simulation, a copy of the outputs of Step 1 is needed, which effectively doubles the memory needed for Step 2. This copy is kept on the CPU (in the tank object) in order to make the most out of the GPUs memory, in turn implying that more CPU-GPU communications are needed for this second step. Every time we inject a slice of fields and particles upstream of the domain, only this amount of data is copied from the CPU to the GPUs, using the injector variable as sketched in Figure 6.
For Step 1, the decaying turbulence ran 10000 iterations, four nvidia V100 GPUs were used with 16GB memory each, corresponding to one complete node of the IDRIS cluster Jean Zay. A total of 2 billion particles (500 million per card) were initialised. The time for the solver on each card reached a bit more than three hours, with a final total coast of about 13 hours of computation time for this simulation, taking into account all four cards, and the variables initialisation and output.
Step 2 was executed on larger V100 32GB cards, providing much more room for the addition of 60000 cometary macro-particles per iteration. During Step 1, 87.3% of the computation time was spent on moments mapping, i.e. steps 0 and 2 in the algorithm of Figure 1, while respectively 2.7% and 0.8% were spent on advancing the particles velocity and position. The computation of the Ohm's and Faraday's laws sums up to 0.5%. 0.9% was utilised for MPI communications of field variables, while only 0.08% was dedicated for particles MPI communications, due to the limited particle transport happening in Step 1.
91% of the total solver computation time is devoted to particles treatment, with 96% of that part spent on particles moment mapping, which might seem a suspiciously large fraction. We note however that such a simulation is characterised by its large number of particles per node, 2000 in our case. 99.6% of the total allocated memory is devoted to particles. The time spent to map the particles on the grid is also remarkably larger than the time spent to update their velocity, despite both operations being based on the same interpolation scheme. However, during the mapping of the particles moments, thousands of particles need to increment the value at particular memory addresses (corresponding to macro-particle density and flux), whereas during the particle velocity advancement, thousands of particles only need to read the value of the same addresses (electric and magnetic field).

First result
We now focus on the result of Step 2, in which cometary ions were steadily added to the turbulent plasma of Step 1, moving at a super-Alfvénic and super-sonic speed. Table 3 lists the physical parameters used for Step 2. After 4000 iterations, the total number of cometary macro-particles in the simulation domain reaches a constant average value: the comet is fully developed and has reached an average "steady" state. From this point, we simulate more than one full injection period, 1500 iterations after which we loop over the domain of the injection tank in Figure 6. As an example, Figure 7 displays the state of the system at iteration 6000, focusing here again on the perpendicular fluctuations of the magnetic field. This time the colour scale is logarithmic, since magnetic field fluctuations are spanning over a much wider range than previously. While being advected through a dense cometary atmosphere, the solar wind magnetic field piles up (augmentation of its amplitude because of the slowing down of the total plasma bulk velocity) and drapes (deformation of its field lines due to the differential pile-up around the density profile of the coma), as first theorised by [34]. This general result was always applied to the global, average magnetic field, and was observed in situ at the various comets visited by space probes.
Without diving very deep in the first results of Menura, we see that the pile-up and the draping of upstream perpendicular magnetic field fluctuations also has an important impact on the tail of the comet, with the creation of large amplitude magnetic field vortices of medium and small size. This phenomenon, together with a deeper exploration of the impact of solar wind turbulence on the physics of a comet, are gathered in a subsequent publication.

Conclusion
This publication introduces a new hybrid PIC plasma solver, Menura, that allows for the first time to investigate the impact of a turbulent plasma u 0 363 km/s (6.66 v A ) Q 5 · 10 26 s −1 ν i 2 · 10 −7 s −1 u 0 1 km/s  flow on an obstacle. For this purpose, a new 2-step simulation approach has been developed which consist in, first, developing a turbulent plasma, and second, injecting it periodically in a box containing an obstacle. The model has been validated with respect to well-known fluid and kinetic plasma phenomena. Menura has also proven to provide the results expected at the output of this first step of the model -namely decaying magnetised plasma turbulence.
Until now, all planetary science oriented simulations have ignored alltogether the turbulent nature of the solar wind plasma flow, in terms of structures and in terms of energy. Menura has been design to fulfill this lack and it will now allow us to explore, for the first time, some fundamental questions that have remained open regarding the impact of the solar wind on different solar system objects, such as: what happens to turbulent magnetic field structures when it impacts on an obstacle such as a magnetosphere? How are the properties of a turbulent plasma flow reset as is crosses a shock, such as the solar wind crossing a planetary bow shock? How do the additional energy, stored in the perpendicular magnetic and velocity field components impact the large-scale structures and dynamics of planetary magnetospheres?
On top of the study of the interaction between the turbulent solar wind and solar system obstacles, we are confident that the new modeling framework developed in this work, in particular its 2-step approach might as well be suitable for the study of energetic solar wind phenomenons, namely Corotating Interaction Regions and Coronal Mass Ejections, which could be similarly simulated first in the absence of an obstacle, to then be used as inputs of a second step including obstacles.

Acknowledgements
E. Behar acknowledges support from Swedish National Research Council, Grant 2019-06289. This research was conducted using computational resources provided by the Swedish National Infrastructure for Computing (SNIC), Project SNIC 2020/5-290 and SNIC 2021/22-24 at the High Performance Computing Center North (HPC2N), Umeå University, Sweden. This work was granted access to the HPC resources of IDRIS under the allocation 2021-AP010412309 made by GENCI. Work at LPC2E and Lagrange was partly funded by CNES.

Appendix A. Normalised equations
In Menura's solver, all variables are normalised using the background magnetic field amplitude B 0 and background density n 0 , or equivalently using the corresponding proton gyrofrequency ω ci0 and Alfvén speed v A0 . We report the table of the background variables definitions in Table A.4. Based on these definitions, one can derive the following main equations of the solver. A normalised variableã is obtained by dividing this variable by its background value,ã = a/a 0 , equivalently a =ã a 0 . We first consider the Faraday's law, and using the background parameters definitions of Table  A.4: In other words, the Faraday's law expressed with normalised variables is unchanged compared to its SI definition. The Ohm's law becomes: with w spec =ñ spec /particle-per-node spec . For the solar wind proton,ñ = 1 and one simply gets w sw = 1/particle-per-node . W (r p ) stands for the shape factor, triangular in our case (in 2 spatial dimensions, one macro-particle affects the density and current of nine grid nodes, with linear weights).