Modeling the compressibility of Saturn's magnetosphere in response to internal and external influences

The location of a planetary magnetopause is principally determined by the balance between solar wind dynamic pressure DP and magnetic and plasma pressures inside the magnetopause boundary. Previous empirical studies assumed that Saturn's magnetopause standoff distance varies as DP−1/α and measured a constant compressibility parameter α corresponding to behavior intermediate between a vacuum dipole appropriate for Earth (α≈6) and a more easily compressible case appropriate for Jupiter (α≈4). In this study we employ a 2‐D force balance model of Saturn's magnetosphere to investigate magnetospheric compressibility in response to changes in DP and global hot plasma content. For hot plasma levels compatible with Saturn observations, we model the magnetosphere at a range of standoff distances and estimate the corresponding DP values by assuming pressure balance across the magnetopause boundary. We find that for “average” hot plasma levels, our estimates of α are not constant with DP but vary from ∼4.8 for high DP conditions, when the magnetosphere is compressed (≤25 RS), to ∼3.5 for low DP conditions. This corresponds to the magnetosphere becoming more easily compressible as it expands. We find that the global hot plasma content influences magnetospheric compressibility even at fixed DP, with α estimates ranging from ∼5.4 to ∼3.3 across the range of our parameterized hot plasma content. We suggest that this behavior is predominantly driven by reconfiguration of the magnetospheric magnetic field into a more disk‐like structure under such conditions. In a broader context, the compressibility of the magnetopause reveals information about global stress balance in the magnetosphere.


Pressure Balance at the Magnetopause
The magnetopause is the magnetic and plasma boundary formed around a magnetized planet, caused by the interaction between the solar wind and the planetary magnetic field. It separates the internal planetary plasma of the magnetosphere from the external shocked solar wind plasma of the magnetosheath. In a steady state system, the magnetopause shape and size is determined by pressure balance across the boundary. The effective pressure exerted by the solar wind on a planetary magnetosphere is principally due to its dynamic pressure D P , defined as m v 2 , where m denotes mass density and v is the solar wind velocity. The dayside magnetosphere is compressed when the upstream solar wind dynamic pressure increases and inflates when the dynamic pressure drops. The magnetopause is therefore in constant motion, with a velocity of order 10 km s −1 for Earth's magnetopause [Berchem and Russell, 1982] and 100 km s −1 for Saturn's [Masters et al., 2011]. However, to first order, its location can be approximated by assuming Newtonian pressure balance across the surface, between the component of D P normal to the magnetopause surface and the total internal magnetospheric pressure just inside the magnetopause. A key source of pressure inside all planetary magnetospheres is the magnetic pressure P B =B 2 ∕2 0 due to the total magnetic field strength B, which comprises the internal planetary field and other sources, such as the field associated with a magnetospheric ring current.
A useful proxy of the overall size scale of the magnetosphere is the "standoff distance," R MP . This is the radial distance of the magnetopause boundary measured from the planet center along the planet-Sun line (i.e., the subsolar nose of the magnetosphere). We can estimate this value to first order by finding the radial distance from the planet center r at which the magnetic pressure associated with the magnetosphere is exactly balanced by the solar wind dynamic pressure incident on the magnetopause at the nose. If we make the assumption B ∝ r − , such that the magnetic pressure P B ∝ r −2 , we can then write where a 1 is a constant that determines the size scale of the system and is the "compressibility parameter," equal to 2 . For a vacuum dipole magnetic field = 3, giving = 6. This is broadly consistent with observations of Earth's magnetosphere [e.g., Shue et al., 1997].
In contrast to the terrestrial system, Saturn's magnetosphere has significant internal plasma sources and also rotates more rapidly, with a period of ∼10.7 h [Desch and Kaiser, 1981]. In particular, it was discovered by the Cassini space mission that the icy moon Enceladus, which orbits Saturn at a distance of 3.95 R S (R S = Saturn's equatorial radius, 60,268 km), ejects plumes of water group molecules into the magnetosphere at around 300 kg s −1 , which are then ionized and form a wide torus of plasma around the planet Tokar et al., 2006;Bagenal and Delamere, 2011]. This dense plasma is accelerated to corotation with the rapidly rotating magnetosphere, producing a significant centrifugal force which is directed radially outward on the plasma in its corotating frame of reference. In order for the magnetic tension force to balance this centrifugal force in Saturn's outer magnetosphere, the magnetic field lines may be pictured as being "stretched" radially outward near the equatorial plane, from a dipole structure into a more "disk-like" structure, supported by a corresponding azimuthal ring current.
This type of magnetic field configuration affects the compressibility parameter , since for a disk-like magnetic field, B varies more slowly with radial distance in the outer magnetosphere than for a dipole. This can be seen in Figure 1, which shows the total magnetic field strength measured by the Cassini magnetometer (MAG) [Dougherty et al., 2004] in Saturn's dayside magnetosphere, taken from three equatorial orbits during Saturn equinox (Revs 120−122, 23 October to 17 December 2009). In the inner and middle magnetosphere the data appear to be well approximated by the relationship B ∝ r −3 . However, in the outer magnetosphere (r ≳ 15 R S ), where the magnetic field associated with the ring current is more significant compared to the internal dipole magnetic field, the data appear better approximated on average by a relationship B ∝ r −2 , i.e., a lower value of , and therefore a lower value of . This means that with a magnetodisk magnetic field structure, the magnetosphere size is expected to be more sensitive to changes in solar wind pressure than for a dipole case, as the index −1∕ in equation (1) is greater.
The Cassini mission also confirmed the existence of a hotter (>3keV) and more highly variable population of plasma originating in Saturn's outer magnetosphere, mainly composed of hydrogen and oxygen ions [e.g., Sergis et al., 2009]. In the middle and outer magnetosphere, this population contributes more to the total plasma pressure than the colder equatorial plasma and also contributes to the formation of the magnetodisk field structure, via an enhancement of the ring current activity [e.g., Sergis et al., 2010].
Values of plasma (the ratio of plasma pressure to magnetic pressure) in the range ∼2-5 have been observed in Saturn's outer magnetosphere [Sergis et al., 2010], meaning that under some conditions the hot plasma population may control compressibility directly. The derivation of equation (1) assumes that the solar wind dynamic pressure is directly balanced by the magnetic pressure inside the magnetosphere; however, for > 1 at the magnetopause boundary, the hot plasma pressure inside the magnetosphere is also significant in controlling pressure balance. Therefore, will be partially determined by how the hot plasma pressure P h varies with radial distance (via the same arguments as laid out for equation (1) above). Gold [1959] used the concept of magnetic flux tubes of plasma expanding isothermally to explain that, in a dipolar magnetic field, one would expect the hot plasma pressure to fall with radial distance according to P h ∝ r −4 . We would thus expect = 4 for a fictitious magnetosphere where compressibility is dominantly controlled by a hot plasma population embedded in a dipole magnetic field and a value even smaller for a magnetic field that varied more slowly than r −3 . We expand on this in section 3.2.
Jupiter's magnetosphere also has a significant plasma population, due to mass loading from the volcanic moon Io, which orbits Jupiter at a distance of 5.9 R J (R J = Jupiter's equatorial radius, 71,492 km). Io ejects material into Jupiter's magnetosphere at a rate of approximately 1000 kg s −1 , and the observed plasma the middle and outer magnetosphere is also significantly higher than that at Saturn, reaching 100 at 45 R J Figure 1. Total magnetic field strength B against radial distance from planet center. Black circles are from Cassini magnetometer (MAG) data, measured on Saturn's dayside, taken from three equatorial orbits during Saturn equinox 23 October to 17 December 2009). Example power law curves shown in green and blue for reference. [Mauk et al., 2004]. The associated centrifugal force is also greater than the Saturnian parallel, mainly due to the larger overall size scale of Jupiter's magnetosphere. This aspect produces a more substantial disk-like magnetic field configuration at Jupiter, which significantly affects the magnetospheric compressibility; has been measured empirically for this system as between ∼4 and ∼5 [Huddleston et al., 1998;Joy et al., 2002;Alexeev and Belenkaya, 2005]. However, these estimates are based on limited spacecraft observations and hence have large uncertainties associated with them.
Saturn's magnetosphere would therefore be expected to show compressibility behavior that is intermediate between that of Jupiter ( ≈ 4) and the Earth ( ≈ 6). The interrelated factors of magnetic field structure, centrifugal force and plasma content that determine the overall stress balance may themselves show behavior that varies with solar wind dynamic pressure (i.e., vary with size of the magnetosphere). It is therefore insightful to investigate Saturn's magnetospheric compressibility in response to changes in both external factors (the incident solar wind dynamic pressure) and internal factors (hot plasma content) in tandem.

Models of Saturn's Magnetopause
Previous studies have mostly employed in situ magnetopause crossing data to create empirical models of the size and shape of Saturn's magnetosphere under different solar wind pressure conditions, in order to make estimates of and a 1 . Slavin et al. [1985] used Pioneer 11 and Voyager 1 and 2 plasma and magnetometer data and modeled the magnetopause surface as a conic section. Arridge et al. [2006] applied the formalism first described in Shue et al. [1997] for the terrestrial system to Saturn using Cassini magnetometer data, estimating the upstream solar wind pressure by assuming pressure balance with the magnetic pressure of the magnetosphere just inside the magnetopause. This model was then augmented by Kanani et al. [2010] who used Cassini plasma data in order to include thermal plasma pressure, as well as magnetic pressure, inside the magnetopause. Further improvements to the treatment of internal plasma pressure sources, and an extension of the crossing data set, were then made by Pilkington et al. [2015]. Previous to these Cassini-based studies, Hansen et al. [2005] used a 3-D global magnetohydrodynamic simulation of Saturn's magnetosphere for the time period of Cassini's initial approach in order to investigate . The relevant parameters calculated in each study are shown in Table 1, where parameters relate to equation (1) with R MP in units of R S and D P in nanopascal (nPa). In the study by Pilkington et al. [2015], they explain that their extended Cassini data set includes more variable and high-crossings than previous studies. A k-means clustering algorithm was used to separate the  Kanani et al. [2010] 17-29 10.3 ± 1.7 5.0 ± 0.8 Pilkington et al. [2015] 14-40 10.8-16.5 5.5 ± 0.2 crossing data set into three groups organized by local plasma ; while the estimates of agreed in each group within the uncertainties, the estimates of a 1 did not and showed a trend of increasing with larger average , from ∼10.8 to ∼16.5.
The estimates of a 1 by Arridge et al. [2006] and Kanani et al. [2010] agree within the quoted uncertainties. In addition, in Pilkington et al. [2015], in order to account for the influence of local plasma on a 1 , these authors repeat their analysis with a modification to equation (1) such that D P is replaced by D P /(1 + ). With this adaptation, Pilkington et al. [2015] measured a 1 = 10.5 ± 0.2, which is consistent with previous results shown in Table 1. The estimates of for each study agree at least within 2 uncertainty level, where appropriate, and in general are consistent with the picture of a magnetosphere that behaves intermediately between the rigid Earth case and the more compressible Jupiter case. However, the role of the global hot plasma content in controlling magnetospheric compressibility is still not fully understood.
The current lack of multipoint simultaneous observations in Saturn's magnetosphere means that the use of in situ data is inherently limited, as the large scale structure of the magnetosphere at the exact time corresponding to one magnetopause crossing cannot be readily obtained. This makes it difficult to interpret the global physical processes that are controlling the compressibility behavior. Instead, empirical studies, such as those referred to above, provide an "average" picture of the magnetopause morphology over varying internal and external conditions. In addition, the empirical studies discussed must assume that the magnetopause is in dynamical equilibrium at the time of crossing observations (i.e., it is not accelerating) in order to estimate the solar wind pressure via Newtonian pressure balance, which is almost never the case at Saturn [e.g., Dougherty et al., 2005;Masters et al., 2011;Pilkington et al., 2015].
In this study we adopt a more theoretical approach, in order to complement these previous observational studies. A 2-D axisymmetric force balance model of Saturn's dayside magnetosphere, first described by Achilleos et al. [2010a], is used to investigate magnetospheric compressibility. This model can be calculated at a chosen range of magnetopause radii. The corresponding estimated upstream solar wind pressure can be readily obtained using calculated plasma and magnetic field information just inside the magnetopause boundary and assuming pressure balance across the magnetopause, which is static in this model. These theoretical estimates of D P can then be used to make estimates of and a 1 . In addition, Achilleos et al. [2010a] referred to a "hot plasma index" in the model, with which the hot plasma content in the magnetosphere can be parameterized. We revisit this concept herein, in order to investigate the influence of the global hot plasma population on magnetospheric compressibility.
The details of the model and its use in this study are outlined in section 2. The analysis of the model outputs are then discussed in section 3, along with possible explanations for the underlying drivers of the magnetospheric behavior. We conclude with a summary of the results and suggestions for future investigations.

The Model
The model used in this study, originally described by Achilleos et al. [2010a], is based on a magnetic field and plasma model originally constructed for the Jovian magnetodisk by Caudal [1986], adapted for the Saturn system. More information can be found in Achilleos et al. [2010aAchilleos et al. [ , 2010b. The model is axisymmetric about the planetary dipole/rotation axis, which are assumed to be parallel. The model is constructed based on the assumption of force balance in the rotating plasma of the magnetosphere such that where J is the current density, B is the magnetic field vector, and is cylindrical radial distance from the axis, with e its unit vector. The plasma properties are isotropic pressure P, temperature T, ion number density n, mean ion mass m i , and angular velocity . This equation represents balance between the magnetic Lorentz body force on the left-hand side (representing both the magnetic pressure and tension forces) and the pressure gradient force and the centrifugal force terms on the right-hand side.
By representing the magnetic field as the gradient of a magnetic Euler potential , Caudal [1986] demonstrated that equation (2) is equivalent to the partial differential equation where = cos , the cosine of colatitude, and g(r, , ) is a source function determined by the distribution of plasma and angular velocity in r, space. This equation can be solved semianalytically using Jacobi polynomials as laid out in detail in Achilleos et al. [2010a, Appendix] to give surfaces of constant , magnetic field lines, in r, space. The model solution also provides a prediction of the azimuthal current density components associated with each of the terms on the right-hand side of equation (2). Since the source function g is itself dependent on , equation (3) must be solved iteratively. First, a solution is obtained for a pure dipole magnetic field. This initial solution is then used to update the source function g, to obtain a successive solution, and this process continues until the maximum difference between the calculated values for successive iterations reaches a prescribed tolerance; we use 0.5% in this study, such that we stop calculations once the difference between iterations is less than 0.5%. In practice, it was found that in order to achieve better convergence, it was necessary, after each iteration n, to take the mean between the function i and i−1 , and use this average rather than i to calculate the next iteration [see Caudal, 1986]. The global plasma properties can then be inferred from the calculated magnetic structure and appropriate boundary conditions as described below.
Caudal [1986] explained that as a consequence of equation (2), with T and constant along magnetic field lines (according to Liouville's theorem and Ferraro's isorotation theorem, respectively), the plasma pressure P is determined by where is the "confinement scalelength" (in ) The subscript 0 means the quantity evaluated at the equatorial crossing point of the magnetic field line. This represents the plasma being confined toward the rotational equatorial plane due to the centrifugal force exerted on it. For the hot plasma population, where the thermal energy associated with the plasma is significantly greater than the centrifugal potential, 2 tends to infinity, such that the hot plasma pressure is not confined to the equator but is constant along magnetic field lines, P h = P h0 . This is supported by observations made using data from the Cassini Magnetosphere Imaging Instrument (MIMI) , such as Krimigis et al. [2007]. These authors observe that the hot plasma population extends to high latitudes, particularly on the dayside, verifying that the plasma can effectively fill their flux tubes due to their high energies. Hence, the full form of equation (4) is only necessary for calculating the cold plasma pressure.
The requisite boundary conditions for the model are the equatorial radial profiles of plasma properties. These were obtained from various studies using results from Cassini plasma instruments, as summarized in Achilleos et al. [2010a] and updated in Achilleos et al. [2010b]. For the cold plasma population, the profiles for m i , , and T were obtained from Wilson et al. [2008] and the flux tube content information from McAndrews et al. [2009], both using data acquired by the Cassini Plasma Spectrometer (CAPS) instrument [Young et al., 2004].
As the hot plasma pressure is assumed uniform along magnetic field lines, the plasma population may be completely characterized by a particular equatorial plasma pressure P h0 and flux tube volume V per unit of magnetic flux, where and ds is an element of arc length along the magnetic field line. The integral limits represent measurement along a field line of total length s B between the southern and northern ionospheric footprints at 1R S . The flux tube volume is therefore dependent on both the shape of magnetic field lines, via ds, and the strength of the field, via B.
Studies using Cassini MIMI data such as Sergis et al. [2007] found that the equatorial pressure associated with the hot plasma population was highly variable with and over time. Achilleos et al. [2010a] combined quantile fits of this data set with a radial profile of V obtained from an empirical magnetic field model [Bunce et al., 2007] in order to show a picture of a highly variable hot plasma population, that follows neither adiabatic (P h0 V 5∕3 = constant) nor isothermal (P h0 V = constant) transport behavior.
In light of these observations the original Achilleos et al. [2010a] model simply parameterized the global hot plasma content as follows. In the middle and outer magnetosphere, beyond 8 R S , the hot plasma pressure was calculated by assuming a profile where P h0 V = K h and K h is a constant, known as the "hot plasma index." The observations discussed above indicate that this index may vary in the range 10 5 − 10 7 Pa m T −1 in this region of the magnetosphere. Inside 8 R S , the hot plasma pressure profile was constructed to decrease linearly to 0 with decreasing , such that P h0 ( ) = P h0 ( = 8 R S ) × ( ∕8). A similar parameterization, though with different values of the constants, was made in Caudal [1986], who argued that for the Jovian system, under the expected conditions of rapid radial diffusion, the hot plasma would be transported isothermally. Further discussion and justification of this parameterization can be found in Achilleos et al. [2010a].
Parameterizing the hot plasma content in this way provides the flexibility to very simply characterize the level of ring current activity in the model, and thus investigate the effect of the hot plasma content on magnetospheric compressibility. We recognize that more detailed studies of hot plasma dynamics may require investigation of different parameterizations in the future. In particular, a more physically realistic hot plasma profile was developed in Achilleos et al. [2010b]; however, the use of such a profile does not affect the basic conclusions of this study.

Solar Wind Dynamic Pressure Estimation
As described in section 1, the magnetopause boundary can be approximated as the location where the effective pressure of the solar wind exerted on the magnetosphere is exactly balanced by the sum of the internal magnetospheric particle and field pressures.
Before impacting on the magnetopause, the solar wind flow is first decelerated via the bow shock that forms upstream of the magnetosphere and is deflected around the obstacle. This acts to reduce the dynamic pressure incident on the magnetopause surface and must be accounted for when assuming pressure balance. Petrinec and Russell [1997] used Bernoulli's equation in combination with the Rankine-Hugoniot jump conditions across the bow shock, assuming adiabatic flow of the solar wind, to show that the relation provides an approximation that is valid across the magnetopause surface, not just at the nose. The subscript MS denotes magnetospheric properties, such that the terms on the left-hand side of equation (7) are the magnetospheric magnetic and plasma pressures, respectively. Ψ is the flaring angle, measured between the upstream flow velocity vector and the normal to the magnetopause, such that Ψ = 0 at the nose. The first term on the right-hand side is associated with the solar wind dynamic pressure, where k is a positive constant ≤1 to account for the aforementioned diversion of flow, and the cos 2 Ψ factor accounts for the reduction in the normal component of dynamic pressure on the flanks and tail of the magnetosphere. The second term on the right-hand side is composed of a "static" pressure P 0 associated with the thermal pressure of the solar wind and a sin 2 Ψ factor to ensure a real (i.e., not imaginary) flow velocity in the subsolar region [see Petrinec and Russell, 1997].
The standoff distance at the magnetopause nose R MP is defined where Ψ = 0. We may therefore use a simplified form of equation (7) to easily estimate the solar wind dynamic pressure D P for a given system size, by using the values of B MS and P MS calculated by the model just inside the magnetopause boundary, at the equator. Note that in the model, P MS = P h +P c , the sum of the hot and cold plasma pressure components. The exact value of k depends on the ratio of specific heats in the solar wind and the upstream sonic Mach number M. For high (≳8) Mach number flow with = 5∕3, k = 0.881 [Spreiter et al., 1966]. These conditions are valid for the solar wind at Saturn's orbit [e.g., Slavin et al., 1985;Achilleos et al., 2006]; hence, we adopt k = 0.881 in this study; however, as it is of order unity. it does not significantly affect our estimates of D P or the conclusions of this study. This value of k is also used in the previous studies discussed in section 1.2.

Model Calculations
In this study a value of K h = 1 × 10 6 Pa m T −1 was initially adopted to parameterize the ring current activity in the magnetosphere, representing broadly average conditions, corresponding to the median quantile fit of the hot plasma pressure from Sergis et al. [2007] [see Achilleos et al., 2010a, Figure 6b]. The model was then calculated at 30 different magnetopause radii R MP , equally spaced over a range of 14−40 R S , in order to match the range in this parameter observed in the data set presented by Pilkington et al. [2015]. For each calculation, the corresponding solar wind dynamic pressure incident on the magnetopause nose was estimated from model magnetospheric parameters as described above. A profile of R MP versus D P for a given K h value was then constructed in order to estimate the model compressibility parameter . The value of K h was then varied over the range 10 5 -10 7 Pa m T −1 , in accordance with observations, and the effect on magnetospheric compressibility was investigated. Figure 2 shows the magnetopause radius R MP for each model calculation as a function of the corresponding solar wind dynamic pressure, calculated as described in section 2.2, on a logarithmic scale, with K h = 1 × 10 6 Pa m T −1 . Linear least squares regression lines have been fitted to two regions of the data as described in the caption.

Magnetospheric Compressibility With Average Hot Plasma Conditions
On such a plot, data that exactly obey equation (1) with constant a 1 and would follow a straight line with slope −1∕ and y intercept a 1 . However, it can clearly be seen that for this simulated "data set," neither a 1 nor is constant with system size. In particular, a distinct shift in behavior can be identified at R MP ≈ 25 R S . We therefore divided the simulated data set into two groups, R MP ≤ 25 R S , which we shall refer to as the compressed regime, and R MP > 25 R S , which we shall refer to as the expanded regime. We then fit each data group with a linear least squares regression line separately, in order to make estimates of relevant parameters and investigate how they differ as the system expands. It should be noted that the value of 25R S , in particular, SORBA ET AL. MODELING MAGNETOSPHERIC COMPRESSIBILITY was selected by eye, and a value up to ≈ 28 R S could be chosen to divide the two regimes and does not significantly affect our conclusions.
The estimates for and a 1 with standard error uncertainties, for the two regimes, are shown in Table 2. The full fitting procedure and calculation of uncertainties are described in Appendix A. It should be noted that such uncertainties should be taken as a guideline only. We do not suggest that the data exactly follow this underlying distribution and have made this simplification to give an overall picture of the behavior, given the limitations of our simulated data set. Real observational data would of course be subject to significant measurement errors, and the resulting parameter uncertainties would be higher than those we have calculated for our simulations.
For the compressed regime, the estimates of a 1 and are broadly consistent with the previous results shown in Table 1 at the 2 level, except the estimate of Pilkington et al. [2015]. For the expanded regime, the estimated values are both significantly smaller than those calculated in any previous studies. The value of = 3.53, in particular, corresponds to the magnetosphere becoming more sensitive to changes in solar wind pressure as it expands, and in fact becoming marginally more compressible than the Jovian system, indicated by the value of < 4 (see discussion in section 1.2).This implies that under appropriate conditions, the compressibility behavior of the Jovian magnetosphere may actually be an intermediate between Earth and Saturn. This is investigated in more detail later in section 3.3.
It was shown in Achilleos et al. [2008], who studied the long-term behavior of the size of Saturn's magnetosphere as measured by Cassini, that the magnetopause radius is well described by a bimodal probability distribution, with local maxima at 22 ± 1.5 R S and 27 ± 1.3 R S (apparently distinct from the typical distributions of the solar wind dynamic pressure). In the earlier empirical studies, the magnetosphere is more often observed in a compressed regime, as shown by the observed ranges in Table 1, due to the more restricted data sets available for use in those studies and the higher average solar activity at the time of observations [e.g., Hathaway, 2015]. Thus, the corresponding estimates of a 1 and are likely to be weighted more toward values typical of such conditions. We also find that our parameter estimates are very sensitive to the choice of hot plasma index K h , which may be a cause of the discrepancy between our results and previous studies. We shall investigate this aspect in depth in the following sections, but first, we will look at which components of the internal structure contribute most to the formation of this "knee" in compressibility behavior at ≈25 R S , as seen in Figure 2. Figure 3 shows the individual contributions from the hot, cold, and magnetic pressure components to the total magnetospheric pressure just inside the model's magnetopause boundary. It can clearly be seen that the magnetic pressure, P B = B 2 ∕2 0 , is the dominant component for all system sizes. Values of plasma have been extracted for the hot and cold plasma populations separately, and the hot plasma beta h and cold plasma beta c are both ≲0.7 across all system sizes, while the total plasma of the collective plasma surpasses unity for the largest system sizes, where R MP ≳ 30 R S . The hot and cold pressure components are comparable for R MP ≳25 R S below, which the cold pressure dominates.
The magnetic pressure profile, in particular, appears to exhibit a change in gradient around 25 R S similar to the solar wind pressure profile and, indeed, has a measured gradient changing from −1∕(4.9 ± 0.3) for the compressed regime to −1∕(4.2 ± 0.2) for the expanded regime, with gradients fitted using the same method as that for the D P profile. These values do not agree exactly with those found for the D P profile regions (see Table 2) due to the minor influence of the plasma pressures, particularly for the expanded regime. However, the magnetic pressure profile shows the same overall trend, and thus suggests that how magnetic pressure varies with system size is a dominant factor in controlling the change in compressibility behavior.
As discussed in section 1, the magnetic field strength (and therefore magnetic pressure) varies more slowly with radial distance for a more disk-like magnetic field structure, i.e., the index in B ∝ r − is smaller than the dipole value. This corresponds to a steeper gradient of the R MP versus magnetic pressure profile on a SORBA ET AL. MODELING MAGNETOSPHERIC COMPRESSIBILITY logarithmic scale, as the gradient of such a profile is −1∕2 . Therefore, the observed behavior of the magnetic pressure profile can be interpreted as the formation of a magnetodisk structure in the magnetosphere, but more so for the more expanded regime. This can be understood theoretically as follows. The magnetodisk forms due to the magnetic tension force increasing, in order to balance the centrifugal force exerted outward on the subcorotating cold plasma in the magnetosphere. This tension force is proportional to B 2 ∕r c , where r c is the radius of curvature of the magnetic field lines and therefore can be increased by a larger magnetic field strength or a smaller radius of curvature. For the compressed regime, B in the outer magnetosphere is still comparatively large and so can maintain a sufficient tension force for a relatively large r c . For the more expanded regime, as B generally decreases, the radius of curvature must also decrease in order to maintain a sufficient tension force, corresponding to the formation of a disk structure. This was also observed empirically by Arridge et al. [2008], who employed Cassini magnetometer data to demonstrate that under strong solar wind pressure conditions, when the magnetosphere was compressed to R MP < 23 R S , the disk structure was effectively destroyed on the dayside, and the magnetic field became quasi-dipolar.
A study by Bunce et al. [2007] used Cassini magnetometer data to construct an empirical model of the ring current in Saturn's middle magnetosphere and used the calculated magnetic fields to estimate the corresponding solar wind dynamic pressures in order to investigate magnetospheric compressibility, similar to our study. They also found a value of that decreased with system size, over a range of R MP = 16 − 30 R S , with overall results in agreement with Arridge et al. [2006]. They explained this in terms of an increase in the ring current magnetic moment for an expanded magnetosphere, which corresponds to an enhancement in the magnetodisk structure. Our preliminary results are therefore in general agreement with this previous study, which suggests that the magnetospheric compressibility is primarily controlled by how the magnetic pressure at the magnetopause boundary varies with system size. However, it is worth noting that the model used by Bunce et al. [2007] assumes a ring current density that falls with cylindrical radial distance as 1∕ ; this is not the case in this study, as the Achilleos et al. [2010a] model calculates azimuthal current profiles directly from the source function g described in section 2.1 [see Caudal, 1986;Achilleos et al., 2010a].
Thus far, we have neglected to include the "pressure" contribution of the centrifugal force exerted on the magnetopause boundary surface and whether this has a role in controlling magnetospheric compressibility. The centrifugal force per unit volume at the magnetopause boundary is given by F c = nm i 2 (see equation (2)). Pilkington et al. [2014] explained that by making the assumption that the magnetopause current layer has a  Achilleos et al. [2008]. The black horizontal bar shows the approximate range of D P typically observed at Saturn (0.008−0.06 nPa) according to Cassini CAPS solar wind data, also from Achilleos et al. [2008]. However, it should be noted that values of D P in the full range shown in this figure (0.001-0.5 nPa) have been observed empirically. comparable density and rotation rate to the plasma just inside the boundary, we can estimate the corresponding pressure contribution by simply multiplying this volume force by the thickness of the magnetopause current layer. By approximating the current layer thickness as 1 R S , following Masters et al. [2011], we find that across all system sizes, the centrifugal term is around an order of magnitude smaller than any other pressure component and therefore is not investigated further in this study.

Influence of Hot Plasma Content
The procedure leading to the results shown in Figure 2, described above, was then repeated using 20 different values of the hot plasma index K h , equally spaced in the range 10 5 -10 7 Pa m T −1 and calculating the model at 20 different sizes in the range R MP = 14-40 R S for each of these K h values. The results are plotted in Figure 4, which shows how magnetopause radius varies with solar wind dynamic pressure on a logarithmic scale, as for Figure 2, with the hot plasma index used in the model represented by the color.
For sufficiently high values of the hot plasma index, the model calculation would not converge to the prescribed 0.5% tolerance for large magnetopause radii, hence the lack of coverage in this region of parameter space. We attempted to mitigate this by, at each iteration, weighting the previous Euler potential solution i−1 (r, ) up to 10 times more heavily than i (r, ) (see section 2.1), such that it would approach a convergent solution more slowly. However, the correspondingly more stringent tolerance required in this case was still not achieved. The reason for this lack of convergence is that in such circumstances an equilibrium force balance solution cannot be found, for a field structure that can contain such a level of simulated hot plasma. It is interesting to note that this only occurs under physical conditions that we would not expect to typically observe in the magnetosphere; for typical values of hot plasma content and solar wind dynamic pressure observed at Saturn, the model calculations are well within the convergence limit, as shown by the black horizontal bar in Figure 4.
The most obvious feature apparent in Figure 4 is the effect of the hot plasma content on the parameter a 1 , reflected in the shift of the y intercept for different compressibility profiles. Estimates of this parameter vary from 10.1 ± 0.2 for K h = 1 × 10 5 Pa m T −1 to 11.3 ± 0.03 for K h = 1 × 10 7 Pa m T −1 , measured by applying a linear least squares regression line to the entire profile for each K h value. This is comparable to the behavior observed empirically by Pilkington et al. [2015] (discussed in section 1.2), who calculated that their estimates of a 1 varied from ≈11 for the data set grouped by low-average plasma ( ≲ 1), to ≈16 for the data set  Table 1, included for comparison and are not associated with any particular K h value.
with high-average ( ≳10). However, it should be noted that, while positively correlated, global hot plasma content K h and local plasma are not interchangeable concepts and have a relationship that is dependent on magnetosphere size, as discussed in detail later.
These observations of variable a 1 correspond to the magnetosphere being observed in a range of sizes for fixed solar wind dynamic pressure, up to 10−15 R S in both Pilkington et al. [2015] and this study. Pilkington et al. [2015] interpreted this observation theoretically as the magnetosphere existing in either a plasma-depleted or plasma-loaded state, transitioning between these states perhaps via Vasyliunas-style reconnection and associated ejection of plasmoids [Vasyliunas, 1983], or interchange events [Mitchell et al., 2015]. If we assume that such transitions can occur at least to some degree independently of the incident solar wind dynamic pressure, then it seems intuitive that there should exist a range of possible system sizes for a fixed D P , but different K h . This can be explained as follows. Consider a magnetopause initially in dynamic equilibrium with the upstream solar wind dynamic pressure and what happens when the plasma pressure inside the magnetopause boundary drops due to some plasma loss process. The pressure across the boundary is now unbalanced, with the incident D P greater than the total plasma and magnetic pressure inside the magnetosphere. The magnetopause location therefore moves inward and the magnetosphere is compressed, causing an enhancement in the magnetic field strength B. The magnetic pressure therefore increases inside the magnetopause boundary, until pressure balance is restored and a new equilibrium magnetopause location is found. We would also expect the remaining magnetospheric plasma pressure to, in general, increase as the magnetopause moves inward. This scenario thus corresponds to the magnetosphere generally being observed at a smaller size when the plasma content or pressure is lower, resulting in a lower estimate for the parameter a 1 .
The variation in compressibility behavior with varying system size and hot plasma content shown in Figure 4 is less intuitive. As with our initial results, there appears to be a shift in behavior at around 25 R S for all profiles, although becoming less pronounced as K h increases. Therefore, for each individual profile at a given K h value, the simulated data were again divided into two groups, a compressed regime and an expanded regime, separated at 25 R S , and we applied linear fits to these regions separately to make estimates of . The resulting estimates are shown in Figure 5. Groups with fewer than three individual data points were not fit; hence, for the expanded regime the estimates do not cover the full range of K h . As with the initial results, it can be seen that across the full range of K h , the expanded regime gives lower estimates of than the compressed regime, corresponding to a magnetosphere that is more sensitive to changes in solar wind pressure as it expands. We suggested in the previous section that this was due to the formation of a significant magnetodisk structure, which is more easily compressible than a dipolar magnetic field structure, only for an expanded magnetosphere. It is also apparent in Figure 5 that there is a general trend of decreasing as hot plasma content increases, corresponding to a more easily compressible magnetosphere. (Note that for K h > 4× 10 6 Pa m T −1 for the expanded regime, results may be affected by applying a linear fit to so few data points.) Achilleos et al. [2010b] found in a previous study using this same model that an increase in hot plasma content did significantly affect the magnetic field configuration, causing a more disk-like field structure and demonstrated that there was reasonable agreement with limited Cassini observations. Therefore, a similar argument as to how the compressibility changes as the magnetosphere expands, can be employed to explain how the compressibility changes as K h increases. This effect can be interpreted as the increased hot plasma pressure augmenting the azimuthal ring current in the magnetosphere, which in turn enhances the associated disk-like magnetic field.
However, this is not the full picture. Figure 6 shows a contour plot of how the value of h at the nose of the magnetosphere, just inside the magnetopause boundary, varies with both magnetopause radius and hot plasma index. Values are extracted from all models where convergence was obtained. For the initial conditions described in this study, K h = 1 × 10 6 Pa m T −1 , we found that h never exceeded ≈0.7 for any system size, and thus, the variation in magnetic pressure with radial distance always controlled the compressibility behavior. However, we can see that in the extremes of allowed parameter space, where either K h or R MP are sufficiently large, h > 1, and therefore, the variation in hot plasma pressure with radial distance also becomes important in controlling the compressibility behavior. These conditions correspond directly to regions where is lower, suggesting that a sufficiently high hot plasma pressure makes the magnetosphere more easily compressible. For values of K h near the upper limits of empirical observations, approaching 10 7 Pa m T −1 , it can be seen in Figure 6 that exceeds unity even for the smallest system sizes, and indeed it is in this region of parameter space that the smallest estimates of are obtained. However, it was found by Pilkington et al. [2015] that such conditions of a compressed magnetosphere with high hot plasma content, were rarely observed empirically, which may partially explain the discrepancy between previous results and our lowest estimates, shown in Figure 5. This can also be seen in Figure 4, which shows that the solar wind dynamic pressure typically observed at Saturn is not sufficiently high to support such compressed magnetospheres with high hot plasma content. (Note that cold plasma beta c was never found to exceed 0.7 anywhere in the allowed parameter space and thus is not discussed further.) Figure 7. The shape of the outermost closed field line in a magnetosphere as the magnetopause radius R MP changes, for three different values of K h corresponding to very quiet, average, and disturbed ring currents, respectively. In each panel the shape of the field line has been normalized in the vertical and horizontal direction by the magnetopause radius R MP , the value of which is shown by the color of the field line according to the color bar. A representative dipolar magnetic field line is shown in black on each plot for comparison.
For a dipolar magnetic field with isothermal plasma transport, it was explained in section 1 that we would expect P h ∝ r −4 , thus giving ≈ 4 for a magnetosphere with compressibility controlled by hot plasma content. However, Figure 5 shows that we find < 4 in the most extreme cases. Indeed, when fitting the profiles of hot plasma pressure specifically for each K h value, it was found that for these model calculations, the behavior varied from P h ∝ r −3.3 for the smallest hot plasma index to ∝ r −2.6 for the greatest, corresponding to a reduction in compressibility parameter . It was also found that, unlike the magnetic pressure profiles, there was no considerable shift in behavior at 25 R S , and so the observed kink in Figure 6 in this region is solely due to change in magnetic pressure, the denominator of h .
This behavior of the hot plasma pressure can be further understood as follows. The parameterization of hot plasma adopted in this model means that P h is fully determined by how the flux tube volume varies with system size, as P h V = K h . The flux tube volume is defined in equation (6) and is thus dependent on both the length of a given field line (ds) and the magnetic field strength (B). For a dipolar field line, B ∝ r −3 and ∫ ds ∝ r, hence V ∝ r 4 and P h ∝ r −4 in our parameterization. However, in this model, as we have discussed, we have found that the magnetic field strength varies more slowly with radial distance than this, particularly in expanded regimes and with high hot plasma content. Indeed, the behavior varies from B ∝ r −2.7 for compressed regimes with low hot plasma content, to ∝ r −1.8 for expanded regimes with high hot plasma content, corresponding to a more significant magnetodisk field. This in turn affects how flux tube volume varies with system size, meaning P h varies more slowly with system size.
In addition, we also observed that for more expanded systems, the length of the outermost magnetic field lines ∫ ds varied more slowly with radial distance than for a dipolar magnetic field, with ∫ ds ∝ r 0.9 for magnetospheres with the lowest K h values, to ∝ r 0.8 for the highest K h values. For more expanded, stretched magnetospheres, it is perhaps intuitive that the magnetic field lines in the outer regions have shorter overall lengths than a corresponding dipolar magnetic field line that crosses the equator at the same radial distance, due to the outward radial stretching of the dipole magnetic field associated with the ring current. Figure 7 illustrates this oblateness of the magnetospheric magnetic field lines observed in our model calculations and how it increases with system size and hot plasma content. Each panel illustrates how the shape of the outermost closed magnetic field line, which crosses the equator at the magnetopause nose, varies with the size of the magnetosphere, for a given K h value. The shape is shown in height above the rotational equator Z and cylindrical radial distance from the planet center , normalized to the magnetopause radius R MP , with a representative dipolar magnetic field line shown in black on each plot for comparison.
It can be seen clearly that as the magnetosphere size increases, the outermost closed field line is confined comparatively much more toward the equator than for the magnetospheres with smaller R MP . This is especially pronounced for the magnetospheres with higher hot plasma content, where the field lines show significant oblateness particularly at lower . This effect can be interpreted theoretically as follows. Considering pressure balance perpendicular to a magnetic field line in the outer magnetosphere, the two dominant forces to consider are the hot plasma pressure gradient and the centrifugal force. The force associated with the plasma pressure gradient acts outward from the center of curvature of the field line, perpendicular to the field line, across the entire field line length. This property arises as a consequence of plasma pressure being assumed constant along a field line (see section 2.1). In contrast, the centrifugal force acts in the direction of increasing . Therefore, the component of centrifugal force perpendicular to the magnetic field line can either act outward away from or inward toward the equatorial plane depending on the geometry, acting inward for smaller and outward beyond the "turning" point of the magnetic field line, when the field line begins to converge back toward the equator. This turning point occurs about halfway along the magnetosphere, at around 0.5 ∕R MP , for the compressed, low K h cases, but as far out as 0.7 ∕R MP for the most expanded and hottest cases. The farther out in radial distance that this turning point occurs, the higher the fraction of the magnetic field line for which the perpendicular component of centrifugal force is acting inward, toward the equatorial plane-and thus effectively balancing the outward hot plasma pressure gradient force. Therefore, in order to maintain global force balance, this turning point moves radially outward as the hot plasma pressure gradient increases, and the field line thus becomes more confined.
We explained in section 2.1 that the parameterization of hot plasma content via the state equation involving the index K h = P h V was a simplifying assumption in light of variable observations of hot plasma pressure. However, it would have also been possible to instead parameterize the hot plasma content by, for example, K h = P h V , with = 5∕3, such that magnetic flux tubes of plasma are considered to expand and contract adiabatically rather than isothermally. The factors of B and ∫ ds that determine the flux tube volume V would be unchanged; however, the previous relationship P h ∝ V −1 would be modified to P h ∝ V − . For a magnetosphere with compressibility controlled by hot plasma content, this corresponds to an increase in the estimate of the compressibility parameter by a factor of and thus a magnetosphere that is less sensitive to changes in solar wind dynamic pressure. This aspect may also contribute to the discrepancy between our results and those of previous empirical studies, shown in Figure 5, as the plasma population may behave intermediately between the two state equations discussed in this study. Indeed, the analysis of Achilleos et al. [2010a] supports this conclusion.

Comparison With the Jovian System
It is insightful to make a comparison between the results presented here and the corresponding calculations for Jupiter's magnetosphere, using our implementation of the model of Caudal [1986] directly. Figure 8 shows the results for Jovian calculations analogous to Figure 2 for the Saturn case. Jovian magnetopause radii in the range R MP = 50 − 100 R J have been used to cover the observed range presented in previous studies [Joy et al., 2002].
It can clearly be seen that, unlike the Saturn case, the Jovian magnetosphere shows a compressibility behavior that is well represented by a uniform across all observed system sizes, demonstrated by a constant gradient of the D P profile. A single linear least squares regression line fit provided an estimate for the compressibility parameter = 3.05 ± 0.02. As for the Saturn case, the data were split into two groups representing a compressed and expanded regime, and a linear least squares regression line was fit to each group separately. However, this did not provide significantly different estimates for , with a variation of only approximately 6% between groups, and so only the fitting for the entire simulated data set is shown and discussed here. This value of = 3.05 ± 0.02 is remarkably smaller than observational studies in the literature, which estimate this value as between ∼4 and ∼5 [Huddleston et al., 1998;Joy et al., 2002;Alexeev and Belenkaya, 2005]. Possible reasons for this discrepancy are discussed at the end of this section. However, it can be seen from a comparison with Figure 5 that this value is consistent with our original theoretical expectation that this value be lower than that of the Saturn system, due to arguments discussed in section 1.
It is worth noting that this difference in behavior between the Saturnian and Jovian magnetospheres is also a consequence of their relative locations in the solar system and the weaker solar wind dynamic pressure experienced at Saturn. If Saturn were to be located closer to the Sun such that it experienced the higher solar wind dynamic pressures typically observed at Jupiter, then the Saturnian magnetosphere would also show a compressibility behavior represented by a uniform , specifically corresponding to its compressed regime. This can be seen by comparing the ranges in solar wind dynamic pressure in Figures 2 and 8. We compared the individual P h , P c and P B components that comprise the D P estimate across all system sizes and found that, as with the initial results for the Saturn case presented in Figure 3, the magnetic pressure component was dominant across all system sizes. Both hot and cold plasma monotonically increased with system size, with c ≤ 0.3 and h ≤ 0.8 across the simulated data set. Unlike the Saturn case, there was no evidence of a shift in gradient of the magnetic pressure (P B ) profile for increased values of R MP ; we measured a constant slope equivalent to P B ∝ r −3.4 , corresponding to B ∝ r −1.7 . This index is smaller than was measured even in the hottest and largest magnetosphere models in our Saturn investigation and corresponds to a very significant disk-like distortion from a dipolar magnetic field for all Jupiter R MP . Comparable behavior is expected theoretically as discussed in section 1, and has also been observed empirically [e.g., Khurana and Kivelson, 1989]. We also found that the hot plasma pressure P h varied as P h ∝ r −2.7 , comparable to the behavior measured for the hottest and largest magnetosphere models in our Saturn investigation. This can readily be explained via the same arguments applied to the Saturn case in the previous section.
As with our initial Saturn investigation, we estimated the pressure contribution from the centrifugal force at the magnetopause boundary, using a very approximate value for the magnetopause boundary layer thickness of 1R J following Delamere and Bagenal [2010]. We found that the centrifugal term was at least an order of magnitude smaller than all other pressure components at the magnetopause boundary across all system sizes, and therefore does not play a significant role in determining compressibility. However, in a study by Nichols [2011], it was suggested that this Jovian model may significantly underestimate the magnitude of centrifugal force in the middle and outer magnetosphere, due to the use of a plasma angular velocity profile that overestimates the breakdown of corotation. This would mean that the model would also overestimate how strongly the centrifugal force falls with radial distance. In combination, these effects may contribute to the discrepancy between our particularly low estimate of the compressibility parameter for the Jovian magnetosphere and previous empirical studies. We intend to investigate this aspect further in future studies, as well as the influence of different hot plasma parameterizations, using more recent spacecraft data sets than the Voyager results employed in the Caudal [1986] model.

Summary and Conclusions
We have employed a 2-D force balance model of Saturn's dayside magnetosphere, first described in Achilleos et al. [2010a], to investigate magnetospheric compressibility, and in particular its response to solar wind conditions and global hot plasma content. We have found that, for average global hot plasma conditions, 10.1002/2016JA023544 the compressibility behavior can be described by equation (1) but with a value of compressibility parameter that decreases with system size, from ∼4.8 for R MP ≤ 25R S to ∼3.5 for R MP > 25R S . This corresponds to the magnetosphere becoming more easily compressible as the upstream solar wind dynamic pressure decreases. We have explained this in terms of the distortion of the magnetic field into a magnetodisk configuration, which is more easily compressible than a dipolar magnetic field and that this distortion becomes significant for more expanded magnetospheres.
We have shown that the global hot plasma content of the magnetosphere, parameterized by the hot plasma index K h = P h V, also plays an important role in determining magnetospheric compressibility. When Saturn's magnetosphere is compressed, a higher value of K h acts to increase the observed magnetospheric compressibility via an enhancement of the aforementioned magnetodisk magnetic field structure, due to its contribution to the associated ring current magnetic field. When the magnetosphere is more expanded, the hot plasma pressure exceeds the magnetic pressure at the magnetopause boundary such that the variation in hot plasma pressure with radial distance becomes significant in controlling magnetospheric compressibility. We have determined that, as the hot plasma pressure P h varies more slowly with radial distance than the magnetic pressure P B , this corresponds to the magnetosphere becoming more easily compressible under such conditions. We also explored the behavior of the P B and P h profiles, and how they are related via the flux tube volume V. Our estimates of are in the range 3.3-5.3 depending on system size and hot plasma content, with lower estimates corresponding to higher values of both of these parameters. However, as we have mentioned, these results are sensitive to our simplifying parameterization of the hot plasma content and, for example, would increase up to a factor of 5∕3 (for regions of parameter space where hot pressure is dominant) if we were instead to parameterize the hot plasma content via K h = P h V .
These results thus suggest that future observational studies of the relationship between R MP and D P may benefit from the assumption that the compressibility parameter is not constant across all observations but varies both with system size and hot plasma content, at comparable levels of significance. In addition, our model analysis demonstrates that centrifugal force at the magnetopause boundary does not contribute significantly to the compressibility of the magnetosphere. On the other hand, global centrifugal force is important in determining the global field structure, and thus compressibility. Improvements are needed in our treatment of the plasma angular velocity profile, in order to make it self-consistent with the changing magnetic field structure so as to further test our initial findings.

Appendix A: Fitting Procedure and Error Estimation
If we have a set of n data points {x i , y i } that we have reason to believe are linearly dependent, we can attempt to find a "best fit line" to describe this relationship,ŷ(x) = mx + b. In our study, this corresponds to the logarithm of equation (1), such thatŷ = log R MP , x = log D P , m = −1∕ , and b = log a 1 . We can then obtain estimates of parameters a 1 and by fitting the logarithms of our simulated data set, {x i = log D P , y i = log R MP } with a best fit line and estimating m and b.
In order to obtain estimates of the parameters m and b, we minimize the sum of the squares of the residuals between the data and the best fit line points, with respect to m and b. This yields the solutions wherex denotes the mean of all x i and similarly forȳ, and all summations are calculated over all n data points. This technique assumes that the errors associated with the data are uniform such that all {x i , y i } are