A Statistical Study of Ionospheric Boundary Wave Formation at Venus

Previous missions to Venus have revealed that encounters with plasma irregularities of atmospheric origin outside the atmosphere are not uncommon. A number of mechanisms have been proposed to discuss their origins as well as their roles in the atmospheric evolution of Venus. One such mechanism involves an ionopause with a wavelike appearance. By utilizing the magnetic field and plasma data from Venus Express, we present the first observational statistical analysis of the ionospheric boundary wave phenomena at Venus using data from 2006 to 2014. Results from the minimum variance analysis of all the photoelectron dropout events in the ionosphere reveal that the ionopause of Venus does not always appear to be smooth but often exhibits a wavelike appearance. In the northern polar region of Venus, the normal directions of the rippled ionospheric boundary crossings lie mainly in the terminator plane with the largest component predominantly along the dawn‐dusk (YVSO) direction. The average estimated wavelength of the boundary wave is 212 ± 12 km, and the average estimated velocity difference across the ionopause is 104 ± 6 km/s. The results suggest that the rippled boundary is a result of Kelvin‐Helmholtz instability. Analysis reveals a correlation between the normal directions and the locations of the boundary wave with respect to Venus. This indicates that the draping of magnetic field lines may play a role in enhancing the plasma flow along the dawn‐dusk direction, which could subsequently set up a velocity shear that favors the excitation of ionospheric boundary wave by the Kelvin‐Helmholtz instability along the dawn‐dusk direction.


Introduction
Due to the absence of an intrinsic magnetic field (Luhmann & Russell, 1997;Russell et al., 1980), the solar wind interaction with Venus is highly dynamic. The Venusian ionopause, which is a boundary separating the shocked solar wind plasma and ionospheric plasma, is subjected to a number of plasma instabilities. The two main instabilities are the Rayleigh-Taylor instability (also known as the interchange instability; Arshukova et al., 2004) and the Kelvin-Helmholtz instability (KHI). The interchange instability only grows when there is a nonmonotonic plasma pressure gradient at the subsolar region (Arshukova et al., 2004). The KHI is a macroinstability that is principally generated by the strong shear flows across a boundary (Chandrasekhar, 1961) and is an important mode of energy transfer at Venus (Futaana et al., 2017). Important factors to the generation of KHI waves include the velocity, density, and temperature gradients (Amerstorfer et al., , 2010Biernat et al., 2007;Ferrari et al., 1982;Huba, 1981;Price, 2008;Wolff et al., 1980). Even though the KHI is considered the more dominant instability for wave excitation, there are occasions when terms such as the magnetic field stress, gravity, and boundary curvature are more significant and can give rise to other instabilities, for example, the Rayleigh-Taylor instability or flute instability (Elphic & Ershkovich, 1984).
The concept of ionospheric boundary waves at Venus has long been studied in a number of model simulations. Wolff et al. (1980) showed that the ionospheric boundary is unstable to KHI and illustrated the formation of flux ropes and atmospheric bubbles as a result of the ionospheric surface wave. Terada et al. (2002) used a two-dimensional global hybrid model to investigate the KHI at Venus and showed that the ionopause in the subsolar region is unstable to KHI. Biernat et al. (2007) studied the growth of the KHI at the Venusian ionopause and found that KHI can evolve regardless of the solar wind conditions. Amerstorfer et al. (2010Amerstorfer et al. ( , 2007 showed that a density increase can influence the growth rate of KHI and characterized the evolution of the KHI into three main phases, that is, linear, nonlinear, and turbulent phases. However, Möstl et al. (2011) showed that the ionopause is not able to reach the nonlinear vortex phase during either low or high solar activity due to the stabilizing density jump across the ionopause. In addition to Venus, the development of the KHI has been studied at other planets, including Mars (Penz et al., 2004), Mercury (Sundberg et al., 2010), Earth (Nykyri & Otto, 2001), and Saturn (Masters et al., 2009).
The continuous process of plasma loss resulting from ionospheric boundary wave events over a prolonged period of time plays a significant role in contributing to atmospheric loss at Venus. However, observational studies are rather limited and only short periods of Pioneer Venus Orbiter (PVO) data (one to two orbits of data) have been utilized. Using PVO data, Luhmann (1990) reported one of the earliest observations of an ionospheric boundary wave, which is interpreted as the terminator wave, based on the significant field behavior change within 15 ∘ of the terminator. The authors suggested a possible mechanism by the reverse orientation of the interplanetary magnetic field (IMF) in the ionosphere and suggested its association with the formation of flux ropes. Brace et al. (1980) observed wavelike structures and interpreted them as PVO passing through the ionospheric surface wave, which can be due to the altitude changes of the Venusian ionopause. Brace et al. (1983) reported transterminator ionospheric waves into the nightside and suggested a potential wave energy by the plasma pressure gradient-driven interchange instabilities or the ion-neutral drag-driven shear instabilities. Walker et al. (2011) andPope et al. (2009) suggested that nonlinear vortex-like structures observed in the magnetosheath region using Venus Express (VEX) magnetic field data were associated with the strong shear flow across the ionopause. However, the altitudes of these observations are not consistent with the nominal ionopause altitude. Chong et al. (2017) presented evidence for the ionospheric boundary exhibiting a wavelike appearance for a single VEX pass along the terminator.
To study the dynamics and characteristic distributions of the ionospheric boundary wave on Venus, a statistical analysis is conducted in this paper. Analysis of the available magnetic field and plasma data from the instruments on board of VEX from 2006 to 2014 reveals that the observations of such phenomena are not uncommon. Investigation of boundary wave formation is fundamental to our understanding of the atmospheric loss mechanisms operating and hence the atmospheric evolution of Venus and unmagnetized planets in general. The paper is structured as follows: The instrumentation is summarized in section 2, the observations and data analysis of the ionospheric boundary wave are presented in section 3, the characteristics of the boundary wave and its possible generation mechanisms are discussed in section 4, and the summary and conclusions are presented in section 5.

Instrumentation: VEX
VEX had an elliptical polar orbit with periapsis ranging from 130 to 463 km at a latitude of about 78 ∘ N. The apoapsis distance was around 66,000 km, and VEX had an orbital period of 24 hr (Svedhem, Titov, McCoy, et al., 2007;Titov et al., 2006). The magnetic field was measured by the VEX Fluxgate Magnetometer (MAG; Zhang et al., 2006). The 1-Hz MAG data are used for the statistical analysis in this paper. These data have been cleaned to remove the dynamic stray fields  and corrected for offset (Leinweber et al., 2008). The data have been rotated into the Venus solar orbital (VSO) coordinate frame, with +X VSO in the Venus-Sun direction, The electron spectrometer (ELS) and ion mass analyzer (IMA) are part of the Analyser of Space Plasmas and Energetic Atoms on board of VEX . ELS provides electron energy spectra in several modes, two of which is used in this study: An electron spectrum between 0.9 eV to 15 keV is generated every 4 s with an approximate energy resolution (ΔE∕E) of 7% (the energy resolution is energy and sector dependent) and an electron spectrum between 9 and 250 eV is generated every 1 s with an approximate energy of 7%. The ion measurements provided by the IMA cover the energy range 0.01-36 keV/q with a sampling time of 192 s and an energy resolution (ΔE∕E) of 7%.
In contrast to PVO, which was able to sample the Venusian ionosphere in the subsolar region over the period of solar maximum, the high-latitude elliptical polar orbit of VEX provides an opportunity to study the dynamics of the Venusian ionosphere in the northern polar region of Venus across nearly a full solar cycle with rather quiet solar activity (Futaana et al., 2017). Figures 1 and 2 show the (a) VEX trajectory, (b) and (d) the electron energy-time spectrogram of differential energy flux, (c) and (f ) the 1-Hz magnetic field magnitude, and (e) average electron energy flux at 22 eV from VEX orbits on 8 November 2011 and 2 October 2011, respectively. Data in (e) are smoothed using a moving average filter of seven data points. The location of the bow shock is highly variable due to the variations in the solar extreme ultraviolet, solar wind Mach number, and IMF orientation (Zhang et al., 2008a). Here the observed altitudes of the bow shock on both orbits are comparably different than the nominal bow shock locations (Zhang et al., 2008a;Zhang et al., 2008b). On the inbound leg of the orbit occurring on 8 November 2011, the bow shock was crossed at around 07:06 UT (at an altitude of 3,395 km) with shocked solar wind appearing in the magnetosheath region around Venus. The broad energy intensity of these electron populations can be seen to become narrower toward the magnetic barrier at around 07:15 UT where the magnetic field magnitude increases until the inbound ionopause was crossed at around 07:18 UT. The ionosphere region is identified by the observation of the ionospheric photoelectron population at 21-24 and at 27 eV (Coates et al., 2008;Cui et al., 2011). This can also be observed as an increase of the electron energy flux (averaged at ∼22 eV) in Figure 1e. These photoelectron populations are mainly due to the photoionization of atmospheric oxygen by solar HeII 30.4-nm photons (Coates et al., 2008). VEX then crossed the ionopause and bow shock (not shown here) on its outbound pass, which can be similarly characterized as described above.

Photoelectron Dropouts
The VEX observations of the bow shock through to the magnetic barrier region on 2 October 2011 are similar to those on 8 November 2011. However, the behavior of the ionosphere on these two orbits is quite different. This is easily seen by comparing Figures 1d and 2d. For instance, the ionospheric photoelectron population on 8 November 2011 is observed continuously within the ionosphere. In contrast, there were 10 separate intervals when the ionospheric photoelectron population disappeared, while VEX was in the ionosphere region on 2 October 2011. These will be termed photoelectron dropout events throughout this paper. These 10 intervals of photoelectron dropouts can also be clearly reflected from the dips in the average electron energy flux at 22 eV in Figure 2e. The photoelectron dropout events on 2 October 2011 are not uncommon. After excluding the orbits where the ELS and/or MAG data are unavailable when VEX was in the ionosphere, as well as orbits when ionospheric photoelectron populations are not observed at all, around 23% (495 orbits) show at least one or more intervals of photoelectron dropouts out of the remaining 2,141 orbits from April 2006 to November 2014.
Throughout this study, the first observation of photoelectrons during an orbit is identified as the inbound ionopause crossing and the last observation of photoelectron population is identified as the crossing of the outbound ionopause. The ionosphere is the region between the inbound and outbound ionopause.
It can be seen from Figure 2 (apart from one interval labeled n/a, which is associated with a gap in magnetic field data and which will be omitted from the following analysis), each of the nine blue-shaded photoelectron dropout intervals on 2 October 2011 corresponds to an increase in the magnetic field magnitude. Compared to the ionospheric regions of low field magnitude, the larger magnetic field magnitude regions are comparable to that observed in the magnetic barrier region just before the first ionopause crossing. In addition, the electron energy intensity during the photoelectron dropout intervals is also comparable to the intensity in the The VEX trajectory plot in R V (Venus radii) that is colored in yellow, red, purple, and green, respectively, to show the different regions (unshocked solar wind, magnetosheath, magnetic barrier, and ionosphere) in which VEX was passing. These regions are also reflected by the colored bar above Figures 1b and 1d. The nominal altitudes of bow shock (BS), induced magnetopause (MP), and ionopause (IP) represented using blue, red, and green dashed lines, respectively (Zhang et al., 2008a;Zhang et al., 2008b) Figure 1, and the ionosphere perturbation is from 05:46:00 to 05:54:30 UT. All 10 photoelectron dropouts are shaded in blue. There are missing magnetic field data in the photoelectron dropout interval labeled n/a. The yellow-shaded regions correspond to the dip-to-peak and peak-to-dip magnetic field fluctuations immediately adjacent to the photoelectron dropouts regions. VEX = Venus Express. magnetic barrier region as seen from Figures 2d and 2e. Note that for the purposes of readability, the plot of average energy flux in Figure 2e is seven-point smoothed. As a result, the presented energy fluxes of a couple of photoelectron dropout intervals appear to be larger than the flux intensity in the magnetic barrier.
The occurrence of the photoelectron dropouts implies that while VEX was in the ionosphere (where ionospheric photoelectron population should be constantly observed), there were periods when VEX detected electron population similar to those in the magnetic barrier region. However, to travel from the magnetic barrier to the ionosphere region or vice versa, VEX would be expected to cross the ionopause.

Ionospheric Boundaries Crossings
To assess if these photoelectron dropout events on 2 October 2011 relate to ionopause crossings, all of the 18 (yellow-shaded) regions adjacent to the periods in which photoelectron dropout events were observed were investigated using minimum variance (MV) analysis (Sonnerup & Scheible, 1998). MV analysis is implemented over the dip-to-peak and peak-to-dip field fluctuations to determine if they are ionopause crossings. In the case that they are boundary crossings, MV analysis is used to find the boundary normal directions. Table 1 shows the components of the MV directions in the VSO coordinate system, the intermediate-to-minimum eigenvalues ratio ( int ∕ min ), the ratio of average magnetic field component along the MV direction (B n ) to the larger field magnitudes on either side of the discontinuity (|B n |∕|B|), the ratio of the change of the field magnitude to the field magnitude (|ΔB|∕|B|), the angle between |B n | and |B|, B n _B , and the number of data points during the intervals analyzed.
The intermediate-to-minimum eigenvalues ( int ∕ min ) of all 18 intervals are greater than 3.5, which implies that the MV direction is well defined. All the |B n |∕|B| < 0.14 (mean value of 0.06) and |ΔB|∕|B| > 0.63 (mean value of 0.86). These values are well within the criteria for tangential discontinuity; |B n |∕|B| < 4 and |ΔB|∕|B| ≥ 0.2, indicating that this boundary represents a tangential discontinuity (Knetter et al., 2004, and references therein). In addition, all the angles between |B n | and |B|, B n _B > 81.8 ∘ (mean value of 86.5 ∘ ), which is approximately 90 ∘ further indicating all of the nine intervals of photoelectron dropouts are bounded by tangential discontinuities, a typical characteristic of the Venusian ionopause (Wolff et al., 1980).
The same approach is applied to all the dip-to-peak and peak-to-dip field fluctuations immediately adjacent to the intervals of photoelectron dropout identified across the full data set. Note that field magnitude in a magnetized ionosphere is similar to that observed in the magnetic barrier region, unlike for the unmagnetized case. Hence, for the orbits when photoelectron dropouts are observed, the dip-to-peak and peak-to-dip field fluctuations in a magnetized ionosphere cannot clearly be identified. This results in only 371 unique orbits (from a total of 495 events) selected for further analysis. In these 371 orbits, 1,043 intervals of photoelectron dropouts are observed, hence 2,086 field fluctuations. Since MV analysis is only valid with three or more data points, only 1,633 field fluctuations (from a total of 2,086), which have six or more data points, are selected. Analysis conducted using more data points would result in smaller data sets, and less data points would result in a higher statistical uncertainty (Sonnerup & Scheible, 1998). The resulting distributions of the boundary normal directions are similar using between three and nine data points. The use of a minimum of six data points is chosen as a compromise between the number of data sets and the statistical uncertainty.
Around 98% (1,603 out of 1,633 intervals) of all the MV directions have |B n |∕|B| < 0.4 and |ΔB|∕|B| ≥ 0.2, indicating boundaries of tangential discontinuity. Similarly, around 95% (1,562 out of 1,633 intervals) have B n _B > 75 ∘ further indicating that the analyzed intervals are tangential discontinuities. This is again consistent with the characteristics of the Venusian ionopause (Wolff et al., 1980). Furthermore, more than 88% (1,446 out of 1,633 intervals) have int ∕ min ≥ 3, which shows that the MV directions are well defined. These results imply that the multiple photoelectron dropout events observed during each of these orbits are due to VEX traversing the ionospheric boundary multiple times. For example on 2 October 2011, VEX traversed through the ionospheric boundary 9 times on a single trajectory.
Moreover, the dip-to-peak and peak-to-dip field fluctuations occur on the timescale of a few seconds. Since the decay of the magnetic field in the ionosphere is on a timescale from minutes to several hours (Luhmann et al., 1984), the possible scenario of VEX traveled through patches of magnetized and unmagnetized regions of the ionosphere consecutively can be ruled out.

Ionospheric Boundary Waves
In addition to the previous analysis, a total of 251 VEX passes with ionospheric boundary crossings that are similar to the case presented on 8 November 2011 (no photoelectron dropout) is collected. The criteria used to select these passes were (1) no observations of photoelectron dropout events, (2) only gradual ionospheric crossings without magnetic intermediary, and (3) clear passes from high field magnetic barrier to low field unmagnetized ionosphere. Results of MV analysis of the ionospheric boundary show that the MV directions are very well defined with more than 96% (242 out of 251) of the directions having int ∕ min ≥ 3. The results also suggest that the crossings are ionospheric boundary with more than 98% (247 out of 251) having B n _B > 75 ∘ and around 100% (250 out of 251) having |B n |∕|B| < 0.4 and |ΔB|∕|B| ≥ 0.2. All of the normal directions for the 238 (out of 251) ionopause crossings with int ∕ min ≥ 3 and number of data points ≥ 6 are binned in a three-dimensional polar statistical histogram with an azimuthal bin size of 7.5 ∘ and elevation bin size of 3.75 ∘ in Figure 3a. The and are the angles between the locational radial vector from the center of the planet and the X-Z VSO and X-Y VSO planes, respectively. The ranges from 0 ∘ to 360 ∘ . While ranges from −90 ∘ (southern polar point) to +90 ∘ (northern polar point). The color bar at the bottom of the histogram is the number of ionopause crossings in each bin. For ease of comparison, all MV directions with negative Z VSO components are rotated into the positive Z VSO direction and the local positions of the boundary crossings are shifted to the northern polar point ( = 90 ∘ ) so that all directions can be visualized and compared in a single directional hemisphere. As seen in Figure 3a, around 91% (217 out of 238) of the boundary normal directions fall in the elevation range of > 67.5 ∘ .
For a spherically smooth Venusian ionospheric boundary in the dayside and the polar regions as illustrated in Figure 3c, the normal directions of boundary crossings should be aligned radially from the center of Venus at their local locations. This is reflected in the results in Figure 3a. As the locations of boundary crossings are shifted to the northern polar region, the close proximity between the boundary normal directions and the radial directions (from the center of Figure 3a) implies that ionospheric boundary is smooth and quasi-spherical in the Y-Z VSO plane (i.e., Figure 3c).
However, this scenario of smooth Venusian ionospheric boundary is not reflected in the results obtained for the photoelectron dropout events. In Figure 3b, all of the 1,446 well-defined boundary normal directions from the photoelectron dropout intervals are binned in a similar three-dimensional polar statistical histogram. Figure 3b shows that the bins are more populated in the azimuthal range of 60 ∘ < < 120 ∘ and 240 ∘ < < 300 ∘ as well as in the elevation range of < 45 ∘ . In comparison to the smooth boundary case presented in Figure 3a, only around 23% (333 out of 1,446) of the boundary normal directions fall in the elevation range of > 67.5 ∘ . This shows that instead of pointing radially outward from the center of Venus, the majority of the normal directions of the boundary crossings lie in the Y-Z VSO plane, with the most dominant component along the Y VSO axis. In contrast to Figure 3a, the results from Figure 3b imply that, for a single orbital trajectory exhibiting multiple crossings of the ionospheric boundary, the ionopause crossings from the photoelectron  (d) an ionospheric boundary that exhibits a wavelike appearance. The green line represents the ionospheric boundary. The blue dashed arrows represent vectors projected radially from the center of Venus through its local locations. The red arrows represent the normal directions of the boundary crossings projected from their local locations, which are denoted by the orange colored dots. Note that the illustration of the symmetric ionospheric boundary is visualized on the basis that the nominal ionopause altitude is not solar zenith angle dependent (Zhang et al., 2008b). VSO = Venus solar orbital. dropout cases do not result from a smooth ionospheric boundary but an ionospheric boundary that can exist as a ripple along the Y-Z VSO plane and which may propagate in a direction dominantly along the Y VSO axis. This is illustrated in Figure 3d. All the MV directions are obtained from a variety of VEX trajectories with a range of azimuthal angles (not shown here); that is, the result of consistent MV directions shown here is not biased toward any particular VEX orbit.

Flux Ropes
Crossings of magnetic flux ropes (Russell, 1990;Russell & Elphic, 1979;Wolff et al., 1980) should be, in theory, in an idealized hydromagnetic state, similar to crossings of the Venusian ionospheric boundary. Both are tangential discontinuity boundaries (Spreiter et al., 1970;Wolff et al., 1979). Based on its unique potato chip-shaped hodogram as a key identifier (Russell, 1990), a total of 132 magnetic flux ropes resulting in 264 dip-to-peak and peak-to-dip field changes is identified in this statistical analysis. The second photoelectron dropout event on 2 October 2011 that corresponds to a dip-to-peak-to-dip field change from 05:48:41 to 05:48:55 is identified as a flux rope. The hodogram (with int ∕ min ≈ 55 ) has a potato chip shape as shown in Figure 4a. In addition, the hodogram of the boundary crossing (with int ∕ min ≈ 55) from 05:47:56 to 05:48:09 (just before 10.1029/2018JA025644 the flux rope is encountered) is shown in Figure 4b. The start of both hodograms are marked with blue circles. In comparison to the potato chip-shaped hodogram in Figure 4a, the hodogram shown in Figure 4b does not indicate field rotation, and thus, this is not a flux rope.
The 177 of these 264 dip-to-peak and peak-to-dip field changes have int ∕ min ≥ 3 and number of data points ≥ 6. Similar to Figure 3b, the MV directions of only the flux rope crossings as well as photoelectron dropout events excluding the flux ropes are binned in a three-dimensional polar histogram, with an azimuthal bin size of 7.5 ∘ and elevation bin size of 3.75 ∘ in Figures 4c and 4d, respectively. Compared to the non-flux rope cases in Figure 4d, the boundary normal directions of the flux ropes in Figure 4c are slightly more randomly orientated. However, the majority of them can still be observed to lie along the Y VSO axis.
In summary, the analysis of all the ionospheric boundary crossings from 2006 to 2014 reveals that the Venusian ionospheric boundary is not always smooth. In particular, 23% (495 orbits) of all the available orbits show that the ionopause can often exhibit a wavelike appearance in the northern polar region, where all of the boundary crossings are observed. Figure 4e presents an illustrative diagram to show how the wavelike characteristic of the Venusian ionospheric boundary can be visualized from the simultaneous observations of photoelectron dropout events (top panels) and the changes in magnetic field magnitude (middle panels).

Discussions
Due to the polar orbit of VEX, the ionosphere is only sampled in the northern polar region; hence, the ionospheric boundary wave events are observed in a rather small range of locations with elevation >39 ∘ and with 90% of the observations made at >61 ∘ . Its periapsis ranges from around 166 to 1,025 km, with 90% of the passes made below 475 km. The boundary wave events are observed everywhere in this small range of locations and do not show any particular preferred location.

KHI as a Boundary Wave Generation Mechanism 4.1.1. Boundary Wave Widths
With the assumption that the boundary wave is stationary with respect to the spacecraft velocity (VEX has velocity ∼9.5 km/s around the periapsis), both the sizes of boundary wave billows and flux ropes can be estimated from the duration of their crossings along with their respective MV orientation. The half-width of the boundary wave, ∕2 as illustrated in Figure 5e, is estimated by the product of the VEX velocity and the time spent between two consecutive boundary crossings. The estimation is made under two criteria: (1) Only two consecutive boundary crossings when mv , the angle between their respective boundary normal directions, is less than 30 ∘ or greater than 150 ∘ are selected to eliminate the crossings of boundary waves that are still in the linear growth phase; that is, only developed waves are selected; and (2) only boundary crossings when vex , the angle between the spacecraft velocity vector and the boundary normal vector is less than 45 ∘ are selected to eliminate the boundary crossings of the near tips of the waves as illustrated in Figure 5e and to eliminate the boundary, which is crossed at a large angle by the VEX. These criteria yield a total of 89 boundary wave events. Note that the data of VEX, a single spacecraft mission, could only allow the estimations of the boundary normal directions but not the shape of the boundary wave. Hence, the ionopause is expressed in a dashed line-shaped question mark in the last (fourth panel) schematic diagram in Figure 5e indicating an undefinable shape (Brace et al., 1982).
A histogram of the estimated widths of the boundary wave of the 89 events with a bin size of 50 km is presented in Figure 5a. The distribution of the boundary wave widths is a single peak positively skewed distribution with a median value of 173 km and a mode class of 100-150 km. Mode class is the most frequent range of values in a distribution. It ranges from 87 to 550 km and has an average width, of 212 ± 12 km. The lower 0.25 and higher 0.75 quantiles are 135 and 255 km, respectively. Additional analysis conducted with criterion vex < 30 ∘ yields only 41 events but results in a similar average width, of 219 ± 17 km.
The diameters of the flux ropes are also estimated in a similar fashion but without the angle criteria mentioned earlier. Both of these criteria can be measured by the magnetic field magnitude along B min of a flux rope crossing. For example, for a spacecraft to cross the exact center of a flux rope, the magnetic field magnitude along B min should be 0. And for the case where the flux rope is not crossed through its center, there should be a finite non zero magnetic field magnitude along B min . To measure the significance of the field magnitude along B min , the |B min |∕|B| value of all the flux ropes is calculated. In fact, the average value of |B min |∕|B| for all the 133 flux ropes is only 0.12 ± 0.01. This indicates that all of the flux ropes observed were crossed at or very close to their center.  The criteria are mv < 30 ∘ or > 150 ∘ and vex < 45 ∘ . The red arrows represent the normal directions of boundary crossings projected from its local locations, which are denoted by the orange colored dots. The ionospheric boundary is represented by green lines. The vex is the angle between the spacecraft velocity vector, and the boundary normal vector. The mv is the angle between the boundary normal directions from two consecutive crossings.
The estimated diameters of a total of 64 (out of 133) flux ropes, which have int ∕ min ≥ 3 and number of data points ≥ 6 are presented in a histogram in Figure 5b with a bin size of 25 km. The distribution of the flux rope diameters is also a single peak positively skewed distribution with a median value of 83 km and a mode class of 500-75 km. It ranges from 48 to 295 km and has an average diameter of 90 ± 6 km. The lower 0.25 and higher 0.75 quantiles are 59 and 102 km, respectively. Additional analysis conducted with criterion vex < 45 ∘ yields a total of 25 flux ropes and an average diameter of 79 ± 6 km. The low value of |B min |∕|B| and the consistency 10.1029/2018JA025644 shown in the estimated diameters regardless if the criterion on vex is applied imply that the majority of the flux ropes are crossed quasi-radially and they do not appear to be ideally circular as depicted in Figure 4e.
Note that if a boundary wave is configured similar to the three phases of the KHI evolution (Amerstorfer et al., 2010), there are many possible ways that VEX can traverse through a boundary wave. Hence, with the limitations of the instruments on board of VEX, which is also a single spacecraft mission, estimations of the exact shape and thus the size of the observed boundary wave and flux ropes are not possible.

Velocity Shear Profile
In addition to the estimation of boundary wave widths, the difference of velocity across the ionopause when boundary waves are observed, |U iono − U Mb |, is also calculated. U iono and U Mb are the average proton velocity in the ionosphere and magnetic barrier regions, respectively. Due to the long IMA sampling time of 192 s, to ensure there are at least two data points, only U iono from orbits that VEX spent longer than 192 s in the ionosphere are considered and U Mb is estimated by taking an average proton velocity 5 min before crossing the inbound ionopause. A histogram of the measured velocity difference of 158 events with a bin size of 25 km/s is presented in Figure 5c. The distribution of the velocity difference is a single peak positively skewed distribution with a median value of 87 km/s and a mode class of 25-75 km. The |U iono − U Mb | ranges from 5 to 341 km/s and has an average velocity of 104 ± 6 km/s. Since the time spent by VEX in the magnetic barrier region is often short, much less than the 192-s-long sampling time of the IMA, the estimation of U Mb by averaging the proton velocity 5 min before crossing the inbound ionopause will often include the plasma populations in the sheath region. These have a much larger magnitude along the X VSO axis due to the main solar wind bulk flow. Therefore, the U Mb presented here is an overestimate and shows bias along the X VSO axis. The magnitude of the actual U Mb is expected to be lower. Hence, the velocity difference profile presented here should be examined with caution.
In general, the results of the estimated ionospheric boundary wave widths and the velocity difference flow across the ionospheric boundary are considerably comparable with the study of KHI wave on Venus. For instance, Wolff et al. (1980) show that with a gyroviscosity coefficient, v L of 250 km 2 /s, a velocity difference of 100 km/s results in a wavelength of ∼ 31 km, while velocity difference of 10 km/s results in a wavelength of ∼305 km. For a typical 30 km-thin ionopause, Elphic and Ershkovich (1984) show that velocity difference of 100 and 200 km/s results in wave growth times of 81 and 32 s, respectively. Considering a density jump (ionosphere to the sheath region) with a ratio of 10, the local magnetohydrodynamic simulation by Amerstorfer et al. (2010) with the lower (30 km) and upper (80 km) limits of ionopause boundary thickness (Elphic et al., 1981) gives a dominant KHI wavelength with a lower range limit of ∼181 km and a higher range limit of ∼483 km. Similarly in Ong and Roderick (1972), for the most dominant mode of KHI, 30 and 80 km of ionopause boundary thickness gives a dominant KHI wavelength with of ∼224 and ∼600 km, respectively. The comparable results shown between the estimated values in this statistical analysis and the simulation results suggest that KHI may act as an excitation seed in inducing the ionospheric boundary wave that is observed in the northern polar region of Venus. However, it is noteworthy to mention that the average estimated ionospheric boundary wave width of 212 ± 12 km in this work lies close to the lower limit of Amerstorfer et al. (2010) and is actually 10 km smaller than the lower limit in Ong and Roderick (1972). This slight inconsistency can be attributed to the different parameters (e.g., gyroviscosity, density gradient, and velocity difference) considered in these mentioned studies, which would result in different KHI wavelength.
Further analysis (not shown here) has been conducted to estimate the velocity difference of all orbits regardless of whether or not boundary waves are observed. The results show that the ionospheric boundary does not always exhibit a wavelike appearance when the velocity difference is large (i.e., |U iono − U Mb | > 150 km/s). This can be due to the following: (1) stabilizing terms (e.g., gravity; Elphic & Ershkovich, 1984), which are more significant and dominate; (2) resolution of VEX data is too low to observe the short boundary wavelength resulting from the large velocity difference; and (3) the ionospheric boundary may exist in a wave, but it is not traversed by VEX.

Orientation of Magnetic Field
Analysis reveals that the magnetic field orientation in the magnetic barrier region is quasi-perpendicular to the Y-Z VSO plane, which is a favorable condition for the excitation of KHI along the along the Y VSO direction. This orientation is evidenced from the bimodal shaped distribution of the magnetic field cone angles in the magnetic barrier region presented in Figure 5d, where the majority of the cone angles are <45 ∘ and >135 ∘ , that is, quasi-parallel to the X VSO axis. The distributions of the histogram are very consistent for all VEX passes regardless of whether boundary wave events are observed. This consistent quasi-parallel orientation of magnetic field lines to the X VSO axis is due to the draping pattern of the magnetic field lines in the northern polar region of Venus. This is discussed further in the next section.

Impacts and Consequences
When the KHI-induced boundary wave becomes nonlinear, the broken wave can detach and form flux ropes, which transfers the shocked solar wind plasma in the ionosphere. The reconnection of the consecutive two troughs or two crests (i.e., ∕2) of a turbulent boundary wave would ideally result in the production of flux ropes with diameters comparable to the half-widths of the boundary wave. Therefore, the results of (1) the similarity in the normal directions between the flux ropes and the boundary wave (they mainly lie along the Y VSO axis) and (2) the similarity between the estimated flux rope diameters and the half-widths ( ∕2) of ionospheric boundary waves suggest that the flux ropes observed are likely to be formed as a result of detached ionospheric boundary waves.
At the same time, atmospheric bubbles that contain plasma with ionospheric origin are also expected to form in a similar fashion to the production of flux ropes. Ideally, atmospheric bubbles can be identified by the observation of ionospheric photoelectron populations outside the ionosphere region. However, the convection of the atmospheric bubbles out of the ionosphere and subsequently downstream with the main solar wind bulk flow can change the magnetization state as well as the characteristic energy signature of the ionospheric photoelectron populations (at 21-24 and at 27 eV; Coates et al., 2008;Cui et al., 2011). In addition, the magnetic barrier region, where atmospheric bubbles populate before they are convected downstream, is highly dynamic. These complications may result in the atmospheric bubbles not being accurately identified in this work. On the other hand, the high-latitude elliptical polar orbit of VEX indicates little opportunity for the spacecraft to encounter the atmospheric bubbles that often populate in the downstream region (−X VSO ) as a result of the convection of solar wind bulk flow.

Wave Propagation along Y VSO : Draping Pattern of Magnetic Field Lines
To assess if there is a preference in plasma velocity direction, the contribution of the average plasma velocity components along the X VSO , Y VSO , and Z VSO axes of the inbound solar wind (U sw ), magnetic barrier (U Mb ), and ionosphere (U iono ) is calculated. They are (91.3, 6.7, and 2.1)% for the inbound solar wind, (75.2, 11.7, and 13.2)% for the magnetic barrier, and (15.1, 45.7, and 39.2)% for the ionosphere. The values are normalized and are measured by taking the ratio of the individual components with respect to their overall magnitudes. They are expressed in percentage, and the median values of their respective distributions are utilized for this calculations. To eliminate the possible orbital dependence of velocity, only quasi-terminator VEX trajectories are considered, that is, when the orbital plane is <30 ∘ to the terminator plane. The U sw , U Mb , and U iono are averaged from a total of 658, 341, and 465 plasma velocity vectors, respectively. Similar to section 4.1.2, only U iono from orbits for which VEX spent longer than 192 s in the ionosphere are considered, while U sw and U Mb are estimated by taking an average of proton velocity 30 min before crossing the inbound bow shock and 5 min before crossing the inbound ionopause, respectively. The results show that the average solar wind plasma velocity outside the bow shock is dominantly along the X VSO axis, while the components along the Y VSO and Z VSO axes are minimal, as expected from the main solar wind bulk flow. However, through the bow shock and into the ionosphere, the contributions of the velocity component along the X VSO axis can be observed to decrease drastically, while the contributions of the velocity components along the Y VSO and Z VSO axes become more significant. In the ionosphere, the plasma velocity is actually dominated by the components along the Y VSO axis (with a contribution of ≈46%) followed by Z VSO (≈39%) and X VSO (≈15%) axes. These results are consistent with Lundin et al.(2011Lundin et al.( , 2013Lundin et al.( , 2014 in which the authors attributed the persistent +Y VSO directed ion flow over the northern polar region to the solar wind aberration.
The domination of the plasma velocity along the Y VSO axis in the ionosphere region further suggests that a velocity shear profile across the ionospheric boundary could be set up and consequently exciting the boundary wave by means of KHI. However, the follow-up question is, what is the possible driving mechanism of the plasma along the Y VSO axis?

Dependence of the Boundary Normals With Respect to their Locations
Next the dependence of the normal directions of the boundary crossings with respect to their locations is studied. The MV x , MV y , and MV z , which are the components of the normalized boundary normal directions along the X VSO , Y VSO , and Z VSO axes relative to the observation location of their associated boundary crossings in the X-Y VSO plane are shown in Figures 6a-6c. The measurement for each normal direction is binned in 0.075 × 0.075 R v bins ranging from −1 R v to 1 R v . Only bins with two or more measurements are shown.
The results presented in Figures 6a-6c show a clear dependence of the boundary normal directions on their respective locations, especially when the azimuthal angle of the locations, Location < |45 ∘ |. These regions are outlined with the black colored dashed line which is just in front of X = 0. For instance, observable from Figure 6a, the MV x components in the region Location < |45 ∘ | are much larger compared to the region Location > |45 ∘ |. This is visible by the stronger red colored patches in the region outlined with the black colored dashed line. In contrast, from Figure 6b, the MV y components in the same region are much smaller, visible from the weaker blue colored patches, while MV z components are rather randomly distributed. Even though the majority of the boundary normal directions are dominantly along the Y VSO axis, the above results suggest that, at the location Location < |45 ∘ |, the boundary normal directions have slightly larger components along the X VSO axis.

Draping of Magnetic Field Lines
In Figure 6d, taking the normal directions of the boundary crossings (represented as purple arrows) across an altitude range from 360 to 370 km as an example, the boundaries (represented as red dotted lines, which are perpendicular to the normal directions) show an alignment that is comparable to a typical pattern of draped magnetic field lines (represented as blue lines) around Venus. These curved patterns are also valid at other altitude ranges (illustrations not shown here). The results from Figure 6d are consistent with Figures 6a-6c. The expected overall draping configuration of magnetic field lines across the cross section of altitude is illustrated in Figure 6e.
Upstream of the bow shock, the IMF may appear in a range of orientations. The IMF moves toward Venus in the direction of the main solar wind bulk flow (−X VSO ). If the magnetic field lines are represented as an ellipse, as the field lines drape around the planet, the eccentricity tends from infinity (upstream of bow shock) toward zero (outside of the ionopause in the polar regions of Venus). This is illustrated in Figures 6e and 6f. Typically, upstream of the bow shock, the solar wind has a higher plasma (defined as the ratio of thermal to magnetic pressure). In this case, the thermal pressure dominates and the magnetic field lines are not frozen into the plasma. The plasma still gyrates around the field lines but has an overall velocity component mainly in the direction of solar wind bulk velocity in the −X VSO direction (represented by the orange arrow in Figure 6f. In contrast, through the magnetosheath the IMF lines start to pile up and drape around Venus more orderly and regularly, with the field lines eventually extending mainly in the direction of the solar wind bulk flow regardless of the IMF orientation; McComas et al., 1986;Masunaga et al., 2011;Tanaka, 1993). At the same time, the plasma is slowly cooled down and slows toward the ionopause, where the thermal pressure is gradually converted into magnetic pressure. This consequently results in a lower where the magnetic pressure dominates and the plasma is now frozen into the magnetic field lines. As a result, just outside of the Venusian ionopause around the polar regions, as illustrated by the green arrows in Figure 6f, the now much lower plasma is forced to move toward the center line of draping (which is along the Y VSO axis) due to the movement of magnetic field lines (which is along the Y VSO axis). This could subsequently set up a velocity shear along the Y VSO axis. In addition, as the density jump is along the Z VSO axis in the northern polar region, the combination together with the velocity shear along the Y VSO axis and the quasi-parallel orientation of magnetic field lines to the X VSO axis should favor the growth of KHI-induced ionospheric boundary waves along the Y-Z VSO plane and with a propagation direction along the Y VSO axis. This proposed mechanism is consistent with the results presented in this work and is illustrated in Figure 6g. Note that the flux rope and the atmospheric bubble are only depicted simply to illustrate their possible presence. Their actual shapes are not known, and the depicted circular shapes are for illustration purposes only.
The dependence of the boundary wave events on the orientations of the IMF has also been assessed. To eliminate the effects when the IMF might have changed between the inbound and outbound bow shock, IMF data are only considered when the angle between the inbound and outbound IMF orientations are <15 ∘ . Results (not shown here) indicate that the IMF orientations for the orbits when ionospheric boundary wave events are observed and the orbits where they are not observed are similar. This implies that the occurrence of boundary wave does not show any dependence on the orientations of the IMF. Figure 6. (a-c) The distribution of MV x , MV y , and MV z , which are the components of the normalized boundary normal directions along the X VSO , Y VSO , and Z VSO projected from their locations in the X-Y VSO plane. (d) An illustration of how all the perpendicular directions of boundary normal directions (from the altitude band 360 km < altitude < 370 km) could be related to field draping. Note that the intended slight solar wind direction aberration is due to the orbital motion of Venus. (e and f ) Illustrates how the draping patterns of magnetic field lines may change as a function of altitude and along X VSO direction for the area in the black square box in (d) (adapted from Chong et al., 2017). (g) Illustration of how field draping might lead to ionospheric boundary wave and the production of flux rope and atmospheric bubble. Illustration is not drawn to scale.

Summary and Conclusions
By utilizing the MAG, ELS, and IMA data on board of VEX, we have conducted the first statistical analysis of the perturbations of the ionospheric boundary at Venus over the period 2006 to 2014. Results from the MV analysis reveals that the Venusian ionospheric boundary does not always appear to be smooth but a rippled boundary that fluctuates mainly along the Y-Z VSO plane and predominantly with boundary normals along the Y VSO axis. Further analysis shows that the estimated widths of the ionospheric boundary wave and the estimated velocity difference flow across the boundaries are consistent to the results from previous simulation studies of the KHI. This leads to the suggestion that the ionospheric boundary wave develops due to the KHI. Furthermore, our analysis suggests that the draping pattern of magnetic field lines plays a principal role in enhancing the plasma flow along the Y VSO axis and subsequently sets up a velocity shear that favors the excitation of ionospheric boundary waves.
When the Venusian ionospheric boundary is excited by the KHI and reaches a nonlinear state, the wave can break off and detach, subsequently result in the formation of atmospheric bubbles and flux ropes. If ionopause waves (Luhmann, 1990) exist all along the terminator region, the draping pattern of magnetic field lines can act to enhance the production of both atmospheric bubbles and flux ropes, particularly in both the northern and southern polar regions of Venus where the magnetic field lines are more tightly draped. This can explain why the majority of the observed flux ropes from this statistical analysis have similar widths and boundary normal directions to the ionospheric boundary wave. In addition, this scenario can also provide an explanation to the detection of atmospheric plasma outflows, which are observed mainly in the polar regions, for example, plasma clouds (Brace et al., 1982) and high-energy O + fluxes (Masunaga et al., 2011). Both studies attributed their respective observations to the limitations of the PVO orbits and the upstream IMF orientations. In addition, rippling ionopause could also potentially lead to the flapping of Venusian magnetotail (Rong et al., 2015) and subsequently the productions of magnetic plasmoids in the magnetotail region (Zhang et al., 2012). However, since the ionosphere region is only sampled in locations limited to the northern polar region due to the polar orbit of VEX, it is not possible to assess and compare the nature of the ionospheric boundary in different regions, that is, subsolar, southern polar, and equatorial flank regions.
Continuous scattering and the subsequent convection of atmospheric bubbles downstream and away from Venus over a prolonged period of time plays an important role in atmospheric loss from Venus. The statistical analysis presented here provides observational evidence of the ionospheric boundary existing in a wavelike appearance. This could have had a significant impact for the planetary evolution of Venus. However, the location of the boundary, which is a standoff boundary between the planetary atmosphere and the incoming solar wind, varies significantly depending on if the planetary body is magnetized. Hence, it is of interest for future work to measure the distributions of the boundary surface wave and assess the role of field draping on the stability of planetary boundaries, particularly for unmagnetized bodies where the boundaries are much closer to the planets. Participants from Sheffield received financial assistance from STFC Consolidated Grant ST/R000697/1. Particle and low-resolution magnetic field data are available from AMDA (http://amda.cdpp.eu/) and from the ESA Planetary System Archive (PSA) at https://www.cosmos.esa.int/ web/psa/venus-express. High-resolution MAG data are available on request from the VEX -MAG PI, TieLong Zhang. The authors would like to thank all the MAG and ASPERA-4 team members for the big effort, which led to the successful operation of the instrument and the calibration of the data. R. Frahm acknowledges contributions from NASA (contract NASW-00003). The authors are very grateful for the insightful discussions with Glyn Collinson.