Magnetotail Reconnection Asymmetries in an Ion-Scale, Earth-Like Magnetosphere

. We use a newly developed global Hall MHD code to investigate how reconnection drives magnetotail asymmetries in small, ion-scale magnetospheres. Here, we consider a magnetosphere with a similar aspect ratio to Earth but with the ion inertial length ( δ i ) artiﬁcially inﬂated by a factor of 70: δ i is set to the length of the planetary radius. This results in a magnetotail width on the order of 30 δ i , slightly smaller than Mercury’s tail and much smaller than Earth’s with respect to δ i . At this small size, we ﬁnd that the Hall effect has signiﬁcant impact on the global ﬂow pattern, changing from a symmetric, Dungey-5 like convection under resistive MHD to an asymmetric pattern similar to that found in previous Hall MHD simulations of Ganymede’s subsonic magnetosphere as well as other simulations of Mercury’s using multi-ﬂuid or embedded kinetic physics. We demonstrate that the Hall effect is sufﬁcient to induce a dawnward asymmetry in observed dipolarization front locations and ﬁnd quasi-periodic global scale dipolarizations under steady, southward solar wind conditions. On average, we ﬁnd a thinner current sheet dawnward; however, the measured thickness oscillates with the dipolarization cycle. During the ﬂux-pileup stage, 10 the dawnward current sheet can be thicker than the duskward sheet. This could be an explanation for recent observations that suggest Mercury’s current sheet is actually thicker on the duskside: a sampling bias due to a longer-lasting “thick” state in the sheet.


Introduction
In the magnetospheres of Mercury and Earth, observations of plasmoids, flux bundles, and dipolarization fronts (DFs) demonstrate a marked asymmetry in their distribution across the magnetotail. At Earth, a number of studies have found magnetotail duskward biases in several magnetic phenomena: flux rope occurrence (Slavin et al., 2005;Imber et al., 2011), dipolarization fronts (Liu et al., 2013), energetic particle injections (Gabrielse et al., 2014), and reconnection (e.g., Asano et al., 2004;Genestreti et al., 2014). Additionally, the current sheet was found to be thinner on the duskside (Artemyev et al., 2011;Vasko et al., 2015). Similarly, at Mercury, Poh et al. (2017b) used MESSENGER data to fit the Harris sheet model to 234 tail current sheet crossings and found a bias towards dusk having thinner current sheets (by ≈ 10 %-30 %). In contrast, however, other MESSENGER studies (Sun et al., 2016;Dewey et al., 2018) found dawnward biases in dipolarization events and reconnection front locations.
The general existence of tail asymmetry is thought to be a result of sub-ion-scale effects (Lu et al., 2018;Liu et al., 2019), though there is still some uncertainty about the exact manifestation and causes of specific asymmetries. It is debated whether Hall electric fields are sufficient to reproduce this or if other ion-/electron-scale physics are required. Although some authors argue that electron-scale physics is required (Chen et al., 2019), we show in this paper that Hall effects are sufficient to cause an asymmetry in some observed features. Furthermore, it is unknown exactly why Mercury and Earth observe different asymmetries; it is hypothesized that system-size effects (relative to the ion inertial length δ i ) play a key role (Lu et al., 2016(Lu et al., , 2018Liu et al., 2019).
Several studies have proposed mechanisms to explain how Hall reconnection induces asymmetry in the magnetotail. Lu et al. (2016Lu et al. ( , 2018 (hereafter Lu+), in studying Earth's magnetotail with global hybrid simulations and localized particle-in-cell (PIC) simulations, showed that the decou-Published by Copernicus Publications on behalf of the European Geosciences Union. 992 C. M. Bard and J. C. Dorelli: Tail asymmetries pling of ions and electrons within the current sheet (the Hall effect; e.g., Sonnerup, 1979) creates an electric field and associated tail current density. The resulting E × B drift is sufficient to create tail asymmetries and indeed may be the primary cause. The duskside magnetic flux is preferentially evacuated via electron transport dawnward, which leads to a smaller normal B z and thinner current sheet on the duskside.
In a similar study, Liu et al. (2019) (hereafter Liu+), using local PIC simulations of embedded, thin current sheets, confirmed that the Hall effect creates electron E × B and diamagnetic drifts which transport magnetic flux dawnward within the current sheet. However, they found that, although the preexisting tail B z initially suppresses the onset of dawnside reconnection, the reconnection B z drives outflows towards dawn and thins out the current sheet on that side. This creates an "active region" of reconnection on the dawnside, which has a thinner current sheet and stronger tail current j y . After analyzing both these studies, Liu+ proposed that, although the Lu+ model provides a explanation for a duskward bias in the initial reconnection onset, the Liu+ active region provides an explanation for dawnward biases within local, inprogress magnetotail reconnection. We test several aspects of this general picture within this paper.
Unfortunately, simulating large magnetospheres such as Earth (few hundred δ i ) while properly resolving the smallscale Hall physics requires grid sizes in the billions of cells. Several strategies have been proposed to evade this constraint; one is to embed regions of detailed kinetic physics within large-scale Hall magnetohydrodynamic (MHD) simulations (Chen et al., 2019). This allows for reproduction of kinetic effects within certain regions of the magnetosphere without having to run an expensive, fully kinetic simulation. However, these simulations assume no kinetic effects outside the embedded regions, which are limited to certain regions in the dayside and/or the tail. It is unclear whether or not this methodology, including boundary handling, affects the local-global feedback dynamics in the magnetosphere. Future studies will eventually be needed to compare magnetospheres from Hall MHD, kinetic, and combined kinetic-Hall simulations to ascertain these effects.
Another strategy suggests that we need only set the Hall scale to some length sufficient to capture the essential physics of Hall reconnection without having to fully resolve the physical length scale. In these simulations, the Hall length is set to ≈ 3 % of the global-scale length (Tóth et al., 2017), which is sufficient to capture the out-of-plane flows and the quadrupolar magnetic field structure induced by the Hall effect. However, recent research in 2D island coalescence (Bard and Dorelli, 2018) suggests that although including the Hall term in MHD simulations is sufficient in itself to generate these signatures of Hall reconnection, the actual reconnection rate depends on resolution and numerical resistivity. Although the Hall term is present, the reconnection itself may be Sweet-Parker-like and slow (unlike fast Hall reconnection). Bard and Dorelli (2018) observed that 20-25 cells per δ i was necessary (within the context of their numerical viscosity) in order to observe fast Hall reconnection. This is much greater than the 5-10 cells per δ i typically used in simulations (Dorelli et al., 2015;Dong et al., 2019;Chen et al., 2019). This suggests that, although artificially inflating δ i allows the Hall effect to emerge and have a global impact, much higher resolution is required to observe the universally fast (∼ 0.1 v A ) reconnection observed in kinetic simulations. Finally, Bard and Dorelli (2018) found qualitatively different behavior for varying ratios of system size to δ i : large systems can produce bursty reconnection (with a low average reconnection rate) even when δ i is sufficiently resolved to produce "fast", instantaneous reconnection. Ultimately, the combined requirements of high resolution and large system size create a computational requirement beyond what is currently possible for magnetospheres. Indeed, we are only setting a resolution appropriate to 5 cells per solar wind δ i in this work, though local density fluctuations in the tail may allow up to 10-20 cells per δ i .
One possible method for dealing with this issue may be to use graphics processing units (GPUs), which have proven to be robust and viable for scientific computing. Indeed, several groups have already utilized GPUs to accelerate plasma simulations throughout heliophysics, astrophysics, and plasma physics (Bard and Dorelli, 2014;Benítez-Llambay and Masset, 2016;Fatemi et al., 2017;Bard and Dorelli, 2018;Schive et al., 2018;Grete et al., 2019;Liska et al., 2019;Wang et al., 2019). GPUs take advantage of parallelism in order to have a higher throughput for floating point operations. Finitevolume schemes are massively parallel: the calculation of how a computational cell evolves from t to t + t is independent of similar calculations for other cells. This makes explicit Hall MHD schemes (such as presented in this paper) quite amenable to GPU acceleration.
In this paper, we undertake a numerical experiment designed to assess the role of the Hall effect on global magnetospheric structure and dynamics within a "small" ionscale magnetosphere, specifically focusing on how it induces asymmetry in the magnetotail. We present a magnetosphere simulation code which accelerates the explicit MHD solver algorithm via GPUs. We simulate an Earth-like analogue magnetosphere which has a similar bow shockmagnetopause distance and magnetotail width as Earth's (relative to the planetary radius); however, δ i is artificially inflated to the planetary radius (R E ). In other words, we are self-similarly scaling Earth from its current size relative to ion scales (magnetotail width ≈ 600δ i ) to a size closer to the ion scale (magnetotail width ≈ 30δ i ). In this "ion-scale Earth", Hall physics plays a greater role in magnetosphere dynamics, and we are able to more readily observe global Hall MHD effects. We view this work as a first step in the study of the system-size dependence of Hall MHD magnetic reconnection in Earth-like magnetospheres; future systemsize studies can be performed by making δ i smaller relative to the planetary radius and increasing the resolution to sufficiently cover the ion scales. This paper is presented as follows: Sect. 2 provides a brief overview of the Hall MHD algorithm as implemented using GPUs; Sect. 3 provides the initial condition and setup of the simulation; and Sect. 4 presents tail asymmetries in the simulation and discusses them in the context of observations and proposed theoretical explanations.

Methods and code
We take a Hall MHD code accelerated by GPUs using the MPI library and the NVIDIA CUDA API (Bard and Dorelli, 2014;Bard, 2016;Bard and Dorelli, 2018) and adapt it to simulate planetary magnetospheres. We review the underlying mathematical equations and algorithms in this section.
Following Tanaka (1994) and Powell et al. (1999), we split the magnetic field vector B into a background component B 0 and a perturbed, evolving component B 1 such that B = B 1 + B 0 . The embedded B 0 is assumed to be static (∂B 0 /∂t = 0), divergence-free (∇ · B 0 = 0), and curl-free (∇ × B 0 = 0). This allows for more accurate handling of the magnetic field, especially near the planet where the dipole field is very strong. In order to preserve the divergence-free constraint on the evolved magnetic field, we solve the Generalized Lagrangian Multiplier (GLM) formulation of MHD (Dedner et al., 2002), with an additional Hall term added via Ohm's law.
The ideal MHD Ohm's law is extended with the Hall term such that the electric field E is given by with c the speed of light, e the elementary charge, n the plasma number density, v the plasma bulk velocity vector, and J = c 4π ∇ × B the current density vector. We note that since the background magnetic field B 0 is curl-free in our formulation, the current density is taken to be the curl of the perturbation B 1 .
We normalize the density (ρ), magnetic field, and length scale to reference values ρ 0 , B w , and L 0 , respectively. v is can be used in place of setting P 0 directly; the conversion between the two is β 0 = 2 P 0 .
This results in the set of equations: ∂ρv ∂t where E = ρv 2 /2 + p/(γ − 1) + B 2 1 /2 is the total energy density, γ is the ratio of specific heats (taken to be 5/3 in all of our simulations), n 0 e 2 is normalized to the reference length such that δ i = δ i /L 0 . The normalized δ i in our simulation is a fixed parameter that can be changed at runtime. We evaluate the normalized current density (J = ∇ × B 1 ) at cell centers and linearly interpolate to the cell edges when needed.
ψ is a scalar function whose evolution is designed to be equivalent to ∇ · B; c h and c p are parameters for the propagation and dissipation of local B divergence errors, respectively. Following Dedner et al. (2002), we set c h as the global maximum wave speed over the individual cells and set c p such that c 2 p /c h = 0.36. Although Dedner et al. (2002) recommended c 2 p /c h = 0.18, and this value works very well to control the magnetic divergence in non-magnetospheric simulations, we find that some level of tweaking is required because of the accumulation of divergence errors at the inner boundary. To ameliorate further complications caused by this issue, we separate the momentum equation into a nonmagnetic flux and a magnetic source term: which prevents divergence errors from inducing a nonphysical acceleration along magnetic field lines (Brackbill and Barnes, 1980) but with some loss of accuracy across shock jumps.
The overall system is evolved via a time-explicit secondorder Runge-Kutta scheme coupled with a simple HLL Riemann solver (Harten et al., 1983;Toro, 1999) and a monotonized central limiter (e.g., Tóth et al., 2008) with the slopelimiting parameter β set to 1.25. For numerical stability, the explicit time step is determined by the global Courant condition: t = C x min /v max . The Courant parameter is C = 0.325, x min = 0.2 R E is the smallest cell length in the simulation, and v max is the maximum wave speed in the simulation. For Hall MHD, v max,i = |v|+|v f |+|v w | is estimated in each grid cell i using the fast magnetosonic (v f ) and whistler wave (v w ) speeds (Huba, 2003;Tóth et al., 2008): with v 2 A = B 2 /ρ the normalized Alfvén speed and v 2 s = γ P /ρ the normalized sound speed. The highest value of v max,i across all cells is used to set the global time step.

Problem initialization
For our simulation setup, we choose normalized solar wind and terrestrial magnetic field parameters such that the magnetopause standoff distance matches that of Earth's magnetosphere (≈ 10 R E ), the bow shock standoff distance is ≈ 3 R E beyond the magnetopause, and δ i is equivalent to the planetary radius.
These relative distances can be controlled by setting four dimensionless parameters 1 : 2. β sw = 8π P sw /B 2 sw , the solar wind plasma beta; 3. B 0 / B sw , the ratio of the dipole field strength at 1 R E to the solar wind magnetic field; 4. δ i , the ion inertial length.
Since only the first three of these parameters control the magnetopause and bow shock standoff distances, we can arbitrarily set δ i to control the relative Hall scale. This allows our theoretical magnetosphere to be simultaneously "Earth-like" (in a relative physical sense) and ion-scale (relative to δ i ). Thus, we choose our reference values as follows: L 0 = R E = 6.371 × 10 8 cm = 6371 km, n 0 = 5 cm −3 , ρ 0 = m i n 0 with the mass of each ion m i = 6.546 × 10 −21 g = 3942.18 amu, and B w = 10 −4 G = 10 nT. From these, we derive v 0 = 97.5 × 10 5 cm s −1 = 97.5 km s −1 , t 0 ≈ 65 s, P 0 = 7.96 × 10 −10 Ba = 0.0796 nPa, and δ i = 6371 km = 1 L 0 = 1 R E .
The solar wind is initialized with values ρ sw = 1 ρ 0 , v sw = 4.09 v 0 = 400 km s −1 , and wind plasma β sw = 0.305 such that P sw = 0.1526 P 0 . The wind magnetic field is initially set to B sw = (-0.174, 0, 0.985)B w for a northward interplanetary magnetic field (IMF) with magnitude B sw = 1 B w ; we later flip the IMF by setting B z = −0.985.
The planetary background magnetic field (B 0 ) is approximated with a magnetic dipole with r the position vector from the center at (0, 0, 0) and the planetary dipole moment m = (M x , M y , M z ) taken as (0, 0, −3000), such that B 0 = 3000 B w = 0.3 G on the magnetic equator at r = 1 R E . This satisfies the requirement that B 0 be both curl-free and divergence-free. As a summary, the dimensionless parameters for this experimental magnetosphere are M A = 4.09; β sw = 0.305; B 0 / B sw = 3000; δ i = 1. For Earth, δ i ≈ 1/70 and the other parameters would be the same. For Mercury, the dimensionless parameters would be, e.g., M A = 6.; β sw = 0.65; B 0 / B sw = 15; δ i = 1/65.
At the inner boundary, we had difficulties with density depletion causing large local Alfvén speeds and small global time steps. We tried several types of inner boundaries, experimenting with different combinations of floating (zerogradient) and fixing the plasma variables. Ultimately, although the following inner boundary conditions are not entirely realistic, they allow for a stable evolution of the magnetosphere in both the dayside and the tail. We set the inner boundary at a radius of 3 R E ; in these ghost cells, we fix the density at 4 ρ 0 , float the pressure, float the radial magnetic field, set the tangential B to zero, and set the velocity to zero. For the divergence cleaning, we find that simply setting the ghost ψ c = 0 works better than having a floating condition. We note that in more realistic magnetospheres, cold plasma from the ionosphere may flow out to the tail and impact the dynamics. We will leave this topic for future studies.
For the outer boundaries, the left edge of the simulation domain fixes the conservative variables to the background solar wind condition; the rest of the box has zero-derivative boundaries for all variables.
The simulation coordinates are defined with −x pointing towards the Sun, z along the planetary magnetic dipole axis, and +y (−y) towards the dusk (dawn) completing the orthogonal set. Although the planet does not rotate, dawn and dusk are used assuming the sun rises in the east. In order to resolve the artificially inflated δ i , we choose 5 cells per δ i , giving a minimum resolution of x, y, z = 0.2 R E . This resolution is set within the range −20 R E < x < 20 R E ; −15 R E < y, z < 15 R E ; beyond this the cell length increases by 7 % with each additional cell up to a maximum of 5 R E or until it hits the boundary. The total size of the grid is 290 × 253 × 253 (just over 18 million cells), covering a domain of We start the simulation in ideal MHD (δ i = 0) with a northward IMF (B sw given above) for 120 t 0 and then flip B z,sw for the southward IMF case and run it for another 120 t 0 . At this point, we turn on the Hall term by setting δ i = 1 and run it for another 12 t 0 in order to allow the perturbations induced by the abrupt change of physics to settle. From this point on, the simulation was run for 45 t 0 under continuous pure southward IMF and with the Hall term on. The main results discussed below come from this final portion of the run with the Hall term enabled.

Hall-induced asymmetry
Prior to turning on the Hall term, the magnetospheric convection is of Dungey-type. Turning the Hall term on, however, induces an out-of-reconnection-plane E × B force which breaks that symmetry and drives convection in a preferred direction (Fig. 1). This convection also modifies the current sheet structure between ideal and Hall MHD (Fig. 2) and causes it to vary in length along the x− direction across the tail y coordinate. For smaller magnetospheres, this effect was first seen in nonideal MHD simulations of Ganymede (Dorelli et al., 2015;Tóth et al., 2016;Wang et al., 2018); this was later seen in 10-moment and embedded kinetic sim-ulations of Mercury (Dong et al., 2019;Chen et al., 2019). Our simulation supports the idea that this Hall-induced drift is sufficient to produce an asymmetry: kinetic effects are not required, but they may manifest different kinds of asymmetry.
Finally, since our choice of min x, y, z = 0.2 R E is based on the normalized δ i = 1 R E , we check here that we are able to correctly resolve δ i in the tail. As Fig. 3 illustrates, the local value of δ i exceeds 1 R E in the magnetotail and may reach up to 5 R E . Thus, we are appropriately resolving the Hall length scales in the tail.      . Best-fit current sheet half-thicknesses (L CS ) derived by fitting Eq. (11) to 5037 cuts of B x along the z direction. These cuts were randomly sampled in the tail xy plane and over the simulation time period (see text). There is a bias towards the current sheet being thinner on the dawnside. However, the dawnside also sees a larger spread in thicknesses: this is a result of temporal effects (see main text for discussion).

Dipolarizations
In our simulation, the Hall electric field induced by tail reconnection accelerates ions towards the duskside and the electrons towards dawn. Since δ i = R E here, the reconnection current sheet spans a significant fraction across the tail; this means that the ions are decoupled from the magnetic field during much of their in-plane convection duskward (green arrows in Fig. 4). The electrons, being coupled to the magnetic field, carry the reconnected, normal B z flux dawnward (blue arrows). Because the reconnected magnetic flux origi-nates over a large region within the tail, there is a significant pileup leading to a reconnecting, active region of plasmoid formation on the dawnside. This pileup-reconnection mechanism may be a general cause of dipolarizations in ion-scale magnetospheres (like Mercury, e.g., Sundberg et al., 2012). Although the drifting electrons themselves do not cause dipolarization events, they can affect where events occur. We speculate that increasing the system size relative to δ i will limit the extent of dawnward flux transport via current sheet electrons and cause the flux pileup to shift closer to tail center (or duskward). In other words, this type of asymmetry we observe here may be more pronounced in δ i -scale magnetospheres and weaker in larger magnetospheres. However, it is unknown how this mechanism will interact with the overall Hall-induced duskward ion convection, which will begin to dominate electron convection at scales δ i .
During the 45 t 0 ≈ 48 min duration of our simulation, there were seven events visually observed on the dawnside (none on the duskside) which followed the general substorm pattern of a buildup/loading phase followed by an unloading (or expansion/relaxation) phase (Rostoker et al., 1980). For each event, we observed pileup of the normal B z magnetic flux over a period of several minutes, followed by a burst of reconnection and the subsequent ejection of plasmoids tailward (Figs. 5 and 6). Three of the seven events produced large plasmoids (on the order of 10 R E = 10δ i ), while the rest resulted in smaller ones (≤ 5 R E ; ≤ 5δ i ). The larger ejecta appeared to build up and release on timescales of around 10 min, while the smaller events had shorter timescales of around 5 min. Most events originated at a down-tail distance ≈ 13-16 R E ; after ejection, their resulting plasmoids traveled to about 30 R E down-tail over several minutes before dissipating.
The observed dawnward bias in dipolarization events for our ion-scale magnetosphere corroborates similar dawnward Figure 9. Cross-comparison of current sheet density magnitude (a, d), current sheet B z flux pileup (b, e; same parameters as Fig. 5) and sampled thicknesses (c, f) during (a, b, c) and after (d, e, f) a global dipolarization event. Current sheet fits are sampled from the area within the wedges (13 R E < R < 17 R E ). The current sheet is thick where the B z flux has piled up and thin where the flux has been unloaded. There is a 2.25 t 0 ≈ 2.4 min time difference between the snapshots of the top and bottom rows. Figure 10. 1D cuts of tail magnetic field B z (a) and current density magnitude J (b) taken from x = 15 R E , −15 R E < y < 15 R E in the z = 0 equatorial plane. J is normalized via J 0 = B w /L 0 ; B z is normalized via B w = 10 −4 G. Colors denote times relative to t = 0 in Fig. 5 (top left panel). Arrows highlight where local pileup of B z on dawnside thickens the current sheet, resulting in a lower current density and impeding local reconnection. 1000 C. M. Bard and J. C. Dorelli: Tail asymmetries biases found in MESSENGER observations (Sun et al., 2016;Dewey et al., 2018) and global simulations of Mercury (Dong et al., 2019;Chen et al., 2019). It is interesting to note that our results are under a steady, southward solar wind condition. As long as there is Hall-driven convection in the tail, the competition between dawnside B z pileup and reconnection drives this cycle. At the moment, it is not clear whether this process is unique to our ion-scale Earth, since its strong planetary dipole field means that flux piles up over a large swath of the tail. It is possible that a similar process may occur at Mercury (which has a weaker dipole field), i.e., that its observed dipolarizations are indeed akin to global substorms (Kepko et al., 2015). Further investigation is needed to determine how varying the magnetospheric parameters (as presented in Sect. 3) affects these observations, especially as system size increases relative to δ i .
We note that, at Earth, there are additional localized (not global) dipolarization fronts resulting from current sheet instabilities or transient reconnection events (e.g., Runov et al., 2009;). We do not see these small-scale fronts in our ion-scale Earth; this may be because we do not have enough down-tail resolution to observe localized current sheet instabilities which form them.

Current sheet thickness
Another test of the active region picture is the predicted thickness asymmetry of the tail current sheet: Liu+ predicted that the sheet would be thinner on the dawnside. We follow Poh et al. (2017a) and estimate the current sheet thickness in our model by using a Harris sheet (Harris, 1962): where B a is the asymptotic lobe field, Z 0 is the current sheet center, L CS is the current sheet half-thickness, and the offset allows for asymmetry between the north and south B x lobes on either side of the current sheet. We take 6000 one-dimensional cuts of B x along the north-south direction between z = ±10 R E in a volume covering the current sheet from 12 R E < x < 16 R E and −15 R E < y < 15 R E , randomly sampled across the box plane and during the final 45 t 0 ≈ 48 min period (example shown in Fig. 7). These cuts are fit to Eq. (11) using the Levenberg-Marquardt leastsquares algorithm in scipy.curvefit (Virtanen et al., 2020); instances that do not fit well (χ 2 > 0.01) or that return nonsensical results (L CS < 0) are rejected. This results in 5037 samples of the current sheet thickness across the magnetotail (Fig. 8). This distribution shows that the dawnward current sheet is thinner on average than the duskward sheet. However, there is a significant scatter in this result; the dawn sheet covers a wider range of thicknesses. This variation is caused by the dawnside pileup-reconnection mechanism.
The current sheet oscillates with the dipolarization cycle (Sect. 4.2) between a "thick state" due to the B z pileup and a "thin state" immediately following the flux unloading and plasmoid ejection. This is demonstrated in Fig. 9, where fitted CS thicknesses during both flux loading and unloading stages are plotted along with snapshots of the B z state. During the loading stage, the piled up flux on the dawnside (5 R E < y < 12 R E ) fattens the current sheet; here, the sampled dawn thicknesses are comparable to and can exceed the dusk thicknesses. However, after the unloading stage, the current sheet on the dawnside is much thinner where the flux has been evacuated (bottom right plot; R > 15 R E ). Interestingly, we can see that where the B z flux remains (R < 15 R E ), the current sheet continues to be thick. Combining all the sample fits over several cycles of loading and unloading results in the picture shown in Fig. 8: a dawnward current sheet moving between thick and thin states depending on the level of flux pileup. Indeed, this is a common pattern throughout the simulation: where there is flux pileup, the current sheet is thicker, and the current density is lower (e.g., Fig. 10).
This cycle may explain the apparent contradiction between the Liu+ prediction of thinner dawnward current sheets in ion-scale magnetospheres and the Poh et al. (2017b) spacecraft observation of thicker dawnward sheets at Mercury. Even though, on average, the current sheet is thinner dawnward (as Liu+ predicts), the sampling of measurements could be producing the opposite result. As shown in Figs. 8 and 9, the sampled sheet thickness can greatly depend on where and when the craft crosses the tail. In our simulation, the current sheet is continuously morphing between thick and thin states; both types of regions exist simultaneously within the dawnside. Most points in the tail preferentially see thicker sheets over time, though some preferentially see thinner sheets. It is possible that these effects combine to produce a sampling bias in time and space towards thicker sheets. We note that this is speculation and will require more investigation with respect to the various solar wind driving conditions and seasons that MESSENGER experiences at Mercury.

Conclusions
We have simulated a small, ion-scale Earth in which the standoff distance and magnetotail width are akin to Earth's as measured in planetary radii but with the solar wind ionscale length δ i set to 1 R E . The resulting tail plasma behavior was more similar to Mercury's magnetosphere than Earth's. Along with Chen et al. (2019) and Dong et al. (2019), our results support the idea that tail asymmetry is a universal consequence of the Hall effect in ion-scale magnetospheres. Essentially, it is the relative size of the magnetosphere compared to δ i , not the absolute size (planetary radii), that controls the importance and influence of Hall-induced consequences. We find that Hall effects are sufficient to generate tail asymmetries in dipolarization, plasmoids, and current sheet thickness. No electron-scale kinetic effects are re-quired, though they may contribute to or modify asymmetries. However, we emphasize that we did not simulate the same magnetosphere as Chen et al. (2019): our ion-scale Earth is smaller relative to δ i than Mercury and has different magnetospheric parameters (Sect. 3). There may be additional effects not being considered, especially with regards to how varying the dimensionless magnetospheric parameters affects the manifestation of tail asymmetries.
In general, our simulation appears to corroborate the Liu+ picture of tail asymmetry in ion-scale magnetospheres; however, the Lu+ finding that the transported tail B z thickens the current sheet is readily manifested here. Although the reconnected B z does drive outflows and thin current sheets on the dawnside, we see that it can pile up and thicken current sheets. There is a continuous cycle between the dawnward transport of B z leading to pileup (which thickens the current sheet) and reconnection (which thins the current sheet); this manifests in an oscillating current sheet thickness. On average, we find the current sheet is thinner on the dawnside, but it can occasionally be thicker in some regions depending on the level of flux pileup.
Further study will be required to confirm or contrast this picture for magnetospheres with system size δ i . Since our simulation is of a experimental magnetosphere, several questions concerning more realistic magnetotails remain: -How does the weaker, offset dipole of Mercury affect the amount of magnetic flux available for transport/pileup and the resulting plasmoid formation/ejection?
-Are the observed dipolarizations at Mercury actually "global", like substorms?
-How does increasing the system size / δ i ratio affect asymmetry formation, tail convection, transport of B z , and plasmoid and DF formation?
-What other effects (e.g., kinetic, ionosphere) cause asymmetries, and how do they interact with the Hall effect and one another?
We look forward to future studies which will investigate these questions in greater detail.
Code and data availability. Observational data were not used, nor created for this research. The model algorithm is described above and in the references, and simulation parameters are given for reproducing the magnetosphere.
Author contributions. CMB developed the code, analyzed the simulation results, wrote the manuscript and produced the figures. JCD edited the manuscript and provided computing resources. Both authors conceptualized the research goals and guided the direction of inquiry.
Competing interests. The contact author has declared that neither they nor their co-author has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Acknowledgements. Christopher M. Bard thanks Alex Glocer for useful discussions concerning simulation techniques. Christopher M. Bard thanks Ryan Dewey for helpful discussions concerning Mercury tail asymmetries and MESSENGER observations. Plots were created using matplotlib (Hunter, 2007). The authors thank Gabor Toth and the anonymous referee for very helpful comments which improved this paper. Observational data were not used, nor created for this research; the model algorithm and simulation parameters are described above.
Financial support. Christopher M. Bard's software development has been partly supported by an appointment to the NASA Postdoctoral Program at NASA Goddard Space Flight Center, administered by the Universities Space Research Association under contract with NASA. This research has been supported by the Goddard Space Flight Center (Internal Scientist Funding Model; grant no. HISFM18-0009).
Review statement. This paper was edited by Minna Palmroth and reviewed by Gabor Toth and one anonymous referee.