We use measurements recorded along traversals through Earth’s magnetotail by the MMS spacecraft over the years 2017–2022. These traversals are identified based on the information provided by the new mission phases and science region of interest of MMS (https://lasp.colorado.edu/galaxy/display/MFDPG/1.3+Mission+Phases+and+Science+Regions+of+Interest) and by visual inspection when no such information is available. In what follows field and particle measurements from MMS1 alone are shown noting that the results returned from the other spacecraft are for our purposes identical. Survey mode magnetic field data from the Flux-gate Magnetometer (FGM) (Russell et al., 2014) and fast-mode electric field data (Ergun et al., 2014; Lindqvist et al., 2014) from the Electric Double Probe (EDP) instrument (Torbert et al., 2014) are employed to evaluate the Poynting flux in GSM coordinates. Measurements from the Fast Plasma Investigation (FPI) instrument (Pollock et al., 2016) in the form of velocity moments evaluated over the energy range above the spacecraft potential and below the 30 keV upper limit of the instrument are employed to derive the kinetic terms. These data are returned every of 4.5 s. This cadence defines the “base” time resolution of the measurements we report. The evaluation of these moments incorporates corrections due to spacecraft charging to eliminate spurious contributions from secondaries and photo-electrons.

From these measurements a database is created by averaging over consecutive intervals defined by advection in the flow corresponding to a distance of 2 π ρ i , with ρ i = m v t h q B being the ion gyro radius and m , q and v t h are respectively the mass, charge, and thermal speed of the protons. This approach eliminates sub-gyro radius variations in the contributions to the moments and in the Poynting fluxes derived, while still quantifying the transport into the outer kinetic range. Notably, variations on electron scales are eliminated by this approach. Those intervals where the peak in the differential ion energy flux spectrum exceeds the 30 keV energy per charge, maximum of the FPI instrument are flagged and estimates recorded at these times are not included in the results shown here. It is worth mentioning that fast flows and hence large bulk kinetic energy fluxes, large enthalpy fluxes, and large heat fluxes that exceed the instrumental limit are not included in the present analysis, so the average values may be somewhat underestimated. While there is no criterion on flow speed, measurements nominally restricted to the plasma sheet by setting the ion number density to N i ≥ 1 c m − 3 , following the approach of Hasegawa (2012). With these criteria a total of 582,414 data points are collected throughout the PS from MMS1, with each of these points composed of an average of measurements over 2 π times an advected thermal ion gyro-radius. Measurements from MMS2, MMS3, and MMS4 produce equivalent results.

The energy flux is decomposed into the electromagnetic energy, or Poynting flux, kinetic fluxes including the enthalpy flux and heat flux, and the bulk kinetic energy flux. The temporal evolution of total energy density W is (e.g., Birn and Hesse, 2005; Liu et al., 2021) d W d t = − ∇ . ρ v 2 2 v + v T r P ⃡ 2 + v ⋅ P ⃡ + Q H + 1 μ 0 E × B (1) Here, we evaluate each term on the RHS, where the first term represents the bulk kinetic energy flux, where ρ = n i m i denotes the ion mass density and v is the ion flow speed. The second term within the parenthesis provides the expression for the enthalpy flux computed from the full pressure tensor, P ⃡ and includes the contribution of the internal energy. The Dual Ion Spectrometer (DIS), part of the FPI suite aboard MMS, provides 3D angular coverage and measures the 3D distribution function of ions in space (Pollock et al., 2016). The zeroth-order moment of which provides the number density, the first-order moment yields the bulk velocity, and the second-order moment determines the pressure tensor, which have been used for flux calculations in this study. Additionally, the third-order moment calculates the heat flux, which contributes to the third term on RHS. To compute the Poynting flux, as denoted by the last term on the RHS, we utilized survey mode magnetic field B data from the FGM and the fast survey mode electric field E data from the EDP instruments. We have used averages over 2 π advected gyro-radii in E and B and report the total Poynting flux as we seek to evaluate the total electromagnetic energy transport. This is consistent with studies by Angelopoulos et al. (2002), Chaston et al. (2012). The FGM and EDP data undergo filtering and down-sampling to coincide with the FPI measurements. Notably, in the plasma sheet Equation 1, the electron bulk kinetic energy flux can be ignored due to the significantly larger ion mass, electron enthalpy flux can also be neglected owing to the ion-to-electron temperature ratio of ∼ 5 or greater in the plasma sheet (Kaufmann et al., 2005). Additionally, all of the contributions to the energy flux are evaluated over an advected spatial scale of 2 π ρ i , so that peak variations in the electron heat flux do not feature prominently in the statistics and on this scale and larger are generally equal to or less than the ion heat flux, while representing a minority fraction of the total transport. A comparison of the ion to electron contributions on this scale are provided in the next section to demonstrate this point.

A representative interval measured from MMS1 on 1 August 2018 is presented in Figure 1 showing time-series fields and plasma parameters and the corresponding derived energy flux components on which the database is based. As shown, the spacecraft journeys from − 15 R E to the more distant magnetotail around − 24 R E . The three components of the ion bulk flow and their absolute value in geocentric solar magnetospheric (GSM) coordinates are presented in Figures 1B, C, respectively. Predominantly, the ion flow is directed earthward, as evidenced by the positive x-component of the ion bulk flow. Four distinct regions of enhanced transport, delineated by colored dashed lines, can be identified as follows: region-1 (magenta; 00:45:00–01:40:00 UT), region-2 (blue; 14:07:00–14:54:00), region-3 (red; 17:53:00–19:25:00), and region-4 (black; 22:47:00–23:59:00). Each comprise intervals when the absolute value of the ion bulk flow exceeds 50 km/s for an extended period and correspond to a rise in differential ion energy flux as illustrated in Figure 1A. In region-1, MMS1 transitions into the plasma sheet, as indicated by a sudden rise in ion temperature (Figure 1E) and an increase in ion density (Figure 1D). Additionally, ion plasma beta (Figure 1F) approached approximately ∼ 1.0 . During this interval, a slow dipolarization in the background magnetic field is observed, while the wave electric field exhibits intensified oscillations as depicted in Figures 1G, H, respectively. The resultant Poynting flux shown in Figure 1I is mostly Earthward (+x-direction), albeit with significant dawn-ward and dusk-ward components (+/− y-direction). The ion enthalpy flux, bulk kinetic energy flux, and heat flux are also predominantly directed earthward as depicted in Panels (j) (k), and (l), respectively with the enthalpy flux dominant. A similar pattern is captured in Regions 2, 3, and 4 for the enthalpy, bulk kinetic and heat flux, while the Poynting flux is predominantly in the y-direction, which is an often repeated pattern of enhanced electromagnetic and kinetic transport in association with fast flows. While this is a well-known feature of magnetotail energy transport (Angelopoulos et al., 1994), there are appreciable variations in the relative contributions of each component that are dependent on location, the field geometry and plasma beta. Indeed, persistent morphological structuring emerges in the examination of the statistical distribution of these energy fluxes as we now describe.

Figures 1M–O presents a comparison of absolute values of enthalpy flux, bulk kinetic energy flux and heat flux for electrons (blue curve) and ions (red curve), respectively. As discussed in Section 2.2, the contributions of ion enthalpy and bulk kinetic energy fluxes are generally higher than those of the electrons, and the electron and ion heat flux contributions are nearly equal over this interval.

Using the procedures outlined in Section 2.1, statistics describing the decomposed energy transport contributions in GSM coordinates were compiled over the spatial range covered by the MMS spacecraft within Earth’s magnetotail. For the purpose of presenting summary distributions in this report, each GSM plane in the magnetotail (i.e., XY, XZ, and YZ planes) is partitioned into bins with a window size of 1 R E × 1 R E . Each bin considered in compiling the statistics contains at least 10 gyro-averaged measurements of the absolute value of each transport quantity. The distribution of events per bin is shown in the Supplementary Figures 1–3. Inspection of the distributions within each bin shows a form symmetric about the mean values we report here. In Figures 2, 3, the color bar indicates the corresponding in situ absolute value of the energy flux (as specified in the figures), measured in m W / m 2 , at a location corresponding to the (abscissa, ordinate) pair depicted in figures. A separation into Earthward and tailward transport is considered in the context of pressure balance and field geometry in Section 3.2 while we reserve a description of the full vector quantities for subsequent work. In this study, the magnetotail region characterized by X GSM < − 15 R E is designated as the distant magnetotail, while the area where − 15 R E ≤ X GSM ≤ − 5 R E is identified as the near-Earth magnetotail region. Moreover, the region with − 12 R E ≤ Y GSM ≤ + 12 R E is demarcated as the central region, and that with | Y GSM | > + 12 R E is considered the edge of the distribution plane.

Figure 2A illustrates the distribution of Poynting flux in the equatorial plane (XY-plane). The Poynting flux exhibits a gradual increase from the distant tail region to the near-Earth region, attaining its maximum value in the near-Earth region. Further, Poynting flux decreases from the central region towards the edges of the XY-plane. The corresponding distribution of enthalpy flux shown in Figure 2B is more irregular or patchy with comparatively larger values than the Poynting flux in the distant magnetotail. Moreover, in contrast to the Poynting flux distribution, the enthalpy flux is somewhat enhanced on the flanks of the near-Earth region, and most prominently in the dusk-ward sector. Additionally, unlike the Poynting flux, where the flux is maximum close to the origin in Y GSM -axis, there are multiple enthalpy flux maxima distributed throughout the X-Y plane. The distribution of bulk kinetic energy flux shown in Figure 2C has relatively larger value in the distant magnetotail central plane and decreases towards the near-Earth region. Its minimum occurs around X ∼ ( − 10 to − 5 R E ) . The significant depression at X ∼ − 15 R E in the dawn-ward sector provides an asymmetric distribution with larger values observed on the dusk-side. Further, the flux value in the central region of the XY-plane is comparatively larger than the edge of the XY-plane in the distant magnetotail region, while it is larger at the edge of the XY-plane than the central portion in the near-Earth region. The heat flux distribution replicates that of the enthalpy flux (Figure 2D).

The distribution of the energy fluxes in the GSM XZ-plane, is depicted in Figures 3A–D. The Poynting flux increases gradually from the distant magnetotail to the near-Earth region, as previously identified (cf. Figure 2A). There is a general depletion at low values of Z GSM with average minima residing at generally small positive Z GSM values. Conversely (Figure 3B), reveals that, unlike the Poynting flux, the peak of the distribution in enthalpy flux is located near the equatorial plane in the GSM system with a progressive shift toward more positive Z GSM values with distance downtail. Significantly, the maximum of the average distribution occupies the region of depletion in the Poynting flux. The bulk kinetic energy and heat flux distributions shown in Figures 2C, D share a similar morphology, albeit the bulk kinetic energy flux is greatly reduced within X GSM > − 10 R E .

The distributions in the YZ-plane are presented in Figures 3E–H. In Figure 3E, the distribution shows that most of the Poynting flux is concentrated within the range − 12 R E ≤ Y GSM ≤ + 12 R E , while occupying the entire range of Z GSM . Furthermore, within this region, the fluxes are accentuated in the south reflecting of the distribution at the inner edge of the plasma sheet shown in Figure 3A. In contrast, the distribution of the enthalpy flux in the YZ-plane (Figure 3F) shows that for a b s ( Y GSM ) > 10 R E , the enthalpy flux is enhanced in the south, while at smaller Y GSM values, the flux is predominantly enhanced in the north. The latter reflects the distribution in the XZ-plane shown in Figure 3B. A similar morphology is found in the distribution of heat flux (Figure 3H) and bulk kinetic energy flux (Figure 3G), albeit with reduced flux.

To provide physical insight into the distribution of these energy fluxes, we resample the statistics into a 2-D coordinate system defined by pressure balance (Baumjohann et al., 1990; Xing et al., 2010) and the magnitude and orientation of the magnetic field in the X-Z plane. Following the formulation outlined by Watanabe et al. (2019) the magnetic field in the lobes is defined as B lobe 2 2 μ 0 = n i k B T i + T e + B 2 2 μ 0 , (2) where, B 2 = B x 2 + B y 2 + B z 2 . The observations are organized according to the normalized fields, s i g n ( B x ) B x 2 + B y 2 / B lobe and B z / B lobe , where s i g n B x = B x / | B x | . The ratio s i g n ( B x ) B x 2 + B y 2 / B lobe varies between − 1 and + 1 , indicating the southern and northern lobes respectively with the neutral sheet located at s i g n ( B x ) B x 2 + B y 2 / B lobe = 0 . Positive and negative values of B z / B lobe signify regions Earthward and tail-ward of a nominal X-line as depicted in Figure 4B. Different regions in the physical XZ-plane and their corresponding counterparts in the transformed system are denoted by A, B, C, D, and E, which are useful in identifying Earthward (B and D), tail-ward (A and C) region, and the reconnection (E) line in the plane. Within this system the statistics are distributed into bins of size 0.1 s i g n ( B x ) B x 2 + B y 2 / B lobe × 0.1 B z / B lobe , and the average values of various energy fluxes and plasma parameters evaluated within each. The x-component of the ion flow and energy fluxes in this system are illustrated in Figure 4H, and Figures 4C–F, respectively. Consistent with this formulation the direction of the flow and energy flux contributions reverse at B z / B lobe = 0 , while s i g n ( B x ) B x 2 + B y 2 / B lobe → 0 along this reversal identifies the x-line where very high plasma beta necessarily occurs (Figure 4G). The distribution of X GSM values over this plane, depicted in Figure 4I, indicates that the statistical location of the x-line is approximately at X = − 20 R E . It is also apparent from Figures 4I, J that with increasing distance down-tail the neutral sheet resides at increasingly positive Z-GSM values thereby accounting for the shift in the flux distributions toward positive Z GSM identified in the X-Z plane in Figure 3.

Several salient features of the energy transport distributions emerge in this coordinate system (Equation 2). Firstly, near the x-line, the fluxes in the tail-ward direction (i.e., negative B z ) without exception are larger than those in the Earthward (positive B z ) direction. This disparity may stem from the formation of a near-Earth neutral line during active times, resulting in more frequent occurrences of tailward flows when fluxes are larger. Secondly, the x-component of Poynting flux, illustrated in Figure 4C, is observed to be smallest near the reconnection point and gradually increases in the Earthward direction in two channels displaced from, and roughly symmetric about the neutral sheet with localized peaks at the inner edge plasma sheet. Based on the X GSM distribution, the location of maximum Poynting flux is estimated to be around X GSM ≥ − 10 R E , falling within the range of the flow braking region (Shiokawa et al., 1997). While our interest here is primarily in the dominant X-directed fluxes, we present in Supplementary Figure 4 accompanying this article the z-component of the Poynting flux. As found in previous studies (Miyashita et al., 2012) this component is convergent on the neutral sheet with values approaching the X-directed fluxes near the lobes. Thirdly, in contrast to the Poynting flux, the x-component of kinetic fluxes (i.e., enthalpy, heat and bulk kinetic energy fluxes) peak mid-tail in the neutral sheet albeit somewhat broadly dispersed in s i g n ( B x ) B x 2 + B y 2 / B lobe . Finally, there is a region of enhanced flow at negative s i g n ( B x ) B x 2 + B y 2 / B lobe corresponding to the southern inner plasma sheet that provides a secondary peak in the kinetic fluxes. The origin of this enhancement is not clear at the time of writing. These features in the electromagnetic and kinetic fluxes complete a picture whereby the kinetic contributions provide the bulk of the energy transport through the central plasma sheet with the Poynting flux providing a comparable fraction of the transport through the high-latitude flanks and dominating at the inner edge of the plasma sheet. The overall pattern is represented schematically in Figures 4A, B representing energy flows divergent from the statistically identified X-line and requiring energy conversion at that location as described by Equation (1).