Magnetic fields and electric currents around the dayside magnetopause as inferred from data-constrained modeling

Based on a new mathematical framework and large multi-year multi-mission data sets


Introduction
The dayside magnetosheath and magnetopause play a principal role in the magnetosphere response to the interplanetary plasma flow.They serve as a main gateway where the first contact occurs between the incoming magnetized solar wind and the geomagnetic field, eventually resulting in a complex chain of magnetospheric processes.Of primary importance here is the mutual orientation of the external IMF and the internal magnetospheric field, defining the reconnection pattern at the boundary.This subject has long been at the center of many studies and extensive debates in the literature, starting from the seminal ideas of Dungey (Dungey, 1961) and followed by a multitude of works, recently summarized in reviews (Trattner et al., 2021;Fuselier et al., 2024).The reconnection geometry has been traditionally addressed in the framework of two basic concepts: the component and antiparallel merging (e.g.(Fuselier et al., 2021), and refs. therein (Qudsi et al., 2023)).A significant contribution to this area was made due to in situ measurements onboard THEMIS (Atz et al., 2022) and MMS missions (e.g.(Trattner et al., 2018;Trattner et al., 2021;Petrinec et al., 2022;Pritchard et al., 2023)).An independent insight into the problem was gained via MHD simulations, revealing in particular the importance of the Earth's dipole tilt angle ψ on the reconnection topology and X-line geometry (Eggington et al., 2020).
Several studies have been made in the past to quantitatively describe the magnetic field in the domain between the bow shock and the magnetopause.Kobel and Flückiger (Kobel and Flückiger, 1994) developed an analytical theoretical model under assumption of a purely potential field and using parabolic approximation for both boundaries.In a much later work by Romashets and Vandas (Romashets and Vandas, 2019), a similar approach was employed, also based mostly on theory.An extended study of the magnetosheath properties was carried out by Zhang et al., 2019, but limited to only a statistical description of plasma and magnetic field parameters throughout the domain, without any numerical model.An in-depth data-based investigation of the IMF-controlled magnetic draping patterns in the magnetosheath was recently performed by Michotte de Welle (Michotte de Welle et al., 2022); however, that study used a direct data-driven approach, without an explicit external input.
This paper presents first results of an effort to implemenent the empirical approach to the problem, based on a formal mathematical framework and a large set of archived space magnetometer and plasma data, collected by several satellite missions over 25 + years of in situ observations.The basic goal of our work is to extend the data-constrained modeling of the magnetosphere, recently reviewed in (Tsyganenko et al., 2021), beyond the dayside magnetopause.An initial step in that direction was described in our previous paper (Tsyganenko et al., 2023) (henceforth TSE23), in which an empirical magnetosheath magnetic field model was developed.In that study, no magnetopause per se was included, and the field inside the magnetospheric boundary was tacitly implied as a smooth inward extrapolation of that in the magnetosheath, such that no Chapman-Ferraro current was included by construction.The present work is based on the same data; a new element is an explicitly defined data-derived magnetopause and an outer magnetosphere module, which makes it possible to study the IMF effects on the magnetopause current patterns.
Figure 1 illustrates the problem geometry with two field line configurations, corresponding to parallel and perpendicular IMF orientations.The plots were obtained using a simplified theoretical model with perfectly shielded vacuum fields on both sides of the magnetopause.In this work we develop a more realistic approach, including a general representation of the magnetosheath magnetic field by toroidal and poloidal components, unconstrained by the current-free assumption and based on a large pool of data.
The paper is organized as follows: Section 2 describes the data set, Section 3 outlines the architecture of the model, starting from a brief overview of its construction logic and followed by a more detailed description of the modular structure.Section 4 presents results of the magnetopause current calculations, Section 5 discusses the model validation and pressure balance issues, and Section 6 summarizes the article.

Data
As in our previous work (TSE23), the original "grand" data base included 1-min average magnetic field and plasma data, obtained over the time period from 1995 to 2022 onboard Geotail, Cluster-1, -3, -4, Themis-A, -B, -C, -D, -E, and MMS-1 missions.The data were restricted to the region sunward from X GSW = − 10 R E (the GSW subscript indicates a system similar to the standard GSM, but with X-axis oriented antiparallel to the solar wind flow).Another limitation is that the data did not extend too far outward from a nominal bow shock, nor too deep inward from a nominal magnetopause.Specifically, only those points were selected whose geocentric distances R i fell within the range 0.7R MP ≤ R i ≤ 1.4R BS , where R BS and R MP were evaluated from the models (Lu et al., 2019;Lin et al., 2010), respectively.Each record was appended with concurrent solar wind and IMF 1-min data from the OMNI resource, whenever available.The total number of records meeting these criteria was 11,120,810; Figure 2 (from Figure 2 of TSE23) illustrates their spatial distribution in three projections.
At the next step, a more accurate selection of magnetosheath and magnetosphere data was performed by means of a procedure similar to that used in TSE23, employing an efficient method first proposed in (Jelínek et al., 2012).Namely, based on the measured magnetic field intensities B and ion densities D, two-dimensional distributions of B/B sw against D/D sw were plotted for each mission, with both B and D normalized by their concurrent values B sw and D sw in the incoming solar wind.Using such diagrams in combination with a special winding algorithm allows one to more or less clearly select data taken in the solar wind (low D/D sw and B/B sw ), magnetosheath (high D/D sw and B/B sw ), and magnetosphere (low D/D sw and high B/B sw ). Figure 3 shows four diagrams for THEMIS, MMS, Cluster, and Geotail, in which the magnetosphere and magnetosheath data areas are demarcated by polygonal boundaries, drawn with light/dark blue dashed lines, respectively.Most of the contour segments were chosen to follow the lines of constant data density around ∼100 records per square 0.1 × 0.1 bins in the B/B sw and D/D sw space, while at some locations they were diverted and closed across the "isthmus" between the two domains.
Inevitably, the adopted selection procedure is somewhat ambiguous and subjective, which is especially evident in the cases of Geotail and Cluster, where the magnetosheath data areas continuously merge with those for the solar wind and magnetosphere, as already discussed in greater detail in TSE23 (Section 4).
In order to further reduce the uncertainties and improve the data discrimination, a second filtering procedure was applied to the data subsets obtained at the first step, based on an independent pair of region identification variables: the plasma bulk speed V b and ion thermal energy W. Both of these parameters suddenly and drastically change on crossing the magnetopause, from which one may expect the same kind of a distinct data grouping into two separate areas.Like B and D in Figure 3, both V b and W were also normalized by the concurrent corresponding quantities V sw and W sw 10.3389/fspas.2024.1425165

FIGURE 1
Illustrating the problem geometry with two IMF draping plots: for 90°(left) and 0°(right) cone angles. in the solar wind, and the binning result is shown in Figure 4 for Themis (left) and MMS (right) data.
Unlike Themis and MMS, the Geotail and Cluster data were not subject to the second V b -W filtering, each for its own reason: the Geotail CPI plasma instrument provided bulk velocity data only in the solar wind, while the Cluster data were obtained mostly in the high-latitude regions with much more diffuse and unstable boundaries, such that only a single wide area could be visualized in the W/W sw -V b /V sw diagram, without any distinct separation between the magnetosheath and magnetosphere.
Table 1displaysbasicstatisticalinformationaboutthefinaldataset, obtained as a result of the above selection.The table's format is similar to that of Table 1 in TSE23, but the numbers are different: first, because the present study concentrates only on the dayside magnetosheath, such that all nightside data with X GSW < 0 have been left out and, second, the outer magnetosphere data are now included as well.Similar to what was done in TSE23, the entire fitting data pool was split into independent training (T) and validation (V) subsets.The splitting method was to divide the data records on the basis of their observation times, such that data belonging to consecutive 30day intervals were alternately placed into the T or V subset.Because of much shorter autocorrelation times of principal interplanetary drivers (e.g.(Marquette et al., 2018)), such a method guarantees that data in the T and V subsets are sufficiently independent and, at the same time, not affected by long-term solar cycle trends.As a result of all the filtering/selection procedures, two separate T and V subsets were created containing 1,515,097 and 1,544,190 records, respectively.In the training subset, the corresponding shares of magnetosheath and magnetospheric data are, respectively, 551,231 and 963,866 records, while in the validation subset they are 557,117 and 987,073.
3 Model architecture

Basic steps
The main problem with creating a joint model of the magnetosheath-magnetosphere interface region lies in drastically different magnetic geometries and field intensities on the opposite sides of the magnetopause, as well as in different timescales associated with these regions.In terms of the model architecture, this prompts to represent the magnetic field with a composite structure consisting of two separate modules and, accordingly, split the model construction into two separate tasks.
At the first step, a model of the outer dayside magnetosphere is constructed, including the derivation of best-fit magnetopause parameters and their dependence on the solar wind and IMF input.In terms of data, our focus on the dayside magnetopause suggests to assemble the modeling set in such a way that only the outermost data are selected, based on the diagrams in Figures 3, 4. Fitting the magnetospheric module to the data allows to derive such basic parameters as the magnetic field magnitude, the magnetopause standoff distance, the curvature of the boundary and its tailward flaring rate, as well as response of all the above quantities to varying solar wind and IMF conditions.
At the next step, the magnetosheath part of the model is derived, similar in its structure to that described in TSE23.More specifically, the magnetic field is represented with a sum of toroidal and poloidal parts, whose generating functions are expanded into Taylor series of the normal distance δ from the magnetopause and spherical functions of the angular coordinates θ and ϕ, shown in Figure 5 (from Figure 1 of TSE23).A more detailed formulation of the model mathematical structure is given below.

Magnetospheric module
A principal element of the magnetospheric module is its boundary, represented here with a simple axisymmetric analytical surface first proposed in (Shue et al., 1997): 10.3389/fspas.2024.1425165 Coordinates used in the model formulation.Cartesian: correspond to GSW system with X-axis antiparallel to solar wind flow.Spherical: radial distance r, cone and clock angles θ = arccos (X/r) and ϕ = arctan (Y/Z), respectively.
where R S is the magnetopause standoff distance and α = 1.4427 ln (R T /R S ) is the parameter, defining the dayside boundary curvature in terms of the ratio of its terminator radius R T to the subsolar distance R S , as well as its tailward flaring rate on the nightside, such that α = 0.5 gives a cylindrical boundary in the asymptotic limit θ → π, while α > 0.5 or α < 0.5 result in a gradually expanding or tapering magnetotail, respectively.
In principle, the magnetic field inside the magnetopause can be represented by a full-scale model, including all principal intramagnetospheric sources such as the Earth's dipole, magnetopause, tail, field-aligned, and ring currents ( (Tsyganenko et al., 2021) or (Tsyganenko, 2013) and refs. therein).This study, however, is focused on a relatively narrow region around the magnetosheathmagnetosphere interface, which suggests to use a simpler single module like that shown in Figure 1, featuring all basic properties of the dayside field near the magnetopause: Here the total magnetospheric field is a sum of the dipole field b (dip) (r, ψ) and its shielding field b (shld) (r, R S , α, ψ), represented by cylindrical harmonic expansions (Tsyganenko, 1995).The sum in Eq. 2 is multiplied by a flexible function F of principal interplanetary drivers (specified in detail in the next paragraph) and, by construction, remains fully confined within a Shue-type magnetopause for any values of the factor F and the dipole tilt angle ψ.In such a setting, the scalable size of the configuration is uniquely defined by the standoff distance R S , while its shape is controlled by the flaring parameter α.In terms of the solar wind/IMF effects, the distance R S must depend on both 1) the ram pressure P d , controlling overall self-similar compression/expansion of the magnetopause, and 2) on the IMF B z , responsible either for the dayside magnetopause erosion or, conversely, the magnetic flux pileup.By contrast, the flaring parameter α is entirely defined by the IMF B z , reflecting the magnetic flux redistribution between the dayside magnetosphere and the tail lobes.This prompts us to treat R S and α as the current ("instantaneous") magnetopause parameters and represent them by separate functions of the interplanetary drivers as follows where five unknown parameters R S 0 , ε, γ, α 0 , and Δα define the magnetopause size and shape and its dependence on P d and IMF B z (normalized for convenience by 2 nPa and 10 nT, respectively).
To further refine the internal field response to changing interplanetary conditions, the shielded dipole field Eq. 2 is allowed to change as a whole by the variable factor F, controlled by interplanetary parameters.Specifically, the magnitude factor was represented as a linear combination parameterized by IMF B z and solar wind ram pressure P d and including four free coefficients a 0 -a 3 .Note here that, formally, the shielded dipole field configuration and magnitude are entirely defined by the dipole moment and the magnetopause size/shape driven by the solar wind pressure and IMF via Eq. 3, such that at first sight the factor F might appear redundant.The reason behind its inclusion is to increase the flexibility of the model by empirically taking into account variations of the subsolar magnetic field due to intra-magnetospheric external currents and a possible interplay between the IMF and P d drivers.In any case, the magnetospheric Eqs 2-4 has only four linear parameters, which definitely rules out any overfitting problems.This conjecture was computationally confirmed by that the obtained values of a 0 were found very close to unity, while the remaining three coefficients a 1 -a 3 turned out rather small (see Table 2

below).
In total, the magnetospheric module has 9 free parameters to be found from data: four coefficients entering in Eq. 4 and five nonlinear parameters in Eq. 3. The fitting was performed using a Nelder-Mead simplex search in 5D parametric space, at each step of which the coefficients were found by a standard SVD algorithm.
In regard to all the above described formalism, a subtle issue should be highlighted: while the magnetospheric Eq. 2 is fitted to only the magnetospheric data (selected in advance by means of the diagrams in Figures 3, 4), it extends throughout the entire space as a smooth curl-free magnetic field and, hence, does not reproduce the required discontinuity of its tangential component at the magnetopause.A natural way to incorporate the magnetopause is to use the coordinate δ = 1 − r S (r, θ) /R S , where r S (r, θ) = r cos 2α (θ/2) (5) as an independent indicator of a data point {r, θ, ϕ} location with respect to the boundary, such that δ > 0 or δ < 0 correspond, respectively, to the inside or outside of the magnetosphere.parameter not known in advance, calculations of δ = δ(r, θ, R S 0 ) and, hence, of the classification factor H(δ) must be made individually for each data record and at each simplex step of the search algorithm.
As confirmed in the fitting experiments, the obtained best-fit magnetopause parameters converge to very reasonable and stable values, for both training (T) and validation (V) subsets (Table 2).Therefore, the developed method can be employed in the future as a instrument for constructing the dayside magnetopause models without using direct crossing data.Table 2 presents basic characteristics of both T and V magnetospheric data subsets: record numbers N and r.m.s.field magnitude ⟨B⟩, as well as the obtained best-fit values of the corresponding module parameters, starting from the residual r.m.s.deviation of the model from data ⟨ΔB⟩.The first thing to note is the strikingly stable values of the average standoff distance R S 0 ≈11.1 R E , equal to each other in T and V variants and in good agreement with many independent estimates (see (Samsonov et al., 2016) and refs.therein).A remarkable fact here is that the calculation used only in situ magnetometer data around the magnetopause, but did not employ any a priori information of the boundary location, nor direct crossing data.The same applies to the parameter ε, quantifying the response of the standoff distance to variations of the solar wind ram pressure: in both cases it is found rather close to its classical theoretical estimate 1/6 ≈0.167 (Mead and Beard, 1964).Here it should be noted in passing that Jelínek et al., 2012 andDušik et al., 2010 found significantly larger exponents, equal to 1/5.5 ≈ 0.18 and 1/4.8 ≈ 0.21, respectively.Such a mismatch is not surprising, because of completely different methods of the magnetopause identification adopted in those works.The next parameter γ, defining the IMF influence on the boundary subsolar distance, is also close to its validation counterpart and reveals a significant sensitivity of R S on the IMF B z , such that negative B z result in decreasing standoff distance.The obtained values of the flaring parameters α 0 and Δα also demonstrate a strong response of the dayside magnetopause shape to the IMF B z orientation associated with the erosion and flux pile-up effects, with negative/positive IMF B z resulting in larger/smaller magnetopause radius R T in the terminator plane.All these aspects will be graphically illustrated in more detail in Section 4.

Magnetosheath module
The magnetosheath magnetic field is represented by a separate mathematical framework, first described in TSE23.It is based on a sum of toroidal and poloidal components, whose generating functions are expanded into triple sums with Taylor series in powers of the normal coordinate δ = 1 − r S (r, θ)/R S and spherical functions Y (mn) (θ, ϕ) of the angular coordinates θ and ϕ (Figure 5), such that and where the summation limits in Eq. 6 are as follows: l = 0, …, 3, n = 1, …, 6, and m = 0, …, n.The coefficients a (lnm) and b (lnm) are further expanded into linear combinations of principal driving parameters and their cross-terms, including IMF B x , B y , B z , a pressure-dependent factor similar to that in Eq. 4, as well as first two powers of the dipole tilt angle ψ.Having imposed appropriate north-south and dawn-dusk symmetry/antisymmetry conditions (addressed in greater detail in Section 2 of TSE23) resulted in a total number of coefficients to be found from data equal to 960.The magnetosheath module also includes the same five magnetopause parameters R S 0 , ε, γ, α 0 , Δα, found earlier from fitting the magnetospheric module and now fixed at the obtained values.They are needed here to properly define the variable reference surface with δ = 0, around which the generating functions Ψ T and Ψ P are expanded in powers of δ.In this case, the expansions Eq. 6 are fitted to only the magnetosheath data subset, also compiled by means of the diagrams in Figures 3, 4. As in the case of the magnetospheric module, an apparent extrapolation problem comes up: namely, inside the magnetopause, the expansions Eqs 6, 7 are no longer constrained by data and, hence, would be inconsistent with magnetospheric observations.The easiest way to resolve this issue is to simply cut off the magnetosheath field inside the magnetosphere using the step function H(−δ).This is justified by the fact that in the present study we do not intend to develop a unified model of the entire magnetosheath-magnetosphere system, but focus on the interface between the two regions, in particular, on the magnetopause electric current.The current surface density is derived from the two independent modules by calculating the total jump of the tangential component of the magnetic field across the boundary.The toroidal component in Eq. 6 is, by construction, tangential to the magnetopause at any location and, hence, can be safely multiplied by H(−δ), which nullifies that part of the magnetosheath field on the magnetospheric side without violating ∇ ⋅ B = 0.As for the poloidal component of Eq. 6, it is normal to the magnetopause and, hence, even though it becomes discontinuous at δ = 0, its elimination inside the boundary does not affect the calculated electric current.
The above described magnetosheath module was fitted to the training (T) and validation (V) subsets, containing, respectively, 351,069 and 355,372 data records.The corresponding r.m.s.field magnitudes were, naturally, significantly lower than those for the magnetosphere: ⟨B⟩ = 20.76nTfor T subset and ⟨B⟩ = 21.27nT for V data.The corresponding residual deviations, however, have nearly the same magnitude (hence, larger relative noise levels around 60%) due to much higher turbulence in the magnetosheath: ⟨ΔB⟩ = 12.44nT and ⟨ΔB⟩ = 13.00nT for T and V subsets, respectively.
Overall, the magnetosheath field magnitude steadily decreases with growing clock angle, reflecting the progressive decrease of the magnetic flux pile-up and increase of the dayside field erosion.In the case of purely southward IMF, two nearly symmetric and rather shallow field compression areas are formed near the magnetopause, centered at ∼10 and ∼14 MLT hours.Given a multitude of factors at play (of physical, observational, and modeling nature), it is hard to offer an unambiguous interpretation; among plausible causes can be strong fluctuations due to intermittent reconnection, and/or partial mixing of data taken in the vicinity of the dayside magnetopause.To get a better idea about the degree of stability of these features, we calculated the model parameters and generated similar plots (Figure 6), based on the independent validation subset.
The magnetic field distributions are quite similar to those in Figure 6; the largest difference is seen in the case of purely southward field (right panel), where the subsolar magnetic field and its radial gradient are substantially weaker, and the off-center compression areas are more localized.As in the previous case, a possible reason can be conjectured as due to a larger field fluctuation level in the subsolar region during the periods of southward IMF.
The next Figure 7 shows in a similar format three equatorial field distributions for more commonly observed IMF orientations, corresponding to Parker spirals with a nonzero radial component: B x = − 4nT and B ⊥ = 7nT.The left, center, and right panels correspond to IMF clock angles ϕ equal to 45°, 90°, and 135°, respectively.The plots also include B-field vectors projected onto the Z GSW plane, demonstrating the magnetic field draping around the magnetopause.Note the strong asymmetry of the field magnitude in the pre-noon and post-noon sectors, due to the different IMF orientation with respect to the magnetopause and bow shock: nearly orthogonal at dawn and parallel at dusk.One also sees a significant decrease of the magnetosheath field magnitude, as the clock angle increases from 0°to 180°.In the discussion Section 5, these issues will be addressed in more detail.The next Figure 9 displays the currents in the same format and for the same IMF orientations, but for the maximum value of the dipole tilt ψ = 30°, revealing in more detail the current configuration in the vicinity of the polar cusp.In this case, the northern cusp shifts to lower GSW latitudes where the data coverage is much denser, which allows to reproduce/visualize the current geometry in more detail.In particular, in the panel (b) with IMF B y = + 4.3nT, B z = + 2.5nT one can see a strong deformation of the northern cusprelated vortex, with a significant dawnward MLT shift of its center by ∼5 R E .An early model-based estimate (Tsyganenko and Usmanov, 1984) predicted the MLT shift of the cusp footpoint around |Δλ| ∼ 4−8°for IMF B y ≶ 0.Here the shift refers to the outermost throat of the cusp, rather than to its footpoint, which is why the effect is much larger.

Electric currents on the magnetopause
Also note that the deformation is not as strong in the panel (c), corresponding to the same IMF B y but negative B z = − 2.5nT.This is due to much stronger azimuthal current, which partially offsets the dawn-dusk asymmetry and is akin to a similar increase of the magnetotail "rigidity" against tilt-related deformations of the magnetosphere (Tsyganenko and Fairfield, 2004;Tsyganenko et al., 2015).Another noteworthy detail is a progressively higher latitude of the current vortex center as the IMF rotates from northward to southward, mostly due to the overall increase of the duskward current density in the subsolar region.

Discussion: validation and pressure balance issues
As already mentioned (Section 4.1; Figures 6, 10) the modeling calculations were tested by reconstructing and comparing magnetic field configurations, based on two independent subsets of data, nearly equal in size.Another commonly used testing method is to generate a model from the training subset and compare its output with the validation data by creating scatterplots of the model field components against their observed values.
Figure 11 shows a result of such a comparison; in the upper three panels, the magnetosheath model field components are plotted against the corresponding data from the training subset, while in the bottom row the model output, calculated using the training parameters, is compared with the validation data.The comparisons reveal a noteworthy difference between the plots for B x and B y , on the one hand, and those for B z , on the other.In the former case, the B x and B y slopes are slightly below unity (as they should be), the validation correlations are reasonably high (0.67 and 0.77) and, as expected, somewhat lower (for B y ) or equal (for B x ) in comparison with the training variants.By contrast, the B z plots reveal abnormally extended green areas of large values above the regression line, implying a large overestimate by the model.Accordingly, the slopes exceed unity (1.07-1.09), the regression lines are shifted upwards from the main red core, and the correlations are substantially lower (0.46).The most plausible cause is the already mentioned ambiguity in separating the magnetosheath and magnetosphere data, discussed in Section 2. As can be seen from Figures 3, 4, that uncertainty applies in greater or lesser extent to data of all missions.As a result, a tangible portion of outer magnetospheric data with large positive B z values creeps into the magnetosheath data selections, manifested in Figure 11 as the broad clouds above the main core of magnetosheath data around the main diagonal.This feature calls for more accurate methods of data selection in the future modeling studies.
The next Figure 12 shows similar scatter plots for the magnetospheric module.Here all correlation coefficients are significantly higher, reflecting more ordered field structure and much lower noise level in the data.At the same time, all the slopes for both T and V subsets significantly exceed unity, which is obviously a result of a somewhat simplified module architecture, in which such important field sources as the ring, tail, and field-aligned currents are missing.Another possible cause is the already noted partial "diffusion" of the magnetosheath data into the magnetospheric subset.In the right panels, an overwhelming majority of data points lie in the first quadrant, due to the largely northward outer magnetospheric field.
Another independent test of the model's consistency deals with the issue of net force balance between the solar wind and the magnetosphere.In particular, according to theory (Spreiter et al., 1966;Petrinec and Russell, 1997), the total pressure along the Sun-Earth line should remain constant.Left and right panels in Figure 13 display the profiles of the model magnetic field variation along the Sun-Earth line, derived from the training and validation subsets, respectively.In both cases, the IMF total magnitude equals 8.0 nT and the cone angle is 120°, while the red and blue lines in each panel correspond to two values of the IMF clock angle: ϕ = 45°(red) and ϕ = 135°(blue) with IMF B z = + 5.0 and −5.0nT, respectively.These two cases correspond to the model configurations shown in the left and right panels of Figure 7.
The field magnitudes jump from 8 to ∼16 − 24nT at the bow shock, roughly consistent with estimates based on Rankine-Hugoniot theory (Petrinec and Russell, 1997), then gradually increase towards the subsolar point to 28 − 32nT, and finally experience the second upward jump ΔB ∼ 18 − 25nT at the magnetopause.In the training variant, the respective geocentric standoff distances for positive and negative IMF B z = ± 5nT were found equal to 11.57 and 10.72 R E , such that the earthward shift of the subsolar point due to the southward IMF reversal ΔR S = 0.85R E ; in the validation variant, the standoff distances are 11.54 and 10.77 R E and ΔR S = 0.77R E .For comparison, in the models by Lin et al. (Lin et al., 2010) and Shue et al. (Shue et al., 1997) the corresponding shifts due to the same IMF B z reversals from +5 to −5 nT are 0.83 R E and 0.51 R E , respectively.
The first thing that draws attention in the plots is that, in both T and V variants, the plots for negative and positive IMF B z virtually merge together inside the magnetopause.This means that the earthward shift of the subsolar point by ∼0.8 R E caused by the IMF B z reversal to south is not accompanied by a compression of the magnetosphere, but represents the classical erosion of the dayside magnetic flux, found as early as half century ago (Aubry et al., 1971) (see also (Tsyganenko and Sibeck, 1994)).
As regards the overall pressure balance, in both T and V cases the subsolar field magnitudes just inside the magnetopause are in the range 50-55 nT, corresponding to the magnetic pressures ∼1.0-1.2 nPa, which is only 50%-60% of the input ram pressure 2 nPa in the upstream solar wind, assumed in the above model calculations.A number of possible reasons of various nature can be envisioned to explain the mismatch.One of them is the unaccounted   Variation of the total magnetic field B along the Sun-Earth line for a typical spiral-type IMF orientation, with northward IMF B z = + 5.0 (red) and southward B z = − 5.0nT (blue).In both cases, the total IMF magnitude is 8 nT and its cone angle equals 120°.
subsolar plasma pressure inside the magnetopause, which can be as high as ∼0.5nPa (Li et al., 2023).Another important factor is the already mentioned still imperfect separation of magnetosheath and magnetosphere data, selected on the basis of Figures 3, 4 diagrams.Next, one cannot ignore the fact that the magnetosheath is an extremely turbulent region and the magnetopause is a dynamic structure, always in incessant motion and only rarely staying in static force equilibrium.An additional factor to also keep in mind is a high fluctuation level in the solar wind and IMF data, due to inevitable propagation errors from the L1 point to the subsolar magnetosphere (Vokhmyanin et al., 2019).All these aspects are still to be explored and understood in future studies.

Summary and outlook
In this paper we presented first results of an empirical databased modeling of the magnetic field and electric current structure at the interface between the dayside magnetosheath and the adjacent outer magnetosphere.Mathematically, the model is constructed as a composite framework consisting of two independent modules: 1) a magnetospheric module reproduced by the field of a dipole with variable strength and tilt, fully shielded within a magnetopause with variable size and shape, and 2) a magnetosheath module, based on flexible toro/poloidal expansions in a coordinate system, specially suited for the magnetopause and bow shock geometry.The experimental database includes multi-year sets of 1-min average magnetic field and plasma data of Themis, Geotail, Cluster, MMS, and OMNI source of interplanetary data.Based on a method by Jelinek et al. (Jelínek et al., 2012), the data are selected into magnetospheric and magnetosheath subsets, occupying distinctly different regions in the 2D space of normalized plasma density and magnetic field.Fitting the magnetospheric module to the data made it possible to derive not only the magnetic field, but also model magnetopause parameters, such as the standoff distance and the curvature/flaring rate of the boundary, as well as their dependence on the solar wind ram pressure and IMF B z .Based on the obtained boundary and complementing the model with a corresponding bestfit magnetosheath module allowed us to derive spatial patterns of dayside Chapman-Ferraro currents for different IMF clock angles and infer their response to the Earth's dipole tilt angle.
This work should be viewed as only the first step in the empirical modeling of the IMF effects on the dayside magnetosphere boundary.Our procedure of data selection and region identification is still far from being perfect, which results in a certain degree of data intermixing between different regions.The employed method based on the diagrams needs to be upgraded by taking into account as much as possible in situ telltale signatures, such as fluctuation levels of the field and plasma parameters.Another improvement area is to refine the outer magnetosphere model by including basic extraterrestrial field sources parameterized by ground-based activity indices, in order to take into account the magnetic flux redistribution between the dayside and tail lobes, which may strongly affect the position of polar cusps on the dayside.One more way to further advance the modeling is to generalize the magnetopause shape by modifying the {δ, θ, ϕ} coordinate system and abandon the assumption of axial symmetry.One should also keep in mind the essentially dynamical nature of the magnetosphere-magnetopause-magnetosheath system, in view of which the adopted parameterization by concurrent 1min averages of the solar wind and IMF observables may be not sufficient.In a more remote perspective, quite an attractive task is to use the already developed technique and large databases for constructing a unified empirical model of the entire system solar wind-magnetosheath-magnetosphere.

FIGURE 2
FIGURE 2Spatial distribution of data records in the preliminary subset.Data of four missions are shown by different colors, as indicated in the inset legend.Only every 50th data points out of the total of 11,120,813 are plotted, limited to planar 4 R E -thick layers centered about equatorial (left), meridional (center), and terminator (right) planes.Solid and dashed contours show, respectively, nominal and extreme positions of the MP (red) and BS (black).

FIGURE 3 Four
FIGURE 3 Four D/D sw vs. B/B sw diagrams for four mission data used in this work: Themis, MMS, Cluster, and Geotail.Light/dark dashed lines encircle the magnetosphere/magnetosheath selection areas.

FIGURE 4 W
FIGURE 4 W/W sw vs. V b /V sw diagrams for Themis (left) and MMS (right).As in the previous figure, light/dark dashed lines encircle the magnetosphere/magnetosheath selection areas.

FIGURE 6
FIGURE 6Equatorial plots of magnetic field magnitude for IMF B = 5 nT, IMF cone angle 90°, and three values of the clock angle: 0°(left), 90°(center), 180°(right).Magnetopause/bow shock and B isointensity contours are shown with white and black lines, resp.Due to the color saturation, the inner magnetospheric field is shown with white.Note the steady growth of the field depression in the subsolar area for larger clock angles.

Figure 8
Figure8shows front views of the electric current distribution on the dayside magnetopause for four values of the IMF clock angle ϕ = 0°, 60°, 120°, and 180°.The current surface density J = n × ΔB/μ 0 was calculated in mA/m from the net field jump ΔB across the boundary with a local unit normal n.The plots correspond to IMF

FIGURE 8
FIGURE 8Color-coded electric current surface density (in mA/m) at the dayside magnetopause as viewed from the Sun, for four values of the IMF clock angle: ϕ = 0°(A), 60°(B), 120°(C), and 180°(D).The J vector projections on Y-Z plane are also shown with arrows.Here and in the next Figure9the training subset was used for the model fitting.

FIGURE 10
FIGURE 10Same as in Figure6but based on the validation data subset.

FIGURE 11
FIGURE 11Scatterplots of the magnetosheath module field components against observations.Panels in the upper and lower rows correspond to the training and validation subsets, respectively.Data point numbers, Pearson correlation coefficients, and best-fit slopes are indicated in legends on each panel.

FIGURE 12
FIGURE 12Same as Figure11but for the magnetospheric module.

TABLE 2
Magnetospheric dataset and module parameters, derived by fitting to the T and V data (see Eq. 3 and text for notations).Difference between the T-and V-based values gives a rough measure of a parameter stability and error.
Multiplying the total model field Eq. 2 by the Heaviside step function H(δ) does not violate the ∇ ⋅ B = 0 condition, since the normal component B MSP ⋅ n is zero by construction.An important thing here is that, since the standoff distance R S 0 is free nonlinear