Oxygen Ion Escape From Venus Is Modulated by Ultra‐Low Frequency Waves

We study the solar wind‐driven, nonthermal escape of O+ ions from Venus in a global hybrid simulation. In the model, a well‐developed ion foreshock forms ahead of the Venusian quasi‐parallel bow shock under nominal upstream conditions. Large‐scale magnetosonic ultra‐low frequency (ULF) waves at 20‐ to 30‐s period are excited and convect downstream along the foreshock with the solar wind. We show that the foreshock ULF waves transmit through the bow shock in the downstream region and interact with the planetary ion acceleration, causing 25% peak‐to‐peak fluctuations in the O+ escape rate. These results demonstrate the importance of upstream plasma waves on the energization and escape of heavy ions from the planetary atmospheres.


Introduction
Our sister planet Venus is extremely dry as compared to the Earth and has likely lost a significant amount of water during the history of the solar system (Greenwood et al., 2018). It is still under debate how the water was lost and how much different atmospheric erosion processes have changed the planet's volatile inventories. Being an unmagnetized body, the upper atmosphere of Venus is subject to the direct, noncollisional solar wind-driven escape of ionized heavy elements, which are gravitationally bound to the atmosphere (Futaana et al., 2017). In the Venus-solar wind interaction, parts of the ionized ionospheric and exospheric particle populations are accelerated to the escape velocity and are lost to space. At Earth the solar wind influence on the atmospheric erosion is mediated by the geomagnetic field (Nilsson et al., 2012;Yau et al., 1985). As the present-day heavy element loss rates from Venus are not very significant on the time scales of planetary evolution, it is essential to quantify all mechanisms for the atmospheric escape to produce a reliable estimate of the volatile erosion history of the planet (Persson et al., 2018). Such results will also be relevant for Mars as well as for any other unmagnetized body in the solar or exoplanetary systems.
The induced magnetosphere of Venus forms when the interplanetary magnetic field (IMF) piles up against the ionosphere creating the magnetic barrier. Similar to other planets, a bow shock forms the boundary between the supermagnetosonic solar wind and the heated and turbulent magnetosheath plasma enveloping the induced magnetospheric boundary or the magnetopause (Russell et al., 1988;Zhang et al., 2008aZhang et al., , 2008b. The foreshock region is magnetically connected to the bow shock, which allows backstreaming of the suprathermal charged particles (Eastwood et al., 2005;Omidi et al., 2017). The interaction between the suprathermal and incident solar wind populations drive a multitude of plasma waves, especially the large-scale ultra-low frequency (ULF) magnetosonic waves, which are typically found under quasi-parallel conditions where the angle between the shock surface normal and the IMF is smaller than about 45 • (Fränz et al., 2017;Keiling et al., 2016). Such waves can have significant effects on the dynamic processes in the planetary space environments especially at unmagnetized bodies, where the bow shock forms close to the planet (Lundin, 2011).
Acceleration mechanisms of ions away from unmagnetized planets include several processes or "escape channels," which can be classified according to the ion energy and region around the planet Dubinin et al., 2011). Significant cold or low energy ion escape is found especially in the induced magnetotail as well as at low altitudes near the ionosphere. At higher altitudes the largest-scale escape channel  Slavin and Holzer (1981) Parameter Symbol Value is the ion pickup, but also, several smaller-scale channels exist. Relative escape rates through different channels can vary as a function of the upstream conditions and observed ion energy (Dong et al., 2017;Fedorov et al., 2011).

Number of macroions
In the high-energy escape, the pickup ions form the heavy ion plume (Jarvinen et al., 2016;Liemohn et al., 2014;Nordström et al., 2013;Wei et al., 2017). In addition to the so-called "north-south" asymmetry, the plume exhibits a significant hemispheric "dawn-dusk" asymmetry in the direction perpendicular to the plane defined by the undisturbed solar wind velocity and convection electric field vectors, when the upstream IMF has a strong flow-aligned component as is typical at the orbit of Venus (Jarvinen et al., 2013;McComas et al., 1986f). Interestingly, the dawn-dusk asymmetry means that the E×B drift turns the escaping heavy ions toward the hemisphere of the quasi-parallel bow shock and the ion foreshock rather than the opposite hemisphere (Jarvinen et al., 2013;Jarvinen & Kallio, 2014). The acceleration of planetary ions by convecting magnetic field fluctuations in the Venus' magnetosheath downstream from the quasi-parellel bow shock, like the foreshock ULF waves, has been suggested in test particle studies (Luhmann et al., 1987) and recently discussed based on observations ( . The three-dimensional field line tracing was started in the upstream region at z = 2,000 km. Small gray spheres give the location of Points P1-P3. Big gray-black sphere has the radius of Venus for context. See Movie S1 for temporal evolution of the parameters. Different aspects of the solar wind-driven ion escape from unmagnetized planets have been studied in self-consistent plasma models including hybrid and magnetohydrodynamic codes (Ledvina et al., 2008).
While several studies have focused on the escape rates and the structure of the induced magnetosphere (e.g., Brain et al., 2010, and references therein), the interaction of the ULF waves and the ion escape has not been analyzed in a self-consistent model. Here we report on a new finding that the foreshock ULF waves have significant effects on the Venusian heavy ion acceleration in the induced magnetosphere and escape.

Model
We simulate the Venus-solar wind interaction using the hybrid model platform RHybrid (Jarvinen et al., 2018(Jarvinen et al., , 2020). In the model, ions of solar wind and planetary origin are treated as macroscopic particle clouds (macroparticles), and their motion is determined by the Lorentz force. Electrons are an isothermal, charge-neutralizing and mass-less fluid. Planetary ions are produced via photoionization of hydrogen and oxygen exosphere coronae and via an upward emission of ionospheric oxygen ions from the model inner boundary. The production rates and profiles of planetary ions are the same as in our earlier Venus works and correspond to solar minimum conditions (Jarvinen et al., 2009(Jarvinen et al., , 2013. The solar wind is injected through the front wall, and macroparticles are removed from the simulation as they encounter simulation boundaries. The simulation run uses nominal, stationary upstream conditions at Venus (Slavin & Holzer, 1981) with the Parker spiral angle of 36 • and, thus, the flow-aligned component of the IMF stronger than the perpendicular component. The simulation setup and the algorithm are similar to our earlier studies of the Venus and Mars space environments (Jarvinen et al., 2013;Kallio et al., 2010), with the exception that the number of grid cells and macroparticles is much higher in the current parallel code compared to the sequential code. See Table 1 for details of the simulation run and Kallio and Janhunen (2003) for further details of the algorithm. We use a planet-centered coordinate system, where the x axis is antiparallel to the incident, undisturbed solar wind flow; the y axis is aligned along the perpendicular IMF component to the undisturbed solar wind flow; and the z axis completes the right-handed coordinate system and, thus, is along the convection electric field in the undisturbed solar wind flow. The hemisphere where the upstream solar wind convection electric field points away from the planet (z > 0) is termed the +E sw hemisphere, and the y < 0 hemisphere is termed the foreshock hemisphere. The radius of Venus (R V = 6,051.8 km) is used as the unit of length in the figures and the text.
Temporal properties of the solar wind and planetary plasma and fields were recorded at every time step between t = 250 … 450 s in grid cells centered at the points P1 (x, y, z) = (0.56, −4.24, 0.01)R V , P2 (x, y, z) = (−2.19, −1.19, 0.01)R V , and P3 (x, y, z) = (0.01, −0.96, 0.86)R V . Figure 1 shows the locations of the points in the simulation domain. Figure 1 shows a snapshot of the Venus induced magnetosphere and the O + bulk flux (see Movie S1 in the supporting information for the dynamics of these parameters). The induced magnetosphere is clearly visible in the B z component with the bow shock lineating the outermost boundary and the magnetosheath separating the induced magnetotail from the upstream solar wind. The foreshock upstream of the bow shock on the y < 0 hemisphere includes large-scale waves visible in B z . Both parameters in the figure and the movie show ongoing variations at different temporal and spatial scales even though the solar wind driving is stationary.  The O + escape rate is analyzed in detail in Figure 2. The O + net fluxes (outward flux-inward flux) integrated over spherical shells at different altitudes show the radial evolution of the escape rate dynamics. The lowest altitudes show little fluctuations in the O + escape, but the fluctuations intensify with increasing altitude. The escape rate through the outer boundaries of the simulation domain shows fluctuations with about 25% peak-to-peak amplitude, with maximum power spectral density at the frequency of 0.03-0.04 Hz (25-33 s) ( Figure S4). Spectral maxima at about the same frequency range can be identified at the spherical shells r ≥ 1.5R V .

Results
The net escape rate increases to about 94% of the value at the domain outer boundary from the r = 1.1R V to 1.5R V shell. This low-altitude increase is caused by a drop in the planetward O + rate and photoion production between the two shells. The planetward O + rate is 3 orders of magnitude smaller than the net escape rate at the r = 2.7R V shell. The O + escape rates through the outer boundaries are 2.9 × 10 24 s −1 for the ionospheric population and 1.9 × 10 24 s −1 for the exospheric photoions. Figure 3 displays the magnetic field time series in the foreshock (P1); in the quasi-parallel equatorial, night-side magnetosheath (P2); and in the low-altitude quasi-parallel terminator region on the +E sw hemisphere (P3). Periodic large-scale waves are evident at the three points with the maximum power spectral density in the ULF frequency range of about 0.03-0.05 Hz (20-33 s) ( Figure S5). In this study, we refer to these waves as the 20-to 30-s waves. At P1, the average B z is nearly zero as expected due to the upstream conditions, whereas at P2 the average B z is slightly negative and at P3 positive because of the IMF piling up and draping around Venus.
The electron density and the magnitude of the magnetic field are positively correlated for the foreshock waves. Minimum variance analysis (MVA) shows that the foreshock waves are left-hand polarized and travel at a small angle (<10 • ) with respect to the magnetic field in the simulation frame. An estimated wave phase speed is below the solar wind bulk velocity projected in the direction of the wave propagation, which implies  that the foreshock waves are propagating upstream and are right-hand polarized in the plasma frame (see details of the MVA and phase speed determination in the model in Jarvinen et al., 2020). Taken together, these imply that the foreshock ULF waves are oblique fast magnetosonic modes excited by the backstreaming solar wind ions in the foreshock. This in agreement with in situ spacecraft observations by Pioneer Venus Orbiter (Luhmann et al., 1983) and Venus Express (Shan et al., 2016) as well as with previous global hybrid models (Omidi et al., 2017).
The O + energization is analyzed in Figure 4. The parameter q ⃗ E · ⃗ U(O + ) gives the average net work by the electric field on an O + ion in each grid cell per unit time. In the foreshock region (P1), the power varies from negative to positive, implying temporal changes from bulk deceleration to acceleration. This is due to a low O + density and statistical variations of the velocity of exospheric photoions sometimes aligned and sometimes anti-aligned with the electric field in the upstream region. In the nightside magnetosheath at P2, the O + flux is also dominated by the exospheric population, but the density is higher, and, thus, the escaping O + flow is more organized along the tail and the electric field than at P1. The O + density is highest at the quasi-parallel terminator (P3) and dominated by the ionospheric population. The average power is positive, indicating net acceleration at both downstream locations. Points P1-P3 show modulation of q ⃗ E · ⃗ U(O + ) with the maximum power spectral density in the ULF frequency range of about 0.02-0.05 Hz (20-50 s) ( Figure S6).

Discussion
Using a global hybrid simulation, we show that the upstream ULF waves interact with the O + ion acceleration and escape from Venus. In the model, under nominal, stationary solar wind conditions, large-scale 20to 30-s magnetosonic foreshock ULF waves are excited in the ion foreshock, and they convect downstream with the solar wind flow and transmit through the quasi-parallel bow shock (Dubinin & Fraenz, 2016;Shan et al., 2014). The waves interact with the O + energization in the upstream, near-equatorial magnetosheath and low-altitude terminator regions on the foreshock hemisphere (Figure 4).
The coupling between the O + acceleration and the ULF waves at the frequency range of the foreshock ULF waves is evident in the upstream region (x = −1.5 … 0.5 R V , y = −3.0 … − 2.0R V ), and in the quasi-parallel magnetosheath (x = −3.0 … −1.5R V , y = −2.0 … − 1.0R V ) (Figure 1 and Movie S2). Furthermore, the coupling of the O + energization and B z is present already at the low-altitude region on the spherical shell at r = 1.29R V where the vantage point P3 is located (Figures 2-4, Movie S03 in the supplementary material).

10.1029/2020GL087462
The ULF waves and the O + modulation at the frequency range of the foreshock ULF waves are clearly visible on the quasi-parallel side of the shell (longitude = 240 … 360 • and latitude = -10 … 60 • ) from t = 100 s onward (Movie S3).
The dynamics of the escaping O + ions in the model can be shown to be consistent with theoretical consideration of an idealized pickup process (Jarvinen & Kallio, 2014): A scatter-free motion of a pickup ion starting at rest in homogeneous electric and magnetic fields includes periods with the z component of the velocity aligned (acceleration) and anti-aligned (deceleration) with the electric field. A time evolution of the energization for an ideal pickup ion is q ⃗ E · ⃗ U(O + ) = qE z V E×B sin(Ω c t), where q is the particle electric charge, E z is the convection electric field, V E×B is the E × B drift velocity, c is the angular gyrofrequency, and t is time. In our case, the upstream conditions give q ⃗ E · ⃗ U(O + ) = 639 eV∕s sin( c t) for a pickup ion, showing that the ideal energization varies from -639 to 639 eV/s, compatible with the values in Figure 4. Even though the energization rates of hundreds of eVs per second are high, planetary heavy ions are observed at tens of keV energies, and such acceleration is available by the electric fields embedded in induced magnetospheres (Futaana et al., 2017;Jarvinen et al., 2018).
It is also important to notice that the ULF waves are not resonant with the gyromotion of the O + ions: In the upstream region, the O + gyroperiod is 104 s, which is well above the foreshock ULF wave period of 20-30 s. An O + ion reaches gyroperiods of ≤ 30 s only when the magnetic field strength is ≥ 34.7 nT, and such strong magnetic fields are limited to the low-altitude dayside magnetic barrier region under nominal upstream conditions at Venus. However, lighter species are more likely to become gyroresonant with the ULF waves (Shimazu et al., 1996).
In order to isolate the effect of the ion foreshock on the O + energization, we performed a test run with a purely perpendicular upstream IMF relative to the solar wind flow. Under perpendicular IMF conditions, the bow shock is quasi-perpendicular throughout the simulation domain, and, consequently, no ion foreshock nor foreshock ULF waves form, and there are no dawn-dusk asymmetries (Jarvinen et al., 2013). In the test run, the O + escape rate does not show the significant fluctuations found in the Parker IMF case; the escape fluctuations had a peak-to-peak amplitudes less than 10% and periods less than 10 s. These weak fluctuations may be associated with turbulence or mirror mode waves in the magnetosheath (Volwerk et al., 2016) or arise as a result of statistical macroparticle noise in the model. Conversely, the test run demonstrates that the statistical macroparticle noise is not the source of the fluctuations in the Parker spiral case.
As the wave modulation of the ion escape may occur also at Mars (Kallio et al., 2006), we will focus future studies on the role of the magnetosheath wave activity and the dynamics of the induced magnetospheres on the heavy ion energization and escape at both Venus and Mars (Dimmock et al., 2018;Futaana et al., 2017;Girazian et al., 2019). Furthermore, foreshock waves upstream of the quasi-parallel bow shock can couple with the low-altitude proton fluxes via the energetic neutral atom (ENA) production at Mars (Fowler et al., 2019), but it is still an open question how the charge exchange, electron impact ionization, and ionospheric photochemistry processes are affected by the ULF waves (Mazelle et al., 2018;Yamauchi et al., 2015); what are the ULF wave properties and their effect on the cold ion escape in or near the ionosphere and further in the tail (Dubinin & Fraenz, 2016;Omidi et al., 2020); or how the ULF modulation of the escape and energization of heavy ions works under different upstream conditions, including flow-aligned IMF cases when the ion pickup is not a significant source of planetary ion acceleration (Luhmann et al., 1993). Resolving the contribution of the ULF waves on the ion escape for different upstream conditions and ionization processes allows us to assess the evolutionary significance of its contribution on the atmospheric erosion at unmagnetized planets.

Conclusions
We analyze the Venus-solar wind interaction using a global hybrid simulation, which demonstrates a strong modulation of the O + ion energization and escape by the foreshock ULF waves. Consistent with the pickup of planetary ions by the solar wind convective electric field (rather than gyroresonance), the O + energization is modulated by the ULF waves in the upstream, magnetosheath, and low-altitude regions leading to the 25% peak-to-peak fluctuations in the global escape rate from the simulation domain. This mechanism is sufficiently effective that it needs to be accounted for in the interpretation of heavy ion observations and possible acceleration of planetary ions by plasma waves at Venus and Mars.