© ESO 2019

1 Introduction and background

The European Space Agency mission Rosetta to comet 67P/Churyumov-Gerasimenko (67P) escorted the comet around the Sun between August 2014 and September 2016, from the approximate heliocentric distance of 4 astronomical units (AU) to perihelion (1.24 AU) and back out again to the outer solar system (Taylor et al. 2017).

During this period, attempts to detect the bow shock were mostly unsuccessful because the spacecraft most of the time orbited in the inner coma (see e.g. Simon Wedlund et al. 2017). Recently, Gunell et al. (2018) reported measurements from the Rosetta Plasma Consortium (RPC) instrument suite onboard Rosetta, which for the first time showed evidence that the cometary bow shock was in a formation process at a heliocentric distance of about 2.3 AU. The boundary crossing showed a simultaneous increase in magnetic field fluctuations and in electron and proton heating, and a decrease in solar wind fluxes downstream of the boundary that was characteristic of a shock-like structure, all of which was also seen in their 3D hybrid plasma simulation results.

Moreover, as shown by Nilsson et al. (2018) with the RPC Ion Composition Analyzer (RPC-ICA), the observed pickup ion energy spectrum from a few eV up to 20 keV has several characteristics that are not associated with ion pickup in an undisturbed solar wind. RPC-ICA onboard Rosetta is a top-hat ion mass spectrometer capable of characterizing the angular and energy dependence of the solar wind and cometary ion particles (Nilsson et al. 2015a,b). For a heliocentric distance of 1.4 AU, Nilsson et al. (2018) showed the water-group ion source density region as a function of energy of the pickup ions, which displays a strong discontinuity in slope occurring around 1 keV. According to Nilsson et al. (2018) and their simple cloud model, this knee in the distribution corresponds to a cometocentric distance of about 4000 km, suggesting the presence of a physical boundary of unknown nature that the high-energy pickup ions, originating from far upstream of the spacecraft, crossed on their way to the inner coma. The nature of this boundary is debated; it could be related to the possible presence of a cometopause-like structure, where cometary ion densities start dominating the solar wind ions, or to the bow shock and associated cometary plasma sheath downstream of it. Models show that at 1.4 AU both structures should already be well developed. Concurrently with the comet’s increasing outgassing, a solar wind ion cavity has also been observed to form in the inner coma (Behar et al. 2017): this cavity corresponds to the disappearance of solar wind ions below about 1.7 AU and confirms the formation of large-scale structures upstream of the spacecraft.

The cometary plasma environment is shaped by the cometary outgassing and the subsequent ionization of the resulting coma. Extending millions of kilometers from the nucleus of the comet, the neutral coma gives rise to mass-loading of the solar wind through momentum transfer and charge exchange processes, which need to be taken into account in global simulations of the cometary plasma environment (Koenders et al. 2013).

Simon Wedlund et al. (2017), subsequently referred to as Paper I, introduced a new quasi-neutral hybrid plasma model developed at Aalto University that is based on earlier approaches at, for example, Mercury, Venus, and Mars (Kallio & Janhunen 2002; Jarvinen et al. 2013; Alho et al. 2015, respectively). We present here a follow-up application of this model to the study of cometary pickup ions. After a summary description of the model, we emphasize two recent developments we have implemented since our first model: an updated resistivity scheme for multiply refined grids, and the addition of an upstream mass-loading condition that takes the extended neutral coma into account.

The goal of this study is to use numerical hybrid plasma modeling to characterize and interpret pickup ion energy spectra that are typically observed in the vicinity of the nucleus where Rosetta was located, in terms of large-scale structures of the cometary plasma environment, using an idealized Haser model for the comet. We explore the parameter space in interplanetary magnetic field (IMF) magnitude, while other parameters are left for further investigation. We perform several simulations with varying solar wind magnetic fields. The kinetic description of ions in the model is then used to extract the velocity distribution function of pickup ions, and is then used to interpret spectral features of the ions. Finally, in this context, we discuss the original observations of the source density (Nilsson et al. 2018) and implications in terms of large-scale structures upstream of the spacecraft.

2 Model description

We used the3D cometary quasi-neutral hybrid (QNH) cloud-in-cell (CIC) architecture developed at Aalto University that is described in detail in Paper I. Only a brief summary is given here.

As described in Paper I, the model treats ions as kinetic macroparticles, whereas electrons are a massless isothermal fluid at temperature T_e moving at bulk velocity U_e. Each macroparticle in the model represents a large number of physical particles, denoted by the statistical weight w_i. The macroparticle weight can vary due to macroparticle splitting, joining, and particle processes (see below). Macroparticle splitting and joining adjust the macroparticle count in a cell toward a set number that varies per cell and per population to limit the shot noise due to low statistics, and, on the other hand, to constrain macroparticle counts with regard to computational requirements. Electric and magnetic fields are solved self-consistently, and the model naturally includes the Hall term. Particle processes such as photoionization (PI), charge exchange (CX), electron ionization (EI), and electron recombination (ER) are included and can be individually parameterized. Ion-neutral collisions are included through a Langevin drag term in the ion equation of motion. The simulation domain is a box, for which a Cartesian coordinate system is chosen to coincide with the cometocentric solar equatorial coordinates (CSEQ), so that the comet nucleus is at the origin, with the X-axis pointing toward the Sun at (+∞, 0, 0), the Z-axis pointing along the solar rotation axis perpendicular to the X-axis, and the Y -axis completing the right-handed reference frame. Cartesian grid refinements in the X, Y, and Z directions can also be used to study a specific region with higher spatial accuracy (see Fig. 1A).

A description of the model inputs and updates is given in the next sections, with emphasis on the salient modifications that we have carried out since the first study in Paper I: a new resistivity scheme, the inclusion of upstream mass-loading on the subsolar walls of the simulation box, and upstream pickup ions. The virtual instrument for producing simulated observations is also discussed.

2.1 Model inputs

Table 1 presents the physical input parameters of the model (solar wind and neutral atmosphere) as well as the simulation parameters (size of the simulation box and the spatial and temporal resolution) shared by our simulations. For this study, we performed four simulation runs with different upstream IMF values using a 90° Parker angle, shown in Table 2. Other parameters such as the outgassing rate were held constant between simulations. We refer to the simulations as Runs 1 through 4, in order of increasing IMF magnitude.

The simulation grid was configured so that the maximum resolution grid covers the shock area, which is determined a priori from preliminary low-resolution simulations. The maximum grid resolution of 250 km is comparable to previous modeling work of cometary (Paper I) and planetary (Jarvinen et al. 2013; Alho et al. 2015) environments that produce bow shocks, and this resolution is taken to be sufficient for modeling of heavy cometary ions in the vicinity of the bow shock.

2.1.1 Cometary outgassing

We considered a neutral atmosphere composed solely of water H₂O, which was the dominating species at comet 67P in April 2015 (Fougere et al. 2016). The cometary neutral density was calculated using a Haser-like profile, which assumes a spherically symmetric neutral outgassing at radial speed v_n for species s (Haser 1957), (1)

with the cometocentric distance r and λ = v₀∕f_PI the characteristic photoionization scale length, a function of the photoionization rate f_PI (Huebner & Mukherjee 2015) and the neutral outgassing speed v₀. In our simulations, v₀ is equal to 700 m s⁻¹, which corresponds to an upper estimate as observed by Hansen et al. (2016) at comet 67P. The density also depends on the water outgassing rate Q₀ (in s⁻¹), given by the pre-perihelion fit from all Rosetta datasets of Hansen et al. (2016) to be (2)

with d_⊙ the heliocentric distance in AU. The water outgassing rate at 1.425 AU is then Q₀ = 4.14 × 10²⁷ s⁻¹.

Neutral outgassing profiles, such as those presented by Hansen et al. (2016), can be used instead of the Haser model. For simplicity, this paper describes the interaction in terms of the symmetric Haser model because our main focus is to arrive at a qualitative description of the spectral features of the pickup ion energy.

Fig. 1 General overview of Run 2. Panel A: magnetosonic Mach number = 2 and cometopause surfaces (orange and red, respectively) and the electron fluid-velocity streamlines, overlaid with the simulation grid displaying the grid refinements. On the X =−1000 km plane, the magnetic field magnitude is displayed in blue. Panel B: magnetic field lines connected to the comet-Sun line, viewed along the direction of the convective electric field. Features of the environment listed in the main text are annotated in the panels. Panel C: values of the magnetic field magnitude B and solar wind H+ bulk values (number density , bulk velocity and thermal velocity vth) relative to upstream values, along the X-axis from X =−1000 to X = 9000 km. Simulation cometopause (red), bow shock standoff (blue) distances, and the extent of the diamagnetic cavity boundary observations(green) of Rosetta up to July 2015 (Goetz et al. 2016) are shown by the shaded areas for reference. Four regions of the environment are numbered in the panels. 1: gradual mass-loading and some draping (panels A–C), 2: a sharp shock-like structure (A–C), 3: magnetic field pile-up approaching the nucleus (B and C), and 4: draping of the magnetic field lines around the inner coma (B).

2.1.2 Solar wind

The chosen solar wind and IMF parameters are also given in Table 1. For this study, we considered the undisturbed solar wind, that is, the solar wind infinitely far upstream of the comet, to be a constant stream of H⁺ and e⁻ traveling at a bulk speed of U_sw = 430 km s⁻¹. For simplicity in the interpretation of the modeling results and in contrast to Paper I, He²⁺ ions are not included in this study. The solar wind proton density and temperature and IMF were given nominal values derived from Slavin & Holzer (1981) at 1 AU (n_sw = 7 × 10⁶ m⁻³, |B| = 6 nT) scaled to the comet’s heliocentric distance. The distribution of H⁺ was set to a drifting Maxwellian with the aforementioned parameters. With a Parker spiral angle ϕ at 1.425 AU of 54.1°, the nominalIMF B_y component has a magnitude of 2.95 nT. To facilitate interpretation, the B_x and B_z components are set to zero throughout this study, which prevents asymmetry from the flow-aligned magnetic field component.

With the aim of characterizing the effect of the upstream magnetic field and to investigate how it affects the pickup ion spectra, four runs with different IMF magnitudes were successively performed, with the values given in Table 2. Run 3 corresponds to the nominal value (B_y = 2.95 nT) at1.4 AU, whereas Runs 2 and 1 have lower magnetic field magnitudes (B_y= 1.96 nT and B_y = 0.98 nT, or 2∕3 and 1∕3 of Run 3 B_y), similar to the magnitudes that are generally expected at 2.2 and 4.5 AU, respectively. For Run 4, the full magnitude of the nominal Parker IMF was taken for the transverse component, so that B_y = 3.64 nT. The change from high to low IMF values corresponds to a weakening coupling between pickup ions and the solar wind flow, as the pickup ion Larmor radius grows to values substantially higher than the effective environment of the comet. As a consequence, from Run 1 to Run 3, and as discussed later in Sect. 3, the bow shock structure will experience a gradual evolution from not quite a bow shock (Run 1, with a sharp but not yet uniformly developed structure; a “caustic”) to a full-fledged bow shock (Run 3, more similar to the bow shock structure at perihelion). The evolution may also be seen to mimic the birth and evolution of such a structure as the comet approaches the Sun.

2.1.3 Processes

Particle processes such as PI, CX, EI, and ER were included in the simulations in terms of reaction rate constants and cross sections between solar wind particles and the neutral atmosphere (H₂O). Consequently, the cometary ion population is composed solely of H₂O⁺ ions produced by PI, CX, and EI processes. For an in-depth discussion of these individual processes, how to implement them in a QNH approach, how they affect the formation of the bow shock and the cometopause, we refer to previous work in Paper I and references therein.

As shown in our previous study, CX is the main driver that regulates the shape and extent of the bow shock, whereas PI and EI processes dominate the plasma dynamics of the inner coma. This conclusion is in accordance with measurements made with Rosetta (Heritier et al. 2017). Energy-dependent CX cross sections for the single capture of H⁺ in H₂O were taken from Simon Wedlund et al. (2017). Electron impact ionization and recombination rates were calculated for an electron temperature T_e = 133 300 K using the parameterizations of Cravens et al. (1987) and Hollenbach et al. (2012), respectively. Compared to the 1.3 AU case chosen in Paper I (calculated for T_e = 300 000 K to mimic heating at the bow shock), the EI rate in our study here is lower by about 65%.

2.2 Hybrid model: new features

In Paper I, upcoming and desirable simulation features were briefly mentioned. We present here two such additions, one for a grid-scale resistivity scheme, and the other for mass-loading upstream conditions. Both are important for the interpretation of the Rosetta observations and will help define the bow shock-like structure and its energy signature.

2.2.1 Local grid-size adapted explicit resistivity

Numerical stability of the model was ensured by a predefined explicit resistivity term η, as in Paper I, but with a modified resistivity scheme. Resistivity is introduced to the model through generalized Ohm’s law: (3)

In earlier versions, such as in Paper I, the resistivity term was constant in the entire domain, leading to relatively high diffusion in multiply refined grid cells. In this version, resistivity is scaled in accordance with the size of the grid cell, so that the magnetic Reynolds number R_m = μ₀UL∕η is constant. The resistivity term η(r) at a point r is thus given by (4)

where we take U = U_sw, that is, the upstream solar wind speed and L(r) = Δx(r), that is, the local grid resolution, and R_m is chosen to ensure stability. This constrains the diffusivity across small grid cells while ensuring numerical stability similar to the coarse grid. In comparison, the similar A.I.K.E.F. model uses a combination of a resistive term in Ohm’s law and a direct smoothing of the magnetic field (Müller et al. 2011).

2.2.2 Extension to upstream mass-loading and pickup ions

The neutral atmosphere of a comet and its associated ion source region expands freely in space and can at perihelion extend millions of kilometers from the comet. As momentum and energy transfer between the solar wind and the cometary ions (mass-loading) will take place over these large scales, their cumulative effect on the solar wind plasma bulk speed and direction may becomeimportant when the inner cometary environment is reached. This is a challenge for any numerical simulation, and especially hybrid and full particle-in-cell models, which typically use a constrained simulation domain to conserve computational resources.

In order to better account for the extended cometary atmosphere and its cumulative effect on solar wind through gradual mass-loading, we have developed an ad hoc extension to the model inflow boundary. Koenders et al. (2013) have employed similar extended boundary conditions to take these upstream effects into account. However, their approach does not capture the finite gyroradius effects of the cometary pickup ions. These are notable sources of asymmetry below the length scales of heavy ion gyromotion for the plasma environment and pickup ion fluxes.

Our 3D upstream mass-loading extension is based on an approximation that takes the finite gyroradius character of the solution into account and is somewhat similar to the works of Behar et al. (2018) and Saillenfest et al. (2018). Our approach is based on two approximations in the upstream mass-loading region outside of the nominal simulation box:

1. Constant upstream solar wind. Upstream of the solar wind inflow boundary, the solar wind parameters are assumed constant in space and time and equal to the undisturbed solar wind. This assumption is used to (a) ascribe exact cycloidal motion to upstream pickup ions, neglecting any scattering processes, and (b) ascribe a ring distribution function f(θ, r)dθdr to the pickup ions at each point r, with θ describing the phase of gyration, thus parameterizing the velocity space part of the distribution. For each point r and gyrophase θ, we can then analytically find the points of ionization for these pickup ions by following their cycloidal trajectory. With a given source distribution, this allows us to evaluate the pickup ion density function at each point (θ, r).

2. Loss of momentum from solar wind is accounted for at the simulation inflow boundary, in accordance with the assumption of constant upstream flow. Using the above evaluation of the pickup ion density function, we can calculate the Lorentz force that is imparted upon the pickup ions from the solar wind flow at each point. For each solar wind parcel, we can then integrate the transfer of momentum to pickup ions from far upstream up to the inflow boundary, with “far upstream” taken to be 10⁷ km. The accumulated change of momentum is then affected upon the solar wind particles at the front wall.

Using this approach, we produced a deflection of the solar wind protons of approximately 5° at the inflow boundary, which is in good agreement with that produced by multifluid magnetohydrodynamics (MHD) modeling (Rubin et al. 2015) at comparable heliocentric distances. Additionally, the produced solar wind deflection and deceleration show spatial variability that is due to the upstream pickup ion distribution. The version of energy and momentum transfer used in our simulations is expected to slightly underestimate the momentum deduction from the solar wind, as can be seen from hybrid simulations with a larger subsolar extent because the solar wind is assumed to be undisturbed and PI is the only upstream source of cometary ions. The development of more kinetically self-consistent upstream mass-loading extension is beyond the scope of this paper and will be developed further in future studies.

In addition to momentum exchange, we also included the upstream pickup ion distribution. At simulation domain boundaries, the flux of pickup ions generated outside of the simulation domain was evaluated, as above, and injected at domain boundaries for each boundary cell. This enabled a spatial variability of the pickup ion fluxes at the boundary grid resolution, which is significant because the domain size is on the order of the H₂O⁺ gyroradius in the solar wind.

2.3 Virtual observations

Following the observations made by the RPC-ICA instrument onboard Rosetta, and to help in their interpretation with our model, we have developed a “virtual instrument” toolkit where pickup ion energy spectra are generated in the vicinity of the nucleus. This consists of three steps:

• 1.

  We define a virtual detector sphere with a radius of 100 km, centered at the origin;

• 2.

  Each pickup ion crossing the surface of the sphere (inwards) is recorded as a data point consisting of statistical weight, velocity, location, and injection and ionization point of the ion;

• 3.

  The pickup ion population data are then collected after the simulation has reached a quasi-stationary state until the end of the simulation (400–600s).

The observations generated by this technique are referred to as virtual observations in the following, to distinguish them from the actual observations made by the Rosetta RPC-ICA instrument. The number of particles collected per run are on the order of 300 000 particles at the 100 km sphere. Estimates for statistical fluctuations, estimated with compound Poisson distributions (Bohm & Zech 2014) to account for the variable particle weights, are shown in Fig. 2D for the total pickup ion fluxes.

The radius of the sphere was chosen to be approximately the distance to the nucleus of the Rosetta spacecraft during the observationspresented by Nilsson et al. (2018). The radius of the sphere was on the order of the maximum grid resolution; therefore nosignificant spatial variations are associated with the model solution at this scale. The modeled convective electric fieldin the vicinity of the nucleus (within the grid cells in direct contact with the origin) is on the order 0.15 mV m⁻¹. This yields an electric potential difference of approximately 30 V over the 200 km diameter of the virtual detector, which is a 30 eV difference in energy for an H₂O⁺ particle. This can be confounding for virtual observations of such low energies, but these particles are dominated by local ionization, whereas the particles pertinent to remote sensing will mostly populate energy regimes on the order of keV. Correspondingly, cometary ions with energies below 100 eV were excluded from the analysis of the virtual observations.

The virtual detector has some limitations. First, the scale of the detector is expected to produce some artifacts for simulated detections of low-energy ions, as described above. Second, even as the detector does accept particles from all directions, the simulation domain back wall boundary is relatively close to the nucleus. Consequently, possible sunward-propagating H₂O⁺ ions from outside of the tailward simulation domain would not be collected, the contribution of which we assume to be negligible based on the following observations. The tail plasma in general flows antisunward, leading to general antisunward motion of newly formed pickup ions. Additionally, backscattering of energetic pickup ions from the tail would most likely consist of cometary ions produced at far larger distances (comparable to the upstream gyroradius of the ions) than those we observe, that is, from a region of neutral density with a negligible relative contribution.

2.4 Overview of simulation features

Figure 1 presents a general overview of Run 2, displaying features of the general solar wind–comet interaction. Starting from the front wall and traveling along the path of a solar wind parcel, from right to left, four main features can be noted:

• 1.

  Initial and gradual mass-loading characterized by deflection and deceleration of the solar wind (Figs. 1A and C) with some draping of the magnetic field;

• 2.

  A sharp shock-like structure (Figs. 1A and C) with a sudden increase in magnetic field (Figs. 1B and C). This increase is associated with the drop in solar wind velocity and heating of solar wind, as we show in detail in Sect. 3;

• 3.

  A progressive pile-up of magnetic field approaching the nucleus (Fig. 1B);

• 4.

  Strong draping of the magnetic field lines (Fig. 1B).

We also note that the diamagnetic cavity boundary is not resolved in the solution because its size would be on the grid cell scale Δx_min = 250 km, as described by Goetz et al. (2016); this applies to other fine structures of the inner coma as well. Because the pickup ions relevant for remote-sensing purposes have a very high kinetic energy compared to convective electric fields at the scales of the close-in coma, we consider this a reasonable approximation of the near-nucleus environment for energetic pickup ions.

We notethat the results given by the chosen neutral profile are quite different from those that would be produced by a more realistic neutral model, such as that by Hansen et al. (2016). Such an asymmetric profile would place a larger portion of the outgassed neutrals sunward from the cometary nucleus, moving the plasma boundaries further upstream from the comet.

3 Results

In this section we present results of our simulations from Runs 1 (B_y= 0.98 nT), 2 (B_y = 1.96 nT), and 3 (B_y = 2.95 nT) and compare them with non-interacting ion pickup control cases, which we use as baseline. The results of Run 4 follow the trend shown by Runs 2 and 3, and are not included in detail for brevity.

3.1 Primer: non-interacting ion pickup

As the primary control cases, we used simulations with non-interacting pickup ions. To produce the idealized control cases, we only considered pickup ion production through PI and forced the electron fluid velocity U_e and magnetic field B to upstream values in the particle propagator. This resulted in the newly formed cometary photoions being accelerated with fixed IMF and convective electric field, forming an ideal pickup ion ring distribution in the velocity space. A control result was produced for each proper simulation, and these results are given as reference pickup ion spectra in Fig. 2D.

3.2 IMF magnitude variation

Figure 2 shows the main results for three different upstream magnetic field B_y values, with columns corresponding to Runs 1, 2, and 3. Panels A–D (rows) present (A) an overview of the plasma boundaries in the X–Z plane, (B) the velocity distribution of collected particles from the 100 km-radius virtual detector sphere around the nucleus, (C) the latitude–longitude maps of inbound cometary ions sorted by energy ranges, and (D) the energy spectrum of the cometary ion flux.

3.2.1 Plasma boundaries

Figure 2A shows a general overview of simulations for Runs 1, 2, and 3. We first calculated the magnetosonic Mach number M_ms as in Paper I:

as the sound and Alfvén velocities, respectively. Here ρ_i is the total mass density of the ions, T_e = 133 300 K and T_i are the electron and ion temperatures (accumulated kinetic temperature per ion population), and γ_e and γ_i are the polytropic indices for electrons and ions (per ion population), respectively. We assumed a polytropic index of 5 ∕ 3 for the calculation of the sound speed v_s for both electrons and each of the ion species. Different choices of γ_e and γ_i, for example, (γ_e, γ_i) = (1, 3), are not found to affect the interpretation qualitatively. We set the plasma velocity U_pl equal to the solar wind proton bulk velocity . A discussion on the choice of bulk velocity and γ values for the calculation of the Mach number is given in Appendix A in the context of solar wind ion kinetics.

Select contour lines of M_ms are displayed in the X–Z plane (Y = 0) in Fig. 2A. The red contour represents a value of M_ms = 2, which indicates the typical position of a mass-loaded bow shock-like structure (hereafter: bow shock) in the upstream solar wind (see Paper I and Galeev & Khabibrakhmanov 1990). The background color displays the magnetic field magnitude in the same plane. The position of the cometopause is overlaid in white, using the definition in Paper I of n_ci = n_sw, where n_ci is the cometary ion (H₂O⁺) number density and n_sw is the local solar wind proton density .

In these plots, the solar wind flows in from the right to the left, and the undisturbed solar wind convective electric field E_conv,sw = −U_e,sw ×B_sw is along the + Z (up) direction. Pickup ion acceleration is thus in the + Z direction, whereas, reciprocally, solar wind protons are deflected toward the − Z direction (see e.g. Kallio & Jarvinen 2012, for details on ion pickup).

We note that the choice of plasma velocity U_pl as is significant only for Run 1 at the position of the bow shock; another possible choice for U_pl, the electron fluid velocity U_e Mach number agrees with the onegiven by from Run 2onward. Furthermore, in Run 1 for the solar wind bulk motion the magnetosonic Mach number is larger than 2 almost everywhere in the domain, leading to the conclusion that Run 1 does not produce a shocked flow.

The contours of M_ms display some additional structure beyond determining a shock boundary, especially in Run 2. Within the cometopause, where the cometary ions dominate the plasma density, the magnetosonic Mach number for increases strongly and is not descriptive of bulk plasma properties: the magnetosonic signal velocity is relatively low (approximately 15 km s⁻¹) around the nucleus, while the H⁺ density decreases strongly (see Fig. 1C). The H⁺ bulk velocity is subsequently calculated from a trace population of relatively fast H⁺ particles reaching the inner coma, giving H⁺ bulk velocities upto 100 km s⁻¹ near the nucleus. In addition to this anomaly, Run 2 shows a complex M_ms = 2 boundary in the Z < 0 km hemisphere, which is partly due to finite gyroradius effects (especially at Z < 4000 km, see Appendix A) and partly due to M_ms being just above 2 in the intermediate region (see Appendix A).

For the simulations displaying shocks (Runs 2–4), the bow shock positions are well defined and the region delimited by M_ms = 2 ± 1 is quite narrow (approximately Δx across), and we take this as an indicator for a sharp boundary. For convenience, we define the stand-off distance for the bow shocks to be the distance between the comet and the M_ms = 2 surface along the X-axis; this is a reasonable compromise between the mass-loaded solar wind flow and the regions of interest with observed pickup ions, as we show in Sect. 4.

In Run 3 (nominal), the stand-off distance along the X-axis is 3000 km. For a smaller IMF magnitude (and increasing upstream M_ms), the stand-off distance shrinks to 2250 km (Run 2); and for increasing IMF magnitude (decreasing upstream M_ms) in Run 4, the bow shock stand-off distance is 3500 km. However, for Run 1, the stand-off distance cannot be defined.

For Runs 2–4, for which it is reasonable to define a bow shock stand-off distance, we show anempirical relationship between the upstream magnetosonic Mach number M_ms and bow shock stand-off distance R_bs in Fig. 3. These runs also display regions of true submagnetosonic flow, which is especially visible in Run 3.

In contrast to the bow shock, the stand-off distance of the cometopause decreases only slightly from 1000 (Run 3) to 800 km (Run 1); this is on the scale of Δx_min and therefore barely significant. Because the neutral profiles between the runs are identical, this is expected. The cometopause is somewhat less strongly extended in the + Z direction than in Runs 2 and 3.

The magnetic field magnitude around the comet presents a behavior similar to the evolution of the bow shock in Fig. 2A, with maximum magnetic field values increasing from 12.3 nT (Run 1) to 18.6 nT (Run 2) and 21.6 nT (Run 3) close to the nucleus in the magnetic pile-up region, while the maximum field strength relative to the IMF decreases. For Runs 2 and 3, an increase in magnetic field magnitude is observed just within the M_ms = 2 boundary location and toward the cometopause boundary and magnetic pile-up region. A strong asymmetry in the ±E_conv,sw direction is seen in all three runs.

This change in behaviour when the IMF magnitude increases corresponds to the expected evolution of a cometary magnetosphere, from a small weakly-interacting structure (Run 1), reminiscent of that described by Gunell et al. (2018), gradually toa more classical collisionless bow shock-like structure, as described in Paper I.

3.2.2 Simulated velocity distributions

The 3D velocity space distributions of particles collected by the virtual instrument are given in Fig. 2B. Each point represents a collected macroparticle; each macroparticle is colored according to its kinetic energy at the moment of collection, that is, as a function of radial distance from the origin at v = [0,0,0] km s⁻¹. The macroparticle weight is not included in the figure. Particles of a given energy are located at fixed radial distance from the origin. The coordinate system is that of the simulation: negative values of v_x represent motion in the antisunward direction, positive values of v_z represent motion in the direction of the (undisturbed) upstream convective electric field, and v_y represents the ions moving in the direction of the (transverse) upstream magnetic field. However, because of the transverse IMF condition, the distribution of v_y generally lies on the v_y = 0 km s⁻¹ plane, so that we chose to show the distribution from the direction of the negative v_y axis. The upstream solar wind is shown at the coordinates v = [-430,0,0] km s⁻¹, which is indicated by the orange triangle on the figures.

A shift can be observed from a distribution resembling an ideal pickup ion ring distribution in Run 1 (as in the non-interacting case discussed in Sect. 3.1) to distributions that are more affected by the environment (Runs 2 and 3), as the pickup ion ring around the solar wind velocity is gradually deformed. Particles with low-to-medium energy up to 1 keV are most affected by the environment, which is due to both their lower kinetic energy and their sensitivity to initial conditions at the ionization site.

Fig. 2 Main results of the transverse interplanetary magnetic field Runs 1–3 (columns). Row A: magnetic field magnitude, normalized to the upstream value for each run, and contours of magnetosonic Mach number (gray and red, see Sect. 3.2.1 for details) and cometopause (green). Row B: velocity space distributions of collected pickup ions; each point is colored according to particle kinetic energy, and the upstream solar wind velocity is indicated by the orange triangle (Sect. 3.2.2). Row C: ovals describe latitude–longitude maps of normalized inbound pickup ion fluxes at the virtual instrument for given energy channels, shown in a Hammer projection, and energy-latitude spectra of the pickup ion flux (Sect. 3.2.3). Row D: pickup ion flux energy spectra, given separately for total flux (black, with 3 σ errors), its constituent populations of PI (purple), CX ions (orange) and EI-produced ions (yellow), and non-interacting control ions (blue; Sect. 3.2.4).

Fig. 3 Bow shock stand-off distances Rbs along the X-axis in Runs 2–4 against upstream magnetosonic Mach number Mms. A least-squares fit is given for Rbs = a∕Mms + b, with associated upper and lower bounds at 1σ confidence. The joint 1σ confidence region for the parameters is shown shaded in light blue.

3.2.3 Latitude–longitude maps of cometary ion fluxes

Figure 2C describes the virtual observations in terms of latitude–longitude maps of inbound cometary ion incidence, between 0.1 and 65 keV in six different logarithmically spaced energy channels: 0.1–0.294, 0.294–0.866, 0.866–2.55, 2.55–7.5, 7.5–22.1, and 22.1–65 keV. The choice is arbitrary, and the channels can be chosen to suit the goals of the data analysis. Additionally, an energy-latitude spectrum is shown for each run, with energy channels marked for reference. In these figures, the particle latitude is defined as arcsin (+v_z∕v), in degrees, and the particle longitude is defined through the four-quadrant inverse tangent atan2(−v_y, −v_x), in degrees.

The maps represent detections of a virtual instrument with a field of view of 4π sr, from the point of view of an observer standing on the comet nucleus (because the collection sphere is contained within the first cell of the simulation), and looking out toward the Sun. Maps are given in the area-preserving Hammer projection and cover the full sphere of the virtual instrument, with a field of view of 180°×360° in latitude–longitude. The angular resolution and binning chosen here is 4.5°. A reference grid with a 30° graticule is overlaid on the latitude–longitude maps. Consequently, the ellipsoid center corresponds to the sunward incidence direction (with the Sun represented as a yellow star), that is, the center of the ellipsoid is usually populated with ions moving antisunward, toward the cometary nucleus. The green dotted circle corresponds to the IMF direction (entrance point out of the plane), 90° off the Sun-comet direction in the east. The direction of the solar wind convective electric field is orthogonal to the Sun-IMF direction, 90° toward the north and depicted by a red cross (exit point into the plane). The color code from blue to yellow describes increasingly high incident differential particle flux in arbitrary units. It is important to note that because of the large span in recorded flux intensities, we had to allow each separate energy channel to have its own normalized flux in log₁₀ scale, whichprevents a direct comparison of the intensities between the energy channels.

Run 1 (Fig. 2C1) shows that at low energies (below 2.5 keV) almost all cometary ions impacting the detector come from the direction of the convective electric field, almost antiparallel to it. The ions display a large angular dispersion, which is connectedto the increasing thermal spread of the ions with diminishing energies and increased draping in the innermost coma. For increasing energies above 4 keV, as the total flux of cometary ions decreases, so does the angular dispersion, with cometary ions entering increasingly from the south, along the convective electric field. For energies above 15 keV, two signals are simultaneously detected, one in the north, one in the south, which correspond to the pickup ion ring distribution expected far upstream the cometary nucleus.

Run 2 (Fig. 2C2) presents a somewhat different case for an IMF value twice as high as that of Run 1. As the interaction region spatially grows, very low-energy particles (<0.2 keV) gain access from the south, whereas low energies (0.2–1.8 keV) mostly come from the northern latitudes (antiparallel to the convective electric field) because of the cycloidal motion in the high-B, low-U_e inner coma. The angular spread is maximum at these energies. The turn-over from northern to the southern hemisphere occurs at 3 keV cometary ion energy, with the flux crossing over through the antisunward direction. For energies between 4 keV and the cutoff around 65 keV, Runs 1 and 2 behave similarly, with the main flux concentrated in the southern hemisphere and the progressive appearance of the pickup ion ring distribution.

Finally, the flux maps obtained from Run 3 (Fig. 2C3) are similar to those from Run 2, but the characteristic evolution of the distribution in energy takes place at higher energies, with the antisunward flux cross-over around 7 keV. The energy-latitude spectra in the bottom of Row C in particular illustrate the tendency of features translating to higher energies with increasing magnetic field.

3.2.4 Energy spectra of pickup ions

Figure 2D shows energy spectra of the omnidirectional particle flux for the virtual observations in units of m⁻² s⁻¹ eV⁻¹ sr⁻¹. The ion spectrometer RPC-ICA on board Rosetta uses 96 energy bins (Nilsson et al. 2007) between 25 and 40 keV. We chose to use for the generation of these spectra 50 logarithmically spaced energy bins from 10 to 100 keV; results below 100 eV are not shown. The resulting relative bin width is approximately 15%. The spectra are shown separately for the non-interacting control case (blue; only the photoion population is included), the total flux (thick black line, ±3σ limits are shown as thin black lines, given by compound Poisson distribution estimation along Bohm & Zech 2014), and fluxes given by particles originating from different source processes (PI, CX, and EI).

For the upstream pickup ion population, only photoions were included in the virtual observations, resulting in underestimated total particle fluxes at energies above 10 keV as CX and EI ion fluxes go to zero at energies corresponding to ions created at the domain boundary. However, the total flux spectra at high energy, when only PI is included, are falling off rather smoothly up to the maximum pickup ion energy (>65 keV, given by a pickup ion moving at twice the solar wind speed), and in accordance with the control runs. In this context, our extended boundary condition model underestimates the effective upstream mass-loading with little loss of physical accuracy, keeping in mind that one of the main motivations for investigating pickup ion energy spectra was to interpret observed features seen by the ICA instrument at 1 keV.

With increased magnetic field magnitude from Run 1 to Run 2, the energy spectra increase in complexity: the spectra from Run 1 more closely resembles the non-interacting control case, while Runs 2 and 3 display several spectral breaks at various energies. Wenote the nearly constant slope of the non-interacting control case with respect to the spectral features, and that the absolute value of the control case is somewhat arbitrary because it only includes photoions. The spectral features are analyzed in detail in thefollowing section.

Fig. 4 Energy channel selection for Runs 1, 2, and 3 (left to right). Panels A: particle flux energy spectra vs. energy; the chosen energy channels of interest are labeled 1–6. Panels B: spectral index α(E) = dlog10f(E)∕dlog10E (blue line) and its derivative in β(E) = dα(E)∕dlog10E. The locations of the six most prominent extrema of the second derivative β are indicated in panels B in electronvolts, and the prominence of the peak is color-coded: blue at the lowest and dark red at the highest values of prominence.

4 Interpretation

To understand the features of the energy spectra, we divided the virtual observations (i) in terms of spectral features to map them to ionization points and plasma regions, and (ii) in terms of the regions of origin of the observed particles to directly describe the contributions of different regions of the spectra. In general, the tendency described in the previoussection from Run 1 to Run 3 is one of increasing complexity of pickup ion spectra by the plasma environment.

4.1 Energy separation analysis

Figure 4A shows the smoothed virtual particle flux energy spectra for cometary pickup ions for Runs 1, 2, and 3, overlaid with energy channels of interest found for each run. These energy channels are defined by the most prominent extrema of the second derivative of the signal. To capture the spectral index and its change in energy, the spectra were first smoothed out using a local linear regression filter (lowess) with a span equal to 10% of the total number of data points in the spectrum. The local linear regression method was chosen for its direct relation to the local slope of the signal, although the choice of the smoothing filter is not important and locations of the extrema vary by approximately 20% depending on the choice of smoothing filter (or no filter at all).

Using the smoothed data, we first calculated the spectral index α(E) = dlog₁₀(f(E))∕dlog₁₀(E). The second derivative of the smoothed signal β(E) = dα(E)∕dlog₁₀(E) was then calculated, and the most significant peaks of β(E) were extracted. As a measure of the significance of the peaks we used peak prominence¹. The six most prominent extrema above 100 eV and below 10 keV were retained to define the boundaries of the energy channels, as shown in Fig. 4B. The lowest-prominence spectral breaks (blue peaks in Fig. 4B) are mostly of little interest and are included only to display the limits of the method. Techniques that use the second derivative to determine the position of inflection points or other features are widely used in analyzing Langmuir probe current-voltage characteristics (e.g. Yang et al. 2016; Johansson et al. 2017).

The energy channels were chosen to characterize the energy of the particles and relate them to physical features at the cometary plasma environment. Seven energy channels are defined by the spectral breaks, and they are labeled 1–7 and are individually color-coded in Fig. 4. Generally, the variation in the locations of the main extrema is consistently around 20% when it is controlled for each of the following: the choice of the smoothing filter (as above), energy bin width (binning with 100 bins), and choice of virtual detector radius (150 km versus the presented 100 km radius).

As described above, Run 1 displays a somewhat featureless energy spectrum that most closely resembles that of the non-interacting pickup ion case. The spectra of Runs 2 and 3 look similar, and both contain a sharp transition on the order of the solar wind proton energy (1 keV), followed by a wide plateau. The start of this feature corresponds to the largest positive maximum of the second derivative encountered when moving from high to lower energies (dark red points in Fig. 4B). For Run 1, no such distinct feature is detected. For Run 2, a drastic change in slope occurs at 1599 eV, while for Run 3, the plateau starts at 2304 eV. The start of the plateau is somewhat more unclear in Run 3, with a wide maximum around 2 to 3 keV. The spectral plateaus for Run 2 and Run 3 (energy channels 5 and 6, respectively) return to the general decreasing trend at energies of 4013 and 8393 eV, respectively. For Run 4 (not shown), the sharp plateau-starting transition is at 3342 eV, and the plateau ends at around 10 keV.

To investigate how this pickup ion energy region maps out to the spatial distribution of cometary ions, we studied the corresponding ionization points of the cometary pickup ions. The ionization points were stored for each particle during the simulation, and these points were collected as the particles reached the virtual detector. Ionization points are shown in Fig. 5 in the X–Z plane. Each colored region corresponds to the seven energy channels of its corresponding non-weighted macroparticles, as defined in Fig. 4. At increasing distances from the nucleus, the effect of splitting macroparticles generated in the coarse grid increases the spread of the ionization-point curve. Controlling the results with no macroparticle splitting for pickup ions in Run 2 did not produce significant differences.

The emergence of a plateau, or a close-zero spectral index in Runs 2 and 3, is marked by a transition from light orange (green) to pink (dark green) for Run 2 (Run 3), as seen in Fig. 4. These breaks in the spectral index map out to the bow shock-like structure marked by the surface M_ms = 2. When we consider the comet-Sun line distance from the nucleus to the surface, these transitions map to bow shock stand-off distances R_bs = 2250 and 3000 km for Runs 2 and 3, respectively. The precision of this mapping is constrained because the distribution of the ionization points is almost tangential to the shock surface, which may explain the wide maximum of β(E) in Run 3. The emergence of the spectral plateaus is discussed further in Sect. 4.2.

In Fig. 5 the cometopause is also shown for reference. The cometopause is found to map to particles in the lowest energies (100 eV and below). Because the energy of these particles is lower than the electric potential across the virtual detector, they represent the actual local fluxes less accurately. In Run 1, the cometopause surface does constrain the lower limit of displayed particles at 100 eV, however, and it is possible that the inner coma environment would also create an observable signature in the spectrum. The scope of this analysis and the resolution of the simulation precludes firm conclusions on this region, however. Additionally, we note that these observations are localized at the vicinity of the nucleus, so that it is possible that repeating the virtual observations at a different location could yield some signal from the asymmetric magnetic field enhancement shown in Fig. 2A1.

Figure 6 shows the view geometry maps of particle fluxes at the energy channels inferred from the breaks in the spectral index. The bottom row of the figure shows the corresponding energy channels in the energy-latitude histogram: the passover of the antisunward direction corresponds closely to the spectral break associated with the shock surface.

From the cometary pickup ion spectral features found for Runs 2 and 3, we conclude that the sharp knee in the energy distribution of the pickup ion fluxes above the solar wind proton bulk energy correlates well with the appearance of a bow shock-like structure in the upstream solar wind. Two broad energy regions in the pickup ion energy distribution start to emerge: one low-energy region inside the shock structure, characterized by high fluxes and a usually complex behavior, and one high-energy region outside of it, characterized by an ion flux that first stabilizes before it decreases, following a power law. In order to unambiguously link the bow shock with this spectral feature, we examine the cometary ion trajectories in more detail in the following.

Fig. 5 Ionization points colored by energy channels for Runs 1, 2, and 3 (left to right). The color code is the same as in Fig. 4. Surfaces for Mms = 2 (red) and cometopause (purple) are included in the figures for context. The visible portion of the surfaces is constrained to Y = ±500 km. In Run 1 a secondary “cometopause” occurs around Z ≈ [2000, 3000] km, showing a thin plasma region, entirely within the Y = ±500 km extent, which is dominated by cometary pickup ions embedded in the solar wind.

Fig. 6 Latitude–longitude maps and energy-latitude histograms of the virtual observations. The format is the same as in panel C of Fig. 2; the energy channels are chosen from the flux energy spectrum of Fig. 4. The limits of the energy channels are shown in the energy-latitude histograms.

4.2 Trajectories of pickup ions

To analyze the mechanisms responsible for the spectral features discussed above, we performed test-particle simulations to map the trajectories of the particles observed within given energy channels. We took the averaged fields U_e and B over a period of time (350–600 s) and traced the motion of H₂O⁺ particles using the Lorentz force. Integration was handled with a Runge–Kutta 2–3 scheme, and the fields were interpolated with a first-order scheme. Theaveraged fields discard the dynamics of the shocked plasma environment, which is highly variable in the dynamic simulation, leading to increased heating and mixing of particles generated within and near the shocked surface. This widens the distribution of injection points shown in Fig. 5 compared to those obtained through test particle simulations shown in this section.

In Fig. 7 we show test particles trajectories, traced back in time from the nucleus; taken from a velocity distribution with energies from 441 to 5 keV (logarithmically spaced with 100 samples; from the lower limit of energy channel 3 shown in Fig. 4 up to energies in channel 6) with ions at each energy launched on the X–Z plane radially outwards, with 1° angular resolution. Only particles that reach a kinetic energy lower than 4 eV are shown; these threshold points are marked in green. Because these particles have been approximately stationary, they are possible pickup ions created at these points. The test particle trajectories are colored with their energy at the nucleus, approximating the virtual observation energy of corresponding particles. The color scheme is such that particles with an initial energy lower than the 1599 eV threshold are blue, and those at higher energies are red-yellow. The sphere at the origin denotes the extent of the virtual detector sphere. The traces continue farther upstream, showing that pickup ions from upstream on this particular trajectory would be detectable at the nucleus as well.

Figure 7 also shows the M_ms = 2 bow shock in magenta, showing the correspondence of the knee energy ion-injection points with the location of the bow shock. Notably, thebehavior of the pickup ions is highly kinetic and almost unmagnetized at these length scales, and the pickup ions move almost parallel to the local convective electric field.

The plateau in pickup ion flux energy spectra is formed through a combination of different complementary factors. Figure 8 shows the estimated energies of cometary pickup ions at their closest approach to the nucleus, that is, approximately corresponding to the collection energy at the virtual detector. The region outlined by the black line segments is such that pickup ion test particles ionized from this region hit the 100 km virtual detector. With the shown energy at closest approach, this shows an effective sampling region for collected pickup ions of a given energy. When we consider the sampling region for a given pickup ion flux at decreasing energies from upstream toward the bow shock, the cross-section S(E) of this sampling volume dV (E) is decreasing.

Furthermore, the relation of the sampling region and the cometary neutral density and pickup ion production is shown in Fig. 9. Up to the shock region, the sampling region on average moves radially outward with increasing particle energy, that is, along the direction of the negative gradient of neutral density n_i. After the bow shock, the pickup ion sampling region moves non-radially with increasing energies, unaligned with the neutral density gradient. Furthermore, the enhancement of particle production P_i both in mass-loaded pre-shock deceleration region and in the shock proper allow the sampling region to sample higher source densities than can be explained by the neutral density alone.

In conclusion, the spectral plateau is formed as a result of the two following factors:

– Geometrical effects: the sampling volume dV (E) =S(E) ⋅ dr(E) varies throughout the acceptance region and has a minimum at the bow shock.

– Ion source density effects: the pickup ion production rate P_i = ν_in_n, with ν_i denoting the ionization frequency due to PI, CX, and EI processes and n_n the underlying neutral density, varies along the acceptance region.

When we take these two factors into account, the virtual observations of a spectral plateau can be qualitatively reproduced starting from the particle ionization locations (and their spread) of the virtual observations (Fig. 5) and the pickup ion production rates given by the simulation.

Fig. 7 Back-propagated ion trajectories for Run 2. The color of the trajectory shows the particle energy in the vicinity of the nucleus (approximately virtual particle collection energy), with the blue-red turnover at the knee energy. The green dots along the trajectories show points on the trajectories where the particle kinetic energy on the trajectory is lower than 4 eV, approximating possible points of pickup ion production. On the Y = 0 plane, the magnetic field magnitude in grayscale and the convective electric field arrows (purple) are shown for context. The Mms = 2 contour is shown in red.

Fig. 8 Close-up of injected test particle behavior in Run 2. The continuous red line shows the extent of the Mms =2 surface. The background color shows a mapping from a test particle created at that location to the energy of the particle at its closest approach to the nucleus. This approximates the energy of the particle as if it were collected by the virtual detector, given that the particle hit the detector. The color scale is centered on the knee energy. The black line segments denote particles with a minimum distance to the nucleus in the range 100–150 km. Traces within the black outlines have a minimum distance to the nucleus smaller than 100 km, and these particles would be collected by the virtual detector (i.e. particles are ionized within the sampling region), and traces outside of the black outlines have a minimum distance to the origin larger than 200 km. The arrows show the direction and magnitude of the local convective electric field; white corresponds to the undisturbed solar wind electric field.

Fig. 9 Close-up of injected test particle behavior in Run 2. The red line shows the extent of the Mms =2 surface, and the black line segments outline the sampling region as in Fig. 8. The purple lines show contour lines of the cometary neutral density, and the background color shows the pickup ion production rate in units of ions m−3 s−1. The value of the production rate does not coincide exactly with the neutral density contour lines, instead, the isocontours of production rates follow the Mms = 2 surface in this region.

4.3 Relation of bow shock distance and primary spectral break

We have presented the correlation between and a causal relationship of the spectral break and plateau at a cometary bow shock above. The IMF magnitude has an effect on the bow shock stand-off distance. This is seen in the simulations, which in addition, describe a relation of the shock and its environment to observable plasma quantities in the vicinity of a cometary nucleus.

Based on the simulation results in Runs 2 through 4, we present a tentative relation between the energy of the knee in the pickup ion spectra to the bow shock stand-off distance in Fig. 10. The relation is not a realistic prediction of the bow shock stand-off with respect to actual observations, noting that a realistic, asymmetric outgassing profile (such as shown in Hansen et al. 2016) would push the environment and the bow shock further sunward. We instead present the relation to show that a cometary bow shock and its responses to the solar wind could be inferred from RPC-ICA observations, for instance, as the knee would be in accordance with our simulations dependent on the physical processes in the vicinity of the bow shock.

It is important to note that the energization of the pickup ions is due to the potential given by the convective electric field E_conv = −U_e ×B. Because varying the upstream magnetic field directly affects the electric field, the energies of pickup ions from a given distance will increaseas well. This accounts for some of the shown energy shift of features. Separating the energy shifts that are due to convective potential differences from effects that are due to bow shock stand-off distance is left for further studies.

To produce the relation in Fig. 10, we took the calculated knee energies (Fig. 4) with their ± 20% uncertainties, and the bow shock stand-off distances along the comet-Sun line as defined by the M_ms = 2 surface. We assumed a grid resolution ± Δx_min as the stand-off distance uncertainty. A line was fit to the data by χ² minimization with respect to both errors, and upper and lower error estimates with 1σ confidence. The shaded region in Fig. 10 shows the joint 1σ confidence region.

Fig. 10 Low-energy limit of the spectral plateau (knee) vs. bow shock stand-off distance for Runs 2–4, with a linear least-squares fit. The ±20% error bars are shown for knee energy and ± Δxmin error bars for stand-off distance. Linear fit parameters are given with associated upper and lower bounds at 1σ confidence. The 1σ joint confidence region is shaded with light blue.

4.4 Relation to observations

A qualitative comparison of the ICA observations (Nilsson et al. 2018) and the simulated spectra is encouraging: the simulated energy spectrum has a slope break at similar energies as in the observations, a following plateau, and a subsequent sloped drop just before the maximum pickup ion energies. A weakly coupled simulation, such as Run 1, does not produce similar features.

Nilsson et al. (2018) noted that an inverse square fall-off with energy was expected because particles with higher energy would be produced farther out in the cometary magnetosphere, which is captured by our model. They noted that a break in the energy spectra at about 1 keV, with an inferred source region density at higher energies, indicated that something happened beyond the distances corresponding to about 1 keV energy (for the case they showed). They suggested that it could be a remote signature of a larger-scale bow shock. We have shown that the 1∕E² fall-off with energy is indeed reproduced in our hybrid simulations, as is a break in the energy spectra at an energy of about 1 keV in Fig. 4. In our model the spectral break is indeed related to the bow shock and to the effect of the electric field on the comet ion trajectories there. This supports the hypothesis that cometary ion spectra as collected by the RPC-ICA instrument can be used to remotely probe the shock region in the escort phase of the Rosetta mission.

5 Conclusions

The cometary ion energy spectrum is found to be sensitive to the surrounding plasma environment. A smooth spectrum corresponds to a weakly coupled cometary plasma (Run 1 with low |B_IMF|; this solution tends toward a spectrum produced by pickup ions in an undisturbed flow), and a correspondingly small cometary environment. In contrast, a strongly coupled system with large-scale features and a bow shock produce several spectral features. This development can be seen to mirror the transitions from an infant bow shock, as discussed by Gunell et al. (2018), toward a proper bow shock.

The (sharp) transition region, that is, the shock-like structure, is shown as the set of upstream points, where the magnetosonic Mach number drops to 2, as discussedfor subcritical mass-loaded shocks in Galeev & Khabibrakhmanov (1990). An application of this Mach number to the system is shown in Appendix A. The low-energy limit of the spectral plateau, described inSect. 4.1, corresponds to ions generated in the vicinity of this surface. As shown in Fig. 10, the relationship between the spectral break energy and the stand-off distance has a linear trend (for the given parameter space). The mechanistic interpretation given in Sect. 4.2 shows that the feature is caused by the characteristics of the mass-loaded shock. Therefore, similar effects can be expected to occur in a more general setting with a similar signature in the pickup ion spectrum observed in the vicinity of the nucleus.

Different cometary ion populations show somewhat different characteristics because the plasma environment affects the reaction rates of the production processes as the plasma density, temperature, and flux participating in charge exchange and electron impact reactions are modified. Although not distinguishable by an in situ measurement, separating the populations through modeling can be used to analyze the remote observations in more detail. Our findings are summarized below.

• 1.

  We showed that the shape of cometary ion spectra can be used to infer large-scale structures, such as a bow shock, and gained some information about the distance to the shock.

• 2.

  We showed that the comet environment changes from essentially no large-scale bow shock to a well-developed bow shock when we increased the solar wind magnetic field strength and therefore the coupling between solar wind and cometary ions.

• 3.

  The energy of the spectral knee is dependent on the upstream magnetic field and the bow shock distance.

• 4.

  The knee observed by Nilsson et al. (2018) is consistent with a bow shock.

Acknowledgements

M.A. acknowledges financial support from the Aino and Kaarlo Tiisala Foundation and the RPC meetings for fruitful discussions. C.S.W. was funded at University of Oslo by the Norwegian Research Council “Rosetta” grant no. 240000. The authors acknowledge ISSI (Bern) and the ISSI International Team “Plasma Environment of comet 67P after Rosetta” for fruitful discussions and collaborations. Aalto Science-IT is recognized for providing computational capacities. The Lawrence Livermore National Laboratories visualization software VisIt was used for the analysis of the simulations and the production of some publication figures. The authors would also like to thank R. Munroe (xkcd) for inspiring us in the use of sophisticated curve-smoothing techniques.

Appendix A Assumptions of Mach numbers compared to ion kinetics

The choice of the Mach number to define shock surfaces in a kinetic regime is non-trivial. Galeev & Khabibrakhmanov (1990), for example, used a single-fluid MHD approximation, with a single plasma-flow velocity. Figure A.1 shows the behavior of Mach numbers for different choices of plasma velocity U_pl and how they relate to solar wind H⁺ trajectories in time-averaged fields. In the calculation of the magnetosonic velocity, γ_e = γ_i = 5∕3 was assumed.

For Run 1, the Mach number of the center-of-mass velocity shows significant influence of the high-mass stationary pickup ions close tothe comet, as the center-of-mass Mach surface extends farther upstream, especially in the heavily mass-loaded + Z hemisphere. For Runs 2 and 3, the center-of-mass M_ms = 2 surface converges toward H⁺ and U_e Mach surfaces.

The velocity of choice in Paper I was the electron fluid velocity U_e, which corresponds to the velocity calculated from for shocks far from the cometopause, as shown in Fig. A.1, confirming that the stand-off distance result of Paper I is robust against the choice of plasma velocity. For Run 1, the electron fluid Mach surface has a discontinuous ridge throughout the minimum region of proton Mach number.

The choice of γ, as noted in Sect. 3.2.1, does not significantly affect the location of the proton M_ms = 2 surfaces in strongly coupled cases (Run 2 and further), but assuming isothermal electrons and 1D motion for ions (γ_e, γ_i) = (1, 3) does allow for a thin approximately 1 Δx region, in which M_ms ≤ 2 for protons at the proton trajectory convergence zone (“caustic”). Otherwise, the robustness of the definition against variations in γ helps to validate the use of ion temperatures in calculating the acoustic mode speed because the ion populations are quite non-Maxwellian.

Figure A.1 also shows solar wind proton trajectories. For Run 1, we see a similar result to that of Behar et al. (2018) and Saillenfest et al. (2018), with the convergence of proton trajectories forming an asymmetric caustic and a solar wind cavity. As the coupling between cometary ions and the solar wind increases through increasing IMF magnitude, we see in detail how the M_ms = 2 surface is formed, and how finite ion gyroradius effects are visible in the sheath as mild ridges in the Mach number values. Notably, from Run 2 onward, the test particles also display flow around the inner coma in the Y direction, which is not observed for Run 1.

Fig. A.1 For Runs 1–3 (left to right), solar wind H+ test particle traces (injected near the front wall with vi =Ue) and Mms = 2 colored with the Y coordinate, i.e. the off-plane distance. The magnetosonic Mach numbers calculated separately for electron bulk fluid velocity Ue (magenta) and center-of-mass velocity vtot (orange; holes in this contour are artifacts from visualization at grid refinenent interfaces) are shown at the Mms = 2 contour, andfor the solar wind proton bulk velocity as background color on the XZ plane. The proton Mach number color scale is chosen so that Mms = 2 is noted by white and true subsonic flow Mms ≤ 1 by yellow to pink for a decreasing Mach number. For Run 1, is larger than 2 almost everywhere. The apparent lack of an SW cavity for Runs 2 and 3 is due to particles moving into and out of the plane of view. This is shown by the color of the traces.

References

All Tables

All Figures

Fig. 1 General overview of Run 2. Panel A: magnetosonic Mach number = 2 and cometopause surfaces (orange and red, respectively) and the electron fluid-velocity streamlines, overlaid with the simulation grid displaying the grid refinements. On the X =−1000 km plane, the magnetic field magnitude is displayed in blue. Panel B: magnetic field lines connected to the comet-Sun line, viewed along the direction of the convective electric field. Features of the environment listed in the main text are annotated in the panels. Panel C: values of the magnetic field magnitude B and solar wind H+ bulk values (number density , bulk velocity and thermal velocity vth) relative to upstream values, along the X-axis from X =−1000 to X = 9000 km. Simulation cometopause (red), bow shock standoff (blue) distances, and the extent of the diamagnetic cavity boundary observations(green) of Rosetta up to July 2015 (Goetz et al. 2016) are shown by the shaded areas for reference. Four regions of the environment are numbered in the panels. 1: gradual mass-loading and some draping (panels A–C), 2: a sharp shock-like structure (A–C), 3: magnetic field pile-up approaching the nucleus (B and C), and 4: draping of the magnetic field lines around the inner coma (B).

In the text

Fig. 2 Main results of the transverse interplanetary magnetic field Runs 1–3 (columns). Row A: magnetic field magnitude, normalized to the upstream value for each run, and contours of magnetosonic Mach number (gray and red, see Sect. 3.2.1 for details) and cometopause (green). Row B: velocity space distributions of collected pickup ions; each point is colored according to particle kinetic energy, and the upstream solar wind velocity is indicated by the orange triangle (Sect. 3.2.2). Row C: ovals describe latitude–longitude maps of normalized inbound pickup ion fluxes at the virtual instrument for given energy channels, shown in a Hammer projection, and energy-latitude spectra of the pickup ion flux (Sect. 3.2.3). Row D: pickup ion flux energy spectra, given separately for total flux (black, with 3 σ errors), its constituent populations of PI (purple), CX ions (orange) and EI-produced ions (yellow), and non-interacting control ions (blue; Sect. 3.2.4).

In the text

	Fig. 3 Bow shock stand-off distances Rbs along the X-axis in Runs 2–4 against upstream magnetosonic Mach number Mms. A least-squares fit is given for Rbs = a∕Mms + b, with associated upper and lower bounds at 1σ confidence. The joint 1σ confidence region for the parameters is shown shaded in light blue.
In the text

Fig. 4 Energy channel selection for Runs 1, 2, and 3 (left to right). Panels A: particle flux energy spectra vs. energy; the chosen energy channels of interest are labeled 1–6. Panels B: spectral index α(E) = dlog10f(E)∕dlog10E (blue line) and its derivative in β(E) = dα(E)∕dlog10E. The locations of the six most prominent extrema of the second derivative β are indicated in panels B in electronvolts, and the prominence of the peak is color-coded: blue at the lowest and dark red at the highest values of prominence.

In the text

Fig. 5 Ionization points colored by energy channels for Runs 1, 2, and 3 (left to right). The color code is the same as in Fig. 4. Surfaces for Mms = 2 (red) and cometopause (purple) are included in the figures for context. The visible portion of the surfaces is constrained to Y = ±500 km. In Run 1 a secondary “cometopause” occurs around Z ≈ [2000, 3000] km, showing a thin plasma region, entirely within the Y = ±500 km extent, which is dominated by cometary pickup ions embedded in the solar wind.

In the text

	Fig. 6 Latitude–longitude maps and energy-latitude histograms of the virtual observations. The format is the same as in panel C of Fig. 2; the energy channels are chosen from the flux energy spectrum of Fig. 4. The limits of the energy channels are shown in the energy-latitude histograms.
In the text

Fig. 7 Back-propagated ion trajectories for Run 2. The color of the trajectory shows the particle energy in the vicinity of the nucleus (approximately virtual particle collection energy), with the blue-red turnover at the knee energy. The green dots along the trajectories show points on the trajectories where the particle kinetic energy on the trajectory is lower than 4 eV, approximating possible points of pickup ion production. On the Y = 0 plane, the magnetic field magnitude in grayscale and the convective electric field arrows (purple) are shown for context. The Mms = 2 contour is shown in red.

In the text

Fig. 8 Close-up of injected test particle behavior in Run 2. The continuous red line shows the extent of the Mms =2 surface. The background color shows a mapping from a test particle created at that location to the energy of the particle at its closest approach to the nucleus. This approximates the energy of the particle as if it were collected by the virtual detector, given that the particle hit the detector. The color scale is centered on the knee energy. The black line segments denote particles with a minimum distance to the nucleus in the range 100–150 km. Traces within the black outlines have a minimum distance to the nucleus smaller than 100 km, and these particles would be collected by the virtual detector (i.e. particles are ionized within the sampling region), and traces outside of the black outlines have a minimum distance to the origin larger than 200 km. The arrows show the direction and magnitude of the local convective electric field; white corresponds to the undisturbed solar wind electric field.

In the text

Fig. 9 Close-up of injected test particle behavior in Run 2. The red line shows the extent of the Mms =2 surface, and the black line segments outline the sampling region as in Fig. 8. The purple lines show contour lines of the cometary neutral density, and the background color shows the pickup ion production rate in units of ions m−3 s−1. The value of the production rate does not coincide exactly with the neutral density contour lines, instead, the isocontours of production rates follow the Mms = 2 surface in this region.

In the text

	Fig. 10 Low-energy limit of the spectral plateau (knee) vs. bow shock stand-off distance for Runs 2–4, with a linear least-squares fit. The ±20% error bars are shown for knee energy and ± Δxmin error bars for stand-off distance. Linear fit parameters are given with associated upper and lower bounds at 1σ confidence. The 1σ joint confidence region is shaded with light blue.
In the text

Fig. A.1 For Runs 1–3 (left to right), solar wind H+ test particle traces (injected near the front wall with vi =Ue) and Mms = 2 colored with the Y coordinate, i.e. the off-plane distance. The magnetosonic Mach numbers calculated separately for electron bulk fluid velocity Ue (magenta) and center-of-mass velocity vtot (orange; holes in this contour are artifacts from visualization at grid refinenent interfaces) are shown at the Mms = 2 contour, andfor the solar wind proton bulk velocity as background color on the XZ plane. The proton Mach number color scale is chosen so that Mms = 2 is noted by white and true subsonic flow Mms ≤ 1 by yellow to pink for a decreasing Mach number. For Run 1, is larger than 2 almost everywhere. The apparent lack of an SW cavity for Runs 2 and 3 is due to particles moving into and out of the plane of view. This is shown by the color of the traces.

In the text