We used electron density measurements from four spacecraft, namely, the International Sun-Earth Explorer 1 (ISEE-1), the Combined Release and Radiation Effects Satellite (CRRES), Polar, and the Imager for Magnetopause-to-Aurora Global Exploration (IMAGE) spacecraft.

Table 1 shows the spacecraft specifications, including the spacecraft orbits and the years and phases of the solar cycle when the electron density measurements were made. The selection of data has good coverage of spatial locations within the inner magnetosphere and different conditions over the solar cycle. Using nonlinear regression (Section 2.3), which takes into account separate dependencies independently, it is possible to combine these measurements to determine the characteristics of the plasmapause density for a variety of conditions (when combining data from different spacecraft, realistic results are not always obtained using simple binning, as demonstrated by Denton et al. (2022)). For instance, some of our formulas depend on the solar EUV index F10.7, which varies with the solar cycle.

All electron density measurements were found from plasma wave measurements, either using the upper hybrid resonance frequency or the lower edge of the continuum frequency within the plasmatrough (LeDocq et al., 1994; Gurnett et al., 1995; Carpenter et al., 2000; Denton et al., 2012). Typical errors in electron density might be roughly 20% (Denton et al., 2022).

The plasma wave instruments used were the Plasma Wave Investigation on ISEE (Gurnett et al., 1979), the Plasma Wave Experiment on CRRES (Anderson et al., 1992), the Plasma Wave Instrument on Polar (Gurnett et al., 1995), and the passive radio data from the Radio Plasma Imager Investigation on IMAGE (Reinisch et al., 2000). Because the electron density is proportional to the square of the plasma frequency, the highest possible electron density that can be measured by one of these instruments is determined by the upper limit of the frequency. The maximum frequency and maximum electron density values are listed for each spacecraft in Table 1. The lowest upper limit of frequency is for ISEE and CRRES, which could measure electron density only up to 1984 cm⁻³ (this explains the cutoff in measurements that is discussed in Section 3.1).

The categorization of data points by region was performed manually using a computer program that generated plots like Figure 1, using data from ISEE observed on 10 November 1997, between 1800 and 2400 UT. The groupings of data points were selected by dragging the mouse over a selection of points. Figure 1 shows data points in the plasmasphere, plasmapause, plasmatrough, and magnetosheath. The slope of the plasmapause was defined using the outermost data point in the plasmasphere and the innermost data point of the plasmatrough, and was required to be at least a factor of 3 within a change in L of 0.4. This corresponds to a slightly smaller slope than that required by Carpenter and Anderson (1992), a factor of 5 within a change in L of 0.5. We used the smaller slope criterion because there were many events for which there was clearly a plasmapause, but the slope was not quite as large as that required by Carpenter and Anderson (1992).

In our database, a plasmapause was defined as a drop in density of at least a factor of 3, but for this study, we use only a subset of the plasmapause observations with a drop in density of a factor of 5. Data points within the plasmapause were within the largest region with the required slope, as shown in Figure 1.

When there was more than one region of L with the required slope, only the one at the lowest L was considered to be the plasmapause for the purpose of this study.

Note that a steep gradient identifying the plasmapause is not always observed, as noted by Denton et al. (2004). This situation can occur when plasmaspheric refilling has occurred for a long enough time so that the entire magnetosphere has density at plasmaspheric levels or when a spacecraft passes radially through a plasmaspheric plume (normally on the dayside). However, the dataset used here only includes those observations for which a plasmapause is observed.

In addition to the condition on the slope, some other conditions were required. For the polar orbiting spacecraft, namely, Polar and IMAGE, only the high-altitude part of the orbit with a geocentric radial distance greater than the minimum L shell sampled by the orbit was used. This criterion was imposed because a plasmapause is not as evident at low altitude and to avoid errors in mapping from low-altitude positions to the equator to determine the L shell.

In order to sample the region of space that is usually within the magnetopause, we limited the L value of the midpoint of the plasmapause, L _pp, to values ≤ 10 , and because of the instrumental limitations discussed below, we also required L _pp ≥ 3.25.

The plasmapause was required to extend over no greater a range in L (from the outermost plasmasphere data point to the innermost plasmatrough data point) than 2, and no greater a range in MLT than 2.5 in order to exclude observations with large data gaps. The slope dL/dR _MLT, where dL is a change in L and dR _MLT = (dMLT)(2 pi L/24) is roughly the change in the distance around the Earth due to a change in MLT of dMLT, was also required to be at least tan 10 ° = 0.176 to reduce the probability of confusing azimuthal structure with that in the radial direction.

Because we found that some formulas depended on the Auroral Electrojet Index AE, we also required that the quality factor for the averaged measurements of AE (defined below) be at least 1, using quality factors analogous to those described by Qin et al. (2007) (basically, a quality factor of 1 means that the observed quantities may be interpolated, but that they need to be within a correlation time of the observed measurements).

With these conditions, we identified 1,276 plasmapause segments, 440 from ISEE, 692 from CRRES, 5 from Polar, and 139 from IMAGE. The CRRES data (at solar maximum) are overrepresented considering the time span of the mission, but our formulas do take into account the geomagnetic activity, so we hope that this will ameliorate that problem.

Figure 2 shows the equatorial distribution of plasmapause segments (between the outermost plasmasphere data point and the innermost plasmatrough data point) in the solar magnetospheric (SM) or dipole equatorial plane. There is a good distribution at all MLT values between L = 3.25 and about L = 8.5.

For the purpose of developing models, the position of the plasmapause was defined as the position of the spacecraft at the midpoint between the outermost plasmasphere data point (data point at the beginning of the steep slope) and the innermost plasmatrough data point (data point at the end of the steep slope), and the electron density at the midpoint of the plasmapause, n _e,pp, was defined as the logarithmic mean of the outermost plasmasphere and innermost plasmatrough density values (see Figure 1).

Other categories of plasma were also defined, as described in the supplementary information. However, the categories described above are sufficient to define the plasmapause for this paper.

In order to model the density at the midpoint of the plasmapause, n _e,pp, we did some preliminary analysis using the Eureqa nonlinear genetic regression program (Schmidt and Lipson, 2009), a proprietary product owned by the DataRobot Inc. company, in a manner similar to that of Denton et al. (2022). However, because of changes in DataRobot Inc.’s pricing structure, we switched to the TuringBot program, which functions similarly; the final formulas that we presented were found using TuringBot (see the Data Availability Statement for links to these software programs).

Both of these programs are nonlinear genetic algorithms that find models using analytical formulas. They search the space of possible formulas, spending more time searching for types of formulas that yield less error. Both programs determine a set of formulas of varying complexity. Our TuringBot runs sampled billions of different formulas, and the best-fitting model for each level of complexity was determined using cross-validation with five groups of randomly sorted data. The TuringBot complexity is described by the mathematical operations: a cost of 1 for introducing a variable or constant or calculating a sum, difference, or product; a cost of 2 for division; and a cost of 4 for all other operations, including logarithms and power laws.

Our modeling of n _e,pp proceeded through several iterations in which we examined the dependence of n _e,pp on a number of different parameters: position; phase of the year; geomagnetic indices like the planetary K index (Kp), The disturbance storm time index (Dst), and AE; solar wind parameters like components of the interplanetary magnetic field (IMF), B _Y and B _Z , the solar wind density, solar wind velocity, V, and solar wind dynamic pressure, Pdyn; the F10.7 index measuring solar radiation at 10.7 cm wavelength; the solar wind convection electric field defined by VB _s , where B _s = −B _Z if B _Z < 0 and B _s = 0 otherwise; and the coupling function of Newell et al. (2007), dΦ_MP/dt.

By far, the strongest dependence that we found was on the L value at the midpoint of the plasmapause, L _pp, but some of our models had a dependence on the magnetic local time, MLT, at the midpoint of the plasmapause, MLT_pp, and some geomagnetic indices. From our early runs using Eureqa, we found that a minimum in n _e,pp with respect to MLT occurred at MLT = 3.43 h. For the final modeling runs using TuringBot, we searched for dependence on L _pp, cos (2π(MLT −3.43)/24), sin (2π(MLT −3.43)/24), cos (4π(MLT −3.43)/24), sin (4π(MLT −3.43)/24), F10.7, Dst, and averages of Kp and AE over the previous 3, 6, 12, 24, 48, 96, and 192 h. Terms averaging up to 96 h appeared in the final models. The TuringBot runs searched for the lowest root mean squared (RMS) error of the difference between the observed and modeled log₁₀ (n _e,pp).

In Figure 4, line segments are drawn from the density of the outermost plasmasphere data point to that of the innermost plasmatrough data point versus L. There is a rough upper boundary for the high-density data points and a rougher lower boundary for the low-density data points. There is a steep L dependence for the density values within the plasmapause.

Figure 4 (unlike the other plots showing our dataset) is plotted without imposing the condition L _pp ≥ 3.25. Note that most of the electron density values are cut off at approximately 2 × 10³ cm⁻³ (see Section 2.1). This cutoff is caused by a high-frequency limit for the instrumentation that made most of our plasma frequency measurements. For that reason, we limited the data that we used for modeling to that above L = 3.25 in order to avoid truncation of the upper portion of some plasmapause segments (which would artificially decrease the corresponding midpoint densities). As can be seen from Figure 4, for L _pp ≥ 3.25, the outermost plasmasphere density (the high-density end of the line segments in Figure 4) is less than 2 × 10³ cm⁻³.

The upper thick solid gold curve in Figure 4 is a best least-squares quadratic polynomial fit to the logarithm of the density of the outermost plasmasphere data points, which is given as follows: log 10 n e = 4.68 − 0.646 L + 0.02632 L 2 . (3) The lower thick solid gold curve in Figure 4 is a best least-squares cubic polynomial fit to the logarithm of the density of the innermost plasmatrough data points, which is given as follows: log 10 n e = 5.37 − 1.58 L + 0.188 L 2 − 0.00865 L 3 . (4) The middle thick solid gold curve is a simplified (only L-dependent) version of our most complicated model (model 14 in Table 2) for the density at the midpoint of the plasmapause. log 10 n e = − 2.042 + 32.67 / 3.981 + L . (5) To derive Eq. (5), we used median values from all of our data for the quantities in model 14 other than L, that is, cosine of the magnetic local time, cMLT’ = cos (2π(MLT −3.43)/24), Kp_96h, F10.7, and AE_6h.

For comparison, the upper dashed yellow curve in Figure 4 represents the Carpenter and Anderson (1992) plasmasphere density (summarized in item 2 of their Section 3), while the lower dashed yellow curve represents the Carpenter and Anderson plasmatrough density given by the average of the two formulas in Carpenter and Anderson’s Eq. (6) (or item 4 of their Section 3). The middle yellow dashed curve is given as follows: n e,pp = 13.4 6.6 / L 4 c m − 3 , (6) which corresponds to model 7 in Table 2. This model is of the same form as Equations (1) and (2) but provides a better fit to our data as described below.

The middle and bottom thick solid gold curves in Figure 4 are fairly close to the middle and bottom dashed yellow curves, showing that the densities of our innermost plasmatrough data points are fairly close to those expected from Carpenter and Anderson’s models and that Eq. (6) is a reasonably good approximation to the L dependence of the more complicated models to be described below, especially at L ≤ 6. The upper thick solid curve in Figure 4, found from the outermost plasmasphere data points, crosses the upper yellow dashed curve from Carpenter and Anderson’s plasmasphere density model at L = 2.7 and L = 9.5 and is lower in the middle range of L values.

The results of our modeling are described in Table 2. Listed there are the model number, N, the TuringBot complexity, Cmp, the formulas for n _e,pp (models 1–8) or log₁₀ (n _e,pp) (models 9–14), the RMS error between the observed and model values of log₁₀ (n _e,pp), and the percentage of data points for which the model value of n _e,pp was between the density of the outermost plasmasphere data point and the innermost plasmatrough data point. In other words, the last column of Table 2 shows the percentage of model values that are within the observed range of plasmapause densities. If a model were to yield 100% in the last column of Table 2, this would mean that the model value of density would always be lower than the plasmasphere value and higher than the plasmatrough value, signifying that the model could always be depended on to separate plasmasphere and plasmatrough populations (at least for positions close to the plasmapause).

The last two columns in Table 2 have three numbers each. The first number is for the data that were used to train the models, the second number is for test data that were not used as input to the models, and the third number is for all the data. The plasmapause events were first sorted by time and then divided into 50 groups. Within each group, the middle one-sixth of the data were separated for test data in a manner similar to that of Denton et al. (2022). The purpose of this procedure was to have test data that were uncorrelated with the training data (Denton et al., 2022). The median time span of the test data intervals was 4.6 days.

The third number of three numbers in the last two columns of Table 2, respectively, showing the RMS error and percentage of events within the plasmapause for all data, including both training and test data, have values that are always intermediate between the first and second numbers. The second number for RMS error is always larger for the test data than for the training data, which is not surprising considering that the models were chosen to minimize the RMS error of the training data. However, the percentage of model values within the observed plasmapauses is sometimes greater for the training data and sometimes greater for the test data. One might reasonably take the largest RMS error and the smallest percentage of events within the plasmapause as an upper-limit measure of error.

Models 1–8 in Table 2 were motivated by the formulas that people have previously used to model the boundary between plasmasphere and plasmatrough data, either constant values in models 1–4 or power law forms in models 5–8. Model 4 is for the optimal constant value for our data set, 48.7 cm⁻³, very close to 50 cm⁻³ (model 2). For model 4, the RMS error of log₁₀ (n _e,pp) for the test data is 0.502, corresponding to a factor of 10^0.498 = 3.2, and 70.2% of the model values occur within the range of density values within the observed plasmapauses.

The use of a power law (models 5–8 in Table 2) results in significantly smaller values for the RMS error and significantly larger percentages of model values within the density range of our plasmapause events. Using the form C (6.6/Lpp)⁴ (models 5–7), the value of C = 15 cm⁻³ used by Denton et al. (2004) has a lower error and a larger percentage of model values within the observed range of plasmapause densities than does the value of C = 10 cm⁻³ used by Sheeley et al. (2001). However, the optimal value for C for our data is 13.4 (model 7), which is between 10 and 15. Based on the numbers in Table 2, there is not much benefit to letting the power law coefficient vary to its optimal value of 4.19 in model 8.

Models 9–14 show solutions from TuringBot. TuringBot finds formulas of varying complexity. We only include more complicated formulas if the TuringBot cross-validation error, our training data set error in Table 2, and our test data set error in Table 2 are less than those for formulas with lower complexity. We do not show the TuringBot cross-validation errors, but they were less than our test data errors in every case.

TuringBot did not find solutions in the form of the power laws because TuringBot’s solution with a complexity of 10 had smaller errors and larger percentages than the power law forms, which also had a complexity of 10 (this could be different if the weights for different operations were altered, but currently TuringBot does not allow that option). TuringBot’s model 10, with a complexity of 8, had comparable errors and percentages to the power law models.

The TuringBot models clearly show that the strongest dependence of n _e,pp is on L (model 9). The more complex models in Table 2 have a dependence on cMLT’ = cos (2π(MLT _pp −3.43)/24), as noted at the bottom of Table 2. The more complex models also have dependencies on F10.7 and averages of Kp and AE. The model with the lowest error in Table 2, that is, model 14, has an RMS error in log₁₀ (n _e,pp) of 0.216 for the test data, corresponding to a factor of 1.64. Model 14 yields model values of n _e,pp that are within the observed plasmapause ranges of density values for 96% of the test data.

All of the models in Table 2 also have model values within the observed range of plasmapause densities for most events. For instance, the power law form in model 7 yields model values that are within the observed range of plasmapause densities for 94% of the test data events. So from the perspective of that number, 96% for model 14, based on the test data, is only a small improvement. On the other hand, considering the number of events for which the model value is not within the observed range of plasmapause densities, the improvement in the model value from model 7 to model 14 is more significant, that is, from 6% to 4%. In other words, although the percent change in correct predictions was insignificant, the percentage of incorrect predictions was affected more significantly.