© The Authors 2023

1 Introduction

The solar wind, a supersonic radially expanding plasma escaping from the Sun, can be described with two levels of complexity. The first level, by considering its average background values, provides a global laminar picture in which most of the seminal studies on the heliosphere and planetary magnetospheres were founded (Biermann 1951; Alf ven 1957; Parker 1958; Dungey 1961). The second level of complexity introduces the turbulent nature of the flow, a combination of chaotic fluctuations within the magnetic field and particle density and velocity, which adds up to their background average values (Bruno & Carbone 2005). In this phenomenology, energy (magnetic and kinetic) cascades from large to small scales, much like described by its neutral fluid analogue (Kolmogorov 1941), corresponding to a spectrum of fluctuations ranging over several decades of temporal and spatial scales. Eventually, this cascade leads to the dissipation of the energy at the smallest scales involved. Turbulence is suggested to play a key role in the acceleration of the solar wind, providing a continuous heating of the plasma from the solar corona and beyond (Cranmer et al. 2015).

The overwhelming majority of our knowledge regarding the formation and dynamics of planetary magnetospheres is based on the laminar description of the solar wind. Specifically, all global numerical simulations of these interactions involve a steady homogeneous plasma flow upstream of the obstacle (Schunk & Nagy 2009; Ma et al. 2008; Kallio et al. 2012). A few recent exceptions can be pointed out, for instance the magnetohydrodynamic (MHD) simulation using time-varying upstream conditions at the Earth (Lakka et al. 2017); two studies of a coronal mass ejection (CME) interacting with the Earth, using either a MHD model (Lakka et al. 2019) or a hybrid particle-in-cell (PIC) model (Moissard et al. 2022); and the hybrid PIC simulation of the effect of a pivoting magnetic field upstream of Mars (Romanelli et al. 2019).

However, a growing interest in the relationship between turbulence and magnetospheres is emerging. In addition to the plethora of publications focusing on turbulence within the Earth’s magnetosheath (Rakhmanova et al. 2021 and references therein), several studies have appeared on the topic of turbulence in the magnetospheres of outer planets and comets (Saur 2021; Ruhunusiri et al. 2020). Turbulence within the geomagnetic tail is also at the centre of many investigations reviewed by Antonova & Stepanova (2021) and El-Alaoui et al. (2021). More specifically, the effect of upstream turbulence on the terrestrial magnetosphere and its dynamics have been investigated for decades, with results reviewed by D’Amicis et al. (2020) and Guio & Pécseli (2021). This sizable quantity of literature outlines one main lacking aspect: a numerical tool for the global simulation of these turbulent interactions. This is where the recently developed code Menura (Behar et al. 2022) positions itself, allowing the injection of a fully turbulent flow upstream of an obstacle.

This publication presents the first application of the code, and focuses on a cometary magnetosphere characterised by a neutral outgassing typical for a heliocentric distance of about two astronomical units (au). Since the dawn of solar and space physics, comets have been emblematic tracers of the solar wind, hinting at its very existence (Biermann 1951; Alfven 1957), while interplanetary sector boundaries as well as CMEs were notably analysed using remote observations based on the intermittent disconnections of comet tails (Niedner & Brandt 1978; Vourlidas et al. 2007). More recently, comets were used once more to trace some of the solar wind’s turbulent parameters (DeForest et al. 2015) as well as its speed (Cheng et al. 2022). In addition to this historic role, a great amount of knowledge on cometary magnetospheres was produced during the last decade in the context of the European mission Rosetta (Goetz 2022). This provides us with a solid background for a first global exploration of such a turbulent interaction.

Because of the multi-scale nature of turbulence, its numerical simulation is intrinsically expensive. For this first study, we made the choice of properly resolving a wide range of scales, from the magnetospheric to the ion scales, below their inertial length. To allow a practical handling of this problem, and specifically to quickly iterate the simulation and its analysis, we chose to work in a 2D spatial domain, with velocities and field components described in a 3D space. This inherently limits the significance of our findings. To this extent, the aims of this first study are not to bring definitive results, but to properly illustrate the capacity of this new numerical approach, to give a first example on how its products can be tackled, and most importantly to highlight new aspects of the interaction, to be later on verified by a 3D approach, in the limit of realistic computing resources. It should not be underestimated, however, that the best resolution achieved by a 2D approach, and therefore the corresponding plasma mechanisms, cannot be verified by a 3D approach with an equal computational power. Therefore a 2D approach cannot be limited to a role of pathfinder; it may very well demonstrate mechanisms otherwise not reproducible.

This first publication focuses on the effect of turbulence on the obstacle itself, looking into scales around and higher than the ion scale, while the characteristics of the turbulent flow within the obstacle is left for future studies.

2 The model

The numerical model used to investigate the interaction of the solar wind (either laminar or turbulent) with a planetary obstacle is described and tested in Behar et al. (2022). It is based on a hybrid PIC model implementation of the Vlasov–Maxwell equations, including a source term for the distribution function: the ions are described as massive particles, the electrons are considered as massless and charge-neutralising, while the electromagnetic fields together with the particle moments are gathered at the nodes of a regular grid. At the core of the model is a generalised Ohm’s law that computes the electric field provided the magnetic field and the particle moments. Menura uses the following formulation: (1)

Here u_i is the ion bulk velocity, n the plasma density, J the charge current, and p_e the electron pressure. Additionally, a term of hyper-resistivity is used to dampen small-scale numerical oscillations, with the coefficient η_h multiplying the Laplacian of the current. Through Faraday’s law, this corresponds to a diffusion term. The value of η_h is taken to be 2.5 × 10⁻⁴ in the entire study. The electron pressure is obtained assuming it results from a polytropic process, with an index of 1 used throughout the study, corresponding to an isothermal process.

The code uses normalised units. Distances are expressed in units of the background proton inertial length d_i0. Time is expressed in units of the inverse of the background ion cyclotron frequency .

The model is based on a two-step procedure. First, a turbulent flow is generated, in the absence of the obstacle, as further described in Sect. 4. Second, this turbulent solar wind is injected in a simulation domain containing an obstacle. Since Step 1 is solved in a fully periodic domain, the injection of Step 1 outputs within the domain of Step 2 can be done periodically as well (see Behar et al. 2022 for more details). For such a medium activity comet, the domain boundaries parallel to the flow are also kept periodic.

During both steps, the equations and all variables are solved and expressed in the solar wind reference frame. It is therefore the obstacle that is moving through the solar wind. The domain is kept centred on the obstacle by regular copies and shifts of the fields and the particles, as illustrated in Behar et al. (2022). To that extent, the solar wind is not really injected in the domain but laid down in front of the moving object. Because the object reference frame is the one used in all planetary simulation codes we have encountered, the vocabulary unavoidably presents some ambiguities between the two frames. In the following sections we describe results and mechanisms in either the solar wind reference frame or the object reference frame, making sure to specify which.

During Step 2, in order to simulate a comet, a collection of cometary ions is added at each time step, as described in the next section.

In order to properly appreciate the influence of the incoming turbulent flow on the interaction, the model can also be used to send a laminar flow on the obstacle, all other parameters kept equal. We describe the laminar case in Sect. 3 before considering the turbulent case in Sects. 4 and 5.

3 The obstacle: a comet

Without gravity, an intrinsic magnetic field, or a solid central body (the size of the nucleus being negligible with respect to the dynamical scales of the system), and without ion collisions, the numerical modelling of such an obstacle is fairly straightforward. Cometary ions are introduced at each iteration according to a given ion production rate (2)

with r the distance from the comet nucleus, v_i the ionisation rate of cometary neutral molecules, n₀ the neutral cometary density, Q the neutral outgassing rate, and u₀ the radial expansion speed of the escaping neutral cometary atmosphere. At the beginning of the run, the number α of particles to be added in the proximity of each grid node is first calculated using Eq. (2) (together with the duration of one iteration). At each time step, we inject a constant number of particles with random positions within the cell surrounding each grid node, a number given by the integer part of α. At each iteration we produce and compare one more random number between 0 and 1 to the remaining part of α to decide whether one additional particle is created or not. In the sub-region surrounding the centre of the comet, where the distribution of neutral molecules shows the highest (radial) derivative, it is necessary to use a finer sub-grid to estimate these α values, and then average them over the main grid nodes closest to the nucleus in order to not underestimate the local creation. We use a sub-grid ten times finer than the main grid, after having verified that no significant change using an even finer sub-grid can be found.

The highest cometary ion density is found close to the nucleus, a region within which noticeable plasma structures appear; this is the interaction region we simulate. At 2 au from the Sun, these plasma structures and boundaries appear over the ion kinetic scales, characterised by the cometary ion gyro-radius, which will be discussed in further detail in Sect. 5. This interaction region is sketched in Fig. 1. Upstream of the nucleus, few newly born cometary ions are picked up by the solar wind electric field and start their cycloid motion. Closer to the nucleus, as their density is much higher, they form a noticeable density structure, which we hereafter refer to as the ‘pick-up plume channel’ (analogous to the pick-up plume at Mars; Dong et al. 2015), the first significant escape channel for cometary ions. This early phase of the gyration is represented in Fig. 1.

As the solar wind permeates the ionised cometary atmosphere, the total bulk velocity of the plasma decreases, while the density of cometary ions increases. The frozen-in magnetic field piles up on the denser coma, and its amplitude increases. Because the coma is denser close to the nucleus than elsewhere, the magnetic field additionally drapes around the nucleus, in the iconic shape established by Alfven (1957). Close to the nucleus, the magnetic field strength and its distortion become so intense that eventually, through the Hall component of the electric field (given by the local curl of the magnetic field under the Darwin approximation ∂_tE << J), the dense inner coma is accelerated downstream, which presents a second escape channel for cometary ions, also indicated in Fig. 1 (Behar et al. 2018b). In the following, we refer to this cometary ion escape channel as the ‘Hall escape channel’.

In the schematic, the cometary pick-up ions are accelerated upwards and, as a result of momentum conservation, solar wind ions are deflected downwards. This kinetic effect results in the formation of a solar wind overdensity highly asymmetric further away from the Sun (Behar et al. 2018a), which transitions to a more symmetric structure at smaller heliocentric distances (Hansen et al. 2007). Together with the structures formed by the pick-up plume and the Hall escape channels, this solar wind over-density is the third main cometary ion structure generated by the solar wind–comet interaction at such a heliocentric distance. Each is discussed in the rest of this study to diagnose the effect of upstream solar wind turbulence on the plasma environment of the obstacle.

Fig. 1 Schematic of a medium activity comet. Shown are the asymmetric dynamics of the solar wind ions. The two major escape channels for cometary ions are indicated.

4 The incoming flow: a turbulent solar wind

During the first of the two simulation steps, a turbulent plasma is obtained by letting the energy of initial perturbations cascade from large to small scales. Sine mode perturbations are initialised at time 0 in the in-plane magnetic field and in the velocity field, in addition to a guiding magnetic field B₀ purely out of the simulation plane. All particles are created with velocities following a Maxwellian distribution, using a thermal speed equal to the Alfvén speed. At time 500 (corresponding to a physical time of 1740 s with the values of Table 1), the turbulence has developed into the omni-directional power spectral density shown in Fig. 2, defined and used in the previous studies of Franci et al. (2015) and Behar et al. (2022). The spectrum displays a Kolmogorov-like scaling over the inertial MHD scales, following the black guide line with slope −5/3, to then adopt a much steeper slope at ion kinetic scales, similarly to the results of Franci et al. (2015) for a very similar simulation. At high spatial frequencies, a flattening of the spectrum is found corresponding to an energy range in which the noise of the particles adds up to the cascading energy. At the highest frequencies, we find a sharp increase in energy due to the noise of the finite differences used by the algorithm, a feature shared with the results of Franci et al. (2015).

The background values of this run are given in Table 1. The simulation domain is 500 d_i0 in size, corresponding to 65 733 km, and a regular grid with 2000 × 2000 nodes is used, corresponding to a grid spacing ∆x of 0.25d_i0. The resolved wave vectors are consequently within . The initial perturbations are injected with wave vectors in the range , and only the remaining non-perturbed spatial scales are shown in Fig. 2. This simulation uses 4000 particles per grid node (equivalently cell). A time step is used during this first step, while a two times finer time resolution is needed to properly solve the physics of the tail during Step 2 (see next sections).

This turbulent plasma has an average out-of-plane magnetic field of precisely < B_z >= 1, but an average total magnitude of < В >= 1.06. Compared to its laminar analogue, which is defined to be out-of-plane with amplitude 1 everywhere in the domain, we find a magnetic energy density (or magnetic pressure) 12% larger in the turbulent plasma, due to its additional in-plane component.

Fig. 2 Final state of the turbulent cascade. The left panel provides a map of the perpendicular (or in-plane) magnetic field fluctuations squared. The right panel shows their power spectral density, providing a guideline with slope −5/3.

5 General comparison

During the first time interval of Step 2, cometary ions steadily increase in number as the magnetosphere develops, before reaching a pseudo steady-state where the total number of particles in the domain evolve around a constant value. After dynamical equilibrium is reached, we further simulate the interaction for an additional .

The densities of the two species are displayed in Fig. 3, showing one snapshot when dynamical equilibrium is reached, for both the laminar and the turbulent upstream conditions. The solar wind overdensity position as well as the pick-up plume taken in the laminar case (left column) are reported in the turbulent results using dashed lines. We find that both the solar wind overdensity and the plume are reduced in size in the turbulent case. On average, the ‘nose’ of the overdensity is 4 d_i0 further upstream with laminar upstream conditions. As for the plume, we can estimate the gyration of the ions using their simple upstream gyroradius (the gyroradius a cometary test particle would have in the upstream wind). For the laminar conditions, with a homogeneous magnetic field, the value is R_laminar = 180 d_i0 everywhere in the domain. However, in the turbulent case the magnetic field now has an additional in-plane component that needs to be accounted for when calculating the gyroradius, which involves the ion velocity component perpendicular to the magnetic field. In addition, the amplitude of the magnetic field is also larger on average in the domain, as described in the previous section. We can compute the local value of the gyroradius, and find an average value over the domain of R_turbulent = 147 d_i0 that is significantly smaller than the laminar value. The lower row of Fig. 3 uses a cycloid of radius 180 (blue dashed line), which shows a good match with the laminar plume, while the turbulent plume is found to be smaller.

Whether we look at this interaction from a kinetic point of view (Behar et al. 2018a), which involves these gyration scales, or a fluid point of view, which involves the upstream magnetic pressure, we find in both cases that the plasma structures around such a comet are expected to be smaller in the turbulent case, which we indeed verified with the present simulations. It should be noted, however, that the definition of the laminar plasma is a choice, arbitrary to some extent, and one may argue that it could very well be defined in such a way that the magnetic energy density or the gyroradius are equal in both the laminar and the turbulent case.

It should be noted that the laminar interaction already shows some level of complexity, with some obvious wave patterns within the magnetosphere in the upper part of the cometary ion density (likely similar to the bi-ion acoustic waves also found by Bagdonat & Motschmann 2002), but also in the lower part of the solar wind density (see also the work of Koenders et al. 2016 on low frequency waves at comet 67P/CG). These fluctuations and their fate when considering a turbulent upstream input, as well as their potential contribution to the inner magnetosphere turbulence, is another important direction to explore in the future.

The Hall escape channel is shown stemming from the comet’s inner region in the lower panels of Fig. 3. In the laminar case (bottom left panel), the homogeneous incoming magnetised solar wind results in a continuous cometary ion acceleration, and in turn in a continuous cometary ion structure, forming the Hall escape channel. On the contrary, in the turbulent case (bottom right panel), the turbulent nature of the incoming magnetised solar wind is responsible for the generation of discontinuities in this same Hall escape channel, which is now found to be composed of discrete high density cometary ion bubbles of similar density to the inner coma.

The modifications of the solar wind turbulent flow are highlighted in Fig. 4, which shows the in-plane magnetic field lines, using a line integral convolution (LIC) representation (Loring et al. 2015). In the upper panel, the in-plane magnetic field lines are coloured using the amplitude of the same in-plane field, while the lower panel uses the density of cometary ions superimposed on the LIC representation. Cometary ion bubbles are found to be enclosed in magnetic islands downstream of the nucleus, which were not present in the upstream turbulent wind. We identify two main regions in this magnetosphere. First the ‘wings’, in which deformed upstream perpendicular field structures can be recognised. And second, the comet tail in which most of the cometary ions are confined, corresponding to very low solar wind densities in Fig. 3. Newly formed magnetic field islands of various sizes can be observed there, labelled ‘loops’ in Fig. 4. In the next two sections, we explore the origin of these two regions.

Fig. 3 Comparison of the turbulent and laminar interactions. The upper row shows the density of the solar wind ions, while the lower row presents the density of the cometary ions. The left column corresponds to laminar upstream conditions, while the right column corresponds to turbulent upstream conditions. The position of the solar wind overdensity in the laminar case is reported in the turbulent case with the dashed red line. Similarly, the estimated gyration of the cometary ions in the plume in the laminar case is shown in the lower row with a blue dashed line.

6 Disconnection of the comet’s head

The pile-up and draping of the magnetic field, which to date in the literature have been discussed in the context of a laminar interaction corresponding to a homogeneous upstream magnetic field, in turn also affect the turbulent structures of the solar wind magnetic field. To our knowledge, the consequences of cometary-induced pile-up and draping of the magnetic field on the incoming turbulent heterogeneous plasma flow have not yet been considered or investigated.

To illustrate the mechanism at the origin of the formation of these high density bubbles and the complex magnetic field structures within the comet’s tail, we designed a dedicated numerical experiment to isolate the interaction between a single upstream perpendicular magnetic field structure and the coma, in an otherwise laminar upstream flow. The laminar version of the interaction is resumed from its steady state shown in Fig. 3, and an ersatz of a magnetic island (i.e. a flux rope in the 3D case) is introduced upstream of the comet. Eventually, the comet meets and interacts with this artificial structure, just as it does with structures of the fully turbulent flow.

The magnetic field island ersatz is defined by adding a perpendicular component to the otherwise homogeneous and out-of-plane background magnetic field. This additional perpendicular component is characterised by circular in-plane field lines, and an amplitude depending on the distance to the centre of these circular lines, with a maximum value at some distance from the centre. We chose a Gaussian profile with a maximum value of 0.5, the background magnetic field at 50 d_i0 from the centre of the structure, and a standard deviation of 25 d_i0. This artificial structure is meant to mimic the large-scale vortices found in the perpendicular magnetic field shown on the left side of Fig. 2, with diameters about 50–100 d_i0, with a larger amplitude found away from their centres. This structure is indicated by a dashed circle in the left panel of Fig. 2. This ersatz is not meant to be a perfect replica of such structures, however. No further considerations are taken into account other than a divergence-free circular structure of about 100 d_i0 in diameter, with a realistic maximum found at some distance from its centre.

This ersatz can be seen still undisturbed upstream of the comet in the uppermost panel of Fig. 5, defining the relative time t = 0 for this experiment. The figure uses the same LIC representation as Fig. 4, with a threshold of 0.01 on the in-plane (perpendicular) magnetic field; if the perpendicular magnetic field is smaller than 0.01, no LIC is shown. The background colourmaps of the first four panels from the top, showing four successive and equally separated times of the experiment, additionally give the amplitude of the perpendicular magnetic field, while for the final time of the experiment, given in the bottom panel, the density of cometary ions is used.

The time t = 0 is chosen as the magnetic island is just starting to interact with the dense inner coma. As a result, the leftmost field lines are found to be deformed, draped inward the island, with a corresponding increase in the field amplitude (pile-up) seen in red tones at the left border of the structure. As the magnetic field amplitude increases and the length scale between the anti-parallel magnetic field lines shortens, a strong current sheet forms around the comet nucleus location (second panel). This strong current sheet is likely unstable to the tearing instability, eventually forming shorter-scale magnetic islands (third panel).

This series of snapshots gives us a very tangible representation of how the comet pierces through an upstream magnetic field structure, similar to a projectile through an obstacle. After the impact, we find two remaining pieces of the initial structure with high perpendicular magnetic field amplitudes looping around two distinct centres, one on each side of the comet’s head. Within these two wings, the upstream information is conserved to some extent, with field structures intensified and deformed. However, closer to the comet’s head and downstream of it, new smaller structures of even higher perpendicular magnetic field amplitude are produced, seen as smaller-scale closed magnetic field loops. The two wings are found to be more or less static in the solar wind reference frame, just as their parent upstream structure, whereas these smaller-scale intense loops have a speed closer to the comet itself. Equivalently, as described in the object reference frame, the wings are advected downstream at more or less the speed of the solar wind, whereas the new loops are transported downstream at a much lower speed.

The difference between these two types of structure originates from the density the upstream parent structure meets during its interaction with the coma. On the one hand, the sides, or wings, of the parent island interact with the comet in regions where the density is still dominated by solar wind protons in terms of number density, and the two wings keep frozen in the solar wind dominated flow, conserving some information from upstream. On the other hand, at the centre of the interaction region, the plasma is largely dominated by the cometary ions. The upstream magnetic field, now with the additional perpendicular component, is piling up in this drastically reduced plasma mean velocity, reaching higher amplitudes than in the wings. More interestingly, this rising magnetic tension through the Hall electric field eventually pulls off the dense inner coma, resulting in one main and two secondary high density bubbles found downstream of the comet’s head at time 56.25 (bottom panel) enclosed in the newly formed intense magnetic field loops. We note that at this time, taken much later than the first four snapshots, the wings are not carrying any significant cometary ion density and are by then long gone downstream of the comet. After this build-up and release phenomenon, and in the absence of additional heterogeneous perpendicular magnetic field structures, the inner coma resumes its laminar continuous Hall escape.

We now take a closer look at this disconnection event. We zoom in on the inner interaction region, and focus on times around t = 6.25.

Fig. 4 Global view of the magnetic field structures within the interaction region. The upper panel provides a line integral convolution of the in-plane magnetic field component, coupled to its colour-coded amplitude. The lower panel shows the same LIC representation, coupled to the cometary ion density along the comet tail. In the lower panel, lower-right corner, the horizontal line length corresponds to 100 di.

Fig. 5 Snapshots of the numerical experiment. All panels use a LIC representation of the in-plane magnetic field lines. The first four rows use the amplitude of the in-plane magnetic field for the colours, while the last panel uses the cometary ion density.

7 Disconnection at the ion scale

Figure 6 shows the comet’s head during the same experiment, zooming in a bit closer on the spatial domain than the previous representation, using a more classic, oriented field line representation superimposed on the cometary ion density. These three snapshots focus around the time when the centre of the upstream magnetic island passes through the inner coma, corresponding to the second panel of Fig. 5. In the upper panel the magnetic island can be identified, already significantly piled up (i.e. strongly compressed along the x-direction). As the plasma is faster above and below the inner coma (red tones on the colourmap), the island is additionally draped, and with the upper and lower parts of the island advected downstream much faster than the central part, an elongation of the structure is also seen. Field lines that were initially circular are now found to describe highly eccentric ellipse-like figures along a main axis, highlighted with a manually added red band describing a bow. Because of this combined compression and elongation along the red band, we find a line separating field lines parallel and of opposite sense. As the comet continues piercing through the structure, the red band gets significantly draped. In the second panel, multiple reconnection points along the red band have occurred, resulting in new, small-scale magnetic islands. At the very tip of the band, where the field keeps piling up, we find one main magnetic reconnection site on which anti-parallel field lines, perpendicular to the flow, keep being compressed against each other.

Under this magnetic tension, the inner ionised coma is slowing down in the plasma frame, or equivalently is pushed downstream in the object frame. This can be easily spotted by looking at the blue crosses in each panel, which indicate the centre of the comet. The tension eventually becomes strong enough to significantly accelerate the inner coma, which starts moving downstream of the comet’s centre (in the comet’s frame). As soon as this pair leaves the central region, the continuous ionisation of neutral molecules starts refilling the inner coma, with most of the magnetic tension released. In the last panel of Fig. 6, we can see the newly detached high density cloud enclosed in closed field lines, while upstream of it a fresh inner coma is getting denser, though still much less dense than the detached cometary plasma bubble.

The same situation can be seen in the fully turbulent interaction. Figure 7 shows one snapshot when a dense cometary ion cloud just got detached from the central region of the comet, enclosed in a magnetic island, while a series of smaller-scale islands (emphasised in the figure by the manually added red band for visualisation purposes) can be seen, which do not capture high densities of cometary ions. Although we did not design the upstream magnetic island to be a perfect replica of some structures of the turbulent flow, we were able to reproduce the density disconnection event in great detail, illustrating how the pile-up and draping of upstream perpendicular magnetic field structures can disrupt the Hall escape channel of the comet, producing cometary ion bubbles within the tail, together with new shorter-scale magnetic field structures not present upstream. Although our ersatz was perfectly circular, the mechanism happens in the same manner with more complex field structures, which do not even need to present closed loops as long as upstream anti-parallel magnetic field lines pile up in front of the inner coma to generate strong unstable current sheets in the same way, resulting in the global picture of Figs. 3 and 4.

Fig. 6 In-plane magnetic field lines, oriented, superimposed on the cometary ion density. The image was taken during the numerical experiment and centre around time t = 6.25 of Fig. 5. The red band was manually added to help follow the main vertical axis of the upstream structure. The blue crosses indicate the centre of the comet.

Fig. 7 In-plane magnetic field lines, oriented to overlap with the cometary ion density, taken during the fully turbulent simulation of the interaction. The manually added red band indicates a reconnecting current sheet, similar to that in Fig. 6. The larger blue cross indicates the centre of the comet, while the smaller cross shows the position of the synthetic probe used in Sect. 8.

8 Draining of the inner coma

This tension-release phenomenon in a fully turbulent flow leads to a pseudo-periodic creation of cometary ion bubbles. The effect of this dynamical flushing is most spectacular when following through time the plasma parameters within the inner coma (i.e. close to the nucleus). A synthetic probe was used within the simulation, static in the comet frame, with a position shown in Fig. 7 as a smaller cross almost superimposed on the centre of the comet. At the probe’s position, all plasma parameters (field values and particle moments) were recorded for a duration of three full periodic injections of the upstream turbulent flow. Figure 8 provides the time series of the ion density in the upper panel, comparing the turbulent and laminar runs, while the lower panel gives the perpendicular and out-of-plane magnetic field for the turbulent case.

The phenomenon previously illustrated in the spatial domain can also be recognised in these time series. Each sudden drop in density (upper red line) corresponding to the time when a bubble is disconnected from the inner coma and moves downstream, leads to a maximum in the perpendicular magnetic field (lower panel) that follows shortly after. Overall, we observe that the two quantities are phased by about π/2, with maxima (respectively minima) of the magnetic field found at a decreasing (respectively increasing) density, crossing its average value. This phasing results from the competing pressures (magnetic on the one hand, dynamic and thermal on the other), never reaching an equilibrium in this dynamic inner region. During this pseudo-periodic draining of the inner coma, the density reaches maximum values that are 147 times higher than the upstream mean density, to then drop to values of around 24. During these drops, the perpendicular magnetic field reaches values of 16 times the background magnetic field, with minimum values around 2, and locally even reaching 0. The out-of-plane magnetic field B_z is found to vary much less, around an average value of −4.

In these time series, recorded during the turbulent interaction, the periodicity corresponding to the cyclic turbulent flow injection is clearly seen. We can deduce from the upper panel density curves that one full injection of the turbulent plasma eventually leads to the formation of three main density bubbles.

The inner coma density recorded at the same point during the laminar interaction shows a much more stable behaviour, around an average value of 120 times the upstream density. The average value in the turbulent case is only 55% of that. Reformulating this number, in the limit of the model exposed above, we find an inner coma that is on average almost twice as thin when we add turbulent fluctuations to the upstream flow, with variations from 20% to 122%, compared to a laminar case with all background parameters equal.

Fig. 8 Time series of the total density (upper panel) and the in-plane and out-of-plane magnetic field components (lower panel), taken in the inner coma, given in normalised units. The horizontal lines of the upper panel correspond to averaged values of the time series.

9 Summary

This work demonstrates for the first time how the global picture of the interaction between an object and a turbulent solar wind, as provided by the dedicated code Menura, can be approached and analysed. On this first application, we have focused on the macroscopic effects of the upstream turbulence on the magnetosphere, especially studying the apparition of high density cometary ion bubbles within the comet tail, not present in the laminar case.

Using an isolated ersatz of a large-scale, upstream magnetic field structure, we have reproduced the two main characteristics found in the fully turbulent interaction. This illustrates that it is the large-scale turbulent structures that carry most of the turbulent energy, which may impact the magnetosphere the most. Another dedicated parameter study will be necessary to explore the topic further, using various initialisations of the turbulence. Based on these early results, however, an important aspect of such a future study would be the role of the amplitude of the initial fluctuations, as well as their range (i.e. including even larger k values for a larger simulation domain).

First, we see the creation of two different regions in the induced cometary magnetosphere, the wings, in which the plasma flow is disturbed to some extent, but keeps some of its upstream information and relative speed, and the tail, in which newly formed structures move downstream more slowly, showing a much more intense magnetic field. Second, we were able with this single structure to reproduce the disconnection of the comet’s head in terms of cometary ion density, leading to high density bubbles correlating to the new magnetic field structures within the tail, just as observed in the fully turbulent run.

The creation of these cometary ion bubbles goes hand in hand with the pseudo-periodic draining of the innermost coma, within which very large fluctuations in terms of density and electric and magnetic field amplitudes are found. These variations in terms of density range within 20% and 122% of the average laminar density. Another striking result is that the turbulent interaction results in a dramatically lower average density in the inner coma, only 55% of that of the laminar interaction.

The temporal evolution of the overdensity position is greater in the turbulent case (not shown in this publication). We could not yet find the right combination of parameters that could account for this precise evolution. The problem is indeed fairly complex, as both the upstream parameters and the inner-coma parameters are highly dynamic, involving different timescales.

The implications of such a reduced average density of the inner coma would be profound for all aspects of its physics. A first example is the balance between the various pressures, magnetic, ram, and thermal, each depending on the magnetic field and the plasma density. Changing so significantly the pressures, within the magnetosphere as well as upstream of it, necessarily and significantly moves the various plasma boundaries from laminar to turbulent upstream conditions. A dynamic input is also expected to introduce more dynamics along these boundaries. Another major topic concerned by this redefined plasma density is the physico-chemistry of the inner coma, with charge-exchange processes heavily dependent on the ion and electron density.

This first work paves the road to a comparative planetary science study for other types of bodies, such as ionospheres (Mars, Venus) and permanent dipoles (Earth, Mercury). The potential effect of upstream turbulence on the escape of planetary ions may introduce a new argument in the topic of the atmosphere evolution.

Though this 2D approach would benefit from a future back up of additional 3D runs, the basic phenomenon illustrated in this publication should remain valid, to some extent: upstream perpendicular magnetic field structures will still pile up and drape, self-consistently forming strong, unstable current sheets at the location where the cometary plasma is densest. Whether the effect on the planetary ions will be that spectacular is not obvious, however. The nature of the upstream turbulence, and especially the large-scale structures relevant for this study, are also quite different when solving the problem with a third dimension (see e.g. Franci et al. 2018, who highlight some differences between 2D and 3D cascades on MHD scales).

Another important direction to explore is the opposite problem, the effect of the magnetosphere on the solar wind turbulence. This is a somewhat more delicate question, which might need a more refined approach.

Acknowledgements

E. Behar acknowledges support from Swedish National Research Council, Grant 2019-06289. 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.

References

All Tables

All Figures

	Fig. 1 Schematic of a medium activity comet. Shown are the asymmetric dynamics of the solar wind ions. The two major escape channels for cometary ions are indicated.
In the text

	Fig. 2 Final state of the turbulent cascade. The left panel provides a map of the perpendicular (or in-plane) magnetic field fluctuations squared. The right panel shows their power spectral density, providing a guideline with slope −5/3.
In the text

Fig. 3 Comparison of the turbulent and laminar interactions. The upper row shows the density of the solar wind ions, while the lower row presents the density of the cometary ions. The left column corresponds to laminar upstream conditions, while the right column corresponds to turbulent upstream conditions. The position of the solar wind overdensity in the laminar case is reported in the turbulent case with the dashed red line. Similarly, the estimated gyration of the cometary ions in the plume in the laminar case is shown in the lower row with a blue dashed line.

In the text

Fig. 4 Global view of the magnetic field structures within the interaction region. The upper panel provides a line integral convolution of the in-plane magnetic field component, coupled to its colour-coded amplitude. The lower panel shows the same LIC representation, coupled to the cometary ion density along the comet tail. In the lower panel, lower-right corner, the horizontal line length corresponds to 100 di.

In the text

	Fig. 5 Snapshots of the numerical experiment. All panels use a LIC representation of the in-plane magnetic field lines. The first four rows use the amplitude of the in-plane magnetic field for the colours, while the last panel uses the cometary ion density.
In the text

	Fig. 6 In-plane magnetic field lines, oriented, superimposed on the cometary ion density. The image was taken during the numerical experiment and centre around time t = 6.25 of Fig. 5. The red band was manually added to help follow the main vertical axis of the upstream structure. The blue crosses indicate the centre of the comet.
In the text

	Fig. 7 In-plane magnetic field lines, oriented to overlap with the cometary ion density, taken during the fully turbulent simulation of the interaction. The manually added red band indicates a reconnecting current sheet, similar to that in Fig. 6. The larger blue cross indicates the centre of the comet, while the smaller cross shows the position of the synthetic probe used in Sect. 8.
In the text

	Fig. 8 Time series of the total density (upper panel) and the in-plane and out-of-plane magnetic field components (lower panel), taken in the inner coma, given in normalised units. The horizontal lines of the upper panel correspond to averaged values of the time series.
In the text