Annales Geophysicae Oxygen ion escape from Venus in a global hybrid simulation : role of the ionospheric O + ions

We study the solar wind induced oxygen ion escape from Venus’ upper atmosphere and the Venus Express observations of the Venus-solar wind interaction by the HYB-Venus hybrid simulation code. We compare the simulation to the magnetic field and ion observations during an orbit of nominal upstream conditions. Further, we study the response of the induced magnetosphere to the emission of planetary ions. The hybrid simulation is found to be able to reproduce the main observed regions of the Venusian plasma environment: the bow shock (both perpendicular and parallel regions), the magnetic barrier, the central tail current sheet, the magnetic tail lobes, the magnetosheath and the planetary wake. The simulation is found to best fit the observations when the planetary O + escape rate is in the range from 3×1024 s−1 to 1.5×1025 s−1. This range was also found to be a limit for a test particle-like behaviour of the planetary ions: the higher escape rates manifest themselves in a different global configuration of the Venusian induced magnetosphere.


Introduction
The unmagnetized Venus planet has a solar wind interaction very different from the Earth.In contrast to the geodipole, the Venusian obstacle to the magnetized solar wind is the highly conducting ionosphere.The field lines of the interplanetary magnetic field (IMF) pile up against the ionosphere Correspondence to: R. Jarvinen (riku.jarvinen@fmi.fi) and drape around the planet.In the dayside the piled field lines increase the field magnitude creating a magnetic barrier.The barrier stops the flow and a bow shock is formed upstream of the planetary plasma environment.Between the bow shock and the region dominated by the planetary ions is the magnetosheath.At the terminator where the field lines slip from the dayside to the nightside the draping results in "a magnetic pole structure" which is connected to the central (or cross) tail current sheet (CTCS).The CTCS separates the two magnetic lobes in the wake of the planet.This kind of a magnetic environment is called an induced magnetosphere (see, for example, Saunders and Russell, 1986).
The solar wind flows about Venus at altitudes close to the top ionosphere and the high exospheric neutral densities.The exospheric neutrals are constantly ionized by the solar EUV radiation, the electron impact ionization and the charge exchange processes.Consequently, the high altitude planetary ions are accelerated by the electric and magnetic fields related to the solar wind flow and the Venus' induced magnetosphere (see, for example, Luhmann et al., 2006, and references therein).These acceleration processes result in a non-thermal escape of the planetary atmosphere.Atomic oxygen (O + ) is the dominant ion species in the Venusian upper atmosphere.Also, atomic hydrogen ion (H + ) densities are considerable at altitudes above 200 km (see Fig. 23 in Brace and Kliore, 1991).Estimates of the total solar wind induced oxygen ion escape from Venus range from the order of 10 24 s −1 to 10 26 s −1 including observational and simulation studies (e.g., Brace et al., 1987;Moore and McComas, 1992;Terada et al., 2004;Lammer et al., 2006;Kallio et al., 2006b;Barabash et al., 2007a;Martinecz et al., 2009).We select to call the value 10 25 s −1 as the nominal O + escape rate.
The Venus-solar wind interaction was first studied in detail by the Pioneer Venus Orbiter (PVO).A recent review of the results is given by Russell et al. (2006).PVO measured, R. Jarvinen et al.: O + escape from Venus among other things, the escaping oxygen ions from Venus.Moreover, PVO observed the changing nature of the Venussolar wind interaction as a function of the solar cycle.The changes are related to the EUV flux and the solar wind conditions.An increase in the EUV radiation modifies the upper atmosphere and enhances the ionization of the exospheric neutral densities.This means that more ions are produced and are available for the solar wind induced escape.The increase of the ion density changes the obstacle to the solar wind and, consequently, can affect to the global Venus-solar wind interaction.For example, the bow shock terminator altitude correlates with the solar cycle being lower for the solar minimum (Russell et al., 1988).
Venus Express (VEX) is an European spacecraft in orbit around Venus.VEX studies the Venusian plasma environment by its particle instrument ASPERA-4 (Barabash et al., 2007b) and the MAG magnetometer (Zhang et al., 2006).Barabash et al. (2007a) reported observations of the O + escape in the near Venus' tail.The statistically presented measurements show the O + , H + and He + flux away from the planet and energization in the direction of the convection electric field typical to the E × B pickup.
Depending on the IMF and the convection electric field, the maximum bulk kinetic energy of a Venusian pickup oxygen ion is about one order of magnitude higher than the energy of a solar wind proton (m O + /m p ≈ 16).The gyroradius of an oxygen ion in the solar wind at Venus is typically about 1-2 planetary radii.
The Venus' solar wind interaction has been examined earlier in several self-consistent simulation studies including MHD and hybrid codes.Brecht and Ferrante (1991) were the first ones to present a global hybrid simulation study of the Venus-solar wind interaction.Bauske et al. (1998) and Kallio et al. (1998) studied the mass loading of the Venusian ions into the solar wind in MHD.Terada et al. (2004) considered in their hybrid code different Venusian ion escape processes.Martinecz et al. (2009) compared their hybrid simulation to the VEX/ASPERA-4 and VEX/MAG observations in an orbit case study concentrating on the boundaries.
In this work we compare our Venus hybrid simulation (HYB-Venus) to the VEX particle and magnetic observations and study how the Venus' plasma environment and the oxygen ion escape change with a wide range of different planetary O + emission rates.We choose here to study only changes related to the emission rates of thermal oxygen ions (O + th ) from the ionosphere.The studied rates are from 10 23 s −1 to 10 27 s −1 .The exospheric sources are kept constant even though in reality a stronger ionosphere would mean different exospheric profiles.However, when we introduce higher ion densities around the planet it may be that their exact origin is not that important.Johansson et al. (2009) considered in a hybrid code an expanding ionospheres of unmagnetized exoplanets by increasing the ion emission velocity upwards from the planet; their study resembles somewhat our extreme parameter cases.
The paper is organized as follows.We begin by describing the simulation and the used Venus Express observations (Sect.2).Next, we compare the simulation to the observations in an orbit case study (Sect.3).In Sect. 4 we present a parameter study in which we change the planetary oxygen emission rates.At the end we discuss the results and summarize our findings (Sects.5 and 6).

The HYB code and the Venus Express observations
HYB is a hybrid code for the planetary plasma interactions developed at the Finnish Meteorological Institute.The code has been used to study the plasma interactions of Mars (e.g., Kallio et al., 2006a), Mercury (e.g., Kallio and Janhunen, 2003), Titan (e.g., Sillanpää et al., 2007), the Moon (e.g., Kallio, 2005) and Venus (e.g., Jarvinen et al., 2008).See Kallio and Janhunen (2003) for the numerical details of the code.Further, Kallio et al. (2006b) introduce the usage of the code in the Venus simulations.
In the hybrid model the ions are treated as particles in the Lorentz force field and the electrons are a massless, chargeneutralizing MHD fluid.The electric and magnetic fields are self-consistently related to the ions and the electron fluid via Maxwell's equations.Recently, we have added to the code the polarization electric field caused by the electron pressure.Next we shortly summarize how the model equations are solved and describe the simulation inner boundary and the planetary ion sources.

Equations
Starting with initial values for the ions (position and velocities) and the magnetic field in the simulation domain, the system is propagated as follows.First, the current density is obtained from the magnetic field by the Amperè's law and the electron charge density comes from the quasineutrality assumption in a grid cell where the summation is over all ion species.The velocity field of the electron fluid is determined as Now all the quantities needed to calculate the electric field from the electron momentum equation are known: where η a is a given (3-D) function for the resistivity and the electron pressure comes from the ideal gas law p e = n e k B T e , where T e is a given constant (isothermal electrons).The resistivity term is included only in the field propagation (Eq.5) to introduce some field diffusion.Further, the electron pressure contributes only to the particle propagation (Eq.6) because the curl of a gradient in Faraday's law is identically zero in the isothermal case.Finally, Faraday's law is used to advance the magnetic field and the particles are moved and accelerated by the Lorentz force: The equations are solved numerically by the leapfrog method and the Boris-Buneman integrator for the Lorentz force (e.g.Birdsall and Langdon, 1991).A Yee lattice type gridding is used to store the magnetic field on grid cell faces to obey the ∇ • B = 0 condition (e.g.Dai and Woodward, 1998).
Table 1 summarizes the main parameters of the simulation runs used in this work.The simulations were run 1100 s in physical time until all the runs had reached a quasistationary solution.The solutions with higher emission rates stabilized later than the nominal and lower emission runs, which got stable at about 400 s.All the visualizations shown are snapshots at t = 1100 s.

Inner boundary and planetary ions
The planetary ionosphere, the obstacle to the solar wind flow, is parametrized as an inner boundary condition in the HYB-Venus simulation.The inner boundary is a spherical surface 300 km above the planet's surface.At the inner boundary the electron velocity and pressure is set to zero.Further, the resistivity is also set to zero at the boundary which means that there is no diffusion of the magnetic field into the obstacle and, therefore, the magnetic field is zero inside the inner boundary.The particles are absorbed (removed from the simulation) at 200 km above the planet's surface, i.e. 100 km below the inner boundary.
Planetary ions are injected into the simulation domain using two different types of sources, the exospheric and ionospheric populations.The exospheric ions are created from neutral corona profiles which are photoionized by constant No other ionization processes are considered explicitly.The applied EUV ionization rates correspond to the solar minimum conditions.See Table 1 for the total production rates and Table 2 for the details of the ionospheric population in the runs.The exospheric ion sources are described in more detail in Kallio et al. (2006b).The O + is dominated by the exospheric photoions in the 10 23 and 10 24 cases.The Venusian ionospheric O + th population is injected into the simulation from a spherical surface above the inner boundary.The emission profile has a cosine dependence on the solar zenith angle (SZA) decreasing towards the terminator on the dayside and a constant nightside (10 per cent of the emission at SZA=0).
The total O + th emission rate is a free parameter since we do not consider a self-consistent ionospheric physics in the code.A self-consistent ionosphere would include chemical processes responsible for creating the neutral and ionized upper atmosphere starting from given (lower) atmospheric profiles and ionization processes.The implementation of the self-consistent ionosphere is complicated by the achievable spatial resolution in the global planetary hybrid simulations.See Table 1 for the radius and the emission rates of the ionospheric population.

Coordinate systems
The Venus Solar Orbital (VSO) system is a planet-centered coordinate system with the x-axis directed towards the Sun, the y-axis pointed opposite to the planet's orbital motion in the orbital plane and the z-axis northward, perpendicular to the orbital plane.
The Venus Solar Electric (VSE) coordinate system is also a planet-centered system with the x-axis directed towards the Sun.The VSE y-axis is directed along the perpendicular component of the IMF to the x-axis and, consequently, the −V SW × B IMF electric field is along the z-axis.The VSE coordinates are related to the VSO coordinates via a rotation around the x-axis.
The simulation coordinate system is the same as the VSO or VSE system with an exception that the negative x-axis Ann.Geophys., 27,[4333][4334][4335][4336][4337][4338][4339][4340][4341][4342][4343][4344][4345][4346][4347][4348]2009 www.ann-geophys.net/27/4333/2009/ is directed along the incident solar wind flow.The discrepancy between the simulation and the VSO/VSE coordinates is caused by the planetary motion around the Sun.The angle between the x-axes of these coordinate systems is about 5 degrees.This causes an aberration of the IMF components less than 1 nT which is not expected to be significant in the current study.
The used value of the planetary radius is R V = 6051.8km.

Venus Express orbit 20 November 2006
Venus Express is observing the Venusian plasma environment in a 24-h elliptical polar orbit.The orbit periapsis is above the geographic north pole of the planet.The orbit on 20 November 2006 was selected for this study from the MAG observations between 2006-2007 since it has a very stable and non-fluctuating IMF close to nominal Venusian values (see the IMF in Fig. 2 and the solar wind in Fig. 3).The stable conditions favor the use of a stationary simulation solution since no dynamics are driven by the upstream changes.Also, the steadiness makes it possible to use a well-defined VSE coordinate system, since the orientation of the convection electric field related to the IMF is almost constant.
The inbound part of the orbit has an IMF x-component of 8.9 nT and a perpendicular (yz) component of 8.3 nT, which is a little higher than the nominal value.In the VSO coordinates the IMF clock angle (the direction of the perpendicular component) is about 160 degrees from the positive y-direction towards the negative z-and y-axis.The IMF spiral angle (the direction of xy-component) is about 41 degrees from the positive x-direction towards the negative y-direction in the VSO coordinates.Thus, Venus is located in the IMF towards sector of the Parker's spiral (see Fig. 4 in Jarvinen et al., 2008).In the VSE coordinates the IMF spiral angle is about 43 degrees from the positive x-axis towards the positive y-axis.
The solar wind consists of H + ions with no distinguishable helium content (see panel 3a) and the dynamic pressure is somewhat weaker than the nominal Venus value.All the conditions are within the 75 per cent confidence intervals of the PVO solar wind and IMF statistics, i.e. the conditions are not in the 25 per cent of the extreme cases of the statistics (Table 2 in Jarvinen et al., 2008).The solar cycle and the EUV flux were close to a minimum on November 2006.The upstream conditions are listed in Table 1.

Simulation results: an orbit case study
In this section we present an orbit case study where we compare the Venus Express magnetic and particle observations made during a single orbit (20 November 2006) to the HYB-Venus simulation.The comparisons provide an extrapolated view at the plasma regions the spacecraft sampled around Venus.We use the run with the nominal O + th emission rate for the comparison.

Magnetic field
Figures 1 and 2 display a comparison of the observed magnetic field and the simulation in the VSE coordinates.Figure 1 shows a 3-D structure of the field magnitude in the simulation and the spacecraft's orbit.Figure 2 compares the components and magnitude of the magnetic field between the VEX/MAG data and the HYB-Venus simulation.Together Figs. 1 and 2 illustrate the regions the spacecraft sampled along its trajectory.
First, VEX crossed the Venusian bow shock in the perpendicular shock region (y < 0 hemisphere) at 06:10 UT close to the z = 0 plane.The MAG data shows a sharp and large gradient, typical to the perpendicular shock wave, and a peak magnitude of 41 nT at the bow shock, whereas the simulation produces a smoother transition from the IMF into the magnetic barrier.The gradient is 500 km long (foot-to-peak) along the spacecraft trajectory in the observations and about 9000 km in the simulation.
After the inbound bow shock crossing the spacecraft moves into the magnetic barrier region where the dynamic pressure of the solar wind flow is converted to the magnetic pressure and the thermal plasma pressure (06:10-06:25 UT).The orbit passes the barrier in the side of a positive B x polarity (see B x > 0 in panel 2a) touching only the edge of the region with the highest field magnitudes (see panel 1c).The observed field magnitude in the barrier is 30-35 nT and this is also compatible with the simulation.
At 06:29 UT, just before the periapsis of the orbit, the spacecraft crosses a region, a current sheet, where B x switches polarity from about 20 nT to −10 nT close to the terminator.In the observations the x-component is zero at the SZA of 88 degrees.The current sheet crossing takes place in the region where the central tail current sheet starts to form from the induced magnetic poles on top of the ionosphere caused by the field line draping (see CTCS in panel 1b).The observed width of the current sheet crossing (along the orbit from a positive B x peak to a negative B x peak) is 800 km whereas in the simulation the width is 16 000 km. Since the current sheet is crossed near the inner boundary of the simulation (see the green parts of the simulation lines in Fig. 2), the timing of the crossing is not exactly correct.In the simulation the B x = 0 point is on the nightside of the planet at the SZA of 128 degrees (06:40 UT).The simulation shows a strong curvature of the magnetic field and the slipping of the field lines from the dayside of the planet to the nightside in this area (figure not shown).
Next, the spacecraft moves into the nightside of the planet (SZA>90 deg) and into a region between the CTCS and a magnetic tail lobe.Here the spacecraft samples the lowest magnitude of the magnetic field for the whole orbit (5 nT) and turbulent wave activity can be seen (06:33-06:43 UT).At about 06:46 UT the spacecraft enters the right magnetic tail lobe with negative B x polarity, i.e.B x is pointing along the tail away from the Sun, and a magnitude of about 10 nT.The observations show a decrease in the field fluctuations inside the lobe (note that there are short data gaps just before and after the lobe).The simulation reproduces the magnitude and the location of the lobe.The maximum SZA (142 degrees) of the orbit is inside the lobe.
The connection of the magnetic field lines changes when moving through the lobe (figure not shown).Before entering the lobe the curvature of the magnetic field close to the current sheet is caused by the magnetic pole structure in the z < 0 hemisphere near the terminator.Inside the lobe the connection changes gradually to the dayside of the planet (near the center of the lobe the field lines trace back into the dayside noon sector) and finally, when leaving the lobe, to the pole region in the z > 0 hemisphere.At about 07:00 UT the spacecraft leaves the lobe and enters the magnetosheath and after that the bow shock.The outbound bow shock crossing takes place in the turbulent parallel shock region (07:30-07:40 UT) near z = 2R V .The jump in the field magnitude at the bow shock is only weak and the data shows a lot of fluctuations.The simulation reproduces the parallel shock wave including its magnitude and location.

Ions
Figure 3 displays a comparison between the observed ASPERA-4 ion energy spectra (for H + and O + separately) and the simulation bulk parameters (bulk flow energy 1 2 mnV 2 , density and temperature).Further, Fig. 4 illustrates the modelled H + and O + densities and bulk energies in the spacecraft's orbital plane (cf. the magnetic field in panel 1c).The peaked modulation in the spectrograms is caused by the electrostatic elevation scan of the ASPERA-4/IMA sensor.Also, note that the instrument collects ions only from a restricted solid angle and with a certain energy threshold (no observations below 10 eV per unit charge) whereas the simulation plots include ions from the total solid angle of 4π and all the energies.Further, there exists some leakage between the H + and O + channels of the instrument that can be seen as, for example, incorrect O + counts in the upstream plasma.When we refer to a certain observed ion species in the spectrograms in this study we have identified the species from the IMA energy-mass spectrum for the discussed time interval (figures not shown).
The bow shock regions are again consistent between the observations and the simulation (see panels a, b and c in Fig. 3).Both show a sharp and clear inbound crossing of the perpendicular shock wave.Decrease of the bulk energy is of the same order of magnitude in the data and in the simulation, approximately from 1 keV (solar wind) to 200 eV (magnetosheath).Further, the temperature of the shocked magnetosheath plasma is about 1.5 orders of magnitude higher than in the solar wind in the simulation, which is approximately consistent with the observed H + energy spectrum (see spreading of the energy distribution at 06:10 UT in panel 3a).The outbound crossing of the bow shock displays a different kind of behaviour.Since the crossing happens in the parallel region of the shock, the bulk energy increases steadily from the magnetosheath into the solar wind showing no abrupt jump in the observations or in the simulation.Further, the temperature decreases also steadily from the magnetosheath to the upstream value.The modelled H + density increases by approximately a factor of two at both bow shock crossings from 14 to about 30 particles per cm 3 .
A low density planetary wake (particle density < 1 cm −3 in the simulation) is seen between the inbound and outbound magnetosheath separated by the proton dropout boundaries (PDB).The PDBs are crossed at about 06:25 and 06:52 UT when the H + counts vanish abruptly and the planetary heavy planetary ions (O + ) start to appear in the spectrograms.In the simulation the oxygen dominates the total ion number density from 06:30 to 06:47 UT (see panel 3d).
The observed O + counts are concentrated in several peaks in the wake (panel 3b).The first O + peak occurs just before the current sheet crossing at 06:28 UT at the altitude of about 500 km.The peak includes an intense (ionospheric) population of low energy counts below and around 10 eV plus a weaker high energy population from about 100 eV to above 1 keV.Accelerated O + and H + is seen in the further peaks also with the energies of several hundreds of electron volts.Since the simulation does not include a realistic ionosphere, the oxygen energies even at low altitudes are always above 100 eV.After the current sheet crossing the spacecraft reaches the orbit's periapsis altitude of 330 km at 06:31 UT.
In the wake the density is dominated by O + (see panel 3d) and the bulk energy of the solar wind H + shows heavy fluctuations (see panel 3c).Also, the planetary ions show heavy fluctuations in energy during the inbound half of the orbit (<06:30 UT).The heavy fluctuations are caused by the particle noise in the low density areas in the simulation.
In the simulation the spacecraft's trajectory intersects the pickup O + region during the outbound half of the orbit (>06:30 UT).The pick-up O + is seen with energies exceeding 10 keV, i.e. up to an order of magnitude higher than the solar wind, (see Figs. 3c, 4b and 4d).However, the pickup O + is not seen in the observed ion spectra.

Simulation results: a parameter study
This section presents a parameter study of the Venus-solar wind interaction expanding the case study in the previous section.Here we explore the response of the Venusian induced magnetosphere to the emission rate of the thermal ionospheric oxygen (O + th ).We have performed four extra runs in addition to the nominal run used in the orbit case study.All runs are identical expect that the O + th emission rate is varied from 10 23 to 10 27 s −1 .The used rates are listed in Table 2 with additional information.In the following we show visualizations of the magnetic field and plasma quantities on the y = 0 plane from three of the runs (10 24 , 10 25 and 10 26 runs).Also, the projection of the VEX orbit is shown.

Magnetic field
Figure 5 displays the magnitude of the magnetic field in the selected three runs.The field magnitude along the orbit in all the five runs is shown and compared to the VEX/MAG observations in Fig. 6.
It is evident from Fig. 5a and b that the magnetic fields in the nominal and lower emission runs are almost identical.The same behaviour is seen in Fig. 6 where the runs 10 23 , 10 24 and 10 25 overlap almost identically.Hence, the configuration of induced magnetosphere is not changed if the O + th emission is decreased by a factor of 100 from the nominal case.Further, we made a run without any planetary ion populations in the simulation and the solution still remained the same as in the nominal emission case.Thus, the planetary ions behave essentially as test particles in the global system (in the sense that their existence does not manifest itself in the global magnetic field structure) when the O + th emission rate is 10 25 s −1 or below.
In the high O + th emission runs (10 26 and 10 27 ) the magnetic configuration of the Venus' system changes when compared to the nominal run.In Fig. 5c the field magnitude is shown in the 10 26 emission case and in Fig. 6 the magnitudes are shown along the orbit also for the 10 26 and 10 27 cases.The increased emission rate creates an expanded magnetosphere around the planet.The bow shock is pushed farther away than in the nominal case.
In the comparison plot (Fig. 6) it is seen that the emission case 10 27 is clearly not consistent with the MAG data.The bow shock crossings (inbound and outbound) and the tail lobe are displaced from the observed locations (see the discussion in the previous section).The 10 26 case has magnetic field closer to the nominal run.The inbound bow shock crossing is displaced but the tail lobe and the outbound bow shock take place near to the locations in the nominal run and in the observations.This implies that the planetary ions start to affect the global structure of the Venusian magnetosphere when the O + th emission rate is between 10 25 −10 26 s −1 .

Ions
In Fig. 3d the O + densities are shown along the orbit for the three cases.The density in the 10 24 run is about 5 times smaller than in the nominal case and the oxygen dominates only very close to the periapsis.In the 10 26 case the oxygen starts to dominate already before the inbound PDB crossing and displays the densities of the order of 10 cm −3 .Figures 7 and 8 show the densities and energies of the H + and O + particles in the selected three simulation runs.Again, the nominal and low emission runs show similar solutions which is a result from the above-mentioned test particle-like behaviour of the planetary ions.The small differences in the H + plots, especially in the wake, are mostly caused by the particle noise.Also, the decreased O + th emission produces naturally lower O + densities (compare panels 8a and 8b).
The rightmost panels in Figs.7 and 8 display the H + and O + properties in the run with the second highest O + th emission.The magnetic field of the solution differs from the nominal case as already seen in Figs. 5 and 6.The increased O + th emission creates an expanded magnetosphere and a high density O + bulk in the wake of the planet.The region of low H + flow energy in the wake is significantly larger than in the nominal run.Also, the increased emission generates a high density and low energy O + buffer around the planet at low altitudes which can be seen in Fig. 8c and f.This buffer diamagnetically shields the planet from the convection electric field and also less solar wind H + impacts the inner boundary.

O + escape
The O + escape is studied in Figs. 9 and 10. Figure 9 shows the morphology of the O + flow in the selected three runs.The vectors are plotted only in the regions with high O + bulk flux (above 10 10 m −2 s −1 ).Most of the escape in the low and nominal runs happens in the pickup ion region in the z > 0 hemisphere where the convection electric field points away from the planet (left and middle panels in Fig. 9).The O + s are energized in the direction of the convection electric field (i.e. the bigger the z-coordinate, the larger the energy).
In the higher emission run most of the escape takes place in the low energy O + bulk in the wake of the planet (rightmost panels in Fig. 9).The O + energization in the bulk is directed along the x-axis.This is caused by the fact that the electric field is small inside the bulk and the thermal motion carries the ions along the tail (the O + th ions are injected with an upward velocity from the inner boundary).The O + s gain bulk energy in the z-direction when they slowly move along the tail.Some high energy pickup ions in the z > 0 hemisphere are also seen in the high emission case.However, the O + buffer shields the lowest altitudes and does not allow the O + acceleration to occur near the regions with highest densities.The resulting pickup flux is weaker than in the nominal case.The dawn-dusk asymmetry of the O + ions on the xy-plane is caused by the non-zero IMF x-component (see middle panels in Fig. 9).The x-component gives rise to the E × B drift in the y-direction.
Figure 10 shows the O + escape rates and the related rate of the kinetic energy escape (the total O + escape power) from the simulation in all the runs.The escape rates and powers are estimated by counting the absorbed O + particles at all outer boundaries of the simulation domain.In the high emission runs it took more time than in the nominal run for the escape rates to stabilize.This was caused by the increasing low energy bulk of O + in the wake.To verify the stationarity of the escape all the rates in Fig. 10 are determined from runs with longer running times (5000 s) and lower resolution than the runs visualized in this study.The runs are otherwise identical.The convergence to same global solutions between the two sets of runs was found to be good.
In the 10 23 and 10 24 cases the O + escape is dominated by the exospheric photoions which are constant sources in all the runs.In the nominal and higher emission runs the escape is dominated by the ionospheric O + th population.The O + escape rate is almost linearly dependent on the emission rate.About 15 to 30 per cent of the injected O + escapes from the simulation in the runs, whereas most of the particles impact onto the inner boundary (if an injected particle does not escape from the box it must impact the inner boundary in a stationary situation).Also the O + escape power from the system depends almost linearly on the emission rate (panel 10b).In the nominal and low emission cases the average kinetic energy per escaping particle is about 17 keV (panel 10c) which corresponds approximately to the solar wind speed for an O + : 1 2 m O + (430 km s −1 ) 2 ≈ 15.3 keV.Thus, most of the O + escape occurs as high energy pickup ions which was already illustrated in Fig. 9.In the 10 26 and 10 27 runs the average energy decreases because of the reduced pickup ion flux and more O + escape through the low energy and high density region (bulk) in the wake of the planet.The ions are energized in the bulk less efficiently than in the lower emission cases because of the decreased electric field.

Discussion
In the orbit case study we ran our hybrid code to make an interpretation of the VEX magnetic and plasma observations.The run gives a 3-D extrapolation of the Venus' induced magnetosphere and plasma environment for the VEX orbit on 20 November 2006 (see Figs. 1,2,3 and 4).In the comparison between the observations and the simulation we identified the perpendicular and parallel bow shocks, the magne-tosheath, the magnetic barrier, the central tail current sheet (or the magnetic pole) and the magnetic lobe with negative B x polarity.In addition, the comparison of the ion quantities display the hot magnetosheath H + flow and the planetary ion environment dominated by O + between the proton dropout boundaries.All the mentioned features are consistent between the simulation and the observations considering the resolution of the simulation.
At the current sheet crossing (see Fig. 3b) the O + spectrogram shows, in addition to a thermal <10 eV population, a ∼ 1 keV oxygen population which is not present in the simulation (Figs.3c and 4b).Similar planetary O + at the solar wind energies is seen in the Martian CTCS and is attributed to the J × B "slingshot" forces (Dubinin et al., 1993).Interestingly, in the simulation we see highly curved magnetic field lines and an intense re-energization of the solar wind H + in this region.A similar re-energization by the slingshot force is discussed in the MHD simulation study of Venus and Mars by Tanaka (1993).During the outbound half of the orbit we see in the simulation a clear population of pickup O + with energies exceeding 10 keV (see Fig. 3c).This is a well collimated beam directed away from the Sun as can be seen in Fig. 9b, e and h.The bulk number flux is more than 10 10 m −2 s −1 in regions where the vectors are plotted in Fig. 9.By comparison, the solar wind flux is 6 × 10 12 m −2 s −1 .However, for an unknown reason no high energy pickup O + counts are seen in the spectrogram in Fig. 3b.This may be a consequence of a limited field of view of the instrument, the attitude of the spacecraft or "smoothing" of the oxygen beam by temporal variations.See th emission rate is a free input parameter in the HYB-Venus simulation because we do not consider currently a self-consistent ionospheric physics.The most realistic emis-sion rate was found to be in the range from 10 25 to 10 26 s −1 for the discussed VEX orbit based on the comparison between the simulation runs and the MAG and ASPERA-4 observations.With the O + th emission rates below 10 25 s −1 the planetary ions do not affect the magnetic structure of the Venus' environment.Also, the oxygen density in the wake drops to 0.1 cm −3 and O + does not dominate the ion composition in the whole planetary wake (Fig. 3d).On the other hand, the emission rate 10 26 s −1 changes the magnetic structure from the nominal case and produces O + densities of the order of 10 cm −3 and the ion composition is changed from H + to O + even before the inbound PDB in the observations.
The O + th emission rates 10 25 s −1 and 10 26 s −1 produce the total O + escape rates of 3 × 10 24 s −1 and 1.5 × 10 25 s −1 , respectively (see Fig. 10).Less than half of the injected oxygen ions escape from the box in any of the runs.Thus, the escape seems not to be limited by the availability of the ions in the system.

The emitted O +
th ions from the ionospheric boundary can be assumed to have a Maxwellian velocity distribution.Then the emission rate F is connected to the density n and the temperature T via the equation for the random particle flux upwards from the ionospheric boundary (Bittencourt, 2004, p. 180): where R is the emission radius, m the particle mass and k B the Boltzmann's constant.Equation ( 8) is valid in thermal equilibrium where the upward and downward particle fluxes are equal.The corresponding emission rates and thermal pressures are listed in Table 2.After the injection the particles are moving in the self-consistent Lorentz force field and the population evolves to a non-Maxwellian velocity distribution.Since less than 50 per cent of the injected O + particles escape, more O + actually impacts the inner boundary than escapes from the simulation box.
The runs with different emission rates show that the planetary particles behave approximately as test particles in the global Venusian induced magnetosphere when the O + th injection rate is 10 25 s −1 or below (see Figs. 5, 6 and 7).With higher rates the planetary particles affect the global configuration of the magnetosphere (see Figs. 5c and 6).This change cannot be explained by the initial thermal pressure of the injected O + th population (see pressures in Table 2).The highest listed pressure is 0.15 nPa which still an order of magnitude lower that the upstream dynamic pressure in the simulation (m p n SW V 2 SW ≈ 4.3 nPa).More likely, the change from the test particle-like behaviour of the planetary ions to the situation where the planetary ions affect the global configuration of the magnetic field happens because inside the magnetosphere the electric and magnetic fields cannot accelerate all the injected particles into the solar wind velocity (about 15 keV for oxygen) as seen in Fig. 10c.Introducing more ions causes higher densities and lower energies (=low electron velocity), a buffer, at low altitudes and in the tail.The low electron velocity means that these regions are diamagnetically shielded from the efficient acceleration processes.
The planetary ions carry away energy from the Venusian system when they are accelerated by the induced solar wind electric field.The average kinetic energy over a gyroperiod of an O + particle accelerated in the solar wind by the E × B pickup for our IMF case is 2|E × B|/B 2 ≈ 14 keV which is less than the found 17 keV average energy per particle in the pickup dominated runs (runs 10 23 −10 25 in Fig. 10c).It is possible that this is caused by more efficient acceleration processes within the Venus' induced magnetosphere than the E × B pickup (cf. the J × B slingshot acceleration of the solar wind H + ) or is an artefact from the counting of the escaping O + particles at the fixed box boundaries (the gyrophases of the particles may not be randomly distributed at the walls).The angle between the E × B drift velocity and the x-axis comes for our IMF case tan −1 (B y /B x ) ≈ 47 degrees which is seen as the dawn-dusk asymmetry of the O + flow in Fig. 9.The asymmetry decreases in the higher emission runs because the pickup is less efficient and the resulting high energy pickup O + flux is smaller (see the lowest panels in Fig. 9).
In the simulation the average energization factor of a planetary particle, i.e. the average kinetic energy of a particle absorbed at the box boundary divided by the average kinetic energy per injected particle, is over several thousands for all the populations.This means that the planetary source of the initial kinetic energy is negligible compared to the total escape power of the energized ions from the box.Practically all the kinetic energy comes from the acceleration caused by the self-consistent Lorentz force in the system.
The total O + escape power of the accelerated planetary particles exceeds 10 10 J s −1 in the runs with the O + th rate 10 26 s −1 and above (see Fig. 10b).On the other hand, the IMF Poynting flux convecting with the solar wind flow into the box through the front boundary is ≈ 5.6 × 10 10 J s −1 .This is the total injection power of the convective magnetic energy from the IMF into the simulation.Furthermore, this is about the same value that is incident into the Venusian bow shock, since the bow shock cross section in the tail is almost the same as the perpendicular box area.Thus, it is interesting to note that this value is of the order of the escape power of the planetary particles when the planetary ions begin to affect the global system.One may ask if the configuration changes because there is no more magnetic energy to accelerate the planetary particles into the solar wind velocity?
That said, there is also bulk flow energy available in the solar wind.The total injection power of the solar wind bulk kinetic energy into the simulation is 1 2 This is two orders of magnitude higher than the Poynting energy flux into the box.Still, the bulk flow energy can be converted to the pickup ions only via the electromagnetic interactions which is a non-trivial matter inside the induced magnetosphere.This conversion and the related energization of the Venusian planetary ions should be a topic for future studies.

Summary
We used the HYB-Venus hybrid code to study the Venus Express ion and magnetic observations of the Venus-solar wind interaction during an orbit of nominal upstream conditions.
In the orbit case study the perpendicular and parallel bow shocks, the magnetic barrier, the central tail current sheet, the magnetic tail lobes, the magnetosheath and the planetary wake were found to be consistent between the observations and the simulation.In the parameter study the simulation was found to best fit the observations with the escape rate of oxygen ions ranging from 3 × 10 24 s −1 to 1.5 × 10 25 s −1 .Further, this range was found to be a limit for a test particlelike behaviour of the planetary ions: the higher rates affect the configuration of the Venusian induced magnetosphere in global scale.

Fig. 1 .
Fig. 1.A 3-D illustration of the magnitude of the magnetic field in the nominal HYB-Venus simulation run for the Venus Express orbit on 20 November 2006.The black line is the orbit and the red dot denotes the inbound phase of the orbit.Panel (a) shows the field magnitude at the y = 0, z = 0 and x = −2R V planes.In panel (b) the z = 0 plane and the spacecraft's approximate orbital plane are shown.Panel (c) shows a normal view to the spacecraft's orbital plane.

Fig. 2 .
Fig. 2.A comparison of the magnetic field observed by VEX/MAG (blue line) and the magnetic field from the nominal HYB-Venus simulation run (red line) along the orbit shown in Fig.1.The lowest panel includes the spacecraft's altitude from the planetary surface (dashed black line) and solar-zenith angle (dashed magenta line) information.The green line in the simulation plot denotes the part of the orbit below 500 km (less than a simulation grid cell) altitude.

Fig. 3 .
Fig. 3.A comparison of the H + and O + time-energy spectra observed by VEX/ASPERA-4 (two uppermost panels) and the bulk flow energy (E = 1 2 mnV 2 ) and density (n) of the H + and O + ions and the total plasma temperature (T ) from the nominal HYB-Venus simulation run (three lowest panels) along the orbit shown in Fig.1.Coloring in the spectral plots denotes particle counts per second.In the simulation plots two different ion populations are shown: the solar wind H + (blue line) and the planetary O + (magenta).In addition, dashed magenta lines in the density panel are the O + densities in the lower (10 24 ) and higher (10 26 ) emission runs.

Fig. 4 .
Fig. 4. Normal views to the spacecraft's orbital plane for the bulk flow energy (upper panels) and density (lower panels) of the solar wind H + (left panels) and the planetary O + (right panels) in the nominal HYB-Venus simulation run.See the orientation of the orbital plane in Fig. 1.

Fig. 5 .Fig. 6 .
Fig. 5.The magnitude of the magnetic field at the y = 0 plane in three HYB-Venus simulation runs with different O + th emission rates.The emission rates are (from the left): 10 24 , 10 25 (nominal) and 10 26 s −1 .The spacecraft's trajectory is projected into the plane.

Fig. 7 .
Fig. 7. Density (upper panels) and bulk flow energy (lower panels) of the solar wind H + at the y = 0 plane in three HYB-Venus simulation runs with different O + th emission rates.The emission rates are (from the left): 10 24 , 10 25 (nominal) and 10 26 s −1 .

Fig. 8 .
Fig. 8. Density (upper panels) and bulk flow energy (lower panels) of the planetary O + at the y = 0 plane in three HYB-Venus simulation runs with different O + th emission rates.The emission rates are (from the left): 10 24 , 10 25 (nominal) and 10 26 s −1 .
Fig. 2a and b in Luhmann et al. (2008) for an example of spectra including the 10 keV O + counts observed by ASPERA-4/IMA.The parameter study explores the Venus' plasma environment with different planetary O + th emission rates.The O +

Fig. 9 .
Fig. 9. Morphology and energy of the O + flow in three HYB-Venus simulation runs with different O + th emission rates.The emission rates are (from the left): 10 24 , 10 25 (nominal) and 10 26 s −1 .Background coloring shows the O + density at the y = 0, z = 0 and x = 0 planes.Topmost planes show a view towards the planet from the positive y-axis, middle row panels from the +z-axis and bottom panels from the +x-axis.Vectors show the O + velocity in regions with high bulk flux of O + (nv > 10 10 m −2 s −1 ).The vectors are normalized to the unit length and the coloring gives the velocity.

Fig. 10 .
Fig. 10.Escape of the planetary O + particles in five HYB-Venus simulation runs with different O + th emission rates.The x-axis is the O + th emission rate.The top panel shows the absolute O + escape rate from the simulation box, the middle panel shows the rate of the O + kinetic energy escape from the box (the total O + escape power) and the bottom panel shows the average kinetic energy per escaping particle.

Table 1 .
Details of the HYB-Venus simulation runs and the used upstream conditions for the Venus Express orbit on 20 November 2006.

Table 2 .
Ionospheric O + th emission rates and thermal pressures from Eq. (8).The emission radius is 400 km and the temperature 2000 K.The O + in the simulation is dominated by the exospheric photoions in the 10 23 and 10 24 cases.