Remote sensing of cometary bow shocks: modelled asymmetric outgassing and pickup ion observations

Despite the long escort by the ESA Rosetta mission, direct observations of a fully developed bow shock around 67P/Churyumov– Gerasimenko have not been reported. Expanding on our previous work on indirect observations of a shock, we model the large-scale features in cometary pickup ions, and compare the results with the ESA Rosetta Plasma Consortium Ion Composition Analyser ion spectrometer measurements over the pre-perihelion portion of the escort phase. Using our hybrid plasma simulation, an empirical, asymmetric outgassing model for 67P, and varied interplanetary magnetic ﬁeld (IMF) clock angles, we model the evolution of the large-scale plasma environment. We ﬁnd that the subsolar bow shock standoff distance is enhanced by asymmetric outgassing with a factor of 2 to 3, reaching up to 18 000 km approaching perihelion. We ﬁnd that distinct spectral features in simulated pickup ion distributions are present for simulations with shock-like structures, with the details of the spectral features depending on shock standoff distance, heliocentric distance, and IMF conﬁguration. Asymmetric outgassing along with IMF clock angle is found to have a strong effect on the location of the spectral features, while the IMF clock angle causes no signiﬁcant effect on the bow shock standoff distance. These dependences further complicate the interpretation of the ion observations made by Rosetta . Our data-model comparison shows that the large-scale cometary plasma environment can be probed by remote sensing the pickup ions, at least when the comet’s activity is comparable to that of 67P, and the solar wind parameters are known.

proposed that the energy spectrum of cometary ions picked up by the solar wind may give an estimate of the ion density in the region of production of these ions. Using the RPC-ICA measurements of pickup ion fluxes, they found that this density followed a 1/r 2 dependence on cometocentric distance r up to an energy of about 1 keV, effectively linking the energy of the pickup ions with the cometocentric distance from which they originate. This technique makes it possible to use in situ observations of the dayside cometary plasma environment to remotely sense the largescale cometary coma with point observations from deep inside the system itself. We anticipate the techniques presented here and in Alho et al. (2019) being useful for a retrospective on the Rosetta escort phase and future missions such as the ESA Comet Interceptor (Snodgrass & Jones 2019).
Here, we present a continuation of our work initiated by Simon Wedlund et al. (2017, henceforth referred to as Paper I) and continued in Alho et al. (2019, Paper II), where we explored this idea with the help of 3D numerical hybrid simulations, and described a method of detecting large-scale structures in the plasma environment of comet 67P, especially focusing on the bow shock, and using the energy spectra of pickup ion fluxes. These signals should potentially be observable by RPC-ICA, even when the Rosetta spacecraft is deep within the coma and close to the nucleus -as it almost invariably was for the duration of the escort phase. Notably, even the extended dayside excursion at subsolar cometocentric distances of about 1500 km in 2015 September-October (just after perihelion) failed to detect a bow shock structure (Edberg et al. 2016). In Paper I, numerical simulations of the perihelion conditions including chargeexchange and ionization processes, even without taking into account the asymmetry of the nucleus outgassing which is expected to push the boundaries further upstream, showed that the bow shock standoff distance should be further than 6000 km from the nucleus for the maximum outgassing rate of the comet 67P. The only, nascent, bow shock-like structure during the mission was found between 2.2-2.4 au by Gunell et al. (2018) using measurements from the RPC magnetometer, electron, and ion spectrometers, and a hybrid simulation to interpret the observations.
In this work, we expand upon our previous study in Paper II with additional hybrid simulations, by including a realistic, empirical outgassing model for the coma (Hansen et al. 2016) over a wide range of heliocentric distances. Specifically, we model a number of cases from 2.0 au to near-perihelion for the inbound phase of the cometary orbit.
From the point of view of our model, the plasma environment of a comet such as 67P is governed by: (i) The impinging solar wind: solar wind ion densities and interplanetary magnetic field (IMF) direction and magnitude (Paper I); (ii) The neutral coma, its extent and evolution, and the corresponding ionization processes [charge exchange (CX), photoionization (PI), and electron ionization (EI)].
Both aspects depend on the solar activity, the location of 67P in the Heliosphere, and in case of the coma, on the orientation of the comet with respect to solar illumination, the shape of the nucleus and its properties (Fougere et al. 2016;Marshall et al. 2019). Any model of the pickup ion contribution to the cometary environment should take these two aspects into account. However, for simplicity, we neglect any contribution of cometary dust, whether neutral or charged. Including dusty plasma effects would likely generate dust-related wave modes and effects in the cometary tail (Mendis & Horányi 2013). Furthermore, while some work on the dust environment of 67P has been done (Gunell et al. 2015;Hornung et al. 2020), there are no models that we are aware of that would describe the sizes, charging, and distribution of dust particles in the cometary environment. While the dust effects are outside the scope of this work, including them in future models is an interesting challenge, as it is well known that dust charging happens and can increase the mass-loading of the comet's environment.
To improve on the description of the neutral coma from that in Papers I and II, we adopt an empirical model of Hansen et al. (2016) for our cometary ion source population, which replaces the spherically symmetric Haser model. The Haser model, previously used by e.g. Rubin et al. (2015) and Koenders et al. (2013), produces bow shock standoff distances of few 1000 km around perihelion for 67P. Further work has extended the estimates for bow shock standoff distances using the Haser model (Paper I and Paper II), while Huang et al. (2016) used their four-fluid model to explore effects of asymmetric neutral outgassing, and found expanded sunward plasma environments combined with increased bow shock standoff distances. In this study, using our hybrid model, we report similar findings.
The paper is organized as follows: In Section 2, we briefly present the hybrid model, which we use to calculate the pickup ion energy spectra and the dayside bow shock structure (see Paper I and Paper II for details). We describe the new addition of a semiempirical analytical model of the asymmetric neutral coma based on the work of Hansen et al. (2016), list the simulated events in detail and introduce the method used to draw conclusions from the pickup ion spectra. Section 3 presents our results applied to 67P, covering eight cases with varying heliocentric distances from 1.86 to 1.25 au and varying IMF configurations. We also present observations from the ROSETTA/ICA instrument and discuss their relevance to our modelling results. Finally, section 4 discusses the shape of the pickup ion energy spectra and the subsolar bow shock standoff distance, in light of measurements from the Rosetta RPC-ICA instrument.

M O D E L A N D M E T H O D S
The applied numerical plasma simulation model is a Cloud-in-Cell, ion-kinetic, quasi-neutral, global hybrid code with a divergencefree face-centred magnetic field solver (Kallio & Janhunen 2002;Sillanpää 2008). Briefly, the model treats ions as kinetic, macroscopic particle clouds (macroparticles), and electrons as a massless, chargeneutralizing fluid. The magnetic field is solved via an upwinding solver and propagated on a Yee-type, divergence-free mesh. The computational mesh is Cartesian, with optional hierarchical grid refinements. The kinetic ion clouds are propagated by the Lorentz force. The electric current and bulk velocity for the charge-neutralizing electrons are obtained from the ion current and total electric current computed from Ampere's law. Generalized Ohm's law is then used to solve for the electric field from the electron bulk velocity (including the Hall term), isothermal electron pressure, and a resistivity term. Faraday's law is then used to propagate the magnetic field, closing the computational loop. The model set-up has been described in detail in Papers I and II, and we refer the reader to those articles for further details.

Coordinate systems
To facilitate the use of the Hansen et al. (2016) neutral model, we added functionality to include auxiliary coordinate systems to the simulation software. We base our coordinate systems on the SPICE kernel framework (Acton et al. 2018), using the Rosetta kernels provided by European Space Agency SPICE Service (2019). Hansen et al. (2016) employ a cometocentric frame aligned with the sunward direction as the principal axis (+X) and the component orthogonal to X of the spin axis of the 67P as the secondary axis (+Z), with the Y-axis completing the right-handed frame. 1 On the other hand, our simulation assumes that the solar wind flow is aligned with the simulation frame X axis, and that the IMF lies on the XY plane. In the case of nominal Parker IMF (Parker 1958) and no solar wind aberration, this corresponds to the Cometocentric Solar Equatorial (CSEQ) reference frame. 2 After the inclusion of the neutral model necessitating an auxiliary coordinate system, it is trivial to include solar wind aberration due to the orbital velocity of the comet. For IMF orientation along the Parker spiral, the aberration accounts for an approximate difference Inbound (Northern autumn) equinox (purple square), perihelion (green diamond), and summer solstice (yellow triangle). For the chosen simulation cases, we show Parker spiral (blue), radial solar wind flow (red arrow), the spin axis of 67P (black arrow), so that the tails of the red and black arrows coincide with the location of 67P at that point in time. The angular component of the outgassing model (Hansen et al. 2016) for each case, equivalent to an antenna beam pattern, is shown as well as a coloured, translucent spheroid, with both colour and distance from the comet at that time showing the value of the angular component. Note the hemispherical flip for the outgassing model during inbound equinox. of 7 per cent to the magnitude of the convective electric field in the cometocentric frame of reference around perihelion. While the solar wind is highly variable, as shown e.g. by MAVEN (Liu et al. 2021), it is not feasible to cover the full solar wind parameter space in the simulations. Instead, we use an essentially constant solar wind profile consistent with its statistical properties, to document trends along the orbit of 67P stemming from the solar wind plasma environment. Appendix A provides more details on aberration effect.
For our simulation coordinate system (SIM), we choose to use the aberrated solar wind velocity and rotate our simulation frame correspondingly, so that the +X axis is antiparallel to the aberrated solar wind flow, the +Z axis is along the convective electric field in the upstream region (assuming away-sector Parker IMF) and the Y axis completes the right-handed coordinate system.

Selection of simulated events
As the parameter space is large, considering the variability of the solar wind and IMF, as well as the temporal and spatial variability of cometary outgassing. Thus, we select a set of representative average cases along the orbit.
For our simulations, we use the inbound leg during 2015 March 27-July 30 with heliocentric distances ranging from 2 au to near perihelion at 1.25 au. Further, we selected nine points equidistant along the orbit of 67P, with a spacing of 0.2685 au. This includes the autumn equinox and approach to perihelion. The solstice just after perihelion and the outbound equinox at 2.6 au are excluded. CO 2 outgassing increases after perihelion and overtakes the water group after the start of 2016 (Hoang et al. 2019), where the wateroutgassing model we use would not be accurate. Post-perihelion carbon-group outgassing might be an issue by itself, with few available cross-sections to describe e.g. the CX processes, but the high outgassing rates after perihelion also limit stable operation of our simulation at the available resolutions. Thus, we limit this study to the inbound phase both to avoid the CO 2 issues and to avoid modelling periods with high outgassing rates, and stay within the range of dates modelled by Hansen et al. (2016). Fig. 1 shows the selected events along the orbit of 67P during the inbound leg, and Table 1 lists the events and their parameters in detail.

Solar wind and IMF
We produce nominal solar wind parameters for each heliocentric distance using Slavin & Holzer (1981), fixing the solar wind parameters at 1 au and scaled according to the relations given. Fixed solar wind values at 1 au are: (i) Solar wind radial velocity 430 km s −1 ; (ii) Proton density 7 cm −3 ; (iii) Proton temperature 80 000 K; (iv) Electron temperature 150 000 K; and (v) IMF strength 6 nT.
For simplicity, we omit the solar wind alpha particle population, and leave the multi-ion solar wind effects on the solar wind dynamic pressure and CX reactions for further studies. Crossing the heliospheric current sheet results in change of IMF polarity, as seen in the nested magnetic draping patterns at comets (Volwerk et al. 2016). We consider both IMF polarities, as they produce asymmetric outgassing profiles. This yields two sets of runs for each event, with 'away' and 'towards' IMF sectors titled according to the IMF vector pointing away or towards the Sun. Table 1 lists the solar wind and IMF parameters for each heliocentric distance along with other simulation parameters. As we limit our study to a steady-state upstream condition, we will not be able to capture nested draping.
To further probe sensitivity of the model to solar wind and IMF parameters, considering the asymmetry of the outgassing profile, we simulate different clock angle cases of the IMF at one heliocentric distance (1.34 au). We generate upstream conditions with the Parker IMF rotated around the CSEQ X-axis at 45 • increments, yielding an additional set of eight runs. The IMF parameters varied from the base run (Run 7, away sector) are listed in Table 2. Here, we employ an unconventional datum for the IMF clock angle, showing the angle as relative to the nominal Parker spiral.

Neutral model
To model the large-scale trends in cometary outgassing, we base our neutral profile on the Hansen et al. (2016) empirical, asymmetric outgassing model, which is, in turn, based on Direct Simulation Monte Carlo (DSCM) modelling (Fougere et al. 2016) and averaged over the rotation period of the nucleus. The Hansen model is given by where (θ , φ, r) are cometocentric coordinates (latitude, longitude, and cometocentric distance in the 67P/C-G SUN SPIN frame), d is the date, whereas f is the empirical angular factor, Q 0 the total neutral production rate, and the radial velocity of outgassing. Setting the angular factor f = 1 recovers the spherically symmetric Haser model (Haser 1957) without the exponential loss term. Although including a full DSMC solution for the outgassing profile could affect the resulting environment, any analysis of diurnal variations is left for further studies. The Hansen et al. (2016) model contains a radial dependence in its angular factor up to the cometocentric distance of 100 km, that is, the end of their model domain. We adopt the angular fit from Hansen et al. (2016) at their maximal cometocentric distance, assume radial gas expansion to larger cometocentric distances to the full domain. That is, we drop the radial dependence from the angular factor, and refrain from extrapolating the radial dependence of the Hansen et al. (2016) model. Fig. 1 shows the angular component of the outgassing profile used for the chosen events.
For the radial part, we include a Haser-like loss term, so that our neutral profile is where r is the cometocentric distance, d is the date, the neutral outgassing velocity, and λ(d) = (d) ( i ν i (d)) −1 is computed from reaction frequencies ν i due to PI, CX, and EI in the undisturbed solar wind. The neutrals are assumed to be H 2 O, with no dissociation processes. We note that the photodissociation of H 2 O would also produce a multispecies coma containing hydrogen and other daughter species as described by the multigenerational Haser-like models (Mendis,  1972;Festou 1981a, b). In our approach, heavy water-group particles (H 2 O, OH, O) are assimilated in to the H 2 O population. The effect of the hydrogen coma to the mass-loading and CX reactions is expected to be small as shown by Simon Wedlund et al. (2019b), so we discard the hydrogen coma altogether in our simulations.

Photoionization, Charge exchange, and recombination
Since the H 2 O population is assumed to be the only cometary population without any dissociation products, we approximate the extent of the H 2 O coma with our neutral model (equation 2), with the neutral sink consisting of the sum of reaction frequencies for PI, CX, and EI, the latter two calculated using upstream solar wind parameters, and PI by scaling photoionization rates by the inverse square of heliocentric distance. We do not take into account proton and hydrogen ENA impact ionization that dominate at energies above the nominal solar wind energies (Simon Wedlund et al. 2019a). Cross-sections and reaction rates are handled as in Paper II, with H + + H 2 O CX cross-sections updated to reflect the results in Simon Wedlund et al. (2019a). The energy-dependent CX cross-sections below the lower limit in Simon Wedlund et al. (2019a) are linearly interpolated between the lowest energy cross-section and zero at zero energy.
EI reaction rate is fitted to Cravens et al. (1987), and electron recombination reaction rates for H 2 O + + e − are taken from Hollenbach et al. (2012) and fitted as in Simon Wedlund et al. (2017), both dependent on the electron temperature, which varies with heliocentric distance. The recombination process for H + + e − is not included.

Simulation domain and grid
The time period we focus on covers a large range of outgassing rates, which leads to a large variation in the simulation domain sizes. To determine the domain size for the simulation set-up, we performed a rough estimation of mass-loading effectiveness, based on the extended upstream boundary condition as in Paper II. We calculate the column density of cometary ions from far upstream and the corresponding momentum density along the X axis, assuming the following conditions: (i) Constant solar wind and a time stationary solution; (ii) Heavy ions from a given neutral profile, with PI and solar wind-caused ion production (CX, EI); where i refers to ion species i, r to cometocentric position and d h the heliocentric distance of the comet; and (iv) Instant acceleration of cometary ions to the upstream solar wind speed U sw .
Further, we compare the momentum density of the produced cometary ions p ci = m i n i U i (m i is the mass, n i is the number density, and U i is the bulk velocity of ion species i with U i assumed equal to U sw ) against the solar wind momentum density p sw , and find the value of x for which p ci /p sw = 1/15. This threshold was found to strike a balance with minimum domain size and reasonable adherence to the above assumptions. Since the outgassing profile is asymmetric with respect to the sunward direction, we numerically search for the maximum X max fulfilling this condition to set the upstream domain extent. Domain extent to other directions is expressed in terms of X max : downstream boundary is set to −0.5X max , and the side boundaries to ±2X max . The same estimation is used to approximate upstream massloading effects for our solar wind inputs. The method was introduced in Paper II and shown to produce qualitative agreement in terms of upstream extent and solar wind deflection with a mass-loaded MHD model by Rubin et al. (2015), who modelled large domains of up to 16 × 10 6 km. The extended upstream condition is formed using the above assumptions: For each upstream boundary cell, the cumulative exchange of momentum from solar wind to pickup ions is calculated and accounted for in terms of corresponding deflection and deceleration of upstream solar wind injected into the simulation. Unlike the extended boundary condition method of Koenders et al. (2013) or the analytical model of Flammer & Mendis (1996), our method does include the finite gyroradius effects.
The base grid resolution is chosen to encompass 30 cubical grid cells from the comet at the origin to the front wall; efforts to increase domain size with coarse upstream grid and close-in mesh refinements were found to produce unstable solar wind behaviour in the coarse upstream cells. t is set so that x/ t = 7500 km s −1 . A resistivity term is introduced to ensure stability using the same method as in Paper II: We set the resistivity term to such a value that the magnetic Reynolds number R m = U sw x/η = 4, where U sw is the upstream solar wind velocity, x is the grid size, and η is the magnetic resistivity.
After the simulation reached a quasi-stationary state of development, simulation particles were collected for a period of time near the nucleus for detailed spectral analysis. The particles were collected when crossing into a virtual detector: a sphere centred at the nucleus with a radius of 1.5 x. See Table 1 for the varying data collection timings, dependent on e.g. the required initialization time.

Spectral break determination
Fig. 2 describes the method of finding the spectral breaks and their respective error estimates. The method for finding the break energies is the same as in Paper II, but we have refined the error estimation routine to account for flat-topped, asymmetric peaks in the peak finding algorithm. As in Paper II, we calculate the spectral index α = log j/log E, where j is the omnidirectional differential particle flux and E the energy, from lowess-smoothed omnidirectional flux of cometary ions j, and find the peaks of the derivative δ(E) = dα(E)/dE.
The breaks are chosen with the following criteria from the local extrema of δ(E): For the ankle, we take the local maximum with the largest prominence within E = ]100, 20 000[ eV (excluding the binning boundary effects at high energies and non-representative low energies). For the knee, we take the largest-prominence local minimum with energy below the ankle energy. For the ankle (knee), we find the upper and lower energy error estimates by finding the closest intersections of the peak with a threshold value of peak maximum (minimum) value minus (plus) half of the peak prominence, which is illustrated in Fig. 2. Intersections are found via linear interpolation.

R E S U LT S
In this section, we present results from the simulation runs. Section 3.1 describes one simulation in detail, while Section 3.2 discusses evolution of the main parameters of interest over the orbit of 67P. The supplementary material presents additional figures of the various runs.

Simulation of 1.34 au in detail
Using our away sector Run 7 as the prototype case, Fig. 3 demonstrates the usual features of our simulation of cometary environments, in the case of a sufficiently strong interaction to create a shock: Mass-loading of upstream solar wind flow, shown by deceleration and deflection before the sharp drop in velocity interpreted as the cometary bow shock, and the shock itself visible here as the sharp velocity drop in solar wind ions and the sharp draping pattern in IMF field lines and the jump in IMF magnitude. The draping pattern within the sheath is consistent with our previous works, showing no diamagnetic cavity and shallow draping after the nucleus due to low resolution.
Figs 4-8 present an overview of results in Run 7, away sector IMF; similar presentation is given for the entire run set in the supplementary materials: Supplement 1 for the away sector, Supplement 2 for the towards sector, and Supplement 3 for the clock angle set. Run 7 is taken here as an example of a simulation producing clear spectral signatures connected with the magnetosonic Mach number M ms = 2 surface used as a shock proxy. The interpretation of M ms = 2 coinciding with the shock is based on Galeev & Khabibrakhmanov (1990), used previously in Papers I and II, and found to be a reasonable indicator of a mass-loaded shock. Fig. 4 shows selected parameters along the X SIM -axis, from front wall on the right, through the nucleus at the origin and down the tail on the left. Solar wind density and IMF values are normalized to upstream values far from the comet. Dimensionless plasma parameters β and M ms,H + are shown as well. Notably, the M ms curve shows a steady decrease approaching the comet from the upstream, with a sharp drop to values below 1 through a narrow region, coincident with enhancements in plasma density and magnetic field strength. This is interpreted as a proper shock in solar wind protons.  Medium grey shading with black outline shows flux from particles ionized from within the region in which M ms < 2. Areas shaded with darker and lighter grey correspond to successive regions with M ms < 1 (dark grey) and M ms < 3 (light grey). The solid vertical line denotes the half-sheath flux, that is the energy at which half of the total pickup ion flux is due to pickup ions ionized from within the sheath region. The horizontal coloured bars, numbered 1-7 annotate different energy regions inferred from spectral breaks, illustrated in the bottom panel. Bottom panel: spectral index α(E) and its derivative δ(E), used to select energy ranges. Energies for spectral knee (3299 eV) and ankle (5762 eV) are shown at the corresponding points of δ(E), whereas other, less prominent, δ(E) extrema are marked with dots. The energy of half-sheath flux is also shown. See Section 2.3.3 for details.
Vertical lines annotate the subsolar extent of the M ms = 2 surface with ± x uncertainty.  Section 2.3.3 describes estimation of the spectral breaks in more detail. The figure presents the total pickup ion spectrum (black). To estimate the connection of the spectral features to the different plasma regions, the shaded regions show the contributions to the spectrum from ions originating from these regions, especially middle grey with black outline corresponding to the M ms < 2 region. It is evident that the spectral knee (the low-energy edge of the first major drop in the flux when moving from lower to higher energies) at 3299 eV and the spectral ankle (the high-energy limit of the steep slope following the knee) at 5762 eV delimit the transition where the contribution of ions from the M ms = 2 region drops substantially. This is in line with results in Paper II. Fig. 6 shows the extent of the M ms < 2 region (grey translucent surface) and the points of ionization for pickup ions collected by the virtual detector (coloured dots). The ionization points are coloured according to the energy of the ion at detection, with colours matching the energy regions shown in Fig. 5. The distribution of the ionization points shows the relationship of the ion energy to the different plasma regions of origin, but there are some regions with mixed contributions (e.g. upstream particles in pink, and highest energy particles from within the M ms < 2 region in pink as well). We note that the points of ionization for pickup ions at these energies lie approximately around a quasi-perpendicular shock. The projections of the ionization points are shown on the principal planes in grey, which shows that the trajectories on the XZ plane fit the cycloidal motion of the pickup ions, and that the XY plane reflects the effect of the Parker angle with a strong asymmetry in the Y direction. Fig. 7 shows the velocity distribution of the collected pickup ions in the CSEQ frame, with the ions coloured according to their energy at detection time and projections on principal planes again in grey. The distribution consists roughly of two parts: A low-energy core population below approximately 4 keV, with antisunward motion, a nearly complete ring beam distribution, and approximate symmetry in y , that is, pickup ions from within the sheath region; and a high-energy population above 4 keV and up to maximum pickup ion energies, with a ring distribution slanted around the IMF in the upstream region.  Fig. 6 for energy ranges), with projections at principal planes in grey. The orange triangle annotates the upstream U SW vector, while the orange arrow shows sunward direction (1, 0, 0) CSEQ , red arrow shows the direction of upstream convective electric field (0, 0, 1) CSEQ , and green arrow shows the direction of upstream IMF (−2.34, 3.17, 0) CSEQ nT. Fig. 8 describes the pickup ion flux in an observation-focused format: The six ellipsoids on the top describe the differential flux in given energy windows (corresponding to energy regions in Fig. 5) in Hammer equal-area projection. The four bottom panels show histograms of the flux of four angular variables against energy: latitude, longitude, cone angle, and clock angle. The orientation of the figure is such that the angular distributions are for inbound ions, i.e. as if the observer is at the origin (or nucleus of the comet) and looking outward. The energy dispersion properties of the angular distributions of the inbound fluxes of high-energy pickup ions are related to the spectral breaks (see Fig. 5 and Section 2.3.3, shown as dashed vertical lines in the bottom panels). As such, the angular distributions provide an independent method for identifying the large-scale plasma structures. Fig. 9 shows the evolution of simulated pickup ion spectra close to the nucleus, over the inbound phase of the orbit from just below 2 au, where we start to see signs of interaction in the spectra, to the perihelion. The shown set of runs assumes nominal away-sector IMF and solar wind parameters from Slavin & Holzer (1981). The inbound equinox takes place between Runs 3 and 4, where the hemispheric asymmetry of the outgassing flips to prefer southern latitudes for the remainder of the run set.

Evolution of pickup ion energy spectra
At 2.0 au (not shown) and at 1.73 au the pickup ion spectra show little effects from the plasma environment, and proper shock-like regions using our working definition of M ms,H + < 2 are not visible. This corresponds to the weakly interacting cases of Paper II, with the spectrum displaying a non-broken power-law shape. Moving closer to the perihelion, we start to see both the M ms < 2 region forming and the pickup ion spectrum showing breaks, which we describe by knee (spectral break around 1000-5000 eV, at the lower energy limit of the steep slope) and ankle (around 3000-9000 eV, at the high-energy end of the steep slope). The common features in the spectra (Fig. 9) from low to higher energies, when the breaks are present, are a flat spectrum below the knee energy (starting from Run 2), followed by a steep spectral slope with a precipitous drop in pickup ion flux (from Run 3), a plateau in the spectrum after the ankle energy (from Run 3), and in some cases even an enhancement in flux (esp. Run 8), and abrupt fall-off around the theoretical maximum energy for pickup ion energy (from Run 5).
For context, Fig. 9 also shows pickup ion fluxes classified by the region of ionization of the particles. Ions are defined as originating from within a region if their point of ionization lies within the volume of the region. We show ions originating from different solar wind Mach number volumes for which M ms,H + < M i , so that M i ∈ [1, 2, 3] as successively lighter grey portions of the spectra, with a black outline at M ms = 2. For Runs 6 to 9, the ion fluxes for regions M ms < 1 and M ms < 2 nearly coincide, which we take as an indication for a strong shock, compared to the ones in Runs 3 to 5.
Comparing to results in Paper II, and the clear correspondence there between the M ms = 2 surface and the spectral break observed, we notice that at large heliocentric distances (>1.51 au), the correspondence is no longer there. This is surprising, as the results of Paper II were taken to be more descriptive of a weak cometary environment, due to the spherically symmetric outgassing profile used. On the other hand, including an asymmetrically sunwardweighted outgassing profile was expected to modify the environment, and we return to these differences in Section 4. Notably, as the correspondence between M ms = 2 and the spectral breaks sets in properly, the model starts to produce a significant portion of ions from the M ms = 1 regions, that is, from the properly subsonic solar wind flow regions. This is not the case for Runs 3-4: In Run 3, we see no M ms < 1 contribution, but some M ms < 1.5 contribution. Fig. 10 shows the relationships between the shock standoff and chosen characteristic energies for the away and towards sector runs. Run 1 is not shown for either IMF sector, as the spectra are too smooth to draw conclusions. The spectral features in Run 2 (both sectors), while very smooth, are still appropriately picked up by the algorithm and thus retained for completeness. We do not draw strong conclusions from these data points in Run 2.

Spectral break analysis
The breaks in the energy spectra (ankles and knees, marked with vertical lines) are clearly identifiable in Runs 3-9. The knee energy corresponds to the low-energy limit of the prominent dip in the pickup ion flux spectrum, the ankle to the high-energy limit and the 'half-sheath flux' is the highest energy at which the simulated pickup Figure 9. Away sector IMF, overview of pickup ion spectra evolution for Runs 2 (top) to 9 (bottom). Each panel describes a simulation at a given heliocentric distance. Each panel shows simulated pickup ion flux energy spectrum close to the nucleus. Thick black line shows the total pick ion flux. Grey shading with a black outline shows flux from particles ionized from within the region in which M ms < 2. Darker and lighter grey shaded areas correspond to successive regions with M ms threshold at values M ms < 1 (dark) and M ms < 3 (light). In Runs 6 to 9, the M ms < 1 and M ms < 2 shadings cover nearly similar areas. Figure 10. Inferred energies versus 'shock' standoff distance, IMF away (top, A) and towards (bottom, B) sector run-sets, excluding Run 1. Run IDs are listed in the figure. Blue: Energy of the spectral ankle (between the dip and the plateau), red: Energy at which half of the flux is from within the M ms = 2 'sheath' and half from upstream, yellow: energy of the spectral knee (low-energy limit of the steep dip). Shock standoff distance increases towards perihelion. Error bars for spectral breaks are represented as upper and lower bounds for the peak width at half prominence (see text for details).
ion flux is produced in equal parts from volumes of M ms > 2 and M ms ≤ 2.
The half-sheath flux connects the spectral features to the shock at the four runs closest to perihelion. The furthest two runs (furthest one not shown) do not show a clearly defined shock, consistent with Paper II Run 1 (weak coupling/large heliocentric distance), and so the energy of half-sheath flux is not defined for these runs. However, in the transition between these regimes, the M ms = 2 surface does not seem to produce a good correspondence between the spectral features nor the 'half-sheath flux'. When compared with the actual injection point curve, the M ms = 2 surface also does not coincide well with the sharp kink in the injection point curve at what could be surmised to be the shock nose.
While the knee and ankle features are properties of the particle energy distribution and this can be observed by spacecraft, the halfsheath flux is a ground truth parameter from the simulation that uses information about the history of the ion. As such, the correspondence of the half-sheath flux to the knee and ankle features serves as validation and a measure of confidence in the correspondence between the spectral features and the tentative bow shock.

Effect of IMF clock angle
Since the outgassing profile is asymmetric with respect to the Xaxis, we decided to explore the effect of IMF clock angle on the system, that is, keeping the cone angle fixed and rotating the IMF around the comet-Sun line. Notably, at certain IMF clock angle orientations, the upstream convective electric field is parallel to the rotation axis of the nucleus and the direction of maximal/minimal outgassing hemisphere, with the potential of influencing the results. The used IMF configurations, as modified from the away sector Run 7, are given in Table 2, and the clock angle set is formed by rotating the IMF vector around the Sun-Comet line in 45 • increments. The axial tilt of the comet at this time is, coincidentally, approximately 45 • , leading to interesting edge cases: clock angle Run 4 having the upstream convective electric field nearly parallel to the maximal outgassing direction (the Southern hemisphere) when projected on the plane perpendicular to the comet-Sun line, and respectively nearly antiparallel for Run 8. For Run 4, this translates to the upstream pickup ions from the maximal outgassing hemisphere being accelerated away from the nucleus, while for Run 8, the ions from the maximally outgassing direction are accelerated towards the nucleus.
Supplement 3 includes the set of IMF clock angle runs based on Run 7, with a varying IMF clock angle, for which we present the pickup ion spectra in Fig. 11, using the same format as in Fig. 9. We note that the spectra vary considerably with varying IMF clock angle, both in shape and spectral break locations. All spectra still exhibit the precipitous drop associated with ions born within the shock surface and in the upstream region, but at varying energies, and with smoother or sharper transitions. We note that the portion of the pickup ions collected at the nucleus and produced in regions M ms < 1 and M ms < 3 varies cyclically: Runs 4 to 7 few or none ions produced in regions with M ms < 3 outside of the shock region, and less ions produced in the region with M ms < 1; and vice versa for Runs 8 and 1 to 3. Fig. 12 shows the effect of IMF clock angle orientation on the spectral break energies obtained from our data reduction pipeline. There indeed is a clear asymmetry with respect to the IMF clock angle and the spectral break energies, especially when the upstream convective electric field is oriented along the spin axis of the nucleus. We also plot the maximum X extent of the M ms = 2 bow shock, which does not appear to be very sensitive to the IMF clock angle. Fig. 12 also shows the dependence of the shock standoff distance with respect to the clock angle. Somewhat unexpectedly, there is only a weak correlation between the IMF clock angle orientation and shock standoff distance, especially when compared to the changes in the spectral break energies that vary significantly. With respect to IMF clock angles, the spectral break energies do not seem to have a clear dependence on the shock standoff distance.

ICA observations
A preliminary survey of ICA cometary ion energy spectra was performed to assess whether there would be more observations 3 that can be interpreted as signatures of large-scale plasma features. As shown in 3.4, different IMF clock angles strongly affect the pickup ion spectra. Therefore, we limited our survey to such 4-h spectra during which the standard deviation of the observed IMF clock angle was less than 15 degrees. 3 See Nilsson et al. (2018) for one previously reported observation. Figure 11. Away sector IMF, overview of pickup ion spectra evolution for varied IMF clock angles. Each panel describes a simulation at a given IMF clock angle, at 1.34 au. Each panel shows simulated pickup ion flux energy spectrum close to the nucleus. Thick black line shows the total pick ion flux. Grey shading with a black outline shows flux from particles ionized from within the region in which M ms < 2. Darker and lighter grey shaded areas correspond to successive regions with M ms threshold at values M ms < 1 (dark) and M ms < 3 (light). In Runs 6 to 9, the M ms < 1 and M ms < 2 shadings cover nearly similar areas.

Figure 12.
Energies of spectral breaks (ankle and knee), the energy of halfsheath flux (defined as before) and bow shock standoff distance (scaled) versus clock angle delta for the clock angle runset. The polar angle corresponds to IMF clock angle relative to nominal Parker spiral, with 0 • on right being the nominal Parker field; radius is in eV for the spectral breaks and in units of 3 km for the shock standoff (scaled to lie around the appropriate energy range), markers show the values from the simulations, connected by cubic splines. The relation of upstream convective electric field with the angle between SIM Y and CSEQ SUN SPIN Z is annotated for the corresponding clock angle deltas. Notably, at 315 • , the observed pickup ions originate from the direction of maximum outgassing (Southern hemisphere), and at 135 • from the direction minimal outgassing (Northern hemisphere). Fig. 13 shows three manually selected cases around 1.5 au (Run 5 being nearest) that exhibit spectral features similar to observed in the simulations, but it is acknowledged that other types of spectra appear as well. The spectra consist of averages of heavy ion fluxes over the available viewing directions. We present these possible matches to highlight that a dedicated survey of the spectra, in conjunction with insights from modelling, can be used to further understand the evolution of the cometary plasma environment.
In terms of this paper, we take a characteristic energy for each spectra from the low-energy limit of the high-energy plateau in the ICA spectra through visual inspection. These correspond to the ankles in the simulated observations (1.5 keV, 6.5 keV, and 4.4 keV for 6th, 7th and 8th of June, respectively). Defining a clear knee is easy for the 2015 June 7 spectrum (at 5 keV), but the other two do not exhibit as clear a knee feature. We note that the June 6th observation has 90 degree difference in clock angle to the June 7th and 8th observations, and correspondingly the spectra are noticeably different with respect to the extent of the high-energy plateau.
In comparison, the ankle energy range at Run 5, at roughly corresponding heliocentric distance are ∼ 2000 eV to ∼ 5300 eV for the away sector and ∼ 4500 eV to ∼ 8000 eV for the towards sector. These energies are not incompatible with the observed spectral features, especially allowing for the factor 2 difference in break energies attributable to clock angle variations.
We note that the ICA observations can include several complicating features when comparing with the simulated spectra. Solar wind and rotational outgassing variability during the 4 h integration period may smear the ICA spectra, and it is possible for field-of-view effects to reduce the observed counts. At low energies below 60 eV, ICA usually observes cometary water-group ion outflow from the nucleus, which is not within the resolution of our simulation, and we do not draw conclusions from energies below 1 keV. However, the more energetic pickup ions observed by ICA display mass-loading effects consistent with solar wind mass-loading Williamson et al. 2020), which gives us confidence that we have correctly reproduced the general shape of the pickup-ion spectra at high energies.

D I S C U S S I O N
We study the spectral features of the cometary pickup ions at the comet 67P using a global hybrid simulation, focusing on how they evolve during the inbound phase and how they respond to changes in the IMF orientation. Several of the features defined as spectral breaks are associated with the boundaries in the cometary plasma environment.
The discrepancy with M ms = 2 surface and spectral features raises the question of the use of the M ms = 2 surface as a robust shock measure. The simulation was checked for grid resolution effects by re-running the simulation with a locally refined grid ( x − → x/2). The similar results in both runs lead us to believe that the effect is not dependent on grid resolution.
We speculate that the discrepancy could be caused by the sunwardextended plasma environment due to enhanced sunward outgassing. Compared to a spherically symmetric Haser model, the sunwardasymmetric outgassing presents a larger portion of the coma to the impinging solar wind for a fixed outgassing rate, possibly leading to a more gradual mass-loading process. This is in contrast with the previously used spherically symmetric Haser neutral models, where the impinging solar wind encounters a more compact obstacle.
Paper I used the value Q 0 = 6.56 × 10 27 s −1 for symmetric outgassing at 1.3 au, comparable with Run 8. With those parameters, they found a shock standoff distance of approximately 6600 km, which is smaller by over a factor of 2 than our comparable Run 8 which yields a standoff distance of ∼ 14 000 km. Further, Paper II with its symmetric outgassing shows strong variation in bow shock standoff distance as function of IMF magnitude, where a maximal standoff of approximately 3500 km was produced for Q 0 = 4.14 × 10 27 s −1 . The heliocentric distance for the runs in Paper II was 1.425 au, approximately corresponding to Run 6 in this study in cometary, solar wind, and IMF parameters, except for IMF Parker angle and asymmetric outgassing used here. The corresponding shock standoff in our Run 6 of some 7500-8000 km is again larger by a factor of 2-3. This is in agreement with the results of Huang et al. (2016) at perihelion, with respect to asymmetric outgassing versus symmetric outgassing.
The spectra between Paper II and this work differ in few aspects. First, the energies at which we identify these spectral features are shifted higher in this work, which is to be expected with the increased shock standoff distance with sunward-asymmetric outgassing: The pickup ions can accumulate more energy along their cycloidal trajectory from their point of ionization up to the detection point, as evident by e.g. the nearly complete pickup ion ring beam at low energies. Secondly, the spectra shown in Paper II contain few more features in the lower end of the energy range, that is, below the break energies identified in Paper II. Some of these features would cause problems with the present algorithm in the correct detection of the knee (being the lower limit of the pickup ion energy spectra associated with the shock transition). We surmise these differences, related to low energies associated with ions within the sheath, were most likely due to higher relative spatial resolution near the nucleus in Paper II.
The angular data shown in Fig. 8 support the conclusion of the shock surface being related to the spectral breaks: At the spectral break energies, we see a shift in the pickup ion angular distributions, from a wide distribution ionized within the bow shock to a narrower ring beam-like distribution. Of these distributions, the former is mainly antisunward in longitude, while the latter shows characteristics of E × B-drifting pickup ions in the Parker IMF. The high-energy ring could potentially be used to diagnose upwind IMF, while the transition in the angular distributions to the intra-shock pickup ions could be an independent measurement of the large-scale structure, in addition to the spectral break. Further, the full virtual distribution function might be used to more accurately classify the collected ions, providing more accurate estimates on the locations of the spectral features. This effort is left for further studies.
The work presented here is highly relevant for the upcoming European Space Agency Comet Interceptor mission (Snodgrass & Jones 2019): a multispacecraft mission to a dynamically new comet, comprising of a three-spacecraft constellation. The planned mission profile consists of a flyby at three different trajectories, sampling the cometary environment simultaneously at several distances, with the constellation assumed to pass within 1000 km from the nucleus. Combining the in situ measurements on the mission trajectories supported with remote mapping and modelling of the environment through methods such as in Paper II and this work could present a highly detailed 3D snapshot of the cometary environment.

C O N C L U S I O N S
The effect of asymmetric outgassing on the plasma environment of 67P/Churyumov-Gerasimenko is a significant factor defining the structure of the cometary environment, expanding the cometary plasma environment in the sunward direction. For a given IMF in the simulations, there is a distinct relationship between the spectral breaks and the size of the plasma environment as the comet approaches the Sun, providing a possible diagnostic for remote sensing purposes with Rosetta plasma data. We note, however, that the strong dependence of the spectral features on the IMF magnitude simulated in Paper II was not readily found in in situ observations at comet 67P. Additionally, as we have shown in this work, the IMF clock angle strongly affects the spectral breaks, which produces another change factor to consider -even if the bow shock standoff distance is not strongly affected. Yet, these simulation results can shed light on the large-scale structure of the cometary plasma environment, which help interpretation of observations made in the vicinity of the cometary nuclei.
For some weak cometary activity conditions and neutral profiles, the modelling here shows that the observed spectral breaks may not have as clear quantitative relation to large-scale plasma features. However, the spectral breaks could still be viewed as a likely indicator of their presence, as the strong spectral features still occur as soon as the large-scale M ms = 2 surface develops in the simulations.
Summary of the simulation results: (i) Asymmetric, sunward outgassing expands the cometary plasma environment sunward by a factor of 2 to 3.
(ii) Shock standoff distance is not strongly affected by the IMF clock angle, even when accounting for asymmetric outgassing around the Sun-comet line.
(iii) Spectral break energies of the ankle and the knee in pickup ion energy spectra vary strongly with IMF clock angle, when accounting for asymmetric outgassing around the Sun-comet line.
(iv) The effect of the IMF sector (away/towards) is similar to the effect of varying IMF clock angle.
(v) For a strong interaction (heliocentric distance up to ∼1.51 au), and a constant IMF clock angle and sector, the breaks in the pickup ion energy spectra correlate with the simulated shock standoff distance over the simulated heliocentric distances.
(vi) For a weak interaction (heliocentric distance over ∼1.51 au), the spectral breaks may be somewhat decoupled from the large-scale plasma boundaries.
(vii) The spectral breaks are accompanied by changes in the angular distribution of pickup ions, which can be used as a potential supporting diagnostic.

AC K N OW L E D G E M E N T S
MA acknowledges the Aino and Kaarlo Tiisala foundation for financial support. The work was supported by the Academy of Finland (Decision No. 310444). We acknowledge the computational resources provided by the Aalto Science-IT project. The work of C. Simon Wedlund is funded by the Austrian Science Fund under project number P 32035-N36. The VisIt 3D visualization software (Childs et al. 2012) was used for data exploration and analysis, and we are thankful for this open-source software, even as no VisIt-produced images are presented in this work.

DATA AVA I L A B I L I T Y
Data available on request. The simulation data underlying this article will be shared on reasonable request to Esa Kallio or Riku Jarvinen. RPC-ICA data are available via the ESA Planetary Science Archive (PSA, psa.esa.int Besse et al. 2018).

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at MNRAS online.

Figure 1.
Overview of the IMF away sector run-set. Figure 2. Overview of the IMF towards sector run-set. Figure 3. Overview of the IMF clock angle set at 1:34 au.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

A P P E N D I X : A B E R R AT I O N
To calculate the effect of solar wind aberration due to the orbital velocity of the comet, we compare the upstream convective electric field E z = (− U sw × B IMF ) z and solar wind velocities U sw between the non-aberrated CSEQ coordinate system and the aberrated SIM frame, which has been rotated to align the upstream solar wind with the X axis. Fig. A1 shows the ratios of these values and the heliocentric distance of the comet 67P. This paper has been typeset from a T E X/L A T E X file prepared by the author.