© The Authors 2025

1 Introduction

Long and structured comet tails lead to stunning displays in the sky, which can even be visible to the naked eye in certain instances. Even before the first measurements of the solar wind (SW) were made, observations of comets suggested the presence of a stream of particles originating from the Sun (Biermann 1951) as well as the work of electromagnetic forces (Alfvén 1957). The partially ionised cometary atmosphere presents an obstacle to the solar wind plasma flow and an induced magnetosphere forms around the nucleus. As the outgassing rates vary enormously, both between different comets and along the comets’ elliptical orbits around the Sun, different regimes of SW interaction can be identified (Goetz et al. 2022). At a high-activity comet a fluid approach captures most of the large-scale features of the magnetosphere (Cravens & Gombosi 2004), while at a low-activity comet kinetic effects play major roles (Nilsson et al. 2018).

The Rosetta mission (Glassmeier et al. 2007) to comet 67P/Churyumov-Gerasimenko put focus on the latter interaction domain, where the typical scales of the magnetosphere are smaller than an ion gyroradius. In the low-activity case, an asymmetric magnetosphere is created (Behar et al. 2016, 2018), with clear boundaries appearing only in one hemisphere. A bow shock is likely not present but the first signs of shock formation are occasionally observed (Gunell et al. 2018). The solar wind ions are more deflected than decelerated when they encounter the cometary plasma. In this regime, both solar wind protons and cometary ions exhibit highly anisotropic distributions, at times resembling partial rings (Moeslinger et al. 2024, 2023b,a; Williamson et al. 2022).

The SW cometary plasma interaction is sensitive not only to the properties of the comet itself, but also to the state of the solar wind. For example, a change in the interplanetary magnetic field due to a crossing of the heliospheric current sheet can lead to a disconnection of the cometary tail (Jia et al. 2007). The impact of a coronal mass ejection (CME) may have the same result (Vourlidas et al. 2007). Furthermore, Edberg et al. (2016a) reported dramatic changes also in the inner coma of comet 67P during a CME event. The compression of the magnetosphere and increased ionisation lead to higher densities, stronger magnetic fields, and the presence of magnetic flux ropes.

Alpha particles are the second most abundant SW species after protons, with typical relative number densities n_α/n_H between 3-5% in the steady-state SW. There is a strong dependence of n_α/n_H with solar activity (Kasper et al. 2007; Song et al. 2022). The lowest n_α/n_H (<1%) occurs during solar minimum in the fast solar wind, while n_α/n_H > 20% is possible during CMEs at solar maximum. Additionally, Yogesh et al. (2023) have shown that n_α/n_H also increases rapidly at the surface of stream interaction regions (SIRs), such as corotating interaction regions (CIRs). Close to the Sun, alphas typically have higher bulk velocities than protons (up to 150 km s⁻¹ higher at 0.3 AU), but this difference decreases with increasing heliocentric distance (e. g. 30km s⁻¹ at 1 AU, 20 km s⁻¹ at 2.5-3 AU) (Reisenfeld et al. 2001; Ďurovcová et al. 2017). This differential flow can cause instabilities, such as electron cyclotron waves (Zhao et al. 2019). CMEs and CIRs were observed by Rosetta at comet 67P and their impact on the plasma environment has been studied by Edberg et al. (2016a) and Edberg et al. (2016b), but the influence of the SW composition has so far not been assessed.

An instructive approach to studying the interaction between the solar wind and the cometary coma is using Poynting’s theorem to focus on the energy transfer between the different particle populations and the electromagnetic fields. Poynting’s theorem expresses the conservation of electromagnetic energy (Olsson & Barger 1987) and can be written as $$-\frac{\mathrm{\partial }}{\mathrm{\partial }t}(\frac{B^2}{2{\mu }_0}+\frac{{\epsilon }_0{E}^2}{2})=(\mathrm{\nabla }\cdot \mathbf{S}\mathbf{+}\mathbf{E}\cdot \mathbf{J}),$$(1)

where B = |B| and E = |E| are the magnitudes of the magnetic and electric fields, respectively, μ₀ is the permeability of vacuum, and ε₀ is the vacuum dielectric constant. Hence, the left hand side is the change over time of the energy density stored in the electromagnetic fields. Often B²/2μ₀ ≫ ε₀E²/2, and ε₀E²/2 can be neglected. To conserve energy the change in electromagnetic energy density is compensated by a flow of field energy to or from the region of interest, given by ∇ ∙ S, where S = E × B/μ₀ is the Poynting flux (Poynting vector). If charged particles are present in the system, energy can also be transferred to or from them through work done by the Lorenz force. The term E ∙ J, where J is the current density, represents this energy transfer. A power density E ∙ J > 0 means that electromagnetic energy is given to the charge carriers, that is, the power density is referred to as a ‘load’. A negative power density (E ∙ J < 0) is a ‘generator’, where energy moves from the charged particles to the electromagnetic fields.

Studies of many different objects in the Solar System make use of Poynting’s theorem for the interpretation of their data in terms of energy conversion. Saunders & Russell (1986) explained measurements made in Venus’ magnetotail by investigating the sign of E ∙ J in a simple model based on the magnetic field observations. Multi-spacecraft constellations made direct observations of the power density easier as the current density can be computed as J = ∇ × B/μ₀. A number of papers make use of this to identify generator and load regions related to the auroral current circuit at Earth and energy conversion in the plasma sheet (Hamrin et al. 2006; Marghitu et al. 2006; Hamrin et al. 2012). In an investigation of the solar wind input to the Earth’s magnetosphere, Lockwood (2019) showed that including the solar wind Poynting flux to the total power input improves the correlation with the am geomagnetic index on timescales of a day and shorter.

It is typically straightforward to extract the data needed to compute the different terms in Poynting’s theorem using numerical models. The power transfer to different ion species in the Martian induced magnetosphere was investigated by Wang et al. (2024) using a hybrid plasma model. These authors showed that both the bow shock and the induced magnetosphere boundary are generator regions, while the planetary ions act as a load and gain energy from the electromagnetic fields. Lindkvist et al. (2018) used a hybrid model and Poynting’s theorem to study the conversion of energy inside the magnetosphere of comet 67P at a heliocentric distance of 1.5 AU for two extreme conditions of solar EUV radiation. They found an asymmetric magnetosphere, where the region of the initial pile-up of the interplanetary magnetic field is a generator and the cometary water ions act as a load, similar to what is seen in the Martian model. In Lindkvist et al. (2018), a single ion solar wind is used, consisting only of protons and electrons. In this paper, we include a second solar wind ion species, alpha particles, in a hybrid model and we investigate the effects a multi-ion SW has on the energy conversion at a comet by simulating two extreme values of the alpha to proton ratio under otherwise constant conditions.

2 Methods

Hybrid models are frequently used to model plasma processes where the kinetic effects of ions are important. The main tool used in this study is the hybrid model code Amitis¹ (Fatemi et al. 2017), which is a GPU-based implementation of the hybrid model that has been successfully used to model numerous objects, such as the plasma environment of Mercury (Fatemi et al. 2018; Aizawa et al. 2021), comets (Gunell et al. 2024), Ganymede inside Jupiter’s magnetosphere (Fatemi et al. 2022). Fatemi et al. (2024) recently presented the first 3D hybrid-kinetic simulations of the solar wind interaction with the entire magnetosphere of Earth at physical scales using Amitis.

In hybrid simulations, ions are modelled as so-called macro particles, where one such macro particle represents a large number of ‘physical’ particles, scaled by the so-called macro particle weight. Position and velocity of each macro particle are computed at every time step using the Lorentz equation of motion (Fatemi et al. 2017). Multiple species of ions (in our case, SW protons, SW alphas, and cometary ions of H₂O⁺ at 18 amu) were simulated. Electrons were implemented as a massless, charge-neutralising fluid, enforcing quasi-neutrality in the entire simulation domain. The simulation domain consists of a cartesian grid of size of 7000 × 8800 × 22 000km (x × y × z), which was split up into cubic cells with a grid resolution of δx = δy = δz = 40 km. The electric and magnetic fields (E and B) were calculated for each grid cell. The x-axis of the simulation coordinate system points towards the sun and the upstream SW velocity, v_SW, is directed along the −x-axis. The interplanetary magnetic field, B_SW, is perpendicular to the SW velocity and is oriented along the y-axis. The convective electric field of the SW (E_SW = −v_SW × B_SW is therefore along the z-axis and splits the simulation domain into a +E (z ≥ 0) and a −E (z ≤ 0) hemisphere.

The ion number and charge density (n_i and ρ_i) and ion current density, J_i, can be calculated from the distribution of macro particles at each time step. The electric current density, J, is calculated via Ampere’s law, ignoring the displacement current (i.e. the Darwin approximation): $$\mathbf{J}=\frac{1}{{\mu }_0}\mathrm{\nabla }\times \mathbf{B}={\mathbf{J}}_i+{\mathbf{J}}_e.$$(2)

This also gives the electron current, J_e, and the electron velocity, v_e. The electric field is given by Ohm’s law, $$\mathbf{E}=-{\mathbf{v}}_i\times \mathbf{B}+\frac{\mathbf{J}}{{\rho }_i}\times \mathbf{B}-\frac{\mathrm{\nabla }{p}_e}{{\rho }_i}+\eta \mathbf{J}=-{\mathbf{v}}_e\times \mathbf{B}-\frac{\mathrm{\nabla }{p}_e}{{\rho }_i}+\eta \mathbf{J}\mathbf{,}$$(3)

with an ion velocity of v_i = J_i/ρ_i, electron pressure of $p_e\propto T_e{n}_i^{\gamma }$, γ = 5/3, T_e = 1.3 × 10⁵K (electron temperature), and plasma resistivity of η = 5.0 × 10³ Ω m. The plasma resistivity helps to dampen numerical oscillations in areas with strong field gradients, while η = 5.0 × 10³ Ωm is a typical value used in hybrid simulations of comets and other planetary bodies. At a weakly outgassing comet, -v_i × B + J × B/ρ_i dominates the electric field E, as expressed in Eq. (3) (Deca et al. 2019). The ions are moved using the non-dissipative electric field E_ND which does not include the resistive term, $${\mathbf{E}}_{\mathrm{N}\mathrm{D}}=\mathbf{E}-\eta \mathbf{J}=-{\mathbf{v}}_i\times \mathbf{B}+\frac{\mathbf{J}}{{\rho }_i}\times \mathbf{B}-\frac{\mathrm{\nabla }{p}_e}{{\rho }_i}$$(4)

Faraday’s law is used to propagate the magnetic field, $$\frac{\mathrm{\partial }\mathbf{B}}{\mathrm{\partial }t}=-\mathrm{\nabla }\times \mathbf{E}.$$(5)

If there are not enough macro-particles in a grid cell (ρ_i → 0), Eq. (5) is approximated by the magnetic diffusion equation, $$\frac{\mathrm{\partial }\mathbf{B}}{\mathrm{\partial }t}=\frac{{\eta }_{\mathrm{v}\mathrm{a}\mathrm{c}}}{{\mu }_0}{\mathrm{\nabla }}^2\mathbf{B}\mathbf{,}$$(6)

where we use the vacuum resistivity, η_vac = 2.5 × 10⁶ Ωm. This also defines the minimum time step for the simulation: Δt ≤ μ₀(Δx)²/ (2η_vac) = 0.4 ms. In our simulations, it is important to include the vacuum resistivity during the startup phase of the simulation. Once the simulation reaches a quasi-steady state, no vacuum regions are present. More information about hybrid models can be found, for instance, in Ledvina et al. (2008). More details on Amitis can be found in Fatemi et al. (2017).

We simulated two different SW compositions and the corresponding simulation runs are subsequently referred to as ‘low-alpha’ and ‘high-alpha’ cases. The major difference between those two runs are the proton and alpha particle densities in the upstream SW. Density values and ratios of the two species are listed in Table 1. The alpha-to-proton ratios were chosen to represent extremes of the SW composition. For a comparable comet-SW interaction the total upstream SW dynamic pressure is conserved between the two simulations. Since both SW species have the same upstream velocity, this also corresponds to a conservation of SW mass density and SW momentum between the two simulations. However, the SW plasma density is lower in the high-alpha case.

Most plasma parameters in our simulations roughly correspond to the plasma environment of comet 67P at a heliocentric distance of 2.5-3 AU. At every time step SW particles are injected into the simulation domain at the upstream boundary at x = +5000 km, with the defined SW parameters following a Maxwell velocity distribution (see Table 1). In both simulations, SW protons and SW alphas have 40 ppc each. The total number of particles (including solar wind and cometary ions) in the simulation is ≈1.8 billion, which are individually tracked at every time step. The cometary neutral density is based on a spherically symmetric outgassing profile with an outgassing rate of Q (Haser 1957) and constant neutral expansion velocity of u_n = 700 ms⁻¹ (Gulkis et al. 2015; Lee et al. 2015). Neutrals are ionised by photoionisation with the photoionisation rate, ν^hν,ioni, and the ion production profile as a function of radial distance, r, is $$P(r)=\frac{Q{\nu }^{h\nu ,\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{i}}}{4\pi u_n{r}^2}$$(7)

with Qν^hν,ioni = 1.08 × 10²⁶ s⁻² (Hansen et al. 2016; Heritier et al. 2018). The cometary ions are produced at every time step based on Eq. (7).

The total electromagnetic power density, E ∙ J, can be directly calculated from the model output variables for each grid point. To calculate the power density of the individual ion species, we first calculated the non-dissipative electric field E_ND used to move the macro particles in the simulation (see Eq. (4)). The power density of a species, n, can then be calculated as E_ND ∙ J_i,n, where J_i,n is the current density of ion species n. The total E ∙ J can be split up into three components: $$\mathbf{E}\cdot \mathbf{J}=({\mathbf{E}}_{\mathrm{N}\mathrm{D}}+\eta \mathbf{J})\cdot ({\mathbf{J}}_i+{\mathbf{J}}_e)={\mathbf{E}}_{\mathrm{N}\mathrm{D}}\cdot {\mathbf{J}}_i+{\mathbf{E}}_{\mathrm{N}\mathrm{D}}\cdot {\mathbf{J}}_e+\eta J^2.$$(8)

The first term on the right-hand side is the summed power density of all ion species and is non-dissipative in the simulation. The second and third term make up the power density of the electron fluid. The second term is the non-dissipative part of the electron power density, while the third term is dissipative and describes processes such as the heating of the electron fluid. We did not analyse the electron power density in this paper; however, we want to emphasise that the presented power densities are not expected to add up to a sum of zero.

3 Results

3.1 Spatial structure

Figures 1 and 2 show a comparison between the low alpha and the high-alpha simulations. The panels in Figure 1 show a slice of the x - z plane, at y = 0 km while the panels in Figure 2 display a slice in the y - z plane at x = −1000 km. The overall shape of the electric and magnetic fields is similar in the two cases. The magnetic field (panels a and b) piles up in front of the comet and drapes around the obstacle consisting of cometary ions. This feature is subsequently referred to as the ‘magnetic field pile-up layer’ and has been observed and modelled in connection with the infant bow shock at comet 67P (Gunell et al. 2018). The maximum magnetic field strength is in the + E-hemisphere right in front of the nucleus and it is approximately the same for both runs (i.e. 6.8B₀, where B₀ is the magnetic field in the upstream solar wind). In the -E-hemisphere (at y = 0), a layer with a steep increase in the magnetic field compared to the upstream value is created. As is seen in Figure 2, in the y - z plane, this layer is a circular structure reaching up into the +E-hemisphere, but with a gap at y = 0 km where the cometary pickup ions dominate. In the high-alpha run, the layer extends about 200 km further upstream. In the low-alpha case, it appears compressed in the +E-hemisphere, thereby forming an oval, rather than a circle. The enhancement of the magnetic field pile-up layer is stronger in the low-alpha case (4B₀), compared to the high-alpha case (3.4 B₀).

The electric field (panels c and d in Figures 1 and 2) is also strongly enhanced in the -E-hemisphere at the location of the magnetic field pile-up layer. The enhancement is stronger in the low-alpha case (3.3E₀) compared to the high-alpha case (2.6E₀). The electric field enhancement is only observed for the lower half (-E hemisphere) of the magnetic field pile-up layer, approximately forming a semi-circle in the y - z plane. The direction of the electric field remains perpendicular to the layer, pointing inward towards the downstream tail. In the downstream tail region behind the nucleus the electric field is significantly decreased compared to the upstream SW. This effect is stronger in the high-alpha case. The downstream + E-hemisphere is characterised by a slight decrease in the electric field to about 0.8 of its upstream value, along with a rotation of the electric field direction towards -x. The flow of electromagnetic energy, represented by the Poynting flux, is shown in Fig. A.1 in the appendix.

Fig. 1 Comparison of the electric and magnetic fields between the low (left column, panels a and c) and high (right column, panels b and d) alpha simulations. All panels show a slice of the x-z plane, at y = 0. Top row (panels a and b): magnetic field, normalised to the upstream value (B0 = 3nT). Bottom row (panels c and d): normalised electric field (E0 = 1.29 mV/m). The background colour shows the magnitude of the electric and magnetic field. The arrows show the strength and direction of the fields in this plane, that is, the x and z components only. The length scales for the arrows are specified above (magnetic field) and below (electric field) the respective colour bars on the right side. The comet is located at (0,0) and marked with a grey ‘x’.

3.2 Energy transfer

To understand the consequences of a change in the SW alpha-to-proton ratio, we focus on the energy transfer from the SW to the cometary plasma. We analysed the energy transfer in three different ways: firstly, the power density for each ion species (E_nd ∙ J), as well as the total power density (E ∙ J), integrated over the entire simulation domain, is summarised in Table 2. Secondly, Figure 3 shows the power density in a plane for two different slices, y = 0 km, and x = −1000 km. Finally, Figure 4 shows the evolution of the energy transfer between the different particle populations and the fields along the x-axis, moving from the upstream SW to the tail region downstream.

Integrated over the entire simulation domain, cometary ions act as a load. The power transfer to the cometary ions (E_ND ∙ J (H₂O⁺)) is comparable for both simulation runs (cf. Table 2). The slightly higher value in the low-alpha case was consistently found for three individual time steps and is expected to be a real (and not a transient) feature and it is not attributed to numerical issues. Similarly, the total E ∙ J is slightly more negative in the low-alpha case. In the case of protons, the power transfer is much larger in the low-alpha case, due to the higher number density of protons in the upstream solar wind. The opposite is true for the alpha particles. The change in electromagnetic energy is small and the total E ∙ J is compensated by the divergence of the Poynting flux (data not shown).

The protons act as a generator in the entire upstream region (see Figure 3 Column 1, rows 1-2). The strongest power density is found in the initial magnetic field pile-up layer (upstreamfacing part). Behind the initial pile-up, the protons gain energy, as seen by the positive power density. This feature is more pronounced in the high-alpha case. In the + E-hemisphere the power density of the protons is mostly negative and appears modulated by wave-like structures. The alphas (column 2) show similar structures as the protons, but the initial generator region is much weaker. Compared to the protons, there are more regions where the alphas act as a load, including large parts of the +E-hemisphere. Upstream of the nucleus, the power density of alphas is close to zero or even slightly positive in the low-alpha case. Comparing the low and high-alpha case (rows 1 and 3 vs rows 2 and 4 in Fig. 3) shows that the structures are similar but the magnitude of the power density of alphas is much smaller in the low-alpha case due to the significantly lower number densities. The cometary ions (Column 3) are a load in the entire simulation domain (E_ND ∙ J > 0). Most of the energy transfer is concentrated near y = 0 in the + E-hemisphere and the tail region. The power density outside the plane at y = 0, as well as in the upstream region, appears in a wave-like pattern. The total E ∙ J (right column, Figure 3) is dominated by SW protons and cometary ions. The effect of alphas is only visible in the high-alpha case, especially in areas where the power density of cometary ions is fairly low (for example at y = ± 1000 km, between z = 0 and z = −800 km, downstream of the magnetic field pile-up layer).

Figure 4 shows two views of the evolution of the power density along x. In the upper panel, we integrate the power density over y and z using the full range of the simulation domain; thereby capturing the energy transfer in the entire simulation box. In the lower panel, only the power density within y = ±200 km is included in the integration and the result is more comparable with the x - z view shown in Figure 3 in the upper two rows. The power density of both solar wind species are shown with the opposite sign (−E_ND ∙ J) for easier comparison with the cometary ions.

Upstream (x ≫ 0), the power density approaches zero with increasing x. There is a net transfer of energy from the SW ions to the cometary ions and the electromagnetic fields. The power density of protons is much higher compared to the alphas; this is true even in the high-alpha case where their mass densities are almost comparable. In the downstream region (x < 0), the energy transfer to the cometary ions as well as the total E ∙ J integrated over the entire y-axis remains almost constant. This is also the case for the protons in the low-alpha simulation run. In this case, the alpha particles act as a load from 200 km upstream of the nucleus to 800 km downstream. In the region further downstream, the alpha population acts as a generator, although its contribution to the total energy transfer is very small compared to the other ion species. For the high-alpha simulation the absolute value of the power density of protons reaches a peak of 400 W/m around x = 0 and decreases in magnitude further downstream. The power density of alphas, on the other hand, starts to increase in magnitude downstream of the nucleus, reaching a value of about −230 W/m at the downstream boundary, comparable to the power density of protons at this point.

The power density close to y = 0 (Figure 4, lower panel) exhibits a clear peak in the absolute value of the power density of protons for both the low- and high-alpha cases (i.e. 53 W/m and 41 W/m, respectively). The peak in the high-alpha case is located approximately 150 km further upstream, consistent with the slightly larger magnetosphere seen in Figure 1). In the downstream region the power density of protons in the low-alpha case is twice that of the high-alpha case. The power density of alpha particles in both cases have a corresponding peak in approximately the same region as the protons. Moving downstream they transition into a load region (power density > 0). Even further downstream, the power density becomes negative again. The power density of cometary ions shows the steepest increase at x = 0 and remains approximately constant after x = −1200 km. The total E ∙ J has a local minimum at the same location as the proton power density, and the values are similar for both simulation runs. Moving further downstream, the total E ∙ J quickly becomes positive, contrary to the results of the full simulation (upper panel). This indicates a transfer of electromagnetic field energy along the y-axis out of the slice around y = ±200 km.

Fig. 2 Magnetic (upper row, panels a and b) and electric (lower row, panels c and d) fields in the y - z plane, for the slice at x = −1000 km. The left side (panels a and c) shows the results of the low-alpha simulation. The right side (panels b and d) shows the results of the high-alpha simulation run. The background shows the magnitude of the electric and magnetic fields, normalised by their upstream values B0 = 3 nT and E0 = 1.29 mV/m. The arrows show the direction of the fields in this plane, that is, the y and z components. The arrow length scales are specified above and below the colour bars for the respective fields.

Fig. 3 Power density, E ∙ J, for the low-alpha (rows 1 & 3) and the high-alpha (rows 2 & 4) cases. The upper two rows show the x - z plane at y = 0 km (same as Figure 1), while the lower two rows show the y - z plane at x = −1000 km (same as Figure 2). Columns 1-3 show the power density of protons, alpha particles, and cometary ions, respectively. The last column shows the total power density E ∙ J. For better visibility, the results of the power density of alphas in the low-alpha case (Column 2, rows 1 and 3) are scaled up by a factor of 20.

Fig. 4 Power density of the individual ion species (protons: red, alphas: blue, and cometary ions: grey) and total E ∙ J (black) along the x-axis. The dotted and dashed lines show the results of the low- and highalpha cases. In the upper panel the power density is integrated along the entire y- and z-axes. For the lower panel, we only integrate between y = ± 200 km along the y-axis. For the SW species the negative values are shown since these species mostly function as generators. The curve for the power density of alpha particles for the low-alpha case (dotted blue line) is scaled by a factor of 20 to improve visibility.

4 Discussion

The bulk properties of the solar wind such as speed, density, and the resulting dynamic pressure are well known to affect its interaction with solar system bodies. For higher pressure the boundaries, such as the bow shock and magnetopause, are found closer to the object (Shue & Song 2002; Edberg et al. 2010). At Earth, a higher solar wind dynamic pressure leads to higher ion escape rates (Schillings et al. 2019). At Mars, a weak inverse relation between the solar wind dynamic pressure and ion escape has been found (Ramstad et al. 2018). In our simulation study, we looked at the effect of varying composition of the solar wind for a fixed dynamic pressure. Our results show that even for a low-activity comet, the composition of the solar wind influences the spatial structure and the mechanisms of energy transfer from the solar wind to cometary ions.

4.1 Spatial structure

The initial formation of the magnetic pile-up layer in the -E hemisphere is mostly driven by the protons. The electron fluid in the model is mass less and charge neutralising so it will immediately follow the ions. The magnetic field in turn is frozen to the electron fluid in our model. In the region dominated by the solar wind ions (i.e. the -E hemisphere), the magnetic field strength distribution will thus be determined by the solar wind ion charge density distribution. Due to different gyroradii, the locations of the density enhancement of protons and alphas are different. In the high-alpha case, the region with enhanced solar wind charge density is broader since the alphas contribute significantly to the total charge density. The relative enhancement with respect to the upstream SW density and, thus, the magnetic field strength in the magnetic field pile-up layer, is therefore lower compared to the low-alpha case (not shown). We may expect a second magnetic pile-up associated with the alpha particles in the high-alpha case. This is not seen in our simulations. The reason appears to be that a magnetic field enhancement has already been caused by the protons and these still dominate in number and charge density, even though the alpha particles carry a similar amount of kinetic energy and momentum flux.

The reason why the magnetic field pile-up and proton density enhancement layer in the -E hemisphere occur further upstream for the high-alpha case can be explained in a similar manner. The amount of mass-loading (cometary ions) is the same as in the low-alpha case, but the mass density of protons is smaller. Thus, there is relatively more mass-loading of the proton population and the boundaries form further upstream. In the high-alpha case, the alpha particles gradually gain importance further downstream from the magnetic field pile-up layer, as discussed below.

4.2 Energy transfer

The energy transfer from the SW to the cometary ions in the upstream region is much more efficient for the protons than the alpha particles. This is a direct consequence of the difference in deflection due to their different gyroradii. As the SW is mass-loaded by the cometary plasma, its individual components react in different ways: the cometary ions are accelerated along the electric field towards + z. To conserve momentum, the SW ions are deflected towards -z, which can also be seen as a rotation or gyration of v_SW around the magnetic field. Upstream of the magnetic field pile-up layer, the electric current density, J, is negligible (cf. Eq. (2)), and the electron velocity, v_e, is approximately the ion velocity, v_i ≈ v_e, where v_i = J_i/ρ_i (J_i ≈ -J_e, ρ_i = -ρ_e). The cometary ions are accelerated along the electric field and the net ion velocity, v_i, and therefore, v_e, differs from the SW flow, v_SW. As a result, the electric field, carried by the electrons (see Eq. (4)), has a component directed against the solar wind velocity. Therefore, the angle between u_SW and the electric field, ∠(υ_SW, E) is > 90°, which also corresponds to a negative E_ND ∙ J of the solar wind ions, namely, a generator region. Equivalently, we can look at the flow direction of solar wind ions with respect to the electrons. A solar wind flow clock-wise of the electron flow indicates a load.

In case of a two-species SW, the SW velocity depends on the velocity and charge density of each individual species. If one SW species dominates the SW mass density, such as the protons in our low-alpha case, the SW velocity will approximately correspond to the bulk velocity of this major species (υ_SW ≈ υ_bulk,major). The deflection of the second, minor species can be approximated as single particle motion in the given electromagnetic fields and depends on its inertia (mass-per-charge m/q). In a typical scenario, the major SW species will be protons, which have the lowest m/q of all ions and, therefore, the highest deflection. The minor species with a larger m/q will be less deflected and ∠(υ_bulk,minor, E) < ∠(υ_buik,major, E). Depending on the m/q of the minor species and the electromagnetic fields, this angle can be above or below 90°, leading to a possible case where part of the solar wind acts as a load (E ∙ J > 0 for the minor species). This can primarily be seen during initial deflection, before particle trajectories get too complex, and it is stronger for ions with larger m/q. As the mass density of the minor species increases, υ_SW and the related electric field direction are more influenced by this species, while ∠(υ_bulk,minor, E) will drop below 90°, ensuring that this species acts as a generator. Our results (see details in Figs. 3 and 4) show that this effect is at work up to 500 km upstream of the nucleus; in the low-alpha case, the power density of alphas is weakly positive in some regions, while in the high-alpha case, the power density is negative, but much smaller in magnitude compared to the protons, despite having 40 % of the proton’s charge density (and therefore current).

Figure 5 illustrates this effect using v_e instead of E. This makes it easier to visualise generator and load regions based on the flow direction of ions and electrons. A region where the ions are rotated clockwise with respect to the electrons is a load region for this ion species. We identified four regions determined by the sign of the power density of protons and alphas, defined as areas ‘1’ through ‘4’ in the insets to the left (based on Figure 3). The illustrated flow lines are calculated from the simulated bulk velocities of protons, alphas, and electrons (note: these flow lines do not necessarily correspond to particle trajectories in the case of the ions if the velocity distributions are not Maxwellian). The upstream case described in the previous paragraph corresponds to area ‘1’ (arrows in the top left). The blue arrows corresponding to the alphas may rotate slower (low alpha, dash-dotted line) or faster (high alpha, solid line) than v_e. In our high-alpha simulation, the alpha bulk velocity is almost aligned with v_e, leading to a negative, but very small (in magnitude) power density. For an efficient energy transfer of alphas to the cometary ions in this region an even higher alpha-to-proton ratio would be required.

Close to the pile-up layer (area ‘2’) the flow direction of the protons (and, to a lesser extent, the alphas) changes more suddenly. The E-field and electron velocity also make a sudden change due to the currents (see J arrows in Fig. 5). The net result is a clockwise rotation of v_e, also seen as a small bump towards +z in the electron flow line of the high-alpha case in Fig. 5. As a consequence, both alphas and protons act as generators in this area. As soon as the current density decreases in magnitude and changes sign downstream of the magnetic field’s pile-up layer, the electric field rotates counter-clockwise again, resulting in E_ND ∙ J > 0 for the alphas (area ‘3’). The structure of this load-generator region near the pile-up layer appears to be similar for both high- and low-alpha ratio cases.

Further downstream (area ‘4’ in Fig. 5), the electric field changes less and the alphas gyrate further, resulting in a negative power density, or a generator region. The deflection of protons decreases as their gyration progresses and they form a load region downstream of the pile-up layer. This feature is seen in both cases and has been previously reported in Lindkvist et al. (2018). They simulate a comet with higher outgassing rate and only include one SW species (protons). Their results show several transitions between generator and load regions in the power density of protons. These striations in the SW proton power density are attributed to the secondary acceleration of protons due to a change in current direction. Compared with our low-alpha case (Figure 3, rows 1 and 3), there are more generator-load transitions in the protons in Lindkvist et al. (2018) due to the broader pile-up layer. The transitions between negative and positive power density, especially the first load region in the protons, are more distinct in our high-alpha case. This is likely due to the two SW species mass-loading each other and affecting the ion dynamics of the other species.

Every cometary ion that is created will eventually be picked up and will attain solar wind speed. Hence, the total transfer of solar wind energy to the cometary ions must be the same regardless of the solar wind alpha-to-proton ratio as long as the newborn ion profiles are identical. Since there can be spatial and temporal variations in the simulations, the power density must be integrated over a sufficiently large volume and averaged in time to cancel out non-negligible fluctuations. Far downstream of the nucleus, the cometary ion production is negligible. The number flux of cometary ions is then conserved, that is, constant through all slices x = x_s of the simulation box, for x_s that is sufficiently far downstream. From Figs. 1c and d, we note that the electric field in the region dominated by cometary ions (downstream of the nucleus and behind the pile-up layer) is not varying much along x. If we assume the field to remain constant along the flow lines of the cometary ions in this region (within the domain covered by the simulation box), the power density of the cometary ions will also be constant if the flux is conserved. This is indeed what we see in the upper panel of Figure 4: downstream of x ≈ −1000 km, E_ND ∙ J for the cometary ions is almost constant for both runs.

The energy transfer from the SW to the cometary ions is affected by the amount of alphas present. In the upper panel of Fig. 4, a clear difference is noted between the high and low alpha runs, especially downstream of the nucleus (x < 0). For the high-alpha run, the decrease in the integrated power density (∫_YZ E_ND ∙ J) of the protons is compensated by the increase in the power density of the alphas. The two solar wind species are coupled to each other and their interaction modifies the spatial distribution of the power density.

In the high-alpha case (where the alpha particles constitute 20% of the total number density) the kinetic energy available in the SW alpha particles is 80% of the energy available in the protons (as the mass of an alpha particle is four times the mass of a proton and the velocities are equal). Still the power density of the solar wind species only reach the same value at the left edge of the simulation box. Within this box, the power transfer from the alphas to the cometary ions and fields is much less than that from the protons (cf. Table 2). This is likely caused by the delayed mass-loading of the alpha particles.

In the spatial domain covered in the simulations, the distance from the pile-up layer to the sun-comet line increases almost linearly with a decreasing x (cf. Figure 1). Provided the magnetic field strength and the thickness of the pile-up layer does not change much the magnetic energy density contained in the layer must increase linearly with -x, as the radius of the almost circular layer increases (cf. Figure 2). This means that the power transfer to the fields must be constant, which is what we see for x < 0 in Figure 4 (upper panel). Close to the nose of the cometary magnetosphere, where there is a significant curvature of the pile-up layer, this will not be true and, in that case, we observe a peak in the magnitude of the power density.

With both the total E ∙ J and power density of the cometary ions being approximately constant downstream of the nucleus, it follows that the sum of the power densities for protons and alpha particles must be constant. In the high-alpha case, it means that when the power density of the alpha population is increasing in this region and there is a corresponding decrease in the power density carried by the protons.

It should also be noted that the magnitude of the total power density is considerably smaller than the power density of the particle populations, with exception of the alpha particles in the low-alpha case. This suggests that mass-loading is the dominant process of energy transfer between SW and cometary ions.

Fig. 5 Illustration of proton, alpha, and electron bulk velocities and derived flow lines in the x - z plane at y = 0 km. In the upper row, arrows show the plasma velocities in different areas (marked 1-4), corresponding to different combinations of generator and load regions for alphas and protons. The deflection angles are exaggerated for illustration purposes. The areas are defined in the two insets to the right, which show the END ∙ J for protons and alphas in the high-alpha case (cf. Fig. 3, row 2). The lower right part of the figure shows the flow lines for the low- and high-alpha case, with the magnitude of the magnetic field (cf. Fig. 1, panels a and b) in greyscale in the background. There, the relative position of the flow lines with respect to the magnetic field pile-up layer can be seen. Solid lines: high-alpha case; dash-dotted lines: low-alpha case. Red, blue, and grey indicate protons, alphas, and electrons, respectively. The grey ‘x’ marker indicates the position of the comet nucleus. The flow line plots for the low-and high-alpha case are separated for improved visibility.

5 Conclusions

In this paper, we used a hybrid model to investigate the effect of different SW alpha-to-proton ratios on the solar wind’s interaction with a comet. We carried out two simulation runs; one with the ratio between the number densities of alphas and protons being 1% and another with a ratio of 20%. The number densities were adjusted so that the dynamic pressure and momentum of the SW were the same in the two cases.

By comparing the two runs, we can conclude that a higher SW alpha-to-proton ratio leads to a larger cometary magnetosphere. The initial pile-up of the magnetic field is governed by protons as well in the high-alpha case. With fewer SW protons available, the relative mass loading of the SW protons is larger, pushing the pile-up layer further away from the nucleus. The electromagnetic fields in the pile-up layer are weaker in this case, as the fraction of the total magnetic field frozen into the electron fluid related to the SW protons is smaller.

The interaction between the cometary atmosphere and the solar wind leads to a transfer of energy from the SW species to the cometary ions. The transfer of energy from the alpha particles to the cometary ions occurs further downstream than for the protons due to the larger inertia of the alpha particles. As a result, there are regions downstream where the power density of alphas is comparable to the power density of protons. The two SW species take turns in acting as loads and generators. In well-defined spatial regions, the sign of the power density is different for protons and the alpha particles. The SW species gain energy in some regions as a result of the alpha and proton populations mass loading each other.

If the SW is comprised of multiple ion species, mass loading and energy transfer to the cometary ions will be most efficient for the species with lowest m/q, typically protons. This is due to the different deflection of the individual species, caused by the different gyroradii. In some situations, it is possible that the SW species with higher inertia has a positive power density, namely, it acts as a load in some regions. When analysing the impacts of CMEs on comets or other unmagnetised bodies, the effects of a potentially increased alpha-to-proton ratio onto the induced magnetospheres shown in this paper should be taken into account.

Acknowledgements

The simulations were enabled by resources provided by the EuroHPC Regular Access at Vega, grant EHPC-REG-2023R03-023, and by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at High Performance Computing Center North (HPC2N) partially funded by the Swedish Research Council through grant agreement no. 2022-06725. Work at the Swedish Institute of Space Physics in Kiruna (IRF) was funded by the Swedish National Space Agency (SNSA) grant 132/19. Work at Umeå university was supported by SNSA grant 2023-00208 (HG) and SNSA grant 2022-00183 (SF).

Appendix A Poynting flux

Figure A.1 shows the Poynting flux of the high-alpha simulations. Overall, the x-component of the Poynting flux dominates and is consistently negative, that is, globally, the Poynting flux points downstream (anti-sunward). The y-component points towards y = 0 km in both hemispheres due to the draping of the magnetic field around the comet. As the magnetic field bends around the comet, there is an inwardly directed flux (towards y = 0 km) of electromagnetic energy. The increasingly negative S_z when moving downstream is due to the increased deflection of the solar wind towards -z.

Fig. A.1 Poynting flux of the high-alpha simulations. The three columns show the x, y, and z component of the Poynting flux (from left to right). The top row is a slice in the x - z plane at y = 0 km, and the bottom row is a slice in the y - z plane at x = −1000 km.

References

All Tables

All Figures

Fig. 1 Comparison of the electric and magnetic fields between the low (left column, panels a and c) and high (right column, panels b and d) alpha simulations. All panels show a slice of the x-z plane, at y = 0. Top row (panels a and b): magnetic field, normalised to the upstream value (B0 = 3nT). Bottom row (panels c and d): normalised electric field (E0 = 1.29 mV/m). The background colour shows the magnitude of the electric and magnetic field. The arrows show the strength and direction of the fields in this plane, that is, the x and z components only. The length scales for the arrows are specified above (magnetic field) and below (electric field) the respective colour bars on the right side. The comet is located at (0,0) and marked with a grey ‘x’.

In the text

Fig. 2 Magnetic (upper row, panels a and b) and electric (lower row, panels c and d) fields in the y - z plane, for the slice at x = −1000 km. The left side (panels a and c) shows the results of the low-alpha simulation. The right side (panels b and d) shows the results of the high-alpha simulation run. The background shows the magnitude of the electric and magnetic fields, normalised by their upstream values B0 = 3 nT and E0 = 1.29 mV/m. The arrows show the direction of the fields in this plane, that is, the y and z components. The arrow length scales are specified above and below the colour bars for the respective fields.

In the text

Fig. 3 Power density, E ∙ J, for the low-alpha (rows 1 & 3) and the high-alpha (rows 2 & 4) cases. The upper two rows show the x - z plane at y = 0 km (same as Figure 1), while the lower two rows show the y - z plane at x = −1000 km (same as Figure 2). Columns 1-3 show the power density of protons, alpha particles, and cometary ions, respectively. The last column shows the total power density E ∙ J. For better visibility, the results of the power density of alphas in the low-alpha case (Column 2, rows 1 and 3) are scaled up by a factor of 20.

In the text

Fig. 4 Power density of the individual ion species (protons: red, alphas: blue, and cometary ions: grey) and total E ∙ J (black) along the x-axis. The dotted and dashed lines show the results of the low- and highalpha cases. In the upper panel the power density is integrated along the entire y- and z-axes. For the lower panel, we only integrate between y = ± 200 km along the y-axis. For the SW species the negative values are shown since these species mostly function as generators. The curve for the power density of alpha particles for the low-alpha case (dotted blue line) is scaled by a factor of 20 to improve visibility.

In the text

Fig. 5 Illustration of proton, alpha, and electron bulk velocities and derived flow lines in the x - z plane at y = 0 km. In the upper row, arrows show the plasma velocities in different areas (marked 1-4), corresponding to different combinations of generator and load regions for alphas and protons. The deflection angles are exaggerated for illustration purposes. The areas are defined in the two insets to the right, which show the END ∙ J for protons and alphas in the high-alpha case (cf. Fig. 3, row 2). The lower right part of the figure shows the flow lines for the low- and high-alpha case, with the magnitude of the magnetic field (cf. Fig. 1, panels a and b) in greyscale in the background. There, the relative position of the flow lines with respect to the magnetic field pile-up layer can be seen. Solid lines: high-alpha case; dash-dotted lines: low-alpha case. Red, blue, and grey indicate protons, alphas, and electrons, respectively. The grey ‘x’ marker indicates the position of the comet nucleus. The flow line plots for the low-and high-alpha case are separated for improved visibility.

In the text

	Fig. A.1 Poynting flux of the high-alpha simulations. The three columns show the x, y, and z component of the Poynting flux (from left to right). The top row is a slice in the x - z plane at y = 0 km, and the bottom row is a slice in the y - z plane at x = −1000 km.
In the text