Comparing Monte Carlo simulations, mean particle theory estimates, and observations of H⁺ and O⁺ outflows at high altitudes and latitudes

We conducted a comparative analysis of the results from Monte Carlo simulations, mean particle theory estimates, and available observational data across different regions of Earth's magnetosphere, including the auroral, polar-wind, central-polar-cap, and cusp regions, focusing on the outflow of H⁺ and O⁺ ions at high latitudes and altitudes. We present altitude profiles for the mean perpendicular energy W_⊥, mean parallel energy W_∥, and mean total energy W_total. The Monte Carlo simulations were carried out using the Barghouthi model (Barghouthi, 2008), while the mean particle theory estimates were derived from Chang et al. (1986), and the observational data were obtained from Nilsson et al. (2013) and Barghouthi et al. (2016). The results of the comparison across different regions reveal the following findings: (1) Monte Carlo simulations and mean particle theory yield similar results in the auroral regions but show discrepancies in the polar-wind region. This discrepancy is attributed to the dominance of wave–particle interactions, which overshadow the effects of external forces (such as gravity, the polarization electric field, the mirror force, and centrifugal acceleration) in the auroral region and compete with them in the polar-wind region. (2) The use of altitude-dependent diffusion coefficients leads to the generation of high-energy particles, which are not consistent with the corresponding observational data at middle and high altitudes. Therefore, we recommend the use of velocity- and altitude-dependent diffusion coefficients. (3) A comparison with observations in the auroral and polar-wind regions demonstrates excellent agreement in the auroral region and good agreement in the polar-wind region, attributed to the implementation of appropriate velocity- and altitude-dependent diffusion coefficients. (4) In the central-polar-cap and cusp regions, both methods and observations show excellent agreement. (5) Based on these comparisons, we conclude that the electromagnetic wave wavelength in the polar-wind and auroral regions is approximately 8 km, and the velocity- and altitude-dependent diffusion coefficients used in the Monte Carlo simulations and mean particle theory are suitable for application in further studies of these regions.

1 Introduction

Numerous studies, including analytical, modeling, and observational approaches, have been conducted to investigate the behavior of ion outflow (e.g., O⁺ and H⁺) from Earth's polar regions into outer space. These studies aim to characterize various properties of the ion outflow, such as velocity distribution, temperature, density, drift velocity, and heat flux. Ongoing research efforts continue to enhance our understanding of the dynamics of these ion flows, with a particular focus on the energy and its components (parallel, perpendicular, and total energy), which constitutes the primary subject of this study.

Various researchers have investigated ion outflow at high altitudes and latitudes. Chang et al. (1986) introduced the mean particle theory, which explains the perpendicular heating of ions in a dipole magnetic field. They suggested that intense broadband electric-field fluctuations observed in the frequency range of 0–100 Hz could induce transverse ion activation through cyclonic resonance heating by left-handed polarized electromagnetic waves. Additionally, by employing a set of equations governing ion motion within the geomagnetic field, they derived expressions for the parallel and perpendicular energies based on the mean particle theory. Retterer et al. (1987) demonstrated how oxygen ions form conic distributions in the auroral zone and utilized the diffusion equation to explain ion velocity distributions obtained through the Monte Carlo method.

Barakat and Barghouthi (1994a, b) enhanced the Monte Carlo simulation and explored the effects of wave–particle interactions (WPIs) on H⁺ and O⁺ ion outflow in the polar wind. Their model incorporated the electrostatic-field, gravity, and geomagnetic-field lines. These studies are considered to be parametric investigations as they used constant values for the quasi-linear velocity diffusion rates along the simulation tube. The velocity distribution function and its corresponding velocity moment profiles were simulated and presented for both ion species. Barghouthi (1997) and Barghouthi et al. (1998) analyzed data from the Plasma Wave Instrument (PWI) aboard the Space Dynamics Explorer 1, calculating the altitude dependence of the velocity diffusion rate. They examined the impact of altitude-dependent wave–particle interactions (WPIs) on H⁺ and O⁺ ion outflow in both the polar-cap and auroral zones using Monte Carlo simulations. Additionally, Barghouthi (1997) compared energy estimates from mean particle theory (Chang et al., 1986) with the corresponding energy results produced by Monte Carlo simulations in the auroral region. Despite the lack of supporting observational data, strong agreement was found between the Monte Carlo simulations and the mean particle theory estimates.

Bouhram et al. (2002, 2003a, b, 2004) developed a two-dimensional Monte Carlo model to investigate ion outflow from the dayside cusp or cleft, focusing on the effects of transverse ion heating. They examined the mechanisms of transverse heating and ion outflow within the cusp or cleft region. Using their model, they interpreted Cluster satellite observations, specifically the saturation of the transverse energization rate, in terms of the influence of finite perpendicular wavelength effects in wave–particle interactions.

Barghouthi and Atout (2006) focused on Monte Carlo simulations of toroidal H⁺ and O⁺ velocity distributions at high altitudes, equatorward of the cusp, by employing an appropriate form for the velocity diffusion coefficient D_⊥. The results of the Monte Carlo simulations, including the toroidal H⁺ and O⁺ velocity distributions, as well as the H⁺ and O⁺ ion temperatures, were compared to the corresponding toroidal H⁺ and O⁺ ion distributions and ion temperatures observed at high altitudes, equatorward of the cusp (Huddleston et al., 2000). These comparisons yielded reasonable agreement.

Barghouthi (2008) utilized Monte Carlo simulations to determine the temperatures and velocity distributions of H⁺ and O⁺ ions at high altitudes in the equatorward portion of the cusp, using various forms of altitude- and velocity-dependent diffusion coefficients $D_{\perp }(r,v_{\perp })$, including the RCC model, Bouhram model, and Barghouthi model. The simulation results were compared with the corresponding observations from Huddleston et al. (2000), and the results from the Barghouthi model showed excellent agreement with the observations. Furthermore, Barghouthi (2008) provided substantial evidence, including comparisons between Monte Carlo simulation results for H⁺ and O⁺ ion outflows obtained using the Barghouthi model and corresponding observational data from multiple sources at various altitudes in the auroral region, further supporting the validity of the Barghouthi model.

Waara et al. (2010) presented a case study of significant heating of outflowing oxygen ions at high altitudes (12R_E) above the polar cap, with energies reaching up to 8 keV perpendicular to the geomagnetic field. The shape of the distribution functions suggests that the majority of the heating occurs locally, within 0.2−0.4R_E altitude. They found that the locally observed wave fields are unlikely to account for the observed ion energization. Additionally, it is improbable that the ions originated from a nearby energizing location and migrated to the observation site. These findings indicate the existence of additional, fundamentally distinct ion energization mechanisms at high altitudes. One possible explanation is that the ions' magnetic moment is not conserved, which would result in slower outflow velocities and an extended period of ion energization.

Waara et al. (2011, 2012) conducted a statistical analysis of ion heating and associated wave activity. They provided average values for coefficients that describe diffusion in ion velocity space at various altitudes, offering a useful framework for studying ion outflow behavior and energy characteristics. Their test particle calculations suggest that the average energies of O⁺ ions correlate with the observed wave activity at high altitudes (8−15R_E) in the cusp and mantle regions. They also found that the electric-to-magnetic-field spectral-density ratios closely match the expected values for Alfvén waves. According to their findings, the diffusion coefficient for O⁺ ions increases with altitude.

Barghouthi et al. (2012) compared simulation results of ion outflow (including ion density, drift velocity, perpendicular and parallel temperatures, and ion velocity distributions at various altitudes) in two distinct regions, the polar-wind and auroral regions, based on the Barghouthi model. They found that wave–particle interactions had a more significant effect in the auroral zone compared to in the polar-wind region and that these interactions had a greater influence on the energization of O⁺ ions than on that of H⁺ ions.

Barghouthi et al. (2016) updated the Monte Carlo model to include the effects of gravity, the ambipolar electric field, centrifugal acceleration, mirror force, and wave–particle interactions in the study of O⁺ and H⁺ ion outflow above the polar cap. They modified several parameters, such as centrifugal acceleration, velocity diffusion coefficients, and boundary conditions at lower altitudes. The results were compared with observational data obtained from various instruments aboard the Cluster spacecraft, and the simulation results were found to be in good agreement with the observed data.

The primary objective of this study is to compare the simulation results (perpendicular energies W_⊥, parallel energies W_∥, and total energies W_total) of O⁺ and H⁺ ions obtained using the Monte Carlo model (Barghouthi model) and mean particle theory with available observational data from various regions of Earth's magnetosphere, including the polar wind, auroral region, cusp, and central polar cap. This comparison seeks to evaluate the importance of including altitude-dependent velocity diffusion coefficients D_⊥(r) or altitude- and velocity-dependent velocity diffusion coefficients $D_{\perp }(r,v_{\perp })$ while also considering the constraints posed by wavelength limitations.

2 Formulations

2.1 Monte Carlo simulation

In the study of space plasma, it is advantageous to describe each constituent species utilizing a distinct velocity distribution function, denoted as $f_s({\mathit{v}}_s,{\mathit{r}}_s,t$). The velocity distribution function is defined such that the expression $f_s({\mathit{v}}_s,{\mathit{r}}_s,t)\mathrm{d}{\mathit{v}}_s\mathrm{d}{\mathit{r}}_s$ quantifies the number of particles of species s that, at time t, possess velocities within the interval v_s and v_s+dv_s and positions between r_s and r_s+dr_s. The temporal evolution of the species' velocity distribution function is governed by the cumulative effects of collisions and interactions, as well as the dynamical movements of the species in phase space under the influence of external forces (Schunk, 1977). This evolution can be mathematically represented by the well-established Boltzmann equation:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial f_s}{\partial t}}+{\mathit{v}}_s\cdot \mathrm{\nabla }{f}_s+\left({\displaystyle \frac{e_s}{m_s}}\right)\left[\mathbf{E}+{\displaystyle \frac{{\mathit{v}}_s\times \mathbf{B}}{c}}\right]\cdot {\mathrm{\nabla }}_{v_s}{f}_s={\displaystyle \frac{\mathit{\delta }{f}_s}{\mathit{\delta }t}}.\end{array}$$

In this equation, the left-hand side represents the evolution of the velocity distribution function $f_s({\mathit{v}}_s,{\mathit{r}}_s,t)$ under the effects of external forces, and the right-hand side represents the Boltzmann collision integral; here, it represents the rate at which f_s changed as a result of wave particle interactions in the region of study. In the above equation, E is the polarization electric field; B is the geomagnetic field; c is the speed of light; ∇ is the coordinate space gradient; ${\mathrm{\nabla }}_{v_s}$ is the velocity space gradient; and e_s and m_s are the charge and mass of species s, respectively. The suitable expression for $(\mathit{\delta }{f}_s/\mathit{\delta }t)$ in the case of wave particle interactions is given by Retterer et al. (1987); they considered the effects of WPI as particle diffusion in the velocity space.

$$\begin{array}{}\text{(2)}& {\left.{\displaystyle \frac{\mathit{\delta }f}{\mathit{\delta }t}}\right|}_{\text{WPI}}=\left({\displaystyle \frac{\mathrm{1}}{v_{\perp }}}\right){\displaystyle \frac{\partial }{\partial v_{\perp }}}\left[D_{\perp }{v}_{\perp }{\displaystyle \frac{\partial f}{\partial v_{\perp }}}\right]\end{array}$$

In the above, D_⊥ is provided by Retterer et al. (1987) and represents the quasi-linear velocity diffusion rate perpendicular to the geomagnetic field:

$$\begin{array}{}\text{(3)}& D_{\perp }=\left({\displaystyle \frac{\mathit{\eta }{q}^{\mathrm{2}}}{\mathrm{4}{m}^{\mathrm{2}}}}\right)\left|E_x\right(\mathit{\omega }=\mathrm{\Omega }){|}^{\mathrm{2}},\end{array}$$

where $\left|E_x\right(\mathit{\omega }){|}^{\mathrm{2}}$ is the measured spectral density of the electromagnetic turbulence, η is the proportion of the measured spectral density of the plasma wave instrument (PWI) on board the Dynamic Explorer 1 (DE-1) spacecraft that corresponds to the left-hand polarized wave, q is the ion's charge, m is the ion's mass, Ω is the ion's gyrofrequency, and ω is the angular frequency of the electromagnetic turbulence.

The expression of the velocity diffusion rate D_⊥ as given in Eq. (3) is independent of velocity and depends on position (altitude) via changes in the ion gyrofrequency Ω along the geomagnetic-field lines. By examining experimental data of the electric-field spectral density obtained by the plasma wave instrument (PWI) on board the DE-1 satellite (i.e., for high-solar-activity conditions), Barghouthi (1997) and Barghouthi et al. (1998) calculated the altitude dependence of D_⊥. They came up with the following expressions for the velocity diffusion coefficient D_⊥ in the polar-wind region (Barghouthi et al., 1998):

$$\begin{array}{}\text{(4)}& D_{\perp }\left(r\right)=\left\{\begin{array}{ll}\mathrm{5.77}\times {\mathrm{10}}^{\mathrm{3}}{\left(\frac{r}{R_{\mathrm{E}}}\right)}^{\mathrm{7.95}}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{3}},& \text{for H}\text{+}\\ \mathrm{9.55}\times {\mathrm{10}}^{\mathrm{2}}{\left(\frac{r}{R_{\mathrm{E}}}\right)}^{\mathrm{13.3}}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{3}},& \text{for O}\text{+}\end{array}\right..\end{array}$$

In the auroral region, D_⊥(r) is given by Barghouthi (1997) as follows:

$$\begin{array}{}\text{(5)}& D_{\perp }\left(r\right)=\left\{\begin{array}{ll}\mathrm{4.45}\times {\mathrm{10}}^{\mathrm{7}}{\left(\frac{r}{R_{\mathrm{E}}}\right)}^{\mathrm{7.95}}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{3}},& \text{for H}\text{+}\\ \mathrm{6.94}\times {\mathrm{10}}^{\mathrm{5}}{\left(\frac{r}{R_{\mathrm{E}}}\right)}^{\mathrm{13.3}}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{3}},& \text{for O}\text{+}\end{array}\right..\end{array}$$

In the central-polar-cap (CPC) and cusp regions, D_⊥(r) is given by Nilsson et al. (2013) as follows:

for the central-polar-cap region,

$$\begin{array}{}\text{(6)}& D_{\perp }\left(r\right)=\left\{\begin{array}{ll}\mathrm{20}{\left(\frac{r}{R_{\mathrm{E}}}\right)}^{\mathrm{9.77}}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{3}},& \text{for H}\text{+}\\ \mathrm{0.5}\times {\mathrm{10}}^{\mathrm{5}}{\left(\frac{r}{R_{\mathrm{E}}}\right)}^{\mathrm{5.5}}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{3}},& \text{for O}\text{+}\end{array}\right.,\end{array}$$

and for the cusp region,

$$\begin{array}{}\text{(7)}& D_{\perp }\left(r\right)=\left\{\begin{array}{ll}\mathrm{1.01}\times {\mathrm{10}}^{\mathrm{6}}{\left(\frac{r}{R_{\mathrm{E}}}\right)}^{\mathrm{5.61}}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{3}},& \text{for H}\text{+}\\ \mathrm{2.5}\times {\mathrm{10}}^{\mathrm{4}}{\left(\frac{r}{R_{\mathrm{E}}}\right)}^{\mathrm{6.4}}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{3}},& \text{for O}\text{+}\end{array}\right..\end{array}$$

The diffusion coefficient was given a new form by Barghouthi (2008), who discovered that it is velocity-dependent in addition to being altitude-dependent.

$$\begin{array}{}\text{(8)}& D_{\perp }(r,v_{\perp })=D_{\perp }\left(r\right)\left\{\begin{array}{ll}\mathrm{1}& \text{for}\left(\frac{k_{\perp }{v}_{\perp }}{{\mathrm{\Omega }}_i}\right)<\mathrm{1}\\ {\left(\frac{k_{\perp }{v}_{\perp }}{{\mathrm{\Omega }}_i}\right)}^{-\mathrm{3}}& \text{for}\left(\frac{k_{\perp }{v}_{\perp }}{{\mathrm{\Omega }}_i}\right)\ge \mathrm{1}\end{array}\right.\end{array}$$

In the above, $D_{\perp }(r,v_{\perp })$ is the quasi-linear velocity diffusion rate perpendicular to the geomagnetic-field lines (altitude- and velocity-dependent), Ω_i is the ion gyrofrequency, and k_⊥ is perpendicular wave number and is related to the characteristic perpendicular wavelength of the electromagnetic turbulence λ_⊥. Equation (8) indicates the diffusion coefficient that is dependent on altitude and velocity. However, as ions are heated and move to higher altitudes, their gyroradius may approach the perpendicular wavelength of the electromagnetic turbulence, and when the ratio ($k_{\perp }{v}_{\perp }/{\mathrm{\Omega }}_i$) exceeds 1, the heating rate becomes self-limited. Bouhram et al. (2004) derived an alternative form for the altitude- and velocity-dependent diffusion coefficient and interpreted their results in terms of finite wavelength effects.

By solving the Boltzmann equation – Eq. (1) – using the Monte Carlo technique, the velocity distribution functions were obtained for each species (in this study, O⁺ and H⁺ ions) and their velocity moments, i.e., density n_s, drift velocity u_s, and parallel T_s∥ and perpendicular T_s⊥ temperatures. The moments considered here are defined as follows (Barghouthi, 1997):

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{n}_s=\int f_s\mathrm{d}{\mathit{v}}_s,\text{(10)}& {\displaystyle }& {\displaystyle }{u}_s={\displaystyle \frac{\mathrm{1}}{n_s}}\int v_{s\parallel }{f}_s\mathrm{d}{\mathit{v}}_s,\text{(11)}& {\displaystyle }& {\displaystyle }{T}_{s\parallel }={\displaystyle \frac{m_s}{n_{s}k}}\int (v_{s\parallel }-u_s{)}^{\mathrm{2}}{f}_s\mathrm{d}{\mathit{v}}_s,\text{(12)}& {\displaystyle }& {\displaystyle }{T}_{s\perp }={\displaystyle \frac{m_s}{\mathrm{2}{n}_{s}k}}\int (v_{s\perp }{)}^{\mathrm{2}}{f}_s\mathrm{d}{\mathit{v}}_s.\end{array}$$

These Monte Carlo results will be used to calculate the mean parallel energy, mean perpendicular energy, and total mean energy as given in the following expressions (Barghouthi, 1997):

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{W}_{s\parallel }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}m{u}_s^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}k{T}_{s\parallel },\text{(14)}& {\displaystyle }& {\displaystyle }{W}_{s\perp }=k{T}_{s\perp },\text{(15)}& {\displaystyle }& {\displaystyle }{W}_s=W_{s\parallel }+W_{s\perp },\end{array}$$

where u_s, T_s∥, and T_s⊥ are given by Eqs. (10–12), respectively, and W_s∥ and W_s⊥ are the mean parallel and perpendicular energies, respectively. W_s is the total mean energy, s denotes the type of ion (O⁺ or H⁺), and k is the Boltzmann constant.

2.2 Barghouthi model

The Barghouthi model was developed to investigate the behavior of ion outflows, specifically for H⁺ and O⁺ ions, at high altitudes and high latitudes. The simulation results generated by this model show excellent agreement with observational data from various regions, including the auroral region (Barghouthi, 2008) and the polar-wind region (Barghouthi et al., 2011). This model accounts for multiple influential factors, including gravity, the polarization electric field, the diverging geomagnetic field, and wave–particle interactions, all of which affect H⁺ and O⁺ ion outflows at elevated altitudes and latitudes. Notably, the model highlights the significant contribution of wave–particle interactions, which are responsible for ion heating. The impact of these wave–particle interactions is characterized by the velocity diffusion coefficient $D_{\perp }(r,v_{\perp })$, which has been formulated as a function of position ($r/R_{\mathrm{E}}$) along Earth's geomagnetic-field lines and the perpendicular velocity of the injected ions (v_⊥). Various functional forms of the velocity diffusion coefficient $D_{\perp }(r,v_{\perp })$ have been employed in Monte Carlo simulations to derive the density, drift velocity, parallel and perpendicular temperatures, heat fluxes, and velocity distribution functions for H⁺ and O⁺ ions at high altitudes and latitudes. In the Monte Carlo simulations of this study, the appropriate velocity diffusion coefficient corresponding to the specific study region, $D_{\perp }(r,v_{\perp })$, was utilized to ascertain the temperatures and velocity distributions of H⁺ and O⁺ ions at these elevated altitudes and latitudes.

2.3 Mean particle theory

Chang et al. (1986) presented a theoretical framework for estimating the mean perpendicular, mean parallel, and total mean energies as functions of geocentric distance. This framework incorporates the average heating rate for each ion into a set of equations that describe the motion of ions along geomagnetic-field lines, as outlined below:

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }{W}_{i\parallel }={\displaystyle \frac{\mathrm{9}{m}_i}{{\mathrm{2}}^{\mathrm{1}/\mathrm{3}}}}{\left[{\displaystyle \frac{r{D}_{\perp }\left(r\right)}{(\mathrm{3}\mathit{\alpha }+\mathrm{1})(\mathrm{6}\mathit{\alpha }+\mathrm{11})}}\right]}^{\mathrm{2}/\mathrm{3}},\text{(17)}& {\displaystyle }& {\displaystyle }{W}_{i\perp }={\displaystyle \frac{(\mathrm{6}\mathit{\alpha }+\mathrm{2})m_i}{{\mathrm{2}}^{\mathrm{1}/\mathrm{3}}}}{\left[{\displaystyle \frac{r{D}_{\perp }\left(r\right)}{(\mathrm{3}\mathit{\alpha }+\mathrm{1})(\mathrm{6}\mathit{\alpha }+\mathrm{11})}}\right]}^{\mathrm{2}/\mathrm{3}},\text{(18)}& {\displaystyle }& {\displaystyle }{W}_i=W_{i\parallel }+W_{i\perp }=(\mathrm{3}\mathit{\alpha }+\mathrm{11}/\mathrm{2}{)}^{\mathrm{1}/\mathrm{3}}{m}_i{\left[{\displaystyle \frac{r{D}_{\perp }\left(r\right)}{(\mathrm{3}\mathit{\alpha }+\mathrm{1})}}\right]}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where W_∥ and W_⊥ are the mean parallel and perpendicular energies, respectively; W_i is the total mean energy; i denotes the type of the ion (H⁺ or O⁺); and α is a fitting parameter. In this theory, the mean energy ratio $W_{\perp }/W_{\parallel }$ asymptotically approaches a constant value.

3 Comparisons

In this section, we will compare Monte Carlo simulations with mean particle theory estimates in the auroral and polar-wind regions (Sect. 3.1). Additionally, we will present a comparison among Monte Carlo simulations, mean particle theory estimates, and relevant observational data in Sect. 3.2.

3.1 Comparison between Monte Carlo simulations and estimates of mean particle theory

By employing the Monte Carlo technique (specifically the Barghouthi model) alongside mean particle theory, we have derived varying altitude profiles for the mean parallel, mean perpendicular, and total mean energies across different regions of Earth's magnetosphere – specifically, the auroral region (Fig. 1) and the polar-wind region (Fig. 2) – for H⁺ (panels a–c) and O⁺ (panels d–f) ion outflow. In both methods, we utilized an appropriate D_⊥(r) as the velocity diffusion coefficient for each region. Overall, our findings in the auroral region (Fig. 1) indicate an excellent agreement between Monte Carlo simulations and estimates of mean particle theory for H⁺ and O⁺ ions. It is important to highlight that, in the mean particle theory, ion heating and energization occur solely due to wave–particle interactions, which primarily depend on the value of the velocity diffusion coefficient. In contrast, the energization process in the Monte Carlo model results from the competition between wave–particle interactions and external forces, such as gravity, the polarization electric field, and the divergence of the geomagnetic field. The close agreement between the two methods can be attributed to the dominant influence of wave–particle interactions, which prevail over the effects of external forces due to the high values of D_⊥(r). In other words, wave–particle interactions are the primary drivers of the energization process for both H⁺ and O⁺ ions in this auroral region. For the polar-wind region and for ${\mathit{\lambda }}_{\perp }\to \mathrm{\infty }$, this means that the velocity diffusion coefficient D_⊥(r) is altitude-dependent. The values of the velocity diffusion coefficient are less than those in the auroral region; see Eqs. (4) and (5). According to Fig. 2, there is no agreement between Monte Carlo simulations and mean particle theory estimates; this is due to the contribution of the external forces in Barghouthi model that competes with the effect of wave particle interaction. However, in mean particle theory, the external forces are not considered, and the heating is due to wave particle interaction. We report here that mean particle theory is not suitable to be used in this region, and Monte Carlo simulations are more appropriate to use, as shown in Barghouthi et al. (2011), where they compared their Monte Carlo results with observations. Also, we have found that, when we use the diffusion coefficient that depends on altitude only, its value becomes very large as altitude increases. Therefore, the values of the particle energies obtained from Eqs. (16) to (18) turn out to be very high, but when the diffusion coefficient becomes velocity- and altitude-dependent according to Eq. (8) (i.e., $k_{\perp }{v}_{\perp }/\mathrm{\Omega }\ge \mathrm{1}$), the produced particle energies turn out to be reasonable, as shown in Fig. 3 (solid blue lines and dashed red lines).

Figure 1 Altitude profiles comparing the mean particle theory estimates (solid blue lines) under auroral conditions with Monte Carlo calculations (dashed red lines). Panels (a)–(c) represent H⁺ ions, while panels (d)–(f) illustrate O⁺ ions. Panels (a) and (d) show the mean perpendicular energy W_⊥, panels (b) and (e) present the mean parallel energy W_∥, and panels (c) and (f) depict the mean total energy W_total.

Figure 2 Comparison of the estimates of the mean particle theory (solid blue lines) for polar-wind conditions with the Monte Carlo calculations (dashed red lines). (d–f) O⁺ ions. (a–c) H⁺ ions. The mean perpendicular energy W_⊥ is represented by (a, d), the mean parallel energy W_∥ is represented by (b, e), and the total energy is represented by (c, f).

Figure 3 Comparison of mean particle theory estimates (solid blue lines) and Monte Carlo calculations (dashed red lines) for auroral conditions (a and c) and polar-wind conditions (b and d), along with observational data represented by maximum (dotted lines), average (dashed-dotted lines), and minimum (black dashed lines) values for O⁺ ions. The top panel shows the mean perpendicular energy, while the bottom panel presents the mean parallel energy. In this analysis, we have considered the wavelength of electromagnetic turbulence to be ${\mathit{\lambda }}_{\perp }=\mathrm{8}\mathrm{km}$.

3.2 Comparison between Monte Carlo simulations, estimates of mean particle theory, and available observations

In this section, we present only the simulation results and the estimates of mean particle theory that have corresponding observations. Barghouthi (2008) compared the Monte Carlo simulation results obtained using the Barghouthi model with the corresponding observations for H⁺ and O⁺ ion outflows in the auroral region at different altitudes in the simulation tube; he obtained excellent agreement, particularly when the typical perpendicular wavelength of the electromagnetic turbulence was 8 km. Also, he observed that there is a broad agreement between the simulation results of the polar wind for this wavelength and the corresponding observations. For these reasons, we chose to use the results of the comparison with corresponding observations when ${\mathit{\lambda }}_{\perp }=\mathrm{8}\mathrm{km}$, i.e., when the velocity diffusion coefficient is altitude- and velocity-dependent. We will compare the outcomes of our Monte Carlo simulations, the estimates of mean particle theory, and available observations obtained from different published articles. Observations of O⁺ ions at various altitudes were obtained for parallel velocity, perpendicular temperature, and parallel temperature for both polar-wind and auroral regions from Nilsson et al. (2013), and observations for parallel velocity and perpendicular temperature for O⁺ ions in the central-polar-cap and cusp regions were obtained from Barghouthi et al. (2016). In this study, we adopt the terms “good” and “excellent” agreement to qualitatively describe how well our simulation and theory curves correspond to observational data. These classifications are based on a visual comparison with the observational ranges, represented by the maximum, average, and minimum values. We define agreement as excellent when the simulation and theory curves largely overlap with the average observational curves and remain within the observed range across all altitudes. Conversely, we describe the agreement as good when the curves generally fall within the envelope formed by the observational maximum and minimum values, even if they deviate partially from the average trend. This qualitative approach is commonly used in comparative studies of ionospheric outflows (e.g., Barghouthi et al., 2011, 2012, 2016), and it acknowledges the significant variability typically present in observational data.

Figure 3 displays the comparison results for the auroral (panels a and c) and polar-wind (panels b and d) regions. It contrasts the Monte Carlo method (dashed red lines) and the estimates from mean particle theory (solid blue lines) at ${\mathit{\lambda }}_{\perp }=\mathrm{8}\mathrm{km}$, alongside observational data represented by maximum (dotted lines), average (dashed-dotted lines), and minimum (dashed black lines) values for O⁺ ions. Across all four panels, it is evident that the results from the Monte Carlo simulations and the mean particle theory estimates are closely aligned, demonstrating excellent agreement at lower altitudes. At higher altitudes, while both methods display similar qualitative behavior and maintain good agreement, we attribute this acceptable and reasonable correlation to the application of a suitable altitude- and velocity-dependent diffusion coefficient in both approaches.

Based on the comparison with observations, it is evident that the simulation results and estimates from mean particle theory fall within the observational range. This indicates that $D_{\perp }(r,v_{\perp })$ is appropriate for use in both the Monte Carlo simulation and mean particle theory. Specifically, the results for perpendicular energy W_⊥ (O⁺) from both the Monte Carlo method and mean particle theory align well with the maximum observed values in the auroral region, and they also exhibit good agreement with the average observed values in the polar-wind region. Regarding the mean parallel energy W_∥ (O⁺), as shown in Fig. 3c and d, the simulation results and estimates from mean particle theory are closely aligned with both the average and maximum observed values in the auroral region, and they demonstrate excellent agreement with the minimum observed values in the polar-wind region.

In the central-polar-cap (Fig. 4a) and cusp (Fig. 4b) regions, we utilized the altitude-dependent diffusion coefficient from Nilsson et al. (2013), specifically when ${\mathit{\lambda }}_{\perp }\to \mathrm{\infty }$. This approach resulted in excellent agreement at all altitudes between the Monte Carlo simulations and the estimates of mean particle theory. Furthermore, we observed a strong correspondence to observational data in both regions; the results closely match the maximum and average values in the central-polar-cap region and align well with the average values from observations in the cusp region.

Figure 4 Comparison of the estimates of mean particle theory (solid blue lines) with the Monte Carlo calculations (dashed red line) for central-polar-cap conditions (a) and cusp conditions (b) in addition to the observations (maximum (dotted lines), average (dotted-dashed lines), and minimum (dashed black lines)) for O⁺ ions. In this analysis, we have considered the wavelength of electromagnetic turbulence ${\mathit{\lambda }}_{\perp }\to \mathrm{\infty }$.

As a result of these comparisons, it is crucial to select the appropriate form of the velocity diffusion coefficient for each region of interest within Earth's magnetosphere when analyzing ion energization and heating processes.

Barghouthi (2008) conducted a comparative study between Monte Carlo simulations – based on the Barghouthi model – and in situ measurements of O⁺ and H⁺ ions at selected altitudes in the auroral region. This work was later extended by Barghouthi et al. (2011) to include the polar-wind region. In both studies, the simulation results sometimes showed good agreement with observational data in terms of energy magnitudes and general trends; however, notable discrepancies were also identified. These deviations were attributed to physical processes not incorporated into the model, such as variations in ion heating or cooling mechanisms, the influence of unmodeled external forces, or other dynamic plasma phenomena. Importantly, these earlier studies did not incorporate mean particle theory, which emphasizes ion heating due to wave–particle interactions.

In the present study, both Monte Carlo simulation results and theoretical estimates based on mean particle theory were found to fall within the observational range, as shown in Figs. 3 and 4. Specifically, we noted that, in Fig. 3a and b, covering altitudes between 6 and 8 Earth radii (R_E), although the simulated ion energies align in terms of magnitude with observed values, the altitude-dependent energy profiles differ significantly in terms of trend. Theoretical and simulated profiles exhibit a steeper gradient with altitude compared to the observational data. Similar behavior can be seen for Fig. 4a, where the theoretical profiles show steeper gradients than those observed, particularly at higher energies. In contrast, Fig. 4b demonstrates a stronger consistency between theoretical predictions and observations.

In our analysis, we examined both absolute energy values and the slopes of the energy profiles. The discrepancies noted in Fig. 3a and b – especially the steeper theoretical and simulated slopes – can be attributed to the use of altitude- and velocity-dependent diffusion coefficients that fully account for processes like ion self-limiting heating (Bouhram et al., 2004) and localized cooling effects. The deviations in Fig. 4a are mainly due to oversimplified assumptions regarding wave mode propagation and the absence of spatial variations in plasma conditions. By contrast, better agreement in Fig. 3c and d and Fig. 4b is achieved due to the implementation of empirically based, altitude- and velocity-dependent diffusion coefficients.

These findings underscore the importance of region-specific parameterizations and highlight the need for future improvements in the modeling framework. In particular, incorporating saturated wave–particle interactions, as well as multidimensional plasma dynamics, would enhance the model's ability to replicate the observed energy profiles more accurately.

4 Conclusions

We conducted a comparative analysis of energy components – mean perpendicular energy W_⊥, mean parallel energy W_∥, and total mean energy W_total – for H⁺ and O⁺ ions using both the Monte Carlo method and estimates from mean particle theory across various regions of Earth's magnetosphere, including the polar-wind, auroral, central-polar-cap, and cusp regions. Utilizing altitude-dependent diffusion coefficients, we found excellent agreement between the two methods in the auroral regions; however, there was a lack of concordance in the polar-wind region. Moreover, we observed that energy values at middle and high altitudes were unrealistically elevated compared to observational data. To mitigate this issue, we implemented velocity- and altitude-dependent diffusion coefficients, which yielded reasonable energy values at low, middle, and high altitudes across both regions. When we compared the results of the Monte Carlo simulations and the mean particle theory estimates with available observational data, we noted good agreement throughout all studied regions – polar wind, auroral, central polar cap, and cusp. Both the simulation results and the mean particle theory estimates aligned well within the observed data range. We conclude that the velocity diffusion coefficients used in both methods, which produced acceptable consistency with the observations, are appropriate for application in these regions.

For future research, it is essential to conduct further observations in various regions of Earth's magnetosphere. This will facilitate additional comparisons and help to identify the most suitable diffusion coefficient, as well as to ascertain which method yields the most accurate results in relation to the corresponding observational data.