An Explanation of the Nightside Ionospheric Structure of Venus

Hybrid simulations of Venus have produced nightside structure. It is believed these are the first simulations to produce the structure reported by spacecraft such as the Pioneer Venus Orbiter (PVO). The simulations reveal that the structure is created by the ▽Pe ambipolar electric field when it aligns with local IMF (interplanetary magnetic field) direction. Thus, the orientation of the IMF in and near the ionosphere controls where on the nightside the tail structure will occur. The structure size and density is found to agree well with that measured by the PVO spacecraft. Further, comparisons with ion energy and velocity distributions also agree well with the summary conclusions from a variety of PVO papers and with some nightside observations by Venus Express. The agreement between the data and simulations lends confidence that the simulation results have merit.

whistler wave spectrum are present in the code. Collisionless shock formation physics is included in this code. Therefore, no assumptions or numerical techniques are needed to capture the shock.
The inner boundary conditions are absorbing with ions entering the planetary region (<100 km) being removed from the calculation. The magnetic field is loaded within the planet as a constant so that no numerical currents are produced on the boundary. We chose to use the IMF (interplanetary magnetic field) for this load. We do not assume that the surface is a conducting boundary.
The HALFSHEL code contains a highly resolved chemistry (including charge-exchange and electron impact ionization), a detailed collisional module, and a very general electric field equation. The ionospheric chemical reactions used in HALFSHEL (Brecht & Ledvina, 2010) are identical to the set used by Ma et al. (2004). The Venus ionization frequencies come from Ma et al. (2013). The chemistry is performed on a high-resolution spherical grid surrounding the planet with radial resolution of 5 km and 0.44° angular resolution (49.8 km at 450 km altitude). Each particle has a density associated with it and the density is modified by the chemistry for all particles within the spherical grid.
The collisional module addresses species-by-species collisions using the same high-resolution grid because both the collisions and the chemistry require that the neutral atmosphere scale height be resolved. The HALFSHEL code treats the ion drag by individually colliding each ion with the local CO 2 and O neutrals using their respective collision cross-sections. The neutral atmosphere included within the code provides the neutral densities and velocities. The energy and momentum loss from each particle is tracked. The neutral profiles are unmodified by the energy and momentum transfer to the neutrals in the HALFSHEL simulations.
The neutral atmosphere complete with neutral winds used in HALFSHEL was provided by Dr. Stephen Bougher of the University of Michigan. The model atmosphere was computed using the VTGCM (Venus Thermosphere Global Circulation Model) a code developed by Dr. Bougher (cf. Bougher et al., 1997. The VTGCM atmosphere used in the HALFSHEL simulation was created when the sun was at solar minimum. The photochemistry module within HALFSHEL creates the ionosphere. A Chapman layer model is included to account for ionization fall off from optical depth effects in altitude and Solar Zenith Angle (SZA). The nightside ionosphere was created by adjusting the photoionization rate until the equilibrium electron densities agreed with the results of Theis and Brace (1993). This mimics ionization by particle BRECHT AND LEDVINA 10.1029/2020JA027779 2 of 17 precipitation. The photochemistry module was run until an equilibrium ionosphere was created and then the plasma dynamics were also activated. Convergence testing of the chemistry package with regard to the required time step was performed before the simulations were begun. The ions created by the photochemistry are produced with a velocity consistent with the local neutral wind velocity.
The generalized electric field equation solved in the code (cf. Brecht et al., 2016) contains the ambipolar term, ∇P e , as well as the Hall and Pedersen conductivities. Equation 1 is the inertialess electron momentum, m e = 0, equation, that is solved for the electric field, E, (1) where the electron current is J e = ∇ x B−J i , the total current is J = ∇ x B, and   is the conductivity tensor containing the Hall and Pedersen conductivity. The requirement of quasineutrality means n e = n i . The Hall and Pedersen terms are taken from a text by Mitchner and Kruger (1973). The electron temperature, T e , needed for the ambipolar term is actively calculated via an energy equation (cf. Brecht et al., 2016) down to 400 km altitude. Thermal conduction is not included in the calculation. Below 400 km, T e , is held constant using a PVO based altitude dependent model. At the start of the simulations, the electron temperature is linearly ramped up from the PVO value at 400 km altitude to the solar wind value at 2 R v . This insures that there are no wild artificial pressure gradients at the start of the simulation. As the simulations are run the solution to the electron energy equation smoothly matches the PVO value at 400 km altitude.
The electron temperature profile below 400 km is taken from Figure 2 of Shinagawa et al. (1987). The Shinagawa et al. (1987) T e profile comes from the Langmuir probe and retarding potential analyzer on PVO and is extended to lower altitudes using the calculations of Cravens et al. (1980). The electron temperature profile used in the simulation falls off with Solar Zenith Angle (SZA) in accordance with the empirical model of Theis et al. (1984). The electron temperature profile and its SZA dependence are for solar maximum conditions which is inconsistent with the neutral atmosphere used in the simulation. However, Brace and Theis (1996) have shown that the electron temperature profile between the altitudes of ∼170-400 km is less for solar maximum than it is for solar minimum and has a gentler slope. Thus, our use of the solar maximum electron temperature profile results in a smaller ∇P e , than would be expected during solar minimum conditions. Figure 2 shows the nightside plasma density and the electron temperature. These profiles reflect the starting conditions of the simulations that drive the instabilities.
There is some ambiguity in the in the meaning of the term "ambipolar electric field" (also known as the polarization electric field). Some consider the ambipolar electric field to be the electric field that decouples the plasma from the frozen in field assumptions. Many in the fusion and space physics communities consider the ∇P e term in Equation 1 to be the ambipolar electric field term (cf. Cravens, 1997;Shunk & Nagy, 2004). The Astrophysics community takes the ambipolar electric field to be due to the ion-neutral collisions represented by the conductivity tensor in Equation 1 (cf. Ballester et al., 2018). Still others consider the ambipolar electric field to be the field that balances out the ∇P e term in an ionosphere (cf. Collinson et al., 2019;Gombosi, 1998)   There were two main simulations and several test simulations performed.
In the first simulation all terms of the hybrid code were turned on, specifically the ambipolar electric field term, (▽P e ), Equation 1. The neutral atmosphere used had no super-rotation (cf. Alexander, 1992;Alexander et al., 1993;Lebonnois et al., 2010;Schubert et al., 2007) in it. Thus, the flow around the planet was essentially symmetric. This is an idealization but makes a good starting point for further study. The second simulation to be discussed uses the same neutral atmosphere and flow patterns but has the ambipolar term set to zero (▽P e = 0).
The 3-D simulations used a solar wind speed of 400 km/s with an ion density of 14 particles/cm 3 . The magnetic field strength was 10 nT (B x = −7.8 nT, B y = 6.3 nT, B z = 0.0 nT) producing a Parker spiral angle of 39°. This Parker spiral angle means that the field is pointing away from the planet at ∼140° SZA on the dawn side; see Figure 3b. The cell size of the Cartesian grid upon which the electromagnetic fields are solved and applied to the particles was 60 km in each direction. The particle count per cell is variable with over 1,000 particles per cell in the ionospheric region, and roughly 80 per cell in the subsolar shock region. The spherical grid where the photochemistry is performed and upon which collisions take place was 0.44° in angular resolution and 5 km in radial resolution. The spherical grid begins at 100 km altitude and extends to 450 km altitude. The overall Cartesian grid covers 3.026 × 10 9 km in each direction (5 R V × 5 R V × 5 R V ) with the Venus at the center of the grid.

Results
All results are reported in a coordinate system where x is directed toward the sun, y is defined for a right-handed system in the ecliptic plane, and z completes the system by pointing perpendicular to the ecliptic plane. Figure 3 displays the magnitude of the magnetic field |B|, via two cuts through the simulation, meridional and equatorial. This figure is shown in order to orient the reader to the magnetic features with respect to the ionospheric densities shown in the following figures. Figure 3a shows the meridional (x-z) cut through the center of the planet. One sees asymmetry in the shock shape. The barrier field and the field in the +z direction are responding to the convection electric field that points in the positive +z direction. One also sees the fine scale electromagnetic EMIC waves in the bow shock region. Figure 3b displays the equatorial (x-y) plane, one notes a more symmetric shock but one that is disrupted by upstream waves propagating along the interplanetary magnetic field with a 39° Parker Spiral. In both cuts through the simulation, one can see asymmetries in the magnetic field on the night side of the planet but overall the night side is relatively quiet. The concentric circles delineate altitudes in 500 km increments.

Nightside Density Structure
The left hand panel of Figure 4 shows an equatorial cut through the simulation with contours of the O + and O 2 + density displayed. It is easy to see that the nightside ionospheric structure is offset from the midnight line. The reason for the offset will be discussed later in the study. The right hand panel of Figure 4 is a meridional cut along the angled line shown in the left hand panel. In the right hand panel of Figure 4 one notes a plume of plasma being lifted off on the dayside in the direction of the convection electric field (+z).
For the purposes of this study, no further discussion of this feature will be provided other than to note a similar feature has been seen in Mars simulations (cf. Brecht et al., 2016 andLedvina et al., 2017) and in the MAVEN data. The nightside structure displayed in the right hand panel of Figure 4 shows the two lobes of density coming off in the nightside separated by the equator. Such structure was reported Brace et al. (1987) BRECHT AND LEDVINA 10.1029/2020JA027779 4 of 17 when examining the PVO data. It is clear that the ions leaving in the tail or "tail rays" are structured and that the structured region is of finite extent. The fact that the structured region is localized means that these features are not always seen as is the actual case with PVO and Venus Express (VEX).
The densities shown in Figure 4 are very comparable to those reported by Brace et al. (1987) Brace et al. (1987) reported (their Figure 3 and Figure 2) maximum ionospheric ion densities ranging from 10 2 to 10 3 particles cm −3 at altitudes above 2,000 km. The structure reported in Brace et al. (1987) ranged from a width of 30° latitude down to a width of a few degrees latitude. At the surface of Venus, 1° represents roughly 105 km. The structure reported ranged from ∼3,000 km down to a few hundred km. This is well within the resolution of the simulation code gridding. The size of the structures seen in the simulations ranged from 120 km on up. There may be a few structures of one cell (60 km) but they are not obvious. Larger structure seen in the simulation could range up to 2,000 km but like the data ( Figure 2 and Figure 3 from Brace et al., 1987) the large-scale structure has smaller scale structure residing within it.
A repeat of the previous simulation with the ambipolar electric field set to zero, ▽P e = 0, was performed. Figure 5 displays O + and O 2 + densities on the same cut as the right hand side of Figure 4. One can see a complete lack of structure formation. Examination of the whole nightside reveals no structure formation at any location. These results strongly suggest the reason for the "tail ray" structure on the night side of Venus is the ambipolar electric field created by the gradient of the electron pressure. It is clear from these plots that the simulations have not been run out for a long enough period of time to assess overall ion loss rates. That was not the purpose of this particular piece of work. The purpose was to assess the cause of the structure seen forming in the simulations after a very short time (5 min of real-time). For example, it is expected that even the run with no ambipolar electric field will eventually show ion loss from the J x B forces that clearly operate on a much slower time scale. In fact, the difference in Figures 4 and 5 give a strong indication of the speed of instability growth compared to magnetohydrodynamic forcing terms. If one only mentioned dependence on the ambipolar electric field seen in the simulation results, there would naturally be some skepticism. If the other aspects of the results aligned/ agreed with observations, then there is more credence to the conclusions drawn from the simulations. In the next section, a more detailed examination of the simulations and comparison to reported data takes place.

Comparison of Mars and Venus
There is another aspect of the ambipolar behavior to be discussed. In Brecht et al. (2017) it was reported that the ambipolar electric field made a significant difference in the loss rate of O + and O 2 + from Mars. Further, using MAVEN driven electron temperatures made a difference in the ratio of the loss rate of the two ion species. However, no evidence was found in the Mars simulations of structure similar to the type reported here. A quick investigation of both the Mars simulations and the Venus simulation reported here suggests an answer. Figure 6 shows the resulting ▽P e generated electric field created at t = 0 s for both the Mars simulation and the Venus simulation. One sees immediately that the peak ambipolar electric field in the Martian simulations is about a factor of 5 smaller than created in the Venus simulation. Further, the Venus simulation shows the peak field to be in the region where the ionospheric density of Venus is still large and below the 400 km transition between the data driven T e profile and the calculated T e . The peak of the ambipolar field in the Martian simulation is above the high-density ionospheric region (>300 km). Figure 6 also shows the Venus ambipolar electric field as the simulation ionosphere evolves. One sees a drop in the maximum strength of the field but it is still considerably above that seen at t = 0 for the Mars simulations. There are also more peaks in the profile found at t = 300 s due to waves in the simulations, but overall the structure is very much the same. Please note: that Figure 6 is a single radial profile, and is not an average of the whole nightside ambipolar electric field.
There is evidence that the Martian ambipolar electric field is reasonable. Akbari et al. (2019) estimate the ambipolar electric field on the dayside of Mars to be 1-2 µV/m based on MAVEN data. This reported value is consistent with the ambipolar electric field strengths found in the code on the nightside of Mars, Figure 6. On the nightside of Mars the density gradient is sharper due to the cooler atmospheric temperature, but the electron temperature is reported to be warmer (Fowler et al., 2015).
At Venus Collinson et al. (2019) make a theoretical estimate of their definition of the ambipolar electric field and find that the potential estimated from VEX is 10 times greater (−9.9 V ± 1.1 V) than their estimate of the potential drop based on the electron pressure gradient (−0.9 V). However, their estimate of the potential drop was focused on the dayside using a draped magnetic field line model. The Venus profile at 300 s, Figure 6, shows changes in the sign of the ambipolar field as the density and electron profiles adjust making an integrated potential drop difficult to estimate. An ambipolar field of roughly 10 µV/m from 200 km to 300 km has a potential drop of 1.0 V which is consistent with estimates of Collinson et al. (2019). Taking an average of −5 µV/m from 300 km to 1,000 km one obtains −3.5 V. Again, in the range of estimates from Collinson et al. (2019). It should be noted that the profiles shown in Figure 6 are simply the fields along the radial line through local midnight. The electric fields vary with space and time especially as large waves move through the system. Thus, Figure 6 is just an illustration of the field magnitudes found within the simulation at any given time.

Earlier Simulations
A reasonable question is: Why have earlier simulations failed to produce this structure? It is unlikely that MHD simulations would produce them because they often do not have the ambipolar term in them (cf. Ma et al., 2013 andTerada et al., 2009).
Hybrid simulations of Venus have been performed but have not, until now, seen this structure. Kallio et al. (2006) and Jarvinen et al. (2008Jarvinen et al. ( , 2013Jarvinen et al. ( , 2015 used a cell size of ∼300 km and had no chemistry or ▽P e . Martinecz et al. (2009aMartinecz et al. ( , 2009b) used a smaller cell size ∼120 km. However, they had no chemistry BRECHT AND LEDVINA 10.1029/2020JA027779 6 of 17 The solid line is the ambipolar electric field on the nightside of the Mars simulations. The dotted line is the ambipolar electric field on the nightside of the Venus simulation after it has evolved for 300 s. Note that the Venus profile weakened but still has a much stronger field than Mars at the lower altitudes of the ionosphere.
or collisional affects but did have a ▽P e term. They used a two electron temperature model, one for the solar wind and one for the ionospheric species and assumed both behave adiabatically. Their study did not make it clear what the actual T e profiles were. Terada et al. (2002) used a 2-D simulation and focused on the ionopause near the terminator. The T e profile comes from Shinagawa et al. (1987) but they do not use a SZA fall off. The smallest grid size is 10 km but it ramps up rapidly and is ∼100 by 500 km altitude. They have a four species ion chemistry and include the ▽P e term. However, the structures seen in the data are three dimensional and in the 2-D simulation, B is perpendicular to the simulation plane.
The simulations reported here have an active photochemistry resolved to 5 km radial cell size along with collisional interactions with the neutral atmosphere also resolved to 5 km radial scale size. Further, the code has the ambipolar electric field and a total electromagnetic resolution of 60 km. Behavior such as ion gyration is included because of the kinetic treatment of the ions. As such, this set of simulations is unlike any performed before and thus it is not surprising that the structure is now manifesting itself. Grebowsky and Curtis (1981) listed six properties of what they called "ionospheric holes." The properties are quoted below:

Nightside Observations
1. "Electron densities in the holes have been observed to drop below 10 cm −1 , i.e., several orders of magnitude lower than in adjacent nightside regions 2. Holes are not typically detected below 200 km 3. Two electron components are detected in the holes -a cold ∼2000°K component and a hot component exceeding 1 eV in temperature 4. The magnetic field in the holes is vertical and enhanced (10-40 nT) compared to adjacent regions 5. The latitudinal extent of the individual holes is 600-1,800 km 6. There are two broad zones (20° in latitudinal extent and ∼3 h in local time extent) where holes are detected. These zones, one north and one south of the equator are symmetrically positioned about ∼10°N latitude and ∼1.5 h in local time while the equatorward boundaries of the two zones are separated by ∼20° in latitude."

Simulation versus Observation
In the next few paragraphs, the six observational features listed in the previous section and discussed in Grebowsky and Curtis (1981) will be compared to the simulations. Features 1 and 2 are found in the simulations results. Figures 7 and 8 show a 3-D rendering of the nightside region where the ionospheric structure manifested itself. These figures are officially called "volume rendering" in 3-D with transparency used to allow one to look through the structure. What this means is that every grid point has a density that is color coded as shown in the color bar.
One is looking at total density in the region with some transparency. The details of these Figures 7 and 8   In Figure 7 one can see a relatively even ionospheric density at altitudes ranging from 100 km altitude (the bottom of the grid) to 200 km. The densities in this region were roughly 10 4 to 10 5 oxygen ions per cm 3 . Above that altitude region, significant structuring of the ionosphere is seen. Although the plot was terminated at 2,000 km, one notes that there is altitude dependence with the structure and of course its density. Figure 8 shows the Figure 7 viewed from the top down. From this view of the 3-d ionosphere one looks into the "holes" in the ionosphere. One can look down into ionospheric density and see levels at lower altitudes. The fact that one does not see regions where the density is 10 4 to 10 5 suggests that the "holes" do not extend down to below 200 km consistent with property number 2. One also sees regions with densities at or below 10 1 again the results are consistent with property number 1.
Property 3, the electron temperature, cannot be compared because of the inability of the code to address a two-component temperature. The code in its current configuration cannot address this observation.
Property number 4 addresses the magnetic field structure. It was stated that the magnetic field enhancements were anticorrelated with the ion density enhancements. The simulations produced results similar to Grebowsky and Curtis (1981). Figures 9 and 10 reveal different aspects of the magnetic field in the structured region. Figure 9 is a "volume rendering" in the same fashion as was Figures 7 and 8 and covers the same spatial region. In Figure 9 there are regions where the local magnetic field strength is in agreement with the data with regard to the enhancement of the field. Also in agreement is the increase of the magnetic field magnitude with altitude, Figure 9. One sees fields up to roughly 40 nT is some regions, which is BRECHT AND LEDVINA 10.1029/2020JA027779 8 of 17   Figure 9. Vectors of the magnetic field at the top layer (2,000 km altitude) of the 3-D volume are displayed, Figure 10. The magnetic field intensity at 2,000 km is plotted as the color contour with the magnetic field. The structured oxygen ion density is shown as a black contour line selected to be 10 particles per cm 3 . The magnetic field is mainly radial to the planet either away from the planet or toward it. 30 nT above the solar wind IMF (10 nT). The density enhancements are not correlated with the magnetic field enhancements, although they are not 100% anticorrelated as the studies referenced above suggested. Please compare Figure 7 with Figure 9 to see that field enhancement resides mainly in a hole. Figure 10 is a color contour plot of the top surface of Figure 9 with field vectors normalized to the solar wind IMF. The dark contour outlines the region of enhance ionospheric density.
One clearly sees in Figure 10 that the magnetic field at 2,000 km is mostly pointing away from the planet. This is in agreement with property 4 from above. Figures 11-13 show the ion density as a surface plot with the outer edges of the density enhancement demarked by black lines. The color contours represent the magnitude of B normalized by the solar wind IMF. The log10 of density is offset by 1 to allow viewing of the surface plot and the color contours. These plots reveal the changing location and magnitude of the structure with altitude as one scans from 1,000 to 2,000 km. Again, the density and magnetic field are not highly correlated but they are not completely anticorrelated as suggested by Brace et al. (1987) and Grebowsky and Curtis (1981).
Property 6 includes among other things the fact that the nightside structured region is localized and does not cover the entire nightside of the planet. Figures 4 and 8 indicate that the simulation results are in agreement with the previous statement. Even the general overall location of the structured region seems to be consistent with the preponderance of structure observations, (Marubashi et al., 1985). It is clear from Figures 11-13, that depending on the altitude of the spacecraft and its trajectory through the structure, one is challenged to compare observation and simulation details. This is particularly true when one considers that the trajectories are not strictly meridional. Marubashi et al. (1985) displays orbital tracks that observed the nightside structure (see his Figure 2). We note that many of the reported orbits were on the dawn side of the midnight longitude. This is consistent with where the structure was found in the reported simulations. However, Marubashi et al. (1985) report passes by PVO that found structure in other parts of the nightside.

Control of the Structure Location
The authors conjectured that the offset of the tail rays is due to the IMF direction. In order to test this conjecture a simulation was begun with the equatorial component of the IMF, By, reversed in sign. Thus, the IMF was changed to (B x = −7.8 nT, B y = −6.3 nT, B z = 0.0 nT). The simulation revealed structure development on the opposite side of the midnight line at the same SZA as the structure shown in Figure 4. Thus, it appears that the location of the structured ion tail depends on the IMF direction. One should note that these simulations were performed with the interplanetary magnetic field loaded within the planet. The purpose of doing this was to avoid having any regions where a current was produced artificially by the boundary conditions (the lower region of the ionosphere, 100 km altitude).
However, there is an additional parameter to consider and that is the orbital velocity of Venus. It is roughly 35 km/s or approximately 10% of the solar wind velocity. This parameter remains unexplored in this study BRECHT AND LEDVINA 10.1029/2020JA027779 9 of 17

Figure 12.
A slice through the 1,500 km altitude on the night side just to the dawn side of midnight. This is the same as Figure 11 but 500 km higher.
but is under active investigation and will be reported soon. However, the overall features found in the simulations seem to have about the right spatial extent and the finding of the nightside "holes" as reported by Marubashi et al. (1985) are explained by changing IMF orientation. Brace et al. (1987) discuss some of the features seen when observing the "tail rays" as they called them. The primary one was the presence of parallel or radial magnetic fields in the regions where the density was low. Grebowsky and Curtis (1981) discuss in their study the idea that the oxygen ions are accelerated in the structured region. The finding that the ambipolar electric field, ▽P e , term, is the controlling term is consistent with their conclusion. The ambipolar electric field does accelerate the ions and drive the instability that filaments the nightside ionosphere. Grebowsky and Curtis (1981) also suggested that this field was the source of the energetic electrons seen in the electron distribution function. Due to the limitations of the hybrid code, we cannot test this theory, but it does seem consistent with the features produced in the simulations. Brace et al. (1987) discussed at length the measurements of the O + and O 2 + energies and currents. They listed a variety of characteristics of the observations. We quote those here:

Observations of ionospheric Ions
1. "The nightside ionosphere becomes increasingly ray like and filamentary with distance downstream. Like cometary tail rays, a central ray is a common feature, with one or more rays often arranged symmetrically on either side 2. At altitude in the range of 2,000-2,500 km, which we have examined in some detail, most of the ions are superthermal O + with energies in the range of 9-16 eV. Virtually no cold thermal ions remain at these altitudes. Super thermal H + is probably also present, but the PVO instruments could not detect them if present 3. A minor, but more energetic ion component (E > 40 eV) is observed in the ionotail, usually well correlated with the 9-16 eV superthermals and determined to be flowing tailward on most occasions. Ions at the intermediate energies are not observable, but the approximate agreement between the superthermal ion density and the electron density is consistent with the absence of the thermal ions and the low density of ions with energies greater than 40 eV. The electron temperature in the tail rays is typically less than 1 eV, not unlike that found in the underlying ionosphere. T e is several times higher outside the rays 4. A strong and steady tailward magnetic field dominates the trough region that surrounds the tail rays.
The magnetic field is weaker in the rays and B x reversals reveal the flow of strong currents in the edges of the rays. An approximate pressure balance exists in the ionotail, with static plasma pressure of rays balanced by magnetic pressure of the surrounding trough regions"

Comparison Between Observations and Simulation Results
Observation number 1 is reproduced in the simulations. Figure 4, Figure 7, and Figures 11-13 all reveal behavior in the simulations similar to that discuss in observation 1. Observation number 2 is also reproduced in the simulations. Figure 14 is a plot of the H + ions entering the nightside ionosphere. These are known to be solar wind protons as we did not include any H + in the ionosphere of the simulation. Note that the H + contours are found where the field is weak. This also agrees with VEX findings (Kollmann et al., 2016). Figures 15 and 16 address both observations 2 and 3. Figure 15 shows the ion energy spectra computed from O + and O 2 + ions found in the region of the tail rays and holes (1). In Figure 16 a comparable region (2) on the other side of the midnight line is shown.
BRECHT AND LEDVINA 10.1029/2020JA027779 10 of 17 Figure 13. A slice through the 2000 km altitude on the night side just to the dawn side of midnight. This is the same as Figure 11 but 1,000 km higher.

Figure 14.
Color contours of the magnitude of B overlaid with a black contour lines of the H + density at 2,000 km altitude.
Region (1) covering the rays and holes extends from 110° to 170° longitude with ±30° latitudinal extent and altitude increments of 500 km starting with 500 km altitude. The second region (2), Figure 16, is the same except that the longitudinal region is 190°-250°.
BRECHT AND LEDVINA 10.1029/2020JA027779 11 of 17 Figure 15. The energy distribution of oxygen ions at different altitude ranges within the "tail ray" region. The thick gray line is the energy distribution of the O 2 + that is only found in the 500-1,000 km altitude range. The volume sampled is the same as for Figure 15. The thick gray line is the O 2 + found in the 500-1,000 km altitude range and no higher.
The thick gray line is the O 2 + spectra in the range of 500-1,000 km above that altitude range little or no O 2 + is found. The dash-dot line representing the 2,000-2,500 shows a steep drop-off of density at energies just below 10 eV (the escape energy from Venus is 9.03 eV) with a plateau just below 20 eV. This is very consistent with the observations listed by Brace et al. Figure 15 shows a density drop off as one samples higher altitude ranges in the tail ray region. However, above 100 eV the spectra match for each altitude. The energy where the plot was stopped, 10 keV, represents an O + velocity of 346 km/s, which is approaching the solar wind speed (400 km/s).
Our speculation is that the ions at the high end of the spectra are ions coming from the dayside of the ionosphere and have been accelerated up to these velocities. Figure 16 (region [2]) shows similar dependence on energy and altitude but the overall density is lower by more than an order of magnitude. The energy cut off seen in the 2,000-2,500 km altitude band and the 1,000-1,500 km altitude band is again consistent with the escape energy from Venus, 9.03 eV. The gravitational escape energy explains why the peak of the energy spectrum move upward as fewer and fewer low energy ions can reach that altitude. Interestingly, the high energy portion of this region has a higher density; again giving credence to the idea that the very energetic ions are coming from the dayside and are not part of the physics of the tail rays.  Figure 18 representing the nontail ray section (2) shows a more mixed response. The 500-1,000 km ranges dominates in the negative velocity region but is the smallest in the positive velocity side again consistent with gravitational forces being applied. However, unlike the tail ray region the highest density region is not the 2,000-2,500 km range but the 1,500-2,000 km ranges. There is a clear shift of the distribution peak from below escape velocity to above escape velocity as one goes higher in altitude indicating a net tailward flow. One notes that the extent of the velocity range of the distribution is larger in the tail ray region (1) than in the nontail ray region (2). This again supports the assertion by Grebowsky and Curtis (1981) that additional acceleration occurs in the tail ray region (1). Brace et al. (1987) observation 4 also discussed magnetic field reversals as a signature of currents in the tail ray area of the magnetotail of Venus. We examine the region around a tail ray and subsequent holes. In Figure 19 the radial magnetic field at 2,000 km altitude is represented by cones. The upward pointing cones point away from the planet. The downward pointing cones are seen as circles. They are unit cones so the size provides no information as to field strength. The color contour plot shown in Figure 19 is the O + and O 2 + density also at 2,000 km altitude. The section of upward pointing cones is just to the side of the high density region. Figure 19 shows regions of high magnetic shear with regard to orientation of the magnetic field. This is consistent with observation number 4 from Brace et al. (1987).
In summary, the results presented in this paper reveals that the driver for the "parallel electric field" discussed by Grebowsky and Curtis (1981) is the ambipolar electric field and the "holes" and "tail rays" come from this driver. The particle spectrum results from the simulations are consistent with the observation reported by Brace et al. (1987) and Kollmann et al. (2016) who reported both tailward and planetward flow of O + . We could not address the pressure balance issue due to the limitations of the codes ability to handle the electron temperature.
BRECHT AND LEDVINA 10.1029/2020JA027779 13 of 17 Figure 18. The oxygen ion radial velocity distribution at differing altitude ranges at a symmetric location on the opposite side of the midnight line. The sampling volume is the same. The thick gray line is the O 2 + found only in the 500-1,000 km altitude range. Figure 19. The radial magnetic field at 2,000 km near a region of tail rays. The cones point radially up or down and are of unity magnitude. The color contour is of the O + and O 2 + ion density. One can see outward pointing cones as well as the circles that are the downward pointing cones.

Issues with Nightside Density
There is an interesting aspect to the simulation and that is the inward flowing H + . Figure 14 reflects the fact that H + deposition is not a uniform. The protons are depositing where the magnetic field is weak in the tail and not where it is strong, consistent with observation from VEX, (Kollmann et al., 2016). One also sees O + ions outflowing in this region as well again consistent with Kollmann et al. (2016).
The plots of the H + coming into the nightside suggest that the simulation of the nightside ionosphere must be addressed. A production rate on the nightside of the simulation was selected via trial and error to reproduce the electron profiles reported by Theis and Brace (1993). There was no pretense to adding physics in an unexplored region. Now that the simulations have been performed, it is thought that two sources of nightside ionosphere production could be added to the code. One is charge exchange of the H + with the nightside atmosphere leading to subsequent chemistry, and charge exchange with inflowing O + and O 2 + . Kollmann et al. (2016) found that there is approximately half as much O + being deposited on the nightside of Venus as leaving Venus altogether.
The second source is addition of electron impact ionization on the nightside. A temperature dependent electron impact ionization equation is already included in the code but the electron temperature profile used is cold consistent with a nightside ionosphere. The finding that there is a hot component suggests that we could improve the code by including a secondary electron temperature in the nightside region consistent with the hot component reported by Grebowsky and Curtis (1981). This would require some complex modeling to determine the actual electron density of this component and handle its spatially localized nature.

Plasma Instabilities
The density and structure size have been examined and they seem to be in quantitative agreement with the structure and plasma density reported by Brace et al. (1987) and Grebowsky and Curtis (1981). The simulations provide an answer as to the cause of the structure, but not what instability is at work. Previous theoretical efforts (Huba, 1992;Huba & Grebowsky, 1993;Linson & Workman, 1970) have suggested that the ambipolar processes drive the instabilities that lead to plasma structure. The work by Huba (1992) and Huba and Grebowsky (1993) address the lower hybrid instability which is a relativity high frequency instability with the most rapidly growing modes having short wavelengths on order 0.1 km. Evidence of this instability has been reported by Huba and Grebowsky (1993). In addition, the conditions on the Venus nightside are very similar to the conditions created in ionospheric barium releases (cf. Haerendel & Lüst, 1968;Haerendel et al., 1967;Wescott et al., 1969) whose structure was determined to be the gradient drift instability (Linson & Workman, 1970). The gradient drift instability is most unstable at wavelengths that are on order of the ion gyroradius. This instability seems to be more in line with the structure being produced in the simulations. The fact that the structure occurs in regions where the radial ambipolar electric field aligns with the solar wind IMF suggests that the cause of the structure might be one of many filamentation instabilities. Additional analysis and simulation will be needed to confirm which instability is a likely and consistent candidate. However, the issue here is not which instability was actually driven, it is the source of the driver and that clearly is the ambipolar electric field, ▽P e .

Conclusion
For the first time a simulation of the solar wind interaction with the ionosphere of Venus has shown nightside structures resembling the nightside ionospheric structures observed by PVO. The creation of these structures was found to be driven by the ambipolar electric field created by the gradients in the electron pressure. The location of the tail ray structure is determined by the IMF orientation; suggesting that the tail rays did not produce radial magnetic fields but were created where radial magnetic fields existed. From an instability standpoint, the instability seems to occur where the ambipolar electric field aligns with the local magnetic field.
Confidence in this assertion is bolstered by the qualitative agreement between the simulated structures and the equivalent structures reported by PVO. The structures become ray like and filamentary with distance from the planet. The longitudinal and latitudinal extent of the structures agrees with the observations reported by Brace et al. (1987). The structure location was found on the dawn region of the ionosphere in agreement with the distribution of the observed ionospheric holes reported by Marubashi et al. (1985) and the orientation of the IMF. The simulated ionospheric holes were found to quantitatively agree with the properties reported by Grebowsky and Curtis (1981). The densities in the holes dropped several orders of magnitude below the adjacent nightside regions. The holes did not extend below 200 km altitude. There is an enhanced magnetic field associated with the holes that is not correlated with plasma density enhancements.
A detailed examination of the ion distributions further reinforces the confidence in the simulation results. The bulk of the ions in the structured regions were found to be accelerated, there was a drop off in the low energy (<10 eV) population as altitudes increased consistent with data and the effect of gravitational acceleration. A subpopulation of ions was found to be moving toward the planet, consistent with observation of currents flowing in both directions. Once again these findings are in agreement with the properties of the observed structures (cf. Brace et al., 1987;Grebowsky & Curtis, 1981;Kollmann et al., 2016).
It is natural to ask why these structures found at Venus and not Mars? The answer lies in the difference in the gradients of the electron densities and pressures of the two planets. The higher gravity field at Venus leads to profiles in the electron density and temperature with much stronger gradients than those found at Mars. This in turn leads to a much stronger ambipolar electric field on the nightside of Venus than at Mars.
The ambipolar electric field was found to be the driver of the nightside structures observed at Venus. However, the plasma instability that leads to the structures has not been identified. Possible candidates include the lower hybrid instability (Huba & Grebowsky, 1993), the gradient drift instability (Linson & Workman, 1970) and one of many filamentation instabilities. Further work is needed to determine what plasma instabilities are occurring.
Research is currently underway to examine the effects of the orbital motion of Venus on the interaction and the energy and momentum transferred to the neutral atmosphere by the solar wind particles and fields. The findings of this continuing research will be presented in a more comprehensive paper now being prepared for submission to the Journal of Geophysical Research.