Statistical Observations of Proton‐Band Electromagnetic Ion Cyclotron Waves in the Outer Magnetosphere: Full Wavevector Determination

Electromagnetic Ion Cyclotron (EMIC) waves mediate energy transfer from the solar wind to the magnetosphere, relativistic electron precipitation, or thermalization of the ring current population, to name a few. How these processes take place depends on the wave properties, such as the wavevector and polarization. However, inferring the wavevector from in‐situ measurements is problematic since one needs to disentangle spatial and time variations. Using 8 years of Magnetospheric Multiscale (MMS) mission observations in the dayside magnetosphere, we present an algorithm to detect proton‐band EMIC waves in the Earth's dayside magnetosphere, and find that they are present roughly 15% of the time. Their normalized frequency presents a dawn‐dusk asymmetry, with waves in the dawn flank magnetosphere having larger frequency than in the dusk, subsolar, and dawn near subsolar region. It is shown that the observations are unstable to the ion cyclotron instability. We obtain the wave polarization and wavevector by comparing Single Value Decomposition and Ampere methods. We observe that for most waves the perpendicular wavenumber (k⊥) is larger than the inverse of the proton gyroradius (ρi), that is, k⊥ρi > 1, while the parallel wavenumber is smaller than the inverse of the ion gyroradius, that is, k‖ρi < 1. Left‐hand polarized waves are associated with small wave normal angles (θBk < 30°), while linearly polarized waves are associated with large wave normal angles (θBk > 30°). This work constitutes, to our knowledge, the first attempt to statistically infer the full wavevector of proton‐band EMIC waves observed in the outer magnetosphere.

Using a two-species cold plasma linear wave theory (protons and electrons), EMIC waves correspond to Left Hand Polarized (LHP) waves with large ellipticity values that exist below the proton cyclotron frequency (ω ci ), with the wavevector aligned to the background magnetic field (André, 1985).They correspond to the same dispersion surface as Alfven waves, when the approximation ω ≪ ω ci does not hold.When cold ions are accounted for, as is generally the case of the Earth's magnetosphere (Toledo-Redondo, André, et al., 2021), the waves can be observed with large growth rate at oblique wave normal angles, with Linear and Right Hand Polarization (LP and RHP, respectively), and also associated with small ellipticity values (J.H. Lee et al., 2019;Omidi et al., 2013;Toledo-Redondo, Lee, et al., 2021).Figure 1 shows the estimated growth rate (a-c) and ellipticity (d-f) for three typical situations in the Earth's magnetosphere, using Waves in Homogeneus Anisotropic Magnetized Plasma (WHAMP) (Roennmark, 1982).The three simulations have B = 40 nT and a hot proton component of T ih = 5 keV, temperature anisotropy T ⊥ /T ‖ = 1.8, and density n ih = 0.5 cm 3 .In addition, simulations 2 and 3 have a cold ion component with temperature T ic = 50 eV, T ⊥ /T ‖ = 1, and densities n ic,2 = 0.5 cm 3 and n ic,3 = 1 cm 3 .There is no relative drift between the two proton components.In the absence of cold ions, positive growth rate occurs for parallel propagation and corresponds to LHP.On the other hand, in the presence of cold ions, positive wave growths extend to oblique wave normal angles, corresponding to LP and RHP.Inferring the k vector of EMIC waves from in-situ spacecraft measurements is a challenging task.Disentangling time variations and space variations cannot be done in the general case if only one-point measurements are available.Using constraints derived from linear wave theory, one can infer the most probable direction of k as a function of frequency from one-point observations of the magnetic field vector using the Single Value Decomposition (SVD) method.If we also have one-point observations of the electric field vector, we can exploit the linearized Faraday's law k × E = ωB to infer the most probable wavevector k as a function of frequency (Oscarsson et al., 1997(Oscarsson et al., , 2001;;Rönnmark & André, 1991;Santolík et al., 2003;Stenberg et al., 2002;Tjulin et al., 2005).If one can measure the current density in the plasma, then it is possible to infer the most probable vector wavevector as a function of frequency exploiting the linearized Ampere's law ik × B = μ 0 J (Bellan, 2016;S. K. Vines et al., 2021).Using four-spacecraft measurements in tetrahedron configuration, one can compare the magnetic field or electric field vectors of the four spacecraft and infer k.These algorithms include for instance the wave telescope (Neubauer & Glassmeier, 1990), k-filtering (Pinçon & Lefeuvre, 1991), phase differencing technique (Pakhotin et al., 2013), or the Dispersion Relation From Timing (DRAFT) methods (W.Zhang et al., 2022).
Toledo-Redondo, Lee, et al. (2021) reported MMS observation of a proton-band EMIC wave with large wave normal angle (θ Bk = 74°) and linear polarization, with a similar number density of anisotropic ring current ions (hot ions) and ionospheric-originating cold ions (n h ∼ n c ∼ 0.5 cm 3 ).Using WHAMP, they showed that the wave was unstable (positive wave growth) for the observed plasma parameters and wave normal angle, and suggested that the presence of cold protons enables wave generation at oblique wave normal angles, which are associated with linear polarization, that is, ɛ close to zero.
Previous statistical studies of EMIC waves conducted in the outer dayside magnetosphere have shown that the wave normal angle of proton-band EMIC waves is typically larger than 30°for all MLT, and above 45°for 03:00 < MLT <09:00 (R. C. Allen et al., 2015;Min et al., 2012;Wang et al., 2017).Min et al. (2012) and R. C. Allen et al. (2015) also showed that while left-hand polarization (ɛ < 0.1) is in general the most probable polarization, a significant portion (varying with MLT and MLAT) of right-hand polarization (ɛ > 0.1) and linear polarization ( 0.1 < ɛ < 0.1) waves is observed, and that typically |ɛ| < 0.4.These results are in contrast with Wang et al. (2017), who found a predominance of RHP and LP waves using MMS data (September 2015-October 2016).They argue that the observed LP and RHP can correspond to deviations from the original LHP due to propagation in the inhomogeneous dipole magnetic field, as shown by Hu and Denton (2009) and Hu et al. (2010) using hybrid simulations.However, they do not rule out the possibility of these waves being generated with LP and RHP.Previous studies using AMPTE mission (equatorial orbit) (B.J. Anderson et al., 1992aAnderson et al., , 1992bAnderson et al., , 1996) ) showed a predominance of LHP in the dusk sector and a predominance of LP in the dawn sector.They also showed that the normalized central frequency of the waves (ω peak /ω fci ) was larger at dawn than at dusk.
In this work, we analyze the Magnetospheric Multiscale (MMS) mission data available to date (September 2015-March 2023) to statistically study the properties of proton-band EMIC waves in the Earth's magnetosphere.We combine the SVD and Ampere methods to infer the full wavevector of the observed waves.Section 2 describes the MMS data set used in this work.Section 3 presents the methods used to detect the EMIC waves and to infer the wavevector properties.Section 4 presents the statistical results obtained from almost 8 years of observations.Section 4 discusses the results and the implications of our findings.Finally, Section 5 presents the conclusions of this work.

Data Set
The MMS mission (Burch et al., 2015) consists of four identical spacecraft flying in tetrahedron formation with spatial separation and time resolution enough to resolve electron scales in near-Earth space plasma environments.It was launched in February 2015 and successfully inserted into a highly elliptical equatorial orbit.It has been gathering high resolution in-situ data at various geocentric distances, from few Earth radii (R E ) to 25 R E .In this work, we make use of data from the Flux Gate Magnetometer (FGM) (Russell et al., 2014) and the Fast Plasma Investigation (FPI) (Pollock et al., 2016) to measure ion density, velocity and temperature moments of the 3D velocity distribution function.We use magnetic field data at a 16 samples/s cadence (fast survey mode) and ion data at a 0.25 samples/s cadence (fast survey mode).
We initially include all MMS data available between the dates 2015-09-01 and 2023-05-31.During that period, the orbit of MMS has undergone modifications, including apogee variations (from 12 R E to 25 R E ), inter-spacecraft distance (from several km to hundreds of km), and shape modifications (tetrahedron, string of pearls, etc).Our science interest is on EMIC wave occurrence and properties in the dayside magnetosphere, and therefore we make a first reduction of the initial data set following these criteria: 1. FGM data are available for all 4 spacecraft and FPI ion data are available on MMS1, both at fast survey cadence (burst cadence is not a requirement for this study).2. MMS1 is at dayside (X GSE > 0) and R GSE < 15R E , that is, the typical region for sampling dayside magnetosphere plasma.
The data available after this first reduction amounts to more than 6,500 hr of observations, divided into 890 continuous time intervals which correspond to Regions Of Interest (ROI) of the MMS mission.A ROI corresponds to the times when MMS operates at high cadence and all instruments are active, their times are designed to maximize the probability of MMS encountering the magnetopause operating at high cadence, and their definition is different on every orbit.Next, we divide the data into segments of 64 s, with each segment overlapping 30% with the previous segment, and obtain 490,582 segments.We select 64 s in order to obtain sufficient time and frequency resolution in Fourier spectrograms.For each of the segments, we impose that the magnetic field direction is constant within 2.5°, to ensure we select segments inside the magnetosphere and discard magnetopause crossings and magnetosheath intervals.We obtain 6,768 segments to which we can apply the method for EMIC detection detailed in Section 3.1, that is, roughly 1,581 hr of data, after taking out the overlapped portions of segments.

Methods
In this section, we describe three key algorithms that will be used to analyze the data.The first algorithm, adapted from Bortnik et al. (2007), checks whether waves in the [0.25-1]ω ci band exist for each given 64 s segment.The second algorithm uses the SVD method (Means, 1972;Samson & Olson, 1980;Santolík et al., 2003) to determine the following wave parameters: planarity, degree of polarization, ellipticity and angle between background magnetic field and wavevector (θ Bk ).Finally, the third algorithm determines k by using cross-correlations between the current density (J) and the magnetic field (B) in the Fourier domain, following Bellan (2016) and S. K. Vines et al. (2021).

Proton-Cyclotron-Frequency Band Wave Detection Algorithm
The basic idea is to search for high wave power in the [0.25-1]ω ci band in a 64 s segment that are well above the median value during the corresponding magnetospheric time of the ROI of the MMS mission, which typically lasts up to several hours.Following (Bortnik et al., 2007), we first compute a covariance matrix (c) of the magnetic field vectors in the frequency domain for each segment (i) inside a region of interest where X(ω), Y(ω), Z(ω) are the three components of the B fast-Fourier transformations, and * indicates complex conjugate.We then combine the off-diagonal terms of c and define C i (ω) represents the squared cross-covariance and has units of T 4 Hz 2 .Using the off-diagonal terms makes C i (ω) more immune to random noise (Means, 1972).We then normalize each frequency coefficient of C i (ω) to its median value in the whole ROI, , where ˜denotes median.
We show an example of the algorithm performance in Figure 2. in log-scale.The requirement of a single peak is intended to avoid wave superposition, which makes k vector estimation unfeasible.For each segment selected, the algorithm extracts the central frequency (ω peak ), bandwidth (Δω), and amplitude of the fluctuations (ΔB).The central frequency (ω peak ) corresponds the frequency where C norm is maximum, see Figure 2b, and the bandwidth is defined as Δω = ω 2 ω 1 , that is, the span in frequency of energies above the threshold (see Figure 2b).

Single Value Decomposition Method
The Single Value Decomposition (SVD) method has been widely used for extracting wave parameters of in-situ observations of waves, including for instance chorus, whistler or EMIC waves (R. C. Allen et al., 2015;Santolík et al., 2014;Taubenschuss & Santolík, 2019).The algorithm we apply follows Santolík et al. (2003), and the basic idea is that, given a wavelet transform of the magnetic field measurement, B(ω), one can find the wavevector direction for each frequency k′ that minimizes B(ω) ⋅ k′.For a single, monochromatic, linear plane wave free of noise measurement, B(ω) ⋅ k′ should reduce to zero and k′ ∝ k.Our implementation of the method only finds the most probable azimuthal (ϕ Bk ) and polar (θ Bk ) angles of k in a semisphere, that is, the wave normal angles.The SVD method can provide full wavectors if both electric and magnetic fields are considered, by means of the linearized Faraday's law, cf.Section 1.However, measuring the small-amplitude electric field fluctuations of EMIC waves in the outer magnetosphere is challenging.Outside the plasmasphere, density is low and spacecraft potential (V sc ) is positive.Here supersonic beams of ionospheric cold (eV) ions are common and can create an extended wake behind the spacecraft.The wake will be filled with electrons and the negative space charge density can create a substantial local electric field near the satellite.This electric field does not correspond to the largescale geophysical electric field relevant for studies of EMIC waves (André et al., 2021;Toledo-Redondo et al., 2019).EMIC waves grow and propagate better when cold ions are present, see Figure 1.Furthermore, ASPOC operations (Torkar et al., 2014) also affect the electrostatic potential structure around MMS and have effects over the electric field measurements.For all these reasons, we decided not to include electric field measurements to infer the statistics of the wavevectors in this work.
For each segment (64 s) where we detect a proton-cyclotron frequency band wave (see Section 3.1), we apply the SVD method to the magnetic field and compute spectrograms of planarity, degree of polarization, ellipticity, θ Bk and ϕ Bk , see Figure 3.The Degree Of Polarization (DOP) and Planarity (P) parameters (Santolík et al., 2003) serve as figures of merit of the solutions found.We consider only solutions with both DOP(t, ω) > 0.85 and P(t, ω) > 0.85, that is, we disregard all data points where DOP or P are below 0.85 in Figures 3d-3f.In addition, we consider only data points for the 25 time bins of the 64 s (each time bin is 2.56 s) in the six closest frequency bins to the wave  central frequency ω peak , selected by the wave detection algorithm (black dashed lines in Figure 3).To consider a SVD result reliable for a segment, we impose that at least half of the surviving points in the ellipticity (ɛ) spectrogram are within ±0.15 of the central value of ellipticity, that is, the ɛ measurement is stable in a 64 s time-scale around the central wave frequency ω peak .In a similar way, we also impose that at least half of θ Bk (t, ω) computed angles (6 frequency bins times 25 time bins) are within ±15°to consider a reliable θ Bk = θ Bk (t, ω).In addition, at least 70% of the 25 time bins must contain one valid θ Bk to accept the result.If the wave polarization is linear, that is, ɛ = 0, SVD has infinite solutions for the direction of the wavevector.Therefore, we consider only angles obtained by SVD if |ɛ| > 0.3.

Full Wavevector Determination
We make use of MMS constellation to extract the full wavevector, that is, direction and magnitude.We choose a method that is based on obtaining the current density vector, J, following Bellan (2016) and S. K. Vines et al. (2021).We approximate the current density using the curlometer technique (Dunlop et al., 1988).The results of this technique strongly depend on two factors: the tetrahedron size, L 4sc , and the tetrahedron quality.With respect to the tetrahedron size, it determines the minimum resolvable spatial scale of J.We consider only results that provide wavelengths λ/2 > L 4sc .For assessing the tetrahedron quality, we use the geometric factor Q RR (Robert et al., 1998), which combines both the planarity and elongation parameters of the tetrahedron.We impose that the geometric factor Q RR < 0.5 to compute the current density using curlometer.We also impose that ∇ ⋅ B/ ∇ × B < 0.7.We chose the 0.7 value because it keeps enough events for the statistics while reducing the amount of events where spatial differentiation yielded incorrect results.Once J is determined, it is possible to obtain the wavevector corresponding to a particular frequency, k(ω peak ), based on the linearized Ampere's law in the Fourier domain The previous equation assumes plane-wave solutions with a single monochromatic propagating packet and neglects any non-linear wave effects.Obviously, this is not necessarily the case in real data and therefore one needs to be cautious when applying this technique to a large data set.Up to now, the technique has provided reliable results using MMS tetrahedron for case studies, (e.g., Gershman et al., 2017;Toledo-Redondo, Lee, et al., 2021).Figure 4 shows an example of wavevector determination following the procedure described by Bellan (2016).We define the wave central frequency where the maximum power spectral density of B is found (Figure 4a).Then, we compute J(ω) × B*(ω), and select its imaginary part (I).We select all frequencies near the central frequency that fulfill I(F J×B ) > I(F J×B ) max /2.For the case in Figure 4, only the central frequency vector is considered (Figure 4b).Finally, we obtain the wavevector following Equation 3. We average among all frequencies with I(F J×B ) > I(F J×B ) max /2 and impose that the wavevectors from selected frequencies agree within 25°, otherwise k cannot be determined.If only one frequency is above the threshold, we accept that value as the wavevector, as in Figure 4c.

Statistical Results
We obtain 18,493 segments of 64 s of duration each where a single proton-band peak is found, following the steps described in Section 3.1.From now on we will refer to them as events.This amounts to roughly 274 hr of protonband EMIC wave observations, after taking out the partial overlaps between segments.This indicates that roughly during 15% of the time spent by MMS inside the dayside magnetosphere, a proton-band EMIC wave was present, although the observation of the wave inside an event can be shorter than 64 s. Figure 5a shows the relative occurrence (total events time divided by dwell time) projected to the GSM equatorial plane.Figure 5b shows the mean normalized central frequency, ω peak /ω ci , for each bin in Figure 5a. Figure 5c shows the mean normalized amplitude of the wave fluctuations, ΔB/B 0 , and the wave relative bandwidth, Δω/ω peak , is shown in Figure 5d.We observe a dawn-dusk asymmetry in the normalized central frequency of the detected waves.
We plot histograms of the normalized central frequency (ω peak /ω ci ), the normalized amplitude of the wave fluctuations ΔB/B 0 , and the wave relative bandwidth (Δω/ω peak ) in Figure 6.All normalization quantities correspond to the parameters observed locally, that is, in the vicinity of the wave observation.We present histograms for all events, dusk events (Y GSE < 0), and dawn events (Y GSE > 0).The observed central frequencies, wave fluctuations and relative bandwidths have been fitted to a lognormal distribution (4) using solid lines.The mean, μ, and standard deviation, σ, are provided in Table 1.The wave fluctuations (6b) and wave bandwidths (6c) follow a log-normal distribution, indicating that the measurements are independent and no particular trends are hidden in the data.Therefore, the mean and standard deviations provided in Table 1 are representative in a statistical sense.On the other hand, the normalized central frequencies (6a) do not have the same properties for dawn and dusk events, a double peak is present for dawn events, and the log-normal fits are not good approximations of the data, indicating that the ensemble of the events do not belong to a single distribution with defined mean and standard deviation, in a statistical sense.Doppler shift corrections are not  Results of these fits are provided in Table 1.
important because the Alfvén velocity associated with the wave detections is typically several hundred to few thousands of km, while the observed plasma velocities are always below 80 km/s.
To further analyze the distribution of wave central frequencies, we arrange the data in the ω ci ω peak parameter space (Figure 7). Figure 7a shows a 2D histogram (counts) of all events.For each bin in Figure 7a, we plot the average Y GSM in Figure 7b.Waves with higher ω peak /ω ci are mainly observed in the dawn region.For dusk observations, the data roughly aligns with ω peak /ω ci = 0.39 (Figure 7c).For dawn observations, two separate classes of waves (1 and 2) are found (Figure 7e), depending on their normalized central frequency and roughly divided by the dashed red line.We perform linear regression analysis separately for classes 1 and 2, obtaining ω peak,1 /ω ci,1 = 0.41 and ω peak,2 /ω ci,2 = 0.67, see Figures 7g and 7h.The dusk waves and the dawn class 2 waves span over a wide range of Y GSE (Figures 7f and  7g), while most of the observations that belong to dawn class 1 are found not far from the subsolar region, that is, moderate Y GSE values (see Figure 7f).
In order to establish which type of wave(s) we are detecting, we plot the events in the proton temperature anisotropy, T ⊥ /T ‖ , versus parallel beta, β ‖ , parameter space.For each wave detection (64 s duration), we computed the mean perpendicular and parallel proton temperatures measured by FPI, and then obtained the mean proton temperature anisotropy.Note that we imposed that the direction of B cannot change more than 2.5°during one 64 s interval, cf.Section 2, and therefore the magnetic field direction is well defined for each wave detection.For computing the parallel plasma beta associated to each wave detection, we computed the mean proton density, n, and parallel proton temperature P i,‖ in the interval from FPI to obtain the parallel proton plasma pressure, P i,‖ = nk b T i,‖ .We also computed the mean magnetic field magnitude to obtain the magnetic pressure, P B = B 2 / (2μ 0 ). Figure 8a shows a 2D histogram of the events.Figure 8a shows the dwell time of MMS1.Bins with less than 15 min of observations have been discarded.The median anisotropy for all segments in our database, which is representative of the outer dayside magnetosphere, cf.Section 2, is T ⊥ /T ‖ = 1.26.Figure 8b shows the relative wave observation time, where bins with less than 5 wave detections have been discarded.Wave occurrence is larger for larger proton temperature anisotropies.Following Verscharen et al. (2016), we computed the growth rate, γ, thresholds for various anisotropy-driven microinstabilities, namely Ion-Cyclotron (IC), mirror, Fast-Magnetosonic, (FM) and firehose instabilities.The 2D histogram indicates that the most probable instability at work is the ion-cyclotron instability (solid black lines).Figure 8c shows the median value of the wave fluctuations, ΔB/B 0 , for each bin in Figure 8a.The largest wave fluctuations are observed for large growth rates of the ion-cyclotron instability, indicating the ion-cyclotron nature of the wave fluctuations.Figure 8d shows the median value of MLAT, which indicates that unstable regions and wave generation are more common in off-equator regions (MLAT 30°), as previously reported (R. Allen et al., 2013;Liu et al., 2012;S. Vines et al., 2019).

Full k Vector Determination
The Ampere method described in Section 3.3, based on Bellan (2016) and S. K. Vines et al. (2021), allows us to obtain the full wavevector, k.However, it assumes that a single, plane, monochromatic wave is present in the data.In addition, it relies on the plasma current density (J) measurement, which is obtained using the curlometer technique.Inferring J from curlometer depends on the spacecraft tetrahedron size, which varies along the MMS data set, and gives an estimate which is representative of the volume filled by the MMS constellation.Our approach to trust the k measurements is as follows.For each event, we compute two independent methods, SVD and Ampere.SVD provides the angle between k and B. If θ Bk from the SVD and Ampere methods agree within 25°, we trust the k provided by the Ampere method.It is also worth noting that the SVD method degenerates (has infinite solutions) for waves with linear polarization (ɛ = 0).Therefore, we only use θ Bk angles from SVD when | ɛ| > 0.3.When |ɛ| < 0.3 from SVD, we accept the θ Bk obtained by the Ampere method without comparison.
Of the initial 18,493 events (64 s each), we find reliable ellipticities ɛ for 5,770 events, that is, DOP and planarity are above 0.85 near the wave frequency, and ɛ near the wave frequency is consistent for most of the interval.However, only 1,598 of them provide a reliable θ Bk .This is either owing to |ɛ| < 0.3 or that θ Bk appears to change in the 64 s period.On the other hand, The Ampere method, see Section 3.3, provides 2,079 wavevectors from the initial 18,493 events.The overlap within the two methods is 740 events.571 out of 740 have |ɛ| < 0.3 and we accept k obtained from the Ampere method.There are 70 events, out of 169, where θ Bk from the two methods agree within 25°and we also accept k obtained from the Ampere method, making a total of 641 accepted k vectors.There are 99 events where the two methods do not provide consistent results, and these cases are left for future investigation.See Table 2 for a summary of these events statistics.
We use the 641 events to produce Figure 9, where a 2D histogram of the occurrences in the k ‖ ρ ik ⊥ ρ i plane is presented (Figure 9a).It is interesting to note that most observations fall in the kinetic regime, that is, the k ⊥ ρ i ≥ 1 (also k ⊥ d i ≥ 1, not shown).It is also noteworthy that most of the observations exhibit rather oblique θ Bk angles.
Figure 9b shows the mean value of the ellipticity for each bin in Figure 9a (bins with only count have been discarded in Figures 9b-9d).Negative ellipticities correspond to LHP waves, and the largest values correspond to small θ Bk angles.For θ Bk > 45°, the ellipticities are closer to 0, that is, linear polarization.Figures 9c and 9d show the ellipticities for LHP and RHP waves only, respectively.

Discussion
Our algorithm finds a proton-band [0.25-1]ω ci EMIC wave on 15% of the checked segments, that is, dayside magnetosphere segments, indicating a large prevalence of these waves in the magnetosphere.We confirm a persistent dawn-dusk asymmetry in terms of the normalized wave frequency (ω peak /ω ci ), previously reported by B. J. Anderson et al. (1992b).Grison et al. (2021) showed that EMIC waves are more common near the magnetopause in the dawn region, making them more sensitive to solar wind pressure pulses.We confirm this tendency in the dawn region, as it can be seen in Figure 5a.Overall, the observations in the dusk magnetosphere present a mean ω peak /ω ci ∼ 0.39, see Figure 7c, while the dawn magnetosphere observations present two distinct classes, 1 and 2, according to normalized central frequency, see Figure 7e, with ω c1 /ω ci ∼ 0.41 and ω c2 /ω ci ∼ 0.67 respectively (Figures 7g and 7h).Class 1 is observed closer to the subsolar point, while class 2 is observed toward  the flank, see Figure 7f.Class 1 has a normalized central frequency consistent with the mean normalized central frequency observed for the waves in the dusk magnetosphere, and they are both observed in the region that is typically influenced by plasmaspheric plumes and plasmaspheric detached material (Toledo-Redondo, André, et al., 2021).The dawn class 2 is not likely to be affected by detached plasmaspheric material due to its location in the dawn flank magnetosphere, but may be more affected by the warm plasma cloak.Identifying the cold ion populations present in the plasma is a challenging task, (e.g., Toledo-Redondo et al., 2019), in particular for the type of data (fast survey mode) used, and it is left for future work to study the ionospheric-originating populations involved for each of the observed wave classes.
When the 18,493 events are plotted in the anisotropy (T ⊥ /T ‖ ) versus parallel plasma beta (β ‖ ), one can see that wave observations are more common for large proton temperature anisotropies (Figure 8b).Above the Ion Cyclotron instability threshold γ > 10 4 , wave observations become much more common than below.This is a clear indication that the nature of the detected waves corresponds to ion cyclotron, and allows us to discard mirror, fast magnetosonic, and firehose instabilities as the cause of the observed waves.Yue et al. (2019) also computed statistics of EMIC waves in the proton temperature anisotropy versus parallel plasma beta parameter space using the Van Allen probes and showed similar relative occurrences.
L. Chen et al. (2019) conducted for the first time a statistical wavenumber analysis of EMIC waves in the magnetosphere.Their study focused in the inner magnetosphere, where they used electric and magnetic field measurements from the Van Allen probes (Mauk et al., 2014) to infer wavevectors, which were then used for compute the minimum resonant energy for electrons.The study was limited to quasi-parallel wavevectors.They obtained mean values for the parallel wavevector k ‖ c/ω peak ≃ 400 800, see Figure 1f from L. Chen et al. (2019), which agree with our results shown in Figure 9, which yield a mean value k ‖ c/ω peak ≃ 640.
The present work constitutes, to the knowledge of the authors, the first wavevector statistics of EMIC waves in the outer magnetosphere, thanks to the capabilities of the MMS constellation.A large data set, including roughly 275 hr of proton-band EMIC wave observations has been used in this work.We chose using two independent methods (SVD and Ampere method, see Sections 3.2 and 3.3) to assess the wavevector.However, only 740 out of 18,493 wave events provided reliable results for both methods.We have been restrictive for accepting the results of each method, imposing that the wave properties remain constant over a 64 s period.In addition, if multiple waves are superposed, the methods are expected to fail.Furthermore, SVD cannot provide wave normal angles when the wave polarization is linear (Santolík et al., 2003), and we did not compute wave normal angles based on SVD whenever |ɛ| < 0.3.Overall, we are left with 541 events of quasi-linear polarization |ɛ| < 0.3, plus 191 events where θ Bk from these two methods can be compared.Finally, 70 out of 169 events provide similar, that is, within 25°results for θ Bk by the two methods.It remains a future goal to further investigate the causes of 99 events providing distinct predictions for θ Bk between the two methods.Finally, 641 events are selected to estimate the full wavevector (Figure 9).It is noteworthy that for most observations k ⊥ ρ i > 1 and k ⊥ d i > 1, while k ‖ ρ i < 1 and k ‖ d i < 1.This result indicates that these waves experience kinetic effects most of the time in the magnetosphere, indicating that they exchange energy with resonant ions most of the time.Finally, we also note that the wave normal angle determines the ellipticity, with larger angles being associated to linear polarization, in accordance with previous studies (R. C. Allen et al., 2015;B. J. Anderson et al., 1992b;H. Chen et al., 2019;Kang et al., 2021;Toledo-Redondo, Lee, et al., 2021).

Conclusions
We conducted a statistical analysis of proton-band EMIC waves using 8 years of dayside magnetosphere observations by MMS, resulting in 274 hr of wave observations, which indicate that these waves are present in the dayside magnetosphere roughly 15% of the time, although some events may contain the wave for durations shorter than 64 s.We confirm that the wave central frequencies are larger in the dawn flank magnetosphere than in the dusk or in the dawn closer to the subsolar region.The waves are consistent with an ion cyclotron instability generated by temperature anisotropy (T ⊥ > T ‖ ) and exhibit larger wave amplitudes for larger predicted wave growths, while the mirror, fast magnetosonic and firehose instabilities can be ruled out as the source of the observed fluctuations.We find that most observations satisfy k ⊥ ρ i > 1 and k ⊥ d i > 1.Finally, we find that small wave normal angles, θ Bk ≤ 30°, are observed with large ellipticities and left-hand polarization, while waves with θ Bk ≥ 30°are associated with linear polarization.
Figure 2a shows a spectrogram of C norm during 8.5 hr in 2015-09-29.Accumulations of energy in the band of interest, that is, between [0.25-1]ω ci (black lines), are observed.The criteria for selecting a segment that contains a proton-band wave is shown in Figure 2b.The low-pass filtered C norm (black line) must present only a single peak above the threshold (red line), which is set to 1 Journal of Geophysical Research: Space Physics 10.1029/2024JA032516 Finally, we define the wave fluctuations ΔB = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ C(ω peak ) 4 √ , and the normalized wave fluctuations ΔB/ B 0 = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ C(ω peak )/ C(0) 4 √ .

Figure 2 .
Figure 2. Automatic detection algorithm of proton-band EMIC waves, based on (Bortnik et al., 2007).(a) Magnetic field covariance matrix (squared) spectrogram (C norm ).Black lines indicate f ci and f ci /4.Red line indicates the corresponding time of second panel.(b) C norm and C norm low-pass filtered.The algorithm detects a proton-band wave above the energy threshold (red line), with central frequency ω peak = 0.4ω ci and bandwidth Δω = 0.2ω ci (red dots).

Figure 4 .
Figure 4. Determination of the wavevector (k).(a) Power spectral density of the magnetic field (4-spacecraft averaged, B 4sc ).Black lines determines the central frequency, based on the maximum of B 4sc .(b) Fourier transform imaginary part of J × B. (c) k vector determination as a function of frequency.

Figure 5 .
Figure 5. Spatial distribution of the proton-band waves relative occurrence, central frequency, fluctuations and bandwidth.(a, e, i) 2D histograms of the relative occurrence of waves projected to the XY, XZ, and YZ GSM planes.Bins containing less than 10 events (64 s each) are discarded.(b, f, j) Median value of the normalized central frequency for each bin.There is a prevalence of larger frequencies in the dawn magnetosphere.(c, g, k) Normalized wave fluctuations.(d, h, l) Normalized wave bandwidth.

Figure 6 .
Figure 6.Distribution functions of wave parameters: (a) central frequency, (b) fluctuation amplitude, and (c) bandwidth.(Blue,red, yellow) represent (all, dusk, dawn) wave detections.Solid lines correspond to log-normal fits of the distribution functions.Results of these fits are provided in Table1.

Figure 7 .
Figure 7. 2D histograms of ω peak versus ω ci .(a) Events in the ω peak -ω ci parameter space, for all MMS locations.(b) Median Y GSE value for each event in ω peak -ω ci parameter space, for all MMS locations.(c) Events in the ω peak -ω ci parameter space, for dusk MMS locations.Solid black line shows linear fit of the data.(d) Median Y GSE value for each event in ω peak -ω ci parameter space, for dusk MMS locations.(e) Events in the ω peak -ω ci parameter space, for dawn MMS locations.Dashed red line divides data into two separate classes.(f) Median Y GSE value for each bin in ω peak -ω ci parameter space, for dusk MMS locations.Dashed red line divides data into two separate clusters (classes 1 and 2).(g) Scatter plot in the ω peak -ω ci parameter space, for class 1 waves in the dawn region.Solid black line shows linear fit of the data.(h) Scatter plot in the ω peak -ω ci parameter space, for class 2 waves in the dawn region.Solid black line shows linear fit of the data.

Figure 8 .
Figure 8. Distribution of proton-band electromagnetic ion cyclotron waves as a function of proton temperature anisotropy, T ⊥ /T ‖ , and parallel plasma beta, β ‖ .(a) Dwell time.Bins with less than 15 min of data are discarded.(b) Relative wave observation time, considering dwell time.Bins with less than 5 wave detections are discarded.(c) Median value of relative wave amplitudes, ΔB/B 0 .(d) Median value of MLAT.Instability growth rates, γ, for Ion-Cyclotron (IC), Mirror, Fast-Magnetosonic (MS) and Firehose instabilities are plotted using black and gray lines, following Verscharen et al. (2016).

Figure 9 .
Figure 9. Statistics of the k vector magnitude and orientation.(a) 2D histogram of wave detections in the k ‖ ρ i k ⊥ ρ i plane.(b) Median ellipticity (E).(c) Median ellipticity (E) for left-hand polarized waves only.(d) Median ellipticity (E) for right-hand polarized waves only.Bins with single counts have been discarded in panels (b-d).