Interaction between a Coronal Mass Ejection and Comet 67P/Churyumov–Gerasimenko

The interaction between a coronal mass ejection (CME) and a comet has been observed several times by in situ observations from the Rosetta Plasma Consortium, which is designed to investigate the cometary magnetosphere of comet 67P/Churyumov–Gerasimenko (CG). Goetz et al. reported a magnetic field of up to 300 nT measured in the inner coma, which is among the largest interplanetary magnetic fields observed in the solar system. They suggested the large magnetic field observations in the inner coma come from magnetic field pileup regions, which are generated by the interaction between a CME and/or corotating interaction region and the cometary magnetosphere. However, the detailed interaction between a CME and the cometary magnetosphere of comet CG in the inner coma has not been investigated by numerical simulations yet. In this paper, we will use a numerical model to simulate the interaction between comet CG and a Halloween class CME and investigate its magnetospheric response to the CME. We find that the plasma structures change significantly during the CME event, and the maximum value of the magnetic field strength is more than 500 nT close to the nucleus. Virtual satellites at similar distances as Rosetta show that the magnetic field strength can be as large as 250 nT, which is slightly less than what Goetz et al. reported.


INTRODUCTION
The Rosetta Plasma Consortium (RPC, Carr et al. (2007)) onboard Rosetta (Glassmeier et al. 2007) was designed to study the physical and chemical properties of comet 67P/Churyumov-Gerasimenko (CG), the evolution of the cometary ionosphere and magnetosphere, as well as examining how the interaction region of the solar wind and the comet is developing.The investigation of the cometary plasma environment and the locations of plasma boundaries for an active comet started with the Giotto flyby of comet 1P/Halley.Different plasma boundaries have been identified, such as a bow shock, which forms upstream of the comet where the solar wind slows down from supersonic to subsonic speeds (Galeev et al. 1985;Koenders et al. 2013).And if the comet is very active with a high neutral gas production rate near perihelion, a diamagnetic cavity (inside which the magnetic field drops to zero) is formed because the outward pointing ion-neutral drag force is large enough to balance the inward-pointing magnetic force (the sum of the magnetic pressure gradient and magnetic curvature forces) (Neubauer et al. 1986;Cravens 1986).Recent studies by Goetz et al. (2016a,b) and Henri et al. (2017) suggested that instabilities may extend the size of the diamagnetic cavity to become larger than expected from numerical simulations or theoretical predictions.Besides, an inner shock, where the cold supersonic cometary ion outflow (from the surface of the nucleus) slows down to subsonic velocities (as it reaches near zero velocity near the contact surface along the sun-comet line), is also proposed (Gombosi 2015).
The Rosetta mission has dramatically improved our understanding of the cometary plasma environment since its arrival at comet CG in 2014.One of the most important discoveries is that the diamagnetic cavity is somewhat different from our previous understanding of a global diamagnetic cavity (Goetz et al. 2016a,b).Many small-scale structures along the diamagnetic surface are reported, suggesting that the diamagnetic cavity of comet CG differs from a traditional, global-scale diamagnetic cavity (like at comet 1P/Halley).The small-scale structures are believed to arise from plasma waves along the global diamagnetic cavity boundary (Goetz et al. 2016a,b).The electron environment near the diamagnetic cavity has been extensively studied (Henri et al. 2017;Odelstad et al. 2018;Engelhardt et al. 2018).Readers can also refer to the review article by Goetz et al. (2022).Huang et al. (2018) provided an alternative explanation that the Hall current along the cavity boundary can introduce magnetic reconnection on the day-side and generate weak magnetic field regions outside the global diamagnetic cavity, similar to the RPC observations.Besides, these shallow magnetic fields were observed much farther from the nucleus than the predicted diamagnetic cavity locations from numerical simulations (Rubin et al. 2015;Koenders et al. 2015;Huang et al. 2016b).Mandt et al. (2016) noticed a 'mystery' boundary separating an inner region and an outer region, which is not predicted by previous studies (Gombosi 2015;Rubin et al. 2015;Koenders et al. 2015;Huang et al. 2016b).They suggested that an ion-neutral collisionopause boundary, inside of which collisions between ions and neutrals are important and outside of which they are not, is responsible for forming the 'mystery' boundary.Recently, Gunell et al. (2018) reported that an infant bow shock was observed several times a few months before and after perihelion, with signatures of an increase in the magnetic field magnitude and fluctuations, an increase in the electron and proton temperatures at the shock and the diminution of the solar wind plasma downstream.
Space weather-like events in comet CG's magnetosphere have also been reported by the Rosetta mission (McKenna-Lawlor et al. 2016;Edberg et al. 2016b,a;Noonan et al. 2018;Goetz et al. 2019).McKenna-Lawlor et al. (2016) studied two CME events at comet CG in September 2014 and noticed a jump in the magnetic field strength as well as ion energy and increased ion count rates associated with the two CME events.Edberg et al. (2016b) investigated the interactions between four Corotating Interaction Regions (CIRs) and the comet and found compressions of the plasma in the inner coma, which can increase fluxes of suprathermal electrons and hence electron-impact ionization rates of the neutrals and elevated plasma densities, as well as the magnetic field strengths.Edberg et al. (2016a) studied a CME impact on comet CG in September 2015 and found a huge compression, which increases the flux of suprathermal electrons, the plasma density, and the magnetic field strength.
They also reported unprecedentedly large magnetic field spikes at a distance of 800 km away from the nucleus and interpreted them as magnetic flux ropes, which are possibly formed by magnetic reconnection in the coma or strong shears causing Kelvin-Helmholtz instabilities in the plasma flow.Goetz et al. (2019) reported that a significant increase in the magnetic field magnitude to about 300 nT was observed on July 3, 2015, which is the strongest magnetic field ever measured in the plasma environment of a comet.They used ENLIL (Odstrcil 2003) to simulate the solar wind conditions during that period at comet CG and concluded that the unusual behavior was caused by the impact of a combination of an interplanetary CME and a CIR.
Numerical modeling has been frequently applied to study the cometary magnetosphere.
There are three major approaches: the fluid approximation (Gombosi et al. 1996;Hansen et al. 2007;Rubin et al. 2014Rubin et al. , 2015;;Huang et al. 2016b;Baranov et al. 2019), the hybrid models (Bagdonat & Motschmann 2002;Koenders et al. 2015;Simon Wedlund et al. 2017;Lindkvist et al. 2018;Alho et al. 2019), and fully kinetic models (Deca et al. 2017;Deca et al. 2019;Gunell et al. 2019;Divin et al. 2020;Beth et al. 2022).In the fluid approach, the plasma is treated as one or multiple fluids, and their motion is described by the magnetohydrodynamic (MHD) equations, while in the hybrid approach, the ions are treated as individual particles, and the electrons are treated as a massless fluid.Recently, fully kinetic (sometimes called Particle-In-Cell, PIC) models, which treat both ions and electrons as particles, have been developed (Deca et al. 2017;Deca et al. 2019;Gunell et al. 2019;Divin et al. 2020) to study the kinetic behaviors of both ions and electrons in the cometary magnetosphere.All these models have successfully modeled the cometary plasma environment, and the simulation results generally agree well with observations.However, there are limited attempts to simulate the impact of space weather events on a comet, mainly due to the high computational cost required for running long-duration simulations to cover the passage of a CME over the cometary magnetosphere.Jia et al. (2007) simulated the crossing process of a typical heliospheric current sheet at a comet with their single fluid MHD model and showed that disconnection events can form in the cometary plasma tail.In a similar approach, Jia et al. (2009) simulated the interaction between a CME and a comet and showed that the CME can trigger tail disconnection events.In this manuscript, we use our state-of-the-art 3-D multi-fluid plasma-neutral interaction model to simulate the interaction between a CME and comet CG and, investigate the cometary magnetospheric response, especially in the inner coma environment, and compare the results with Rosetta observations.

METHODOLOGY
The 3-D multi-fluid plasma-neutral interaction model has been used to simulate the global plasma boundaries of comet CG (Huang et al. 2016b) as well as understanding the diamagnetic field regions (Huang et al. 2016a(Huang et al. , 2018)).Rubin et al. (2014) have shown that multi-fluid plasma models can capture the effect of gyration of different ion fluids, and the simulation results, within the limitations of the fluid approach, show good agreement with a hybrid model (Müller et al. 2011).
In addition, our multi-fluid MHD model has been applied in simulating the plasma environment of Mars (Najib et al. 2011;Bougher et al. 2015;Ma et al. 2011;Dong et al. 2018b,a) and Europa (Rubin et al. 2015;Jia et al. 2018;Harris et al. 2021).These studies demonstrated that the model accurately describes the plasma environment, as confirmed by comparisons with in-situ observations.
Here, we briefly describe our state-of-the-art 3-D multi-fluid plasma-neutral interaction model.A more detailed discussion can be found in Huang et al. (2016bHuang et al. ( ,a, 2018)).Four different fluids are simulated in the model: one for the cometary neutral gas, two for ions (cometary water ions and solar wind protons), and the last one for electrons.The cometary neutral gas fluid is described by the Euler equations as it does not interact with the magnetic and electric fields.In contrast, the plasma (ions and electrons) are represented by the multi-ion MHD equations (Rubin et al. 2014;Huang et al. 2016b).The source and loss of the ions and electrons are adequately described by elastic and inelastic collisions between different fluid particles, various chemical reactions (e.g., photo-ionization, electron impact ionization of the cometary neutral gas), charge exchange between neutrals and ions, as well as ion-electron recombination.The magnetic field is obtained from Faraday's law, neglecting the Hall velocity and the electron pressure gradient term in the induction equation.
The Euler and multi-ion MHD equations can describe the behavior and interactions of the cometary neutral gas, the cometary ions, the solar wind protons, and the electrons self-consistently.These equations are solved by the BATS-R-US (Block-Adaptive Tree Solarwind Roe-type Upwind Scheme) code (Powell et al. 1999;Tóth et al. 2012) within the Space Weather Modeling Framework (SWMF, Tóth et al. (2005Tóth et al. ( , 2012)); Gombosi et al. (2021)) on a 3D block adaptive grid.The adaptive grid can resolve different length scales in a global system, which is critical for studying the interactions between solar wind and comet CG.The nucleus in the model can be either an idealized comet with a spherical shape or the realistic shape model of comet CG (Preusker et al. 2015).At the inner boundary, the neutral gas outflow is illumination-driven based on the empirical relation suggested by Davidsson et al. (2007); Tenishev et al. (2008).The ions, electrons, and magnetic field follow a reflected boundary condition, meaning they cannot penetrate the body.At the outer boundaries, a floating boundary condition (zero gradient) is applied, except for the upstream boundary, where the solar wind plasma and the interplanetary magnetic field (IMF) conditions are prescribed.This work focuses on the situation where comet CG is near its perihelion at 1.3 AU.We use the same photo-ionization rate (6×10 −7 s −1 ) as in Huang et al. (2016b).We apply a spherical body with the radius of 2 km, and the neutral gas outflow is illumination-driven to take into account the day-night asymmetry.The computational domain extends ±1.05 × 10 6 km in the X direction, and ±5.24 × 10 5 km in both Y and Z directions, with 18 levels of refinements, increasing the resolution by a factor of two at every level.The smallest cell is 0.125 km near the nucleus, while the largest cell is about 16,384 km far away from the comet.There are about 26,000 grid blocks with 1.64 million grid cells.We use the body-centered Solar EQuatorial (CSEQ) coordinate system, so the upstream boundary condition is specified in the +X direction.
In order to simulate the interaction between a CME and comet CG, the upstream solar wind and IMF conditions need to be prescribed as input.Unfortunately, RPC did not provide solar wind measurements because the Rosetta spacecraft remained within the cometary magnetosphere while in orbit around the comet and was typically pointed in a very different direction than the Sun.
One of the motivations of the study is to investigate if an extreme CME can produce the unusually high magnetic fields (about 300 nT) in the inner coma reported in Goetz et al. (2019), so we use the solar wind observation during the 2003 Halloween event, which is one of the strongest CMEs observed in recent decades and caused one of the strongest geomagnetic storms (Pulkkinen et al. 2013).To save computational cost, we only simulate the interaction of the comet with the solar wind and IMF measured at Earth between 2003 October 29, 04:00 UT and October 30, 00:00 UT, corresponding to a total duration of 20 hours.The solar wind data is directly obtained from CCMC, which is reconstructed from several instruments on ACE and Geotail to best represent the solar wind conditions due to a large portion of missing observational data from the ACE solar wind instruments (Pulkkinen et al. 2013).As a crude estimation, B y and B z can be scaled as 1/r, while the proton density and Bx can be scaled as 1/r 2 .For a more accurate description, a model (e.g., MSWIM2D by Keebler et al. (2022)) is needed to propagate the solar wind from 1 to 1.3 au.As this study focuses on the cometary magnetospheric response to an extreme CME, it is not critical to scale the solar wind observations to the exact same heliocentric distance.Figure 1 shows the upstream solar wind boundary conditions applied in our model.the CME arrives at the bow shock (Panel (c)): the bow shock is strongly compressed, moving from 10,000 km to around 2,000 km and the proton number density is much higher (> 100 cm −3 ) than during quiet time (< 22 cm −3 ).

The Inner Coma
Figures 3 and 4 provide the time evolution of the impact of the CME on the plasma environment in two perpendicular planes (Y =0 and Z=0).Both figures show the cometary ion (H 2 O + ) and solar wind proton (H + ) densities in the inner coma.Furthermore, the magnetic field magnitude (B) with magnetic field lines are given.We can see that the cometary ion density in the radial distances between 50 km and 200 km in Panels (c) and (d) is smaller than in Panels (a) and (b), but the proton density is higher, which indicates that the cometary ion density is anti-correlated with the proton density.This anti-correlation can be explained by the compression of the cometary ion density, which is caused by the higher dynamic pressure of the solar wind.This compression is confirmed by the higher cometary ion density within 50 km in Panels (c) and (d), as indicated by the streamlines of the cometary ions.
The magnetic field is also significantly compressed as the CME arrives at comet CG: the maximum magnitude just outside the diamagnetic cavity jumps from around 70 nT to more than 500 nT in Panel (c) of Figure 3 and more than 250 nT in Panel (d) of Figure 4.The diamagnetic cavity almost disappears in Panel (c), which is caused by the extremely large dynamic pressure of the CME.Gombosi (2015) provided a rough approximation that the location of the diamagnetic cavity (in the +X direction) is inversely proportional to the solar wind dynamic pressure ρ sw u 2 sw .The maximum solar wind speed during the CME event is about 2200 km/s and the corresponding number density is about 13 cm −3 , as shown in Figure 1.Because we use the same parameters for comet CG (the heliocentric distance, the total neutral gas outflow, and the collision rates between different fluids) as Huang et al. (2016b), the predicted location of the diamagnetic cavity is about 8 times smaller, which is about 12 km in the +X direction.The size of our simulated diamagnetic cavity is consistent with the theoretical expectation.
Asymmetries can also be noticed in the inner coma due to different upstream solar wind flows and IMFs.In Figure 3, we can see that the magnetic field lines are slightly draped in the −Z

Virtual Satellites
We placed five virtual satellites in the inner coma along the three coordinate axes at x = 180 km, y = ±180 km and z = ±180 km, respectively.Extracting the model solutions at these virtual satellites provides predictions of what Rosetta would have measured at those locations during the event.Figure 5 shows the simulated cometary ion density, proton density, magnetic field cone and clock angles (based on the formulas in Goetz et al. (2019)), and magnetic fields at the virtual satellite locations.Fluctuations of varying amplitude in all quantities can be readily noticed, which indicates that there is no stationary state of the plasma environment if a dynamic solar wind condition is applied.The cometary ion and proton densities show similar trends at all virtual satellite locations and the anti-correlation of solar wind and cometary ion densities is also found, which is consistent with the discussion in Section 3.2.Besides, the differences between the solar wind proton densities at y = +180 km and y = −180 km (and similar for z = ±180 km) are associated with the deflation of the solar wind, which is the gyration around the charge-averaged ion velocity, due to different magnetic fields, as indicated by the cone and clock angles.We also notice a cometary ion density increase at around 8 hours of simulation time, which may be associated with the increased electron impact ionization due to the high solar wind temperature, as suggested in Edberg et al. (2016a,b).
The B y and B z components show similar trends at those locations during the event, while the

SUMMARY AND DISCUSSIONS
This paper simulates the interaction between a Halloween-class CME and comet CG.We apply the solar wind data for the Halloween event between 2003 October 29 04:00 UT and October 30 00:00 UT from CCMC and use the 3-D multi-fluid plasma-neutral gas interaction model (Huang et al. 2016b).
Comet CG is treated as an idealized comet with a spherical shape, and the neutral gas outflow is driven by solar illumination.Our results show that the cometary bow shock can be significantly compressed, and its standoff distance is reduced from around 10,000 km to about 2,000 km when the Halloween class CME hits comet CG.We also discover that the cometary ion density in the region between 50 km and 200 km is anti-correlated with the solar wind proton density due to the compression caused by the solar wind.
The maximum magnitude of the magnetic field as observed at three virtual satellites at x/y/z = 180 km is about 250 nT, which is slightly less than the 300 nT observation reported by Goetz et al. (2019).The Halloween class CME is one of the strongest CMEs observed in the last few decades, with propagation speed exceeding 2,200 km/s.The standoff distance of the diamagnetic cavity is inversely proportional to the solar wind dynamic pressure (n sw u 2 sw ), as suggested by Gombosi (2015).A CME with a similar magnetic field strength but smaller dynamic pressure will cause less compression in the inner coma, leading to a larger standoff distance of the diamagnetic cavity.In such a case, the interaction may produce the 300 nT value observed by Goetz et al. (2019) because the most significant magnetic field magnitude is more than 500 nT just outside the diamagnetic cavity in our simulation (See Panel (c) in Figure 3).Nearly all previous numerical simulations used a steady solar wind flow (Hansen et al. 2007;Rubin et al. 2014Rubin et al. , 2015;;Huang et al. 2016b;Simon Wedlund et al. 2017;Lindkvist et al. 2018;Alho et al. 2019).Jia et al. (2007) applied approximated current sheet configurations, while Jia et al. (2009) used a rotation of the IMF vector as a simplified flux rope.This is the first time realistic solar wind data has been used to simulate the solar wind-comet interaction.Our results show ubiquitous fluctuations in the cometary ion density, proton density, and magnetic field, which suggest that the plasma environment is not steady.Furthermore, our simulation assumed an idealized comet with a spherical body and illumination-driven neutral gas outflow to save computational costs.A more realistic and dynamic plasma environment is expected if the real shape (e.g., SHAP5.1 (Preusker et al. 2015), which was previously employed in Huang et al. (2016b)), and the comet's rotation were included in the simulation.
We find that the plasma structures are asymmetric during the event due to the time-varying upstream solar wind flow and IMF.Huang et al. (2018) showed that plasma structures are highly asymmetric due to the Hall current along the diamagnetic cavity boundary.They showed that the Hall current can cause magnetic reconnection on the day-side, which generates weak magnetic field regions outside the global diamagnetic cavity.We expect that the Hall effect would further enhance the asymmetry and cause more complex plasma structures in the inner coma if it is taken into account.
However, a Hall MHD simulation is more than 50 times more expensive than a simulation without the Hall effect.With limited resources after the Rosetta mission was concluded, it is impossible to carry out such a costly simulation for the CME interaction with comet CG.However, we plan to

Figure 1 .
Figure 1.This figure shows the solar wind and IMF conditions for the Halloween event.

Figure 2
Figure 2 shows the simulated bow shock at four different time stamps.It shows that the bow shock changes significantly during the CME event under different solar wind conditions.Huang et al. (2016b) suggested that the bow shock is at about 10,000 km for a fixed solar wind condition with the illumination-driven neutral gas outflow, which is similar to the bow shock location in Panel (a) in this simulation.However, the shapes are different: the bow shock in Huang et al. (2016b) is titled towards -Z direction while the bow shock in Panel (a) is titled towards +Z direction.This is caused by the different solar wind flows and IMFs: only the V x and B y component are non-zero (V x =-400 km/s and B y = 4.8 nT) in Huang et al. (2016b) while all three components of the solar wind velocity and IMF are non-zero (v = [-593, -88, 7.00] km/s, B = [7.65,-5.06, -2.59] nT) at the beginning of the current simulation.In Panel (b), the bow shock is tilted slightly towards the -Z direction, similar to the bow shock in Huang et al. (2016b).Koenders et al. (2013) also investigated how the bow shock distance changed due to different upstream solar wind parameters based on a symmetric neutral outflow.Their predicted distances were smaller than the values from our model because we use an asymmetric neutral outflow driven by solar illumination.Besides, they also discussed the physical mechanics responsible for the discrepancy of bow shock distances obtained from the hybrid model and MHD model under the symmetric neutral background.Significant compression can be seen when

Figure 2 .
Figure 2.This figure shows the simulated H + density at different simulation times (labeled at the top left corner of each subplot).The color range is adjusted each time for better visualization of the compression.

Figure 3 .
Figure 3.This figure shows the simulated H 2 O + and H + densities with the associated velocity streamlines, and the magnetic field magnitude with magnetic field lines in the inner coma in the Y = 0 plane, within 200 km of the nucleus at different simulation times during a Halloween class CME event at comet CG.The simulation time is shown in the bottom of each column.The range of the magnetic field is varied for better visualization of the compression.

Figure 4 .
Figure 4.This figure shows the same physical quantities as Figure 3 but for the the Z = 0 plane.

B
x component varies among different locations.The magnetic field lines are bent to surround the diamagnetic cavity, so B x should be 0 along the x-axis in the upstream direction of the diamagnetic cavity, which explains why B x is almost 0 at x = 180 km.The distinct time evolution of the B x component at y = +180 km and y = −180 km (and similar for z = ±180 km) is caused by the different upstream IMF and how those field lines are bent around the diamagnetic cavity, as discussed in Section 3.2.For the same reason, B x is asymmetric around y = ±180 km and z = ±180 km, as shown in Figure 5.The reversal of B x and B y is also noticed by Goetz et al. (2019); however, B z shows different behavior, which may come from the different magnetic field structures, as they suggested their observations were caused by a CME combined with a CIR.The maximum magnitude of the magnetic field is around 250 nT, which is slightly less than what Goetz et al. (2019) observed.

Figure 5 .
Figure 5.This figure shows the simulated H 2 O + , H + densities, cone and clock angels, as well as magnetic fields (B x , B y , B z and the magnitude B) at the virtual satellite locations.The left panels are associated with the virtual satellites at x = 180 km, while the middle and right panels correspond to the y = ±180 km and z = ±180 km, respectively.