Spatially high-resolved solar-wind-induced magnetic field on Venus

The current work investigates the Venusian solar-wind-induced magnetosphere at a high spatial resolution using all Venus Express (VEX) magnetic observations through an unbiased statistical method. We first evaluate the predictability of the interplanetary magnetic field (IMF) during VEX's magnetospheric transits, and then map the induced field in a cylindrical coordinate system under different IMF conditions. Our high-resolution mapping enables resolving structures on various scales, ranging from the thin ionopause and the associated electric currents to the classical global-scale draped IMF. Our mapping also resolves two recently-reported structures, a low ionospheric magnetization over the terminator and a global"looping"structure in the near magnetotail, both of which are not depicted in the classical draping configuration. In contrast to the reported IMF-independent cylindrical magnetic field of both structures, our results illustrate their IMF dependence. In both structures, the cylindrical magnetic component is stronger in the hemisphere with an upward solar wind electric field ($E^{SW}$) than in the opposite hemisphere. Under downward $E^{SW}$, the"looping"structure even breaks, which is attributable to an additional draped magnetic field structure wrapping toward $-E^{SW}$. In addition, our results suggest that these two structures are spatially not overlapping with each other. The low ionospheric structure occurs in a very narrow region, at about 87--95$^\circ$ solar zenith angle and 190--210~km altitude, implying that future simulation to reproduce the structure entails at least a spatial resolution of about 10 km. We discuss this narrow structure in terms of a Cowling channel.

the solar wind (SW), inducing electric currents and shaping comet-tail-like magnetospheres (Bertucci et al. 2011). Among these bodies, Venus does not have Mars-like crustal fields, is relatively close to the sun and surrounded by dense SW, and has been observed using artificial orbiters. Therefore, the Venusian magnetosphere is controlled sensitively by the SW, and relevant observations are relatively abundant, making the Venusian magnetosphere a natural laboratory to study the interactions between the SW and unmagnetized bodies. The interactions between SW and Venus plasma environment have been investigated observationally (e.g., Russell et al. 2006;Bertucci et al. 2011), and essential processes and structures have been reproduced in global hybrid simulations (e.g., Jarvinen et al. 2013;Jarvinen & Kallio 2014). Here, we briefly summarize the processes and structures, and for details, readers are referred to the relevant reviews (e.g., Brace & Kliore 1991;Baumjohann et al. 2010;Bertucci et al. 2011; "T between VEX polar crossing and the closest bow shock crossing (min) 2007; Dubinin et al. 2020). The SW plasma surrounding Venus is supersonic and super-Alfvenic (the sound speed and Alfven speed are significantly lower than the solar wind bulk velocity) with the frozen-in interplanetary magnetic field (IMF, B IM F ). When the SW plasma encounters Venus, the SW motional electric field (E SW ) drives electric currents in the highly conducting ionosphere and the surrounding plasma. The magnetic fields associated with the currents decelerate the SW flow, forming a bow shock, magnetosheath, and magnetic barrier or magnetic pileup region. At the bow shock, SW plasma is decelerated to subsonic. The deceleration thermalizes the plasma and excites intense waves and turbulence. The shocked SW plasma remains collisionless (collision frequencies much smaller than gyro frequencies), which populates the magnetosheath but cannot easily penetrate the magnetic barrier and is mostly deflected around the planet. The deflection gives rise to a wake with low-density plasma downstream, forming the magnetotail, where happen reconnection events and particle acceleration and loss. The boundary of the low-density plasma region is known as the induced magnetosphere boundary (IMB) or sometimes magnetopause. The induced currents distribute most intensely on a thin layer, where the ionospheric thermal pressure holds off the magnetic barrier's magnetic pressure. The layer is the bottom boundary of the magnetic barrier and is known as the ionopause. The ionopause currents are closed through currents in the magnetic barrier or on IMB. According to, e.g., Russell et al. (2007) and Brace & Kliore (1991), the magnetic barrier is referred to as the most inner part of the magnetosheath, and both border to the dayside ionosphere at the ionopause, whereas according to, e.g., Zhang et al. (2008) and Baumjohann et al. (2010), the magnetic barrier is defined as an independent region between the magnetosheath and the ionopause. The boundary between the magnetosheath and the magnetic barrier is characterized by a discontinuity of magnetic distortion and magnetosheath wave terminations, which is referred to as the dayside IMB or magnetopause (Zhang et al. 2008). Referring to Zhang et al. (2008), we sketch the structures involved in the current study in Figure 1. The magnetic field configuration of the induced magnetosphere is characterized by a draping configuration, elongating along the Sun-Venus axis (e.g., Luhmann 1986). In the plane normal to the Sun-Venus axis, the orientation of the induced magnetic field is mainly determined by E SW and the cross-flow IMF component, which is usually described in a coordinate system defined referring to E SW , known as the Venus-Sun electric field coordinate system (VSE for short, also called SW velocity and magnetic field coordinate system, e.g., Dubinin et al. 2014a;Du et al. 2013). The x-axis of the VSE coordinate system is antiparallel to the upstream SW flow velocity (v SW ), the y-axis is parallel to the IMF cross-flow component, and the z-axis is parallel to E SW , therefore also termed as the E SW axis. In VSE coordinates, the classical induced magnetic field is overall anti-symmetric between the ±y hemispheres, and symmetric between the ±z hemispheres (namely ±E SW hemisphere, e.g., Bertucci et al. 2011;He et al. 2016). The ±E SW symmetry was reported broking in both the magnetic barrier and the magnetotail: the magnetic field is stronger in the +E SW hemisphere than that in the −E SW . These features were investigated in statistical studies (Du et al. 2013;Bertucci et al. 2003Bertucci et al. , 2005Zhang et al. 1991Zhang et al. , 2008Zhang et al. , 2010 and reproduced in simulations (Brecht 1990;Brecht & Ferrante 1991;Jarvinen et al. 2013), and therefore have been well understood.
The above-mentioned magnetospheric structures are on planetary scales. Recently, two relatively small-scale magnetospheric structures were observed in the low-altitude ionosphere (below 300 km altitude, Zhang et al. 2012Zhang et al. , 2015 and in the magnetotail (Chai et al. 2016), referred to as the ±E SW asymmetrical low-ionospheric magnetization (e.g., Dubinin et al. 2014a) and "looping" magnetic field (Chai et al. 2016), as indicated by the red and pink symbols in Figure 1, respectively. The lowionospheric asymmetry is characterized by a preference of induced magnetic field pointing towards +y in both ±E SW hemispheres, namely parallel to the cross-flow IMF in the +E SW hemisphere but antiparallel in the −E SW hemisphere, which was also described, alternatively, as a dawnward preference over the Venusian north pole (Zhang et al. 2015). Various mechanisms were proposed, including giant flux ropes formed in the magnetotail (Zhang et al. 2012), Hall currents induced regionally in the low ionosphere (Dubinin et al. 2014a), and antisunward transports of low-altitude magnetic belts build-up on the dayside (Villarreal et al. 2015). The "looping" magnetic field is characterized by a planetary-scale cylindrical symmetry of cylindrical magnetic field in the near-Venus magnetotail, which was explained in terms of a global induced current system distributed on double cylindrical layers (Chai et al. 2016). Most studies on these two structures are observational, based on either manually selected cases (e.g., Dubinin et al. 2014a;Zhang et al. 2015) or observations without case selections (Chai et al. 2016). Studies with selected cases resolved the structures at relatively high spatial resolutions but were potentially subject to prior knowledge in the selection, whereas unbiased studies, without any guiding selection, presented the structures at typically low spatial resolutions. Consequently, the spatial distributions of these structures are not determinedly resolved. For example, it is not known whether the low-ionospheric asymmetry overlaps with the "looping" structure. To deal with this problem, we present unbiased statistical results of the near-Venus magnetic field at a high spatial resolution up to about 50 km. Our results suggest that, as sketched in Figure 1, the low-ionospheric asymmetry was not covered by a recent relevant magnetohydrodynamic simulation in Villarreal et al. (2015), and the low-ionospheric asymmetry is not overlapped with the "looping" field. We further compare the spatially highly resolved depictions under different IMF conditions, demonstrating that the "looping" structure is not cylindrically symmetric and breaks under some IMF conditions. Below, after explaining our data analysis in Section 2, we demonstrate the capability of our method in resolving small-scale magnetospheric structures in Section 3 and investigate the "looping" field and low-ionospheric asymmetry in Section 4.  Figure 3. (a) The correlation between 20-min averaged B IM F x at time t and that at t±30 min at Venus. The color code represents the IMF disturbance σ total 20 min defined in Section 2.1. The red points display the cases collected during the 15% most disturbed states characterized by σ total 20 min >3.7 nT. Indicated in the top left corner are the correlation coefficient r ± ∆r, and p-value, estimated from a bootstrapping analysis detailed in Section 2.1. For a better visibility, shown is only a portion of randomly selected pairs. (b and c) Same as Panel (a) but for B IM F y and B IM F z . Read Section 2.1 for details.

DATA ANALYSIS
Venus Express (VEX) operated in a polar orbit. The orbit had a periapsis at 170-350 km altitude at about 78 • N latitude in the Venus Solar Orbital coordinate system (VSO, the x-axis points from Venus towards the sun, the z-axis is normal to the Venus orbital plane, pointing northward, and the y-axis completes the right-handed set, e.g., Svedhem et al. 2007;Titov et al. 2006). Due to the high orbital eccentricity, VEX repeatedly encountered SW, magnetosheath, and magnetosphere. VEX was not accompanied by any independent spacecraft for monitoring the IMF condition in VEX's magnetospheric transit. Therefore, the IMF has to be estimated. Here, after introducing and evaluating our IMF prediction in Section 2.1, we introduce a cylindrical coordinate system, our IMF binning, and a high spatial resolution mapping approach in Sections 2.2, 2.3, and 2.4, respectively. We use all available 4-s-averaged magnetic field observations from about 3000 orbits during VEX's eight Venusian years of operation (Zhang et al. 2006). Readers are referred to He et al. (2016) for the distribution of the orbits as functions of local time, solar zenith angle, and altitude.

IMF predictability for VEX's magnetospheric transit
A popular IMF estimation approach compares the IMF conditions immediately before the inbound bow shock crossing and after the outbound crossing (e.g., Zhang et al. 2015). However, the current study focuses on the magnetosphere over the northern polar cap, which is typically much closer to one bow shock crossing, either inbound or outbound, than the other in any given orbit. This situation is similar to that illustrated by He et al. (2017) where the authors estimated IMF for MESSENGER's Mercurian magnetospheric transit using the actual IMF observation from the nearest SW encounter. Following the IMF estimation approach in He et al. (2017), we estimate the IMF conditions as the average of actual IMF observations in a remote 20-min-wide interval of VEX's nearest SW encounter at the closest bow shock crossing, either inbound or outbound. The bow shock crossings are identified manually and detailed in Rong et al. (2014). In each orbit, we first identify the VEX's maximum latitude arriving as a reference point of the magnetospheric transit, and calculate the temporal separation between the reference point and the corresponding closest bow shock crossing. The temporal separations of all orbits are distributed in Figure 2, which peaks at about 17 min with an average of 20.5 min. By average, the reference point is separated away from the center of the SW interval by ∆t= δt/2+20.5 min =30.5 min, where δt=20 min is the sampling window width. Approximately IMF is estimated as the actual IMF observations 30 min preceding or succeeding. Below, we evaluate the estimation using actual IMF observations through a correlation analysis.
For the correlation analysis, we first select the magnetic field measured in the SW and average them within discrete 20-min-wide intervals. We denote B IM F (t) as the average in the interval from t − 10min to t + 10min. Then, B IM F (t ± 30min) are 'estimations' of B IM F (t). In total, the VEX SW observations enable about N pair =6.5× 10 5 independent pairs of [B IM F (t), B IM F (t ± 30min)], the three components of which are scatter-plotted in the three Panels in Figure 3. The points are color-coded according to the IMF variability defined as σ total 20 min = σ 2 Bx + σ 2 By + σ 2 Bz . Here σ Bx , σ By , and σ Bz are the standard deviation of the three magnetic components within the 20-min sampling window centering at t ± 30min. According to ascending σ total 20 min and following He et al. (2017), we sort the N pair =6.5× 10 5 pairs into ten groups evenly, and calculate the correlation coefficient r between B IM F (t) and B IM F (t ± 30min) in each group, illustrating that r decreases with increasing σ total 20 min (not shown). For a robust analysis (following Section 2.5, He et al. 2017), we exclude the VEX orbits associated with disturbed IMF states, characterized by σ total 20 min >3.7 nT, from the following analyses. In the end, we get N o =2945 orbits for which bow shock crossings are identified when IMF states are not disturbed.
After removing the 15% most disturbed IMF pairs from the N pair =6.5× 10 5 pairs, we carry out a bootstrapping analysis as follows. First, we select randomly N o =2945 samples from the reduced set of 6.5×10 5 pairs with replacement to calculate the correlation coefficient r and p-value. The sampling is repeated for N s =5000 times, yielding N s values of r and p. Here, N s =5000> N o is selected subjectively. The average and standard deviation of r are present in each panel of Figure 3, indicating the correlations are significant for all IMF components and demonstrating that the IMF estimation is reasonable. The correlation coefficients of the three components are close to each other, suggesting that their predictability is comparable. This is different from the IMF predictability at Mercury (He et al. 2017) where different components' predictability is significantly different.

A cylindrical coordinate system for cataloging IMF conditions
Venus does not have a significant intrinsic magnetic field, neither global dipole field nor regional crustal field, which is magnetically spherically symmetrical. The interaction between SW and Venusian ionosphere induces a global magnetic field B V and a magnetosphere. The interaction is largely magnetically rotationally symmetrical with respect to the axis along v SW . Below, we describe the symmetry mathematically and use it to discuss the VSE system and introduce a new coordinate system.  In each panel, the red axis points from the Venusian center to the sun while the green axis is normal to the Venus orbital plane and points to the north. Read Section 2.2 for details.

The rotational symmetry in VSO
The large-scale topology of B V is characterized by a draped configuration, determined mainly by the upstream SW conditions, mainly its magnetic field B IM F (e.g., He et al. 2016). Mathematically, Here, B IM F is assumed to be homogeneous on the planetary scale at Venus. Therefore, it is not a function of r but a function of time B IM F (t).
Investigating the details of Equation 1 is a sort of an essential task of SW and Venus interactions.  magnetic variance near Venus (as quantified in Figure 10 in He et al. 2016), and therefore is often neglected for a first-order approximation.) Accordingly, Equation 1 is often simplified as, B V 's response to the IMF is characterized by rotational symmetry. Mathematically, for any angle ω, Here, a product of a vector with the rotation matrix R (ω) :=

In VSE
Specially, we define the IMF clock angle For an arbitrary vector a, its product with R −θ IM F defines a coordinate transformation from VSO to VSE: a V SE = R −θ IM F a V SO (see Figure 4a). Therefore, Equation 4 writes has four independent scale variables, one less than in the VSO representation in Equation 2, which can therefore facilitate data processing, analysis and interpretation. However, the VSE position vector r V SE is determined using the IMF estimation. Therefore, the r V SE determination is subject to the error in the IMF estimation and r V SE cannot be determined at a singular point B IM F ⊥ ≈ 0. To conquer these difficulties, in the following subsection, we develop a cylindrical coordinate system.

A cylindrical coordinate system
The previous subsection uses the IMF clock angle θ IM F to define the rotation matrix R −θ IM F . We can also define R −θ V EX using the clock angle θ V EX := arctan r V EX ] is VEX's position vector. Then, the product R −θ V EX a V SO defines a coordinate transformation of the vector a V SO from VSO to a cylindrical coordinate system along the VSO x-axis. The three components of R −θ V EX a V SO correspond, in order, to the axial, azimuth, and radial directions, which could be arranged into the conventional order (ρ, φ, and ζ) of cylindrical coordinate systems by multiplying with a permutation matrix U =: Figures   4b and 4c, the ζ-axis is the cylindrical axis pointing from the sun to Venus, the ρ-axis is the radial direction parallel to the position vector r V EX , and the φ-axis is azimuth coordinate normal to the ρ-ζ plane. This coordinate system is called hereafter as Solar-Venus-cylindrical (SVC) coordinate system. The product a SV C : ≡ 0 denotes that the VEX is on the plane determined by the sun, Venus, and VEX.) Therefore, all VEX observations distribute two-dimensionally on the ρ-ζ plane in the SVC system. Substitute Distributing all on the ρ-ζ plane, the VEX observations allow only investigating the induced field B V only on the ρ-ζ plane. Accordingly, Equation 6 writes Similar to Equation 5 in VSE, Equation 7 also comprises four independent variables and describes the response of B V to B IM F ρ and B IM F φ on the ρ-ζ plane. The SVC position vector is not subject to the IMF estimation and therefore is capable of handling the VSE singular point at B IM F ⊥ ≈ 0. In the following subsections, we construct B V (r ρ , r ζ ) under five different B IM F ⊥ SV C conditions.

Five IMF conditions
Figure 5 displays the IMF estimations for all VEX magnetospheric transits in SVC ρ-ζ plane. Each point corresponds to one vector of 4-s-resolved B V observed below the altitude h <4400 km and within the range of solar zenith angle (SZA) between 65 • and 150 • , which is the targeting region of the current work. The arc-shaped traces reflect cylindrical rotations of VEX with respect to B IM F ⊥ . One arc represents once targeting region transit. The portion colored in cyan corresponds to the 15% lowest magnitude B IM F ⊥ , while the rest portions, in magenta, blue, green, and red denote four 60 • -wide ranges of the clock angle θ IM F . The five colors represent approximately the five SW/IMF conditions, namely, Under each of these conditions, we combine the VEX observations to construct B V in ρ-ζ depiction, using the method detailed in the following subsection.

An non-uniform high spatial resolution method for mapping B V (ρ, ζ)
We distribute the B V observations under the 0B IM F ⊥ condition, colored in cyan in Figure 5, equally into fourteen altitude levels, then into fourteen SZA bins at each altitude level, resulting in 14×14 spatial bins with approximately identical sampling count, e.g., N sample =478-479 in all bins. In each bin, the median of the position vector (r ρ , r ζ ) of the 478-479 samples is shown as a point in the ρ-ζ plane in Figure 7a and is color-coded with the median of the component B ρ . On each point, the black cross is the error bar that represents the upper and lower quartiles of (r ρ , r ζ ) of the 478-479 samples. Similarly, B V φ and B V ζ components, and the magnitude B V are also constructed in Figures 7b, 7c, and 7d, respectively. For other four IMF conditions, B V (r ρ , r ζ ) are also constructed and displayed in Figures 7e-h, 7i-l, 7m-p, and 7q-t. In addition, to check the details, we zoom and project B V φ and B V , namely, Figures 7b, 7d, 7f, 7h, 7j, 7l, 7n, 7p, 7r, and 7t, into the SZA-altitude depiction in Figure 8.
The IMF dependence of Venusian magnetospheric B V distribution has been observationally investigated in both cases (e.g., Zhang et al. 2012) and statistical studies (e.g., Dubinin et al. 2014a;Zhang et al. 2015), which are mostly based on manually-selected observations and therefore are subject to potential prior expectations. Some studies have also constructed unbiased-statistical B V under different SW/IMF conditions (Chai et al. 2016;Zhang et al. 2010), using data without any guiding selection. These unbiased studies focused mainly on large structures at planetary scales, because they used a popular data averaging approach. The approach divides the VEX samplings spatially uniformly into cubes, and averages data in each cube separately (e.g., Du et al. 2013;Zhang et al. 2010). The resultant maps have uniform spatial resolutions but non-uniform standard error of the mean (SEM) that is proportional to 1/ N sample . Here, N sample counts observations in a bin, which is highly spatially non-uniform ( Figure 6a, also cf, Du et al. 2013;Chai et al. 2016). In the current study, we map VEX magnetic observations with uniform SEM by averaging the same number of samplings N sample in each bin, which yields a uniform significance with non-uniform spatial resolution, different from the uniformly spatially resolved works (e.g., Du et al. 2013).
Our resolution is about 1-1.5 • in SZA and up to about 50 km in altitude over 90 • SZA at 200 km altitude, at least one order of magnitude higher than the typical resolution of results on the uniform spacing grid (e.g., about 600×600 km in Du et al. 2013;Chai et al. 2016).     Figure 7d). In Panels in the second column, the arrow, cross, and dot symbols at SZA=75 • denote the orientation of the ionopause current.

CLASSICAL IMF DRAPING CONFIGURATION AT A HIGH SPATIAL RESOLUTION
The current section illustrates that our mapping approach can resolve structures at various spatial scales, including the planetary-scale IMF draping configuration and the thin ionopause. For convenience, hereafter, we denote the induced B V on the in the ρ-ζ plane as B{ The anti-symmetries in B ρ and B ζ reveal the draped IMF: the magnetic field is bending from ±ρ direction in upstream SW to ±ζ direction in the magnetosphere (for a three-dimensional schematic draping pattern, readers are referred to He et al. 2016).
Although above the three magnetic components exhibit dependence on SW/IMF conditions, B exhibits similarities under different SW/IMF conditions. The first similarity is the distinctive daynight difference: on the dayside, e.g., at SZA=70 • , B is stronger than that on the nightside at SZA> 90 • , and it maximizes vertically at about 400-800 km altitude (see Figures 8b, 8d, 8f, 8h, and 8j, zoomed versions of Figures 7d, 7h, 7l, 7p, and 7t) which corresponds to the magnetic barrier. For example, under +E SW in Figures 7h and 8d, B in the magnetic barrier is a few times stronger than the IMF magnitude: B barrier > 40 nT vs. B IM F = 7.3 nT. Since the B in Figures 7  and 8 share similarities between different IMF conditions, in Figure 9a we depict B as a function of h using all VEX's observations at SZA=70-75 • without considering IMF conditions.

Venus
-j ζ (a) viewing from sun cross-tail current sheet (b) viewing from +y Figure 10. Sketches of the ionopause current in the planes (a) r x =Rv/4, and (b) r y =0 in the VSE coordinate system, viewing from the sun and the +y flank, respectively.
In Figure 9, B maximizes at 30-35 nT between 400 and 900 km altitude. In this altitude range, the probability density of B exhibits two maxima, at 25-40 and 0-6 nT, as illustrated in Figure 9d. The 25-40-nT maximum corresponds to the magnetic barrier, whereas the 0-6-nT maximum reveals the unmagnetized ionosphere which exists below the ionopause. We display the probability p of B <6 nT as a function of altitude in Figure 9c to use the weak B as an indicator of the unmagnetized ionosphere. At 400-900 km, for about p <10%, the ionosphere is unmagnetized. (Note that this indicator is imperfect. Otherwise, (1) p should maximize at 100% at a low altitude and decrease monotonously with h, and (2) its vertical gradient should be proportional to the probability density of in-situ occurrences of ionopause.) In Figure 9a, below and above the 400-900 km region, B decreases with decreasing and increasing h, respectively. The vertical gradients are signatures of meridional currents flowing oppositely on the ionopause and IMB, as sketched in Figure 19 in Baumjohann et al. (2010), which is quantified in Figure 9b. In Figure 9b, the B gradient maximizes at h=320 km, with a full width at half maximum (FWHM) of about 90 km, extending from 280 to 370 km. Note that the gradient peak at h=280-370 km describes the altitude of the maximum probability density of ionopause occurrences, rather than the complete occurring range or thickness of the ionopause. The ionopause at every instant is a thin layer without a vertical extension. The altitude h=280-370 km is above the exobase (150-200 km, see Table 2 in Hodges 2000), suggesting that during VEX's operation, the dayside (SZA<70-75 • ) bottom ionosphere is unmagnetized for most of the time. The low-ionospheric magnetization reported in VEX observations (e.g., Dubinin et al. 2014a;Zhang et al. 2012) occurs mostly around the terminator (as detailed below in Section 4.2), since most VEX's low-ionospheric magnetic observations are distributed not far from the terminator due to the constraints of VEX's obit as displayed in Figure 6 and also pointed out by Dubinin et al. (2014a). At h < 200 km, 75.7% of the observations are distributed at 85 • ≤SZA≤95 • (Figure 6b). The B similarities among Figures 7d, 7h, 7l, 7p, and 7t are associated with different magnetic components under different SW/IMF conditions. Under ±E SW (among Figures 7e-g and 7i-k) the gradient is sharper in B φ than that in the other two components, whereas under ±B IM F ⊥ (among Figures 7m-o and 7q-s) the gradient in B ζ is the sharpest. The gradients of these components suggest that the dominant components of the ionopause current in the ±E SW and ±B IM F ⊥ hemispheres are in the directions of ±ζ-and ±φ-axis, respectively, as displayed by the arrows, cross, and point in Figures 8d, 8f, 8h, and 8j. Figures 10a and 10b sketch the topology of the ionopause current in VSE coordinates, in the plane r x =Rv/4 and r y =0, respectively. In the current work, Rv=6052 km denotes Venus radius. The currents exhibit symmetry between ±B IM F ⊥ hemispheres and anti-symmetry between ±E SW hemispheres. On the dayside, ±j φ in the ±B IM F ⊥ hemispheres closes +j ζ in the +E SW hemisphere with −j ζ in the −E SW hemisphere. Over the poles of the ±E SW hemispheres, ±j ζ flows across the terminator, and then closes through the SW or through currents on the IMB or in the magnetic barrier. Further, the vertical B φ gradient in, e.g., Figures 7f and 8c, allows estimating the ionopause current density: the drop by more than 20 nT from the magnetic barrier maximum to low altitude is associated with a current at an intensity more than 15 A/km under the assumption of sheet-like current distribution.

±E SW asymmetries
Although Section 3.1 illustrates the similarity of B between the ±E SW hemispheres, the B intensity is stronger in +E SW hemisphere than that in the other. Referring B =25 nT as a boundary, the magnetic barrier extends to about 85 • SZA under the +E SW condition (Figure 8d) but about 75 • SZA under the −E SW (Figure 8f). To quantify the ±E SW difference, in Figure 11a we display B{+E SW } − B{−E SW } , namely, the difference between Figures 7h and 7l. The main contributing component to the asymmetric B is B φ which exhibits ±E SW anti-symmetric polarity (Figures 7h and 7l) as discussed in Section 3.1. To quantify the break of anti-symmetry between Figures 7h and 7l, we Figure 11b. The ±E SW difference reflects the asymmetry of the IMF wrapping between the ±E SW poles, as reported in Zhang et al. (2010).
In Figure 11a, the B difference exhibits three peaks, on the dayside in the magnetic barrier, on the nightside in the very near magnetotail, and exactly over the terminator at low altitude (in the ionosphere, see the contour line of 13 nT). The barrier and terminator B peaks are associated with peaks in the B φ difference in Figure 11b (above 13 nT, see the contour lines). The barrier and magnetotail peaks extend into global scales, whereas the terminator one is restricted quite regionally. The global-scale asymmetries are known as magnetic ±E SW asymmetry (investigated in both observations and simulations, Saunders & Russell 1986;Zhang et al. 1991;Jarvinen et al. 2013), which also exists at Mars (Vennerstrom et al. 2003). The underlying mechanism is still under debate. Finite Larmor radii effects of pickup ions were suggested to account for the asymmetry (Phillips et al. 1987), but hybrid simulations (Brecht 1990) could reproduce the phenomenon without planetary ions, and indicated that the Hall effect contributes significantly to the asymmetry. In addition, global MHD simulations with multifluid treatment at Mars (Najib et al. 2011) reproduced the ±E SW asymmetry, suggesting that the decoupling of separate ion fluids could serve as an alternative explanation.
Compared with the global-scale asymmetries, the regional one over the terminator is relatively less understood and will be discussed below in Section 4.2.

NEW LIGHTS ON THE INDUCED FIELD
The previous section demonstrates the capability of our mapping approach in resolving small-scale strictures. In the current section, we investigate two structures that have not been depicted in the classical draping configuration.

The "looping" magnetic field
A cylindrical symmetrical preference of positive B φ was reported on a cylinder in the magnetotail, which was called the induced global magnetic field "looping" (Chai et al. 2016   structure is denoted by the pink region in Figures 1 and 12a. For a comparison with Chai et al. (2016), in Figure 11c we average B φ under the SW/IMF conditions ±B IM F ⊥ and ±E SW (Figures  7f, 7j, 7n, and 7r). As a reference, the black box in Figure 11c sketches the distribution of the B φ symmetry according to Chai et al. (2016). In Figure 11c, the tail-like B φ maximum on the nightside is largely consistent with the black box and wrapping slightly towards low latitude.
The tail-like B φ maximum exhibits a quantitative difference between the different IMF conditions (among Figures 7f, 7j, 7n, and 7r). The B φ maximum is strongest under +E SW (Figure 7f) and weaker under ±B IM F ⊥ (Figures 7n and 7r). Under −E SW (in the black rectangle in Figure 7j), B φ even partially reverses to negative. This IMF dependence suggests a cylindrical asymmetry, as sketched in Figure 12b. Chai et al. (2016) explained the symmetry of +B φ preference in terms of currents on two cylinders with the identical longitudinal axis pointing towards ζ and −ζ (as sketched in Figure 4b in Chai et al. 2016, and here in Figure 12a by the red symbols). The breaking of the symmetry of +B φ preference under the −E SW condition reported here suggests that the current system is not cylindrically complete in the near tail. A correction of the "looping" structure and two-cylinder currents is sketched in Figure 12b.
This corrected "looping" structure could be explained as the superposition of the classical magnetic draping structure with an additional draping configuration. This superposition theory was proposed for explaining a similar structure on Mars (Dubinin et al. 2019). Similar to Venus, Mars does neither have an internal dipole magnetic field, and its interaction with solar winds produces a similar "looping" structure (Chai et al. 2019). Dubinin et al. (2019) attributed this "looping" structure to an asymmetrical pileup of IMF and an associated additional draping configuration in which magnetic fields bend toward the −E SW direction in the planes normal to the SW velocity. The authors also reproduced this structure in hybrid simulations. The additional draping can produce closed loops over the −E SW hemisphere and weaken the magnetic field and potentially also magnetic reconnection (also see Rong et al. 2014;Ramstad et al. 2020). Another morphological character of this "looping" in Figure 11c is that it does not overlap with, but is isolated from, the low-ionospheric preference of B φ over the terminator as detailed in the following subsection.
4.2. Ionospheric ±E SW asymmetry over the terminator Section 3.2 mentioned a ±E SW asymmetry over the terminator (highlighted in the black circle in Figure 11b). In this small region, the polarity of the induced magnetic component B φ is unresponsive to the upstream B IM F φ polarity and prefers to be positive (see the circle in Figure 7j and the green isoline in Figure 8e). Otherwise, in all other regions, the B φ polarity is consistent with the IMF polarity of B IM F φ . The asymmetry has been investigated using selected cases (e.g., Dubinin et al. 2014a; Zhang et al. 2015) and was termed as the ±E SW asymmetrical response of low-ionospheric magnetization to B IM F ⊥ polarity (e.g., Dubinin et al. 2014a). Here, Figures 11b, 8c and 8e present the unbiased statistic at high spatial resolution.
Not surprisingly, our unbiased statistical results are significantly different from those based on selected cases. For example, the maximum B φ is about 13 nT at +E SW in Figure 8c and about 4 nT at −E SW in Figure 8e, which are significantly lower than the averaged magnetic intensity, 45 nT, of the 77 selected events in Zhang et al. (2015). This significant difference between +E SW and −E SW is also inconsistent with the conclusion in Zhang et al. (2015) that the asymmetry does not show a preference for any particular IMF orientation. Further, our results illustrate that the asymmetry occurs in a very narrow region, about solar zenith angle 85 • <SZA< 95 • and altitude h <230 km, referring to B φ >4 nT in Figure 8e. To specify further the extensions in h and SZA, we construct the one-dimensional distribution of B φ as a function of h and SZA in Figure 13 under ±E SW conditions. The vertical dashed line in each panel displays the half-maximum of B φ , referring to which B φ peaks at about h=185-210 km, h=185-250 km, SZA=88-95 • and SZA=88-93 • in Figures 13a, 13b, 13c, and 13d, respectively. (Note that B φ maxima will be smeared out if using broader h or SZA windows. Therefore, the estimation of the B φ peak width is subject to the sampling window widths.) Among the identifications of these four peaks, the peak identification in Figure 13b is most vulnerable to disturbances of the half-maximum of B φ . For example, if the half-maximum of B φ in Figure 13b increases by 2 nT, the peak width will shrink from h=185-250 km to 200-230 km. Therefore, below we will describe the altitude variations referring mainly to Figure 13a. The regional distribution implicates that simulations to reproduce the structure entail a high spatial resolution, at about 1 • in SZA and 10 km in altitude. Note that this SZA range was not included in the simulation in Figure  7b by Villarreal et al. (2015) where the simulation presents the magnetic distribution in the VSO x-y plane at z=1 Rv, as sketched by the yellow line in Figure 1 and the cyan line in 7h and Figures 11a.  Figure 15. The expectation of the magnetic component B φ induced by (a) j ζ (r) = δ(r) in unit of Ampere, namely, a unit line current along ζ direction centering at the origin of coordinate r = 0, and (b) j ζ (r) = j C δ(r − r j ), namely, a unit line current along ζ direction centering at a random position r j . Here, position vectors r and r j are two-dimensional defined in the r φ -r ρ plane, δ denotes a Dirac delta function, and r j ∼ N 2 (r 0 , d 2 ) denotes a two-dimensional Gaussian random variable with a mean of r 0 := [r 0φ , r 0ρ ] and variance of d 2 as illustrated in (c). (c) The probability distribution of a two-dimensional Gaussian random variable r j ∼ N 2 (0, d 2 ). (d) A least-square fit of VEX observations to the model in (b). In (b) the green and white contours display the same distribution as in (c) but centering at r 0 . In (d), the magenta line presents the B φ along the magenta line in (b) when j C =191.1 A, h 0 =179 km, and d=10 km, which are fitted using values of the black crosses, namely, the medians displayed as the black dots in Figure 13a. Read Section 4.3.1 for details.

The Cowling channel at Venus
Several efforts were made to interpret the asymmetry detailed in Section 4.2, e.g., in terms of the giant flux rope, magnetic belt, intrinsic field, or induced field, and antisunward transports of low-altitude magnetic belts (Zhang et al. 2012;Dubinin et al. 2014a;Villarreal et al. 2015). The potential interpretations were discussed qualitatively in Dubinin et al. (2014a) and Dubinin et al. (2014b), which suggested the most promising explanation is a Cowling current. As sketched in Figure A.1 and detailed in Appendix A, a Cowling channel entails a particular electromagnetic configuration, including a narrow band of high Hall conductivity σ H elongating perpendicular to an ambient magnetic field B 0 , and a primary electric field E 0 , perpendicular both to the gradient of the σ H band and to B 0 . A Cowling channel might occur over the poles of ±E SW hemispheres with the orientation as sketched in Figure 14. In the daytime Venusian ionosphere, the Hall conductivity σ H maximizes vertically nearby the maximum vertical peak of the ionospheric electron density (V2 layer, cf, Dubinin et al. 2014a;Pätzold et al. 2007). The peak is thin vertically and extends broadly in horizontal directions, which is typically under ionopause and therefore is unmagnetized. However, as suggested by Figures 7 and 8, the ionopause does not cross the terminator but below about 85 • SZA at 200-400 km altitude. Therefore, SW magnetic field B IM F ⊥ might penetrate into the ionosphere over the ±E SW poles (Dubinin et al. 2014a). Over the poles, the V2 layer is normal to ρ-axis and σ H gradient, while the B IM F ⊥ is parallel to the φ-axis and might serve as B 0 . Under this configuration, an antisunward primary electric field E 0ζ might induce a Cowling current along the ζ-axis, j Cζ . The magnetic field associated with j Cζ might account for the low-ionospheric ±E SW asymmetry. However, few issues about the Cowling current are still open. Here we discuss the specific distribution of the current and the primary field E 0ζ . 4.3.1. The distribution of the Cowling current Dubinin et al. (2014a) have not specified how does the Cowling current j Cζ distribute. The simplest distribution model is an infinite straight line current with a constant intensity. A line current would induce a magnetic field inversely proportional to the distance to the line, and therefore, the induced magnetic intensity would monotonically decrease with increasing distance. The associated cylindrical component B φ would also decrease monotonically with increasing altitude above the current, as illustrated by the magenta line in Figure 15a. Inconsistent with this monotonic altitude dependence in the current model, the B φ component of the VEX observation (Figures 13a-b) maximizes at about 200 km. To address this inconsistency, we propose a minor correction to the infinite line current model.
We assume that the line current does not occur at a fixed position but at a random position, e.g., r j ∼ N 2 (r 0 , d 2 ), a two-dimensional Gaussian random variable with a mean of r 0 and variance of d 2 . The probability density function of r j for the particular case of r 0 :=[r φ0 ,r ρ0 ]=[0 0] is displayed in Figure 15c. A line current with constant intensity j C at r j will induce a distribution of magnetic B φ as depicted in Figure 15b. Above the current, B φ exhibits a peak, similar to that in Figure 13a. Further, we use values at the black points in Figure 13a for a least-square fit to the B φ along the magenta line in Figure 15b, resulting in j C =191 A, d=10 km, and r ρ0 =1Rv+179 km. The fitted results, as displayed in Figure 15d, suggest that a current locating at altitude 179±10km with an intensity of 191 A could account for the B φ variation in the VEX observation (Figure 13a). Dubinin et al. (2014a) conjectured that the primary electric field E 0ζ might arise from downward drifting electrons. In the Cowling channel, electrons are supposed to be decoupled from ions due to their different ratios of the gyro-frequency over the collision frequency (κ i << 1 but κ e >> 1, see Appendix A.1). The downward electron drift is supposed to be at the order of magnitude of 100 m/s and being driven by SW motional electric field E SW mapping along SW magnetic field lines. However, this mechanism still misses details. For example, the downward electron drift is associated with an upward Pedersen current j P 0 , accumulating potentially polarization charges on the upper and lower ionospheric boundaries. The polarization charges and the induced polarization electric field would superimpose to E SW , resist the downward electron drift and terminate E 0ζ . Another detail is that the electron drift mechanism cannot explain why the asymmetry occurs at a very narrow SZA range over exactly the terminator since the Cowling channel and electron drift might occur in a broader SZA range. An alternative explanation of E 0ζ is a horizontal polarization electric field. Over the terminator, maximizes the solar radiation gradient in ζ direction, and therefore maximize the gradients of the ionospheric density ∇ ζ n e and conductivity ∇ ζ σ P . It is also over the terminator in the ±E SW hemisphere where the ionopause current merges the cross-tail current. By preventing a divergence of currents, polarization charges and electric field E ζ might be created and serve as the primary electric field E 0ζ in the Cowling channel. Note that the polarization mechanism might be equivalent to the electron drift mechanism but might describe the equilibrium from a different perspective. Evaluating mechanisms and quantitatively explaining the equilibrium entail simulations at a spatial resolution higher than the order of magnitude 10 km due to the narrow vertical extension.

SUMMARY
The current work uses all VEX magnetic observations to map the near-Venus induced magnetic field under different IMF conditions, at a high spatial resolution through an unbiased statistical method. Using VEX SW observations, we predict the IMF for VEX's magnetospheric transit and evaluate the prediction. For the mapping, we define the Solar-Venus-cylindrical (SVC) coordinate system. Similar to the VSE system, SVC also makes use of the cylindrically symmetrical response of the induced magnetic field to IMF, but is superior to VSE because the determination of SVC axes is not subject to the error of the IMF estimation and SVC can deal with the near-zero IMF condition. Our mapping resolution is spatially not uniform, maximizing in the ionosphere where VEX collects most observations. The high resolution enables depicting the small-scale structures, such as the ionopause and the associated, as well as the planetary-scale structures, such as the classical draping configuration and the magnetic pileup region. At 70-75 • solar zenith angle (SZA), the ionopause occurs with the highest probability density between 280 and 370 km altitude, above the expected exobase (150-200 km), suggesting that the bottom ionosphere at SZA<75 • is unmagnetized for most of the time during VEX's operation. Our mapping reveals that the magnetic pileup region extends to about 85 • SZA over the +E SW pole but about 75 • SZA over the −E SW pole. In addition, our mapping resolves the ±E SW asymmetry in the low-ionospheric magnetization and the planetaryscale "looping" magnetic field. The "looping" structure is cylindrically asymmetrical: much stronger under +E SW than under −E SW . Under −E SW , the "looping" breaks, which can be attributed to the presence of an additional IMF draping configuration occurring in the planes perpendicular to the SW velocity. The "looping" occurs in the magnetotail, separated from the low-ionospheric ±E SW asymmetry occurring in the narrow region, at about 88-95 • SZA and 185-220 km altitude. Our least-square fit of VEX observations to a line current model at a Gaussian random position suggests that the altitude variation of the low-ionospheric asymmetry is attributable to an antisunward line current with an intensity of 191.1 A at 179±10 km altitude. We explain this line current as in terms of a Cowling channel.
Our results also implicate that to reproduce the Venusian low-ionospheric magnetization through simulations entails an altitude resolution of about 10 km. As artificial satellites orbit extraterrestrial planets mostly with high orbital eccentricities, their observations are typically highly spatially uneven distributed. Our mapping method exploits these unevenly distributed datasets to resolve structures on various scales, with higher resolutions in regions with denser observations. Further works could implement the method on observations at Mars and Titan where similar low-ionospheric magnetization and the planetary-scale "looping" structures are expected or observed.
VEX magnetic data are available in the ESA's Planetary Science Archive. MH and ED acknowledge the supports from the Deutsche Forschungsgemeinschaft through grants HE6915/1-1 and TE664/4-1, respectively. APPENDIX A. COWLING CHANNEL A.1. Atmospheric electrical conductivity and its altitude dependence In magnetized plasma, Lorenz force drives charged particles drifting gyratorily and periodically perpendicular to the ambient magnetic field B 0 . As a result, the particles are bound to magnetic field lines. The gyratory drift might be interrupted by collisions with other particles, either charged or neutral. In an extreme collisional case, the gyratory is relatively negligible, and particles are freed from the magnetic field and move freely at the bulk velocity of the ambient particles. A measure of the relative importance between the gyratory drift and collision is the ratio of the gyrofrequency to collision frequency κ = Ω ν . The ratio κ quantifies the expectation of the number of gyratory drifting circles that one particle could achieve in-between twice collisions with other particles. The collision frequency is proportional to particle densities. Therefore, in planetary atmospheres, the collision frequency decreases exponentially with increasing altitude as the neutral density is decreasing, whereas the gyro frequency is relatively less dependent on altitude. Accordingly, κ is overall increasing with altitude. At low altitude, collisions are superior to gyratory drifts (κ << 1) whereas at high altitude, the collisions are rare and gyratory drifts are more important (κ >> 1). In-between the low and high altitudes, there is a transitional altitude h κ≈1 , where κ ≈ 1. However, the vertical increasing rate dκ dh is different for different species, steeper for electrons than for ions dκe dh > dκ i dh , yielding higher transitional altitudes for ions than for electrons h κ≈1 i > h κ≈1 e . Here, i and e in the subscripts denote ions and electrons, respectively.
As a result, the atmospheric electrical conductivity and its isotropy are altitude-dependent. Above h κ≈1 i when plasma is collision-less (κ i >> 1, κ e >> 1), in response to an applied electric field E 0 , both ions and electrons drift towards the direction of E 0 × B 0 on the macroscopic scale, and no net current is generated. Under h κ≈1 e in the collisional case (κ i << 1, κ e << 1), the applied field E 0 drives ions and electrons flowing along E 0 oppositely, generating a current j P = σ P E 0 parallel to E 0 , independent of B 0 . The current j P is known as Pedersen current and the coefficient σ P is known as Pedersen conductivity. Between h κ≈1 i and h κ≈1 e in the intermediate case ( κ i << 1, κ e >> 1), electrons are already bound to B 0 , whereas ions could still essentially move with ambient winds. The applied field E 0 decouples electrons from ions: electrons drift in the direction E 0 × B 0 , whereas the ions prefer following towards E 0 . The caused current comprises a component j P following along E 0 and a component j H = σ H E 0 parallel to E 0 × B 0 . The current j H is known as Hall current and the coefficient σ H is known as Hall conductivity.

A.2. Cowling channel: a model and two instances at Earth
In a particular configuration of conductivity and magnetic field, the electric current in the direction of E 0 might be much intenser than j P due to a superposition of j H . The geometry of configuration is characterized by a band of σ P elongating infinitely in one dimension but restricted in a perpendicular direction. Figure  through three steps. In the first step, a primary electric field E 0 , parallel to the elongating direction and within the ambient magnetic field B 0 that is normal to the plane spanned by ∇σ P and E 0 , generates Hall and Pedersen current j H and j P . Secondly, on the boundary where ∇σ P >0, polarization charges and electric field E are created by preventing a divergence of j H . In the last step, the secondary field E further generates secondary Hall and Pedersen currents j H and j P . j H flows along j P , and therefore, the total current parallel to E 0 equals to j C =j P +j H , intenser than j P . The reinforcement of the current in E 0 direction is called Cowling effect (e.g., Yoshikawa et al. 2013), the reinforced current j C is known as Cowling current (e.g., Tang et al. 2011), σ C := j C /E 0 is Cowling conductivity, and the whole current system is called the Cowling channel. Note that the horizontal separations of elements shown in Figure A.1a represent only the causative relations in the explanation but not the spatial distribution. In space, the elements coexist under an electrodynamic equilibrium, as sketched in Figure A.1b. The three steps illustrated in Figure A.1a are not a unique explanation for the equilibrium. In the end, only the equilibrium matters.
For detailed discussions, readers are referred to the literature in the context of two instances of Cowling channels in Earth's ionospheric E-layer, over the equator and auroral zone, respectively (e.g., Kelley 2009;Fujii et al. 2011, and references therein). The corresponding Cowling currents in these two instances are well known as the equatorial electrojet (EEJ, Forbes 1981, and references thereafter) and auroral electrojet (AEJ, Boström 1964, and references thereafter). The orientation of these two Cowling channels are sketched in Figures A.1c and A.1d. In both cases, the Cowling channels are zonally-elongated and the geomagnetic field serves as B 0 . In the case of EEJ, B 0 is northward, the gradient ∇σ P is vertical, and the eastward component of the electric field induced by the tidal wind serves as the primary electric field E 0 . For the boreal westward AEJ, B 0 is downward, ∇σ P is in the magnetic meridional direction, and the zonal electric field mapped from the magnetosphere serves as E 0 .