Plasma Density and Magnetic Field Fluctuations in the Ion Gyro-Frequency Range Near the Diamagnetic Cavity of Comet 67P

We report the detection of large-amplitude, quasi-harmonic density fluctuations with associated magnetic field oscillations in the region surrounding the diamagnetic cavity of comet 67P. Typical frequencies are ∼ 0.1 Hz, corresponding to ∼ 10 times the water and ≲ 0 . 5 times the proton gyro-frequencies, respectively. Magnetic field oscillations are not always clearly observed in association with these density fluctuations, but when they are, they consistently have wave vectors perpendicular to the background magnetic field, with the principal axis of polarization close to field-aligned and with a ∼ 90 ◦ phase shift with respect to the density fluctuations. The fluctuations are observed in association with asymmetric plasma density and magnetic field enhancements previously found in the region surrounding the diamagnetic cavity, occurring predominantly on their descending slopes. This is a new type of waves not previously observed at comets. They are likely ion Bernstein waves, and we propose that they are excited by unstable ring, ring-beam, or spherical shell distributions of cometary ions just outside the cavity boundary. These waves may play an important role in redistributing energy between different particle populations and reshape the plasma environment of the comet.


Introduction
The plasma environments of active comets are dominated by the interaction of the solar wind (hereafter SW) with newly born cometary heavy ions. These are mainly water group ions H 2 O + and H 3 O + , produced by ionization (predominantly by solar EUV radiation, but also charge exchange and electron impact reactions with the SW and high-energy electrons) of cometary neutral volatiles (mostly H 2 O) over large distances (∼10 5 -10 6 km) in the extensive and diffuse cometary coma. The resulting vast comet-solar wind interaction region hosts an abundance of plasma instabilities, waves and turbulent phenomena, and thus constitutes a formidable natural laboratory for studying such processes. Waves are important in determining many of the properties of the cometary plasma environment. They can, for example, heat or cool plasma populations, produce supra-thermal electrons, reduce plasma anisotropies and gradients, couple different plasma species, and provide anomalous resistivity.
For at least the first three of these comets, the dominant magnetic wave phenomenon in the SW interaction region was found to be very strong hydromagnetic turbulence in the ultralow-frequency (ULF) range, 10.1029/2020JA028592 f < 1 Hz in the spacecraft (S/C) frame, with maximum power near the local water ion cyclotron frequency, in all cases about 10 −2 Hz, and a ∼f −2 power law drop-off at higher frequencies typical for turbulent cascade processes Tsurutani et al., 1995;Tsurutani & Smith, 1986a, 1986b. (Detailed analysis of the spectral properties of the turbulence observed by DS1 at 19P/Borrelly does not appear in the surveyed literature, but Richter et al., 2011, suggests similar turbulence also at this comet.) The prominent "pump wave" near the ion cyclotron frequency has been attributed to instabilities caused by the highly anisotropic velocity distribution of the newly born cometary ions in the SW frame, where they essentially form a ring, beam or combined "ring-beam" distribution in velocity space depending on the angle between the interplanetary magnetic field and the SW velocity. These distributions can lead to the generation of a multitude of ULF instabilities (see, e.g., Tsurutani, 1991, for a review of this topic). For large (∼90 • ), there is the ion cyclotron instability, a parallel-propagating non-oscillatory mode and a fluid mirror instability (Hasegawa, 1969). For small (≲ 70 • ) there is a right-hand resonant helical beam instability (Wu & Davidson, 1972) and a fluid nonresonant or firehose instability. The resulting waves act back on the particle distribution, isotropizing the pitch angle distributions and thus play an important role in the process of incorporating the newly picked up cometary ions into the SW flow (Coates, 2004).
At both GZ and Halley, the pump wave at the water cyclotron frequency was clearly present only for quasi-parallel ( ∼ 0 • ) magnetic field. For quasi-perpendicular ( ∼ 90 • ) field, left-hand waves generated by the ion cyclotron instability were expected but absent (Glassmeier et al., 1989;Tsurutani, 1991). This is in accordance with results by Richardson et al. (1988) that there was little or no pitch angle scattering of the pick-up ions in the quasi-perpendicular regime. At GZ this regime instead featured detections of single-cycle magnetic pulses (solitary waves) with durations close to the local proton cyclotron period and identified as proton cyclotron waves associated with the pick-up of cometary protons resulting from the dissociation of water molecules . It is not clear whether these waves were also present in the quasi-parallel regime and just drowned out there by the heavy ULF turbulence of the water group ions, or if they were really confined to the quasi-perpendicular regime only. Also, some solitary waves close to the water group gyro-period were reported in Tsurutani et al. (1990). These latter waves were later found near the Earth's bow shock (Schwartz et al., 1992) and were given the name "SLAMS." While the spectral characteristics of the ULF turbulence was similar for the three comets, polarizations and wave forms varied. At comet GZ, where the turbulence was observed at least up to 7 ⋅ 10 5 km from the nucleus (∼7 times the bow shock stand-off distance of 10 5 km), wave polarization changed from essentially elliptical at large distances to nearly linear close to the bow shock, where also significant phase steepening of the waves was observed. The waves were identified as right-hand magnetosonic (fast MHD) waves, propagating roughly parallel to the magnetic field and in the sunward direction (in the SW frame), in good agreement with wave generation by the right-hand resonant helical beam instability (Tsurutani et al., 1987). In the case of phase-steepened waves near the bow shock, these were found to be lead by large-amplitude, parallel-propagating ion-scale whistler packets at frequencies around 0.3 Hz (∼30 times the water cyclotron frequency and about twice the proton cyclotron frequency), possibly resulting from generation of dispersive whistlers, pick-up of heavy ions and protons at the distorted steepened wave fronts, or trapping of heavy ions by the whistler wave train (Tsurutani, 1991).
At comet Halley, where the ULF turbulence was observed at least up to 2 ⋅ 10 6 km from the nucleus (∼2 times the bow shock stand-off distance of 10 6 km), several higher harmonics of the fundamental "pump" wave were additionally present in the spectra. The fundamental mode generally had linear polarization, whereas the higher harmonics typically exhibited elliptical polarizations. These waves were identified as Alfvénic type fluctuations, the wave mode (fast or slow MHD) apparently varying, possibly due to the plasma varying between > 1 and < 1, which would change the nature of the Alfvén wave between fast and slow modes (Glassmeier et al., 1989). Unlike at GZ, the waves had no obvious structure (Tsurutani et al., 1995), and there was no sign of phase steepening or development of leading whistler packets, not even close to the bow shock. Also these waves have been suggested to be generated by the right-hand resonant helical beam instability, which has been shown to exhibit linear polarization under certain conditions (Gary & Winske, 1986).
At comet GS, the ULF wave signatures were observed out to a distance of ∼6 ⋅ 10 5 km from the nucleus (∼30 times the observed bow shock/wave distance of ∼2 ⋅ 10 4 km) in the form of discontinuous wave packets, and more or less continuously inside of 2.6 ⋅ 10 5 km. This is in spite of the fact that Giotto's flyby occurred during quasi-perpendicular ( ≳ 50 • ) conditions, for which no such turbulence was observed at GZ or Halley. The At intermediate distances (3-18 ⋅ 10 4 km) from comet Halley, inside the magnetosheath, the VEGA spacecraft observed small-scale magnetic field depressions with a thickness of about one water ion gyro diameter (∼800 km) in conjunction with corresponding increases in ion density (i.e., out of phase) (Russell et al., 1991). These were linearly polarized fluctuations at an oblique angle to the background magnetic field, propagating at an angle of about 70 • . They were interpreted as slow magnetosonic waves or mirror mode waves generated by the mirror mode instability. These fluctuations disappeared closer to the nucleus where ion densities exceeded 1,000 cm −3 . Similar, though less prominent structures were also detected by Giotto and ICE, by the latter also at comet GZ (Tsurutani, 1991).
At GZ and Halley, there were also plenty of plasma waves detected at higher frequencies, in the ELF (10-1,500 Hz) and VLF (10 3 -10 6 Hz) frequency ranges (Scarf, 1989). (No instrument on Giotto had the ability to detect waves at these frequencies, so no such observations are available from GS.) At GZ, the short electric antenna of the ICE plasma wave instrument observed bursts of strong waves in the ion acoustic frequency range (0.6 ≲ ≲ 10 kHz) almost continuously within about 2⋅10 6 km of the nucleus. These bursts occurred preferentially under quasi-parallel conditions ( ≲ 60 • ) and have been attributed to a beam-type instability excited by the pick-up photoelectron population Richardson et al., 1989). The long electric antenna and search coil magnetometer detected electromagnetic waves at frequencies characteristic of the electron-scale whistler mode ( ≲ 100 Hz, corresponding to one fourth to one half of the electron plasma frequency) and near the hydrogen lower-hybrid (LH) frequency (6-12 Hz) . The latter have been proposed to be generated by an ion-loss cone instability due to the pick-up of water group cometary ions into a perpendicular ring distribution ) (thus, yet another instability driven by the pick-up ions, this time in the ELF range; Hartle et al., 1973). Just upstream of the bow shock, electron plasma oscillations were additionally detected, and the shock crossing and downstream region featured broadband impulsive turbulence (Scarf, 1989;Tsurutani, 1991). Electron-scale whistlers (∼300 Hz) were also prominent at Halley (out to at least 1.3 ⋅ 10 5 km), and waves around 1 kHz, likely ion acoustic waves, were also observed (Savin et al., 1987). Emissions near the hydrogen LH frequency (∼30 Hz) were detected closer to the comet (within 5-7 ⋅ 10 4 km) (Grard et al., 1985). Galeev (1987) suggested that the whistler waves could in fact be exited by supra-thermal electrons accelerated by the LH waves.

Plasma Waves Observed by Rosetta at 67P
The European Space Agency's Rosetta spacecraft accompanied the comet 67P/Churyumov-Gerasimenko (hereafter 67P) in its orbit around the Sun from August 2014 (at 3.6 au from the Sun) through perihelion in August 2015 (at 1.24 au) until the end of September 2016 (3.8 au). This provided the instruments in the Rosetta Plasma Consortium (RPC) (Carr et al., 2007) with an unprecedented long-term view of the near-nucleus cometary plasma environment of an intermediately active comet. During this time, the production rate of 67P varied from ∼4 ⋅ 10 25 to 3.5 ⋅ 10 28 s −1 Heritier, Henri, et al. 2017) (i.e., comparable to GZ) and the spacecraft generally stayed in close to terminator orbit within about 400 km of the nucleus, with the exception of a month-long sunward ("dayside") excursion out to ∼1,500 km in September-October 2015, and a tailward ("nightside") excursion out to ∼1,000 km for just over 2 weeks in March-April 2016.
The observations revealed a highly variable and dynamic plasma environment, and several different types of plasma waves have been observed at 67P. Richter et al. (2015Richter et al. ( , 2016 reported on low-frequency, large-amplitude ( B/B ∼ 1) compressional magnetic field oscillations at ∼20-50 mHz (a.k.a. "singing comet ODELSTAD ET AL. 3 of 26 waves") in the early and late low-activity phases of the mission, but disappearing during the high-activity phase between March 2015 and Spring 2016 (see also Breuillard et al., 2019, andGoetz et al., 2020, for further investigations into the properties and circumstances of these waves). While generally close to the local proton gyro-frequency, the variations in peak frequency of the waves did not correlate with observed variations in the ambient magnetic field magnitude, so it was argued that the waves were in fact not in proton cyclotron resonance. Unlike previous cometary encounters, Rosetta's prolonged stay at 67P was generally characterized by a gyro-radius of newborn cometary ions much larger than the scale size of the innermost interaction region, where the S/C spent almost all of its time. Thus, the ring-beam type pick-up distributions characteristic of previous encounters should not develop here, certainly not during the low-activity phases of the mission where the "singing comet waves" were detected, and were indeed not observed during this time (Behar, Lindkvist, et al. 2016;Berčič et al., 2018). Instead, the cometary plasma environment featured essentially unmagnetized cometary ions and magnetized electrons, a configuration resulting in an electric current perpendicular to the magnetic field. Meier et al. (2016), based on a linear homogeneous dispersion analysis assuming a cold, three-component plasma consisting of magnetized electrons, magnetized solar wind protons and a beam of unmagnetized cometary water ions, suggested a modified ion-Weibel instability driven by the cross-field current as a generation mechanism for the waves. However, it should be noted that neither the assumption of homogeneity, nor that of magnetized solar wind protons, are really supported by observations. For the latter, the absence of a clear ring distribution of solar wind protons is particularly pertinent .
During the high-activity phase close to perihelion, Volwerk et al. (2016) observed strong quasiperiodic dips in the magnetic field strength (relative peak-to-peak amplitude 2ΔB∕B ≳ 1), typically anti-correlated with variations in the plasma density and with minimum and maximum variance directions perpendicular and parallel to the background magnetic field, respectively. This led them to be identified as mirror mode waves which, as discussed above in the context of the previous comet encounters, can be generated by the unstable ring-beam type pick-up distributions. Volwerk et al. (2016) thus inferred such distributions to have developed, at least intermittently, in the heavily mass-loaded plasma and piled-up magnetic field in the inner coma close to perihelion. Observations of energy-angle dispersion of accelerated heavy ions by Nicolaou et al. (2017) indeed suggest that pick-up ion distributions are, at least sometimes, influenced by the effects of ion gyro-motion, although this gyro-motion would be more complex than for the classical ring or partial-ring distributions since the plasma here exhibits significant inhomogeneities on scales comparable to local ion gyro-radii.
One of the most significant findings of the plasma instruments onboard Rosetta was the diamagnetic cavity, a magnetic field-free region in the innermost part of the coma into which the interplanetary magnetic field cannot reach (Goetz, Koenders, Hansen, et al., 2016;Goetz, Koenders, Richter, et al., 2016). First predicted theoretically by Biermann et al. (1967), it was also observed by Giotto at Halley . The cavity at 67P was observed in the form of intermittent magnetic field drop outs ranging in duration from 8 s up to 40 min. The low velocity of Rosetta with respect to (w.r.t.) the comet (ôŔřA¸≲ 1 m/s) suggests that these highly transient events were the result of the cavity expanding and contracting over Rosetta's position, rather than resulting from the spacecraft moving into and out of a stationary cavity. Another possibility is blobs of unmagnetized plasma detaching from the main cavity structure and convecting past the spacecraft (Odelstad et al., 2018, and references therein). Further background and context on the diamagnetic cavity and the surrounding region will be given together with the observations presented in section 3. Gunell, Goetz, et al. (2017) and  reported on plasma density oscillations associated with ion acoustic waves at frequencies ∼200 Hz, both in the magnetized plasma in the early low-activity phase and in the unmagnetized plasma inside the diamagnetic cavity during the high-activity phase close to perihelion. The generation mechanism of these waves is still unknown, though a current-driven instability has been proposed, at least for the waves inside the cavity. Karlsson et al. (2017) and André et al. (2017) reported observations of electric field oscillations in the range of the local H 2 O + LH frequency from October and November 2015, close to peak activity of the comet, and attributed them to a LH drift instability caused by gradients associated with observed local density fluctuations. Madsen et al. (2018) found electrostatic waves of similar frequency also inside the diamagnetic cavity and suggested that they were ion acoustic waves resulting from oscillations of the cavity boundary at the LH frequency triggered by LH waves in the magnetized plasma outside the cavity.

RPC-LAP
The Rosetta Langmuir probe instrument (RPC-LAP) (Eriksson et al., 2007) (hereafter LAP) consists of two spherical Langmuir probes (LAP1 and LAP2) with radii of 2.5 cm and surface coating of titanium nitride (TiN), mounted on booms of 2.24 and 1.6 m lengths, respectively, protruding from the spacecraft main body. LAP has capability for three basic modes of operation: current measurements at fixed bias potential, potential measurements at fixed bias current (or with a floating probe, i.e., disconnected from the biasing circuitry) and Langmuir probe bias potential sweeps. In the first mode, which is the one used in this paper, the bias voltage is held at a constant value (with respect to the spacecraft, which is floating ground for the measurements) while the probe current is sampled continuously at sample rates ranging from 0.5 to 60 Hz, depending on the available telemetry rate. The bias voltage is typically about 30 V positive or negative, for sampling of plasma electrons or ions, respectively.
In the orbit motion limited (OML) regime, where the Debye length D is much larger than the radius of the probe, the current due to collection of ions by a spherical probe at a negative potential with respect to the ambient plasma can be related to the ambient plasma parameters by the following formula (Fahleson, 1967): where I i0 and i are the random thermal current and normalized potential, given by Here, a is the probe radius, V p the probe potential with respect to the ambient plasma and n i , q i , T i , m i , and u i are, respectively, the ion number density, charge, temperature (in kelvin), mass and bulk drift velocity. LAP uses the spacecraft as electrical ground, thus V p is related to the controlled bias potential U B as V p = U B + V S/C , where V S/C is the spacecraft potential. When the probe is at a fixed negative bias potential V p , the probe current is directly proportional to the ion number density n i , with proportionality constant depending on the temperature and drift velocity of the ions. Provided that changes in T i and u i are small compared to variations in density on relevant timescales, variations in probe current can be attributed to density fluctuations in the ambient plasma. However, variations in the spacecraft potential, which are to be expected if the density fluctuations are large since V S/C depends heavily on n i , can greatly affect such measurements since the probe bias potential w.r.t. the ambient plasma is then not fixed. For the ion current used here, the effect would be to amplify the probe current fluctuations since the associated density enhancements would drive the spacecraft potential more negative , increasing the effective probe bias voltage w.r.t. the ambient plasma. The amplitude of the probe current variations may therefore overestimate the magnitude of the inferred density fluctuations, but the sign and phase will be correct. Here, we will generally rely on the Mutual Impedance Probe (section 2.3) to gauge the magnitude of the density fluctuations.

RPC-MAG
The Rosetta fluxgate magnetometer experiment (RPC-MAG) (Glassmeier et al., 2007) (hereafter MAG) comprises two triaxial fluxgate magnetometer sensors, inboard (IB) and outboard (OB), mounted 15 cm apart 1.5 m out from the spacecraft main body on the same boom as LAP2. Three orthogonal components of the magnetic field are sampled at a resolution of 31 pT in a range of ±16 μT at a frequency of 20 Hz, although sometimes this is downsampled to 1 Hz onboard due to telemetry constraints. The magnetic field measurements are subject to disturbances from the spacecraft and the other instruments onboard; however, the most prominent of these (e.g., the influence from the reaction wheels) lie in the frequency band 2-10 Hz (Goetz, Koenders, Hansen, et al., 2016) and are above the frequency range of interest in this study. There is also an unknown offset depending on sensor temperature as well as spacecraft influences. The unmagnetized nature of the plasma inside the diamagnetic cavity allows for good calibration of the sensors during the cavity crossings, which can be combined with a temperature model to obtain high-quality measurements in the surrounding region. Such data was produced by Goetz, Koenders, Hansen, et al. (2016) for 1 Hz data from the OB sensor and this is used in this study to provide accurate magnitude and direction of the background magnetic field. However, the detailed wave analysis is based on full-resolution 20 Hz data (also from the OB sensor) that has been calibrated only with the temperature model, but this should not be a problem for the frequency range of interest here.

RPC-MIP
The Mutual Impedance Probe (RPC-MIP) (Trotignon et al., 2007) (hereafter MIP) consists of two pairs of dipole antennas, which can be used to obtain the plasma density from characteristic signatures that appear in the mutual impedance spectra at or near the plasma frequency. MIP is used here to provide accurate background values of the plasma density, since the LAP fixed-bias current cannot reliably be used for absolute density measurements (cf. section 2.1). However, the limited time resolution of MIP of ≳ 2.5 s is not suitable for detailed wave analysis in the frequency range of interest here, thus for this we rely on the LAP fixed-bias current measurements, as described above.

RPC-ICA
The Rosetta Ion Composition Analyzer (RPC-ICA) (Nilsson et al., 2007) (hereafter ICA) measures three-dimensional distribution functions of positive ions with an electrostatic analyzer (ESA) of a spherical top-hat configuration. It also has a magnetic momentum filter that resolves the major ion species such as protons, helium and water group ions. ICA nominally has an energy range of 5 eV/q to 40 keV/q with an energy resolution dE∕E = 0.07, effectively changing up to dE∕E = 0.30 for low-energy (<30 eV) ions (due to pre-acceleration into the ESA), a time resolution of 192 s and a field of view (FOV) of 90 • × 360 • . In response to the highly variable and dynamic cometary plasma environment encountered by Rosetta at 67P, a high time resolution mode was implemented and used intermittently throughout the mission . Here, the energy and FOV were reduced to 5-97 eV/q and 5 • × 360 • , respectively, allowing for an improved time resolution of 4 s.
The lower cut-off energy in ICA ion spectra is a good proxy for (the negative of) the spacecraft potential. ICA ion energies are subject to an uncertainty in the absolute level of the energy, which after corrections through comparison with LAP data has been reduced to a few eV , but should reliably capture variations in V S/C . This is the main purpose for its use in this paper, since LAP cannot do this at sufficiently high time resolution while in fixed-bias current mode.

The Primary Event on 20 November 2015
Figures 1a and 1b show data from a ∼2.5 h time period on 20 November 2015, at a cometocentric distance of about 130 km, during which the spacecraft entered the diamagnetic cavity a number of times (indicated by purple patches in Figures 1a and 1b). LAP1 probe current (black line, to be read off left-hand axis) and MIP plasma density (red line, to be read off right-hand axis) are shown in Panel a. The magnitude and the x, y, and z components of the MAG magnetic field are shown in Panel b (black, red, green, and blue lines, respectively). Here, the x, y, and z components are given in the Comet-centered Solar Equatorial (CSEQ) coordinate system, where the x axis points toward the Sun, the z axis is the component of the solar north pole that is orthogonal to the x axis, and the y axis completes the right-handed coordinate system. Also shown in Panel b is the polar angle between the magnetic field and the comet-S/C ("radial") direction (yellow line, to be read off the right-hand y axis). We note that we are in a region of substantial magnetic field pile-up, where the field is generally (close to) perpendicular to the radial direction. The region surrounding the cavity is characterized by plasma density (Figure 1a, right-hand y axis) and magnetic field ( Figure 1b) enhancements, as previously reported by Goetz, Koenders, Richter, et al. (2016) and Henri et al. (2017). These large ( n/n > 1 and B/B > 1) compressive features are highly asymmetric: The density and magnetic field increase much more rapidly up to its peak value than the rate at which they decrease afterward (see also Hajra et al., 2018). An investigation into the nature of these magnetized structures is beyond the scope of this study, but this should be further examined in future work. We focus here on the significant fluctuations that are intermittently observed in the LAP1 fixed-bias ion current (Figure 1a, left-hand y axis), predominantly on the descending slopes of these steepened magnetized structures, and which coincide with similarly large scatter in the MIP density measurements.
A zoom-in of the data in Panel a during the interval 05:55-06:20 (marked by a black dash-dotted rectangle in Figures 1a and 1b) is shown in Figure 1c, with corresponding magnetic field data in Figure 1d. Here, it can be seen that these fluctuations are in fact large-amplitude ( I/I ∼ 1), quasi-harmonic oscillations, with a frequency on the order of 0.1 Hz. The coincident scatter observed in MIP plasma density measurements is clearly also attributable to these oscillations and we can confirm that their relative amplitude in terms of density is n∕n ≳ 1. No associated signatures of comparable prominence can be discerned in the magnetic field data in Panel d. In order to better quantify the frequency of the oscillations and look for fainter signatures of them in the magnetic field, we compute wavelet scalograms over the frequency range 0.01-3 Hz 10.1029/2020JA028592 (for details, see section A2). The results, which represent the power spectral density at each frequency as a function of time, are shown in Figures 1e and 1f. The contours show selected harmonics of the water ion cyclotron frequency. The oscillations are clearly seen in the scalograms of the LAP1 current (Panel e) at about 05:57-06:07 and 06:15-06:18, and at frequencies between about four and twelve times the water ion cyclotron frequency, with peak power generally around the 6th to 8th harmonics. In this time-frequency representation, we do in fact observe clear signatures of the oscillations also in the magnetic field (Panel f, which shows the sum of the power in the three components). Panel g shows the result of calculating the coherence between the LAP1 current and magnetic field by means of the cross-wavelet transform (for details, see section A2). The coherence has been calculated in this way between the LAP1 current and each of the three magnetic field components; the result shown in Panel g is the weighted average of these individual coherences, with the weights being the relative power of the respective components. This quantity is independent of the specific coordinate system used to represent the magnetic field data and is a natural generalization of the coherence concept to the combination of scalar and vector-valued time series. The coherence in Panel g exhibits clear peaks at the times and frequencies of the density and magnetic field oscillations, thus showing that they are indeed related. No similar coordinate-independent generalization can be obtained for the relative phase of the LAP1 current and magnetic field; this will instead be addressed in the context of minimum variance analysis and wavelet-based spectral polarization analysis below.  Figures 1c and 1d), for minimum variance analysis (Sonnerup & Scheible, 1998) of the magnetic field measurements. We first apply a bandpass-filter between 0.05 and 0.25 Hz (forward-backward filtering with a third-order elliptic filter) and in each subinterval normalize by the maximum component standard deviation (STD). Figures 1h, 1k, and 1n show the resulting principal components of the magnetic field. B 1 (red lines) is the maximum variance component, corresponding to the presumed principal axis of polarization (Sonnerup & Scheible, 1998). B 2 (green lines) is the component with intermediate variance and B 3 (blue lines) is the minimum variance component, corresponding to the presumed direction of propagation. (The sign ambiguity of the principal components has been resolved by requiring B 1 and B 3 to have positive components along the background magnetic field, with B 2 completing the right-handed system). Also shown are similarly bandpass-filtered LAP1 current (black lines), also normalized by its STD over each subinterval. The oscillations are clearly observable in the first two principal components (B 1 and B 2 ). The phase difference between these is quite variable, sometimes 90 • out of phase, indicating elliptical polarization, sometimes in (anti-)phase, suggesting linear polarization, and sometimes somewhere in between. It is possible that the polarization axes of the wave change during the span of each subinterval, so that the minimum variance components, which are calculated for each subinterval in its entirety, are not accurate representations of the polarization axes throughout the interval. Another alternative is that the wave polarization really is that sporadic. We note that B 1 generally lags the (negative of the) LAP1 current (i.e., the density) fluctuations by ∼90 • in all three cases.

Minimum Variance Analysis
Figures 1i, 1l, and 1o show the background magnetic field (hereafter B 0 ) (the 1 Hz high-quality data of Goetz, Koenders, Hansen, et al., 2016, cf. section A3), its magnitude (black line) and components in a coordinate system whose basis consists of the principal component directions of the wave magnetic field; this in order to expose its direction w.r.t. these components. The component of the background field along the maximum variance (principal axis) direction clearly dominates over the other two, suggesting that B 0 is close to parallel to B 1 and close to perpendicular to B 3 (the wave vector, or propagation, direction ±k Sonnerup & Scheible, 1998). Figures 1j, 1m, and 1p display the polar angles between B 0 and each principal component of the wave magnetic field. For all three subintervals, B 3 is close to perpendicular to B 0 , indicating cross-field propagation. Also for all three subintervals, B 1 is almost parallel to B 0 (the angle between them is ≲20 • ). Also shown in Figures 1j, 1m, and 1p are the angles between the comet radial direction r, and B 0 , and B 3 , respectively. The former is consistently close to 90 • , as previously noted, but the latter varies substantially between the three subintervals, indicating that the direction of the wave normal vector varies in the plane perpendicular to B 0 .

Frequency Domain Polarization Analysis
The principal component variances from the minimum variance analysis above are shown below Figures 1j, 1m, and 1p, respectively, for the three subintervals. The ratio of maximum to minimum variance provides an indicator of how well defined the maximum and minimum variance directions are  (Sonnerup & Scheible, 1998). In the first case this is quite high (∼20), while in the other two cases this ratio is much lower (∼7-8). We are close to the limit where the wave parameters change on timescales comparable to the wave period, so we also perform a wavelet-based polarization analysis in the time-frequency domain. Our analysis is based on the approach of Santolík et al. (2003), applying a Singular Value Decomposition (SVD) to an enlarged real matrix obtained by taking the real part of the (complex) spectral matrix and vertically concatenating the imaginary part to it. The spectral matrix is obtained from wavelet rather than Fourier transforms of the magnetic field components (for details, see section A2), facilitating a dynamic time-frequency analysis of the polarization parameters. Figure 2 shows the results of the polarization analysis for the time interval shown in Figures 1c-1g.
Panel a shows the same magnetic field scalogram as Figure 1f to provide context. Dash-dotted rectangles (black or white, whichever give better contrast) have been inscribed to indicate the time-frequency domains corresponding to the minimum variance analyses in Figures 1h-1p. Additional intervals of wave activity, for which minimum variance analyses were not shown in Figure 1 for lack of space, have been similarly accentuated using ellipses. Panel b shows the ellipticity, that is, the ratio of the lengths of the two axes of the polarization ellipse (Santolík et al., 2003) (analogous to the ratio 2 / 1 in the minimum variance analysis above). Panel c-d show the cosines of the angles between the background magnetic field and the wave vector and the major axis of polarization, respectively (analogous to ∠B 0 B 3 and ∠B 0 B 1 , in Figures 1j,  1m, and 1p). Panel e shows the 3-D degree of polarization (DOP), that is, the ratio of the power carried by the polarized component to the total power in the field (Fowler et al., 1967;Samson, 1973). Panel f shows the planarity of polarization, which is based on the ratio of the length of the smallest and largest axes of the polarization ellipse such that it ranges from 0-1, where 0 corresponds to completely isotropic polarization without any preferred direction and 1 corresponds to an ideal plane wave (Santolík et al., 2003). Panel g shows the coherence between the LAP probe current and the component of the wave magnetic field along the principal axis of polarization (the latter is analogous to B 1 in the minimum variance analysis). Panel h shows the phase of this component w.r.t. the (negative of the) LAP1 probe current (i.e., the density). Negative values thus imply that the principal component of the wave magnetic field lags the density fluctuations.
Turning our attention first to Panels e and f, we note that the degree of polarization and planarity generally stand out as being quite large for the outlined sectors, with values typically ≳ 0.7. This suggests that the magnetic field fluctuations accompanying the density fluctuations are indeed quite consistent with plane coherent waves. In the corresponding sectors in Panels c and d we observe clearly that the wave vectors of these waves are perpendicular to the background magnetic field, with the major axis of polarization closely aligned with B 0 . Thus, the waves have associated variations in the magnitude of magnetic field, not just its direction. The ellipticity in Panel b comes out at ∼ ±0.5, indicating the presence of a second minor component of the wave field, perpendicular and phase-shifted 90 • w.r.t. the major component, thus in phase/antiphase with the density fluctuations, but most often clearly smaller than the principal component (i.e., elliptic rather than circular polarization). The sign of the ellipticity is predominantly positive, indicating right-handed polarization (about B 0 ). However, this is not entirely consistent between and within the highlighted sectors; in at least two cases, left-handed polarization is indicated, at least for some frequency components in the relevant range. We note that the definition of sense of polarization w.r.t. B 0 is somewhat ambiguous when B 0 lies in or close to the plane of the polarization ellipse, as is the case here. Thus, the determination of the sense of polarization of the waves can be expected to come with great uncertainty and difficulty of interpretation. Finally, the coherence and phase plots in Panels g and h show that the wave component along the principal axis of polarization generally exhibits strong coherence with the density fluctuations, with a consistent phase lag of ∼90 • . These results are generally in good agreement with the minimum variance analyses in Figures 1h-1p. We note that in a few cases, the polarization analysis is not entirely consistent within an individual highlighted sector (as was noted for the ellipticity above). For example, the first of the two ellipses in Panels c-g exhibit lower DOP, planarity and coherence near their lower edges, where also the wave vector direction and principal axis of polarization deviate from the general results. This could be interpreted as the sought-after waves not coming through at those lower frequencies in these cases due to higher levels of other magnetic field fluctuations at low frequencies. We also note that the DOP of the last rectangle in Panel e is conspicuously low, indicative of poor coherence between the different components of the magnetic field, in spite of the good coherence between the major component and the LAP1 probe current. Possibly, this can be attributed to elevated levels of incoherent noise in the other magnetic field components, as can perhaps be surmised from the bandpass-filtered time series in Figure 1n.

Adding More Events: The Precarious Use of LAP2
The events shown so far are from the period around 20 November 2015, during one of most prominent clusters of cavity crossings. Other prominent clusters of cavity crossings occurred around 25 and 30 November 2015. During this period, LAP1 was generally in unsuitable operational modes to observe these waves (not sampling ion current), but LAP2 was occasionally run in appropriate mode(s). However, LAP2 is mounted in a less favorable position on the spacecraft for sampling the cometary plasma, being more prone to wake effects and influence from the electrostatic field of the negatively charged spacecraft (Odelstad et al., 2018). It was also more prone to being being shadowed from the sun (by the spacecraft main body). Furthermore, it exhibited strange behavior during large parts of the mission, possibly attributable to surface contamination during the long (≳2 years) spacecraft hibernation en route to the comet. During this time LAP2 was consistently in shadow and hence presumably very cold and prone to condensation of materials of spacecraft origin (Schläppi et al., 2010). Thus, there is ample reason not to expect its electrical coupling to the ambient plasma to be the same as that of LAP1. Nevertheless, similar density oscillations to those presented above from LAP1 data in November 2015 are also present in LAP2 data from 30 July, when another prominent cluster of cavity crossings occurred (Goetz, Koenders, Hansen, et al. 2016;Odelstad et al., 2018). Here, the fluctuations were not generally as prominent or prevalent.
In order to determine if these general differences in character of the oscillations between these two different mission phases represent actual physical differences or just different measurement premises, we make use of the fact that LAP2 was run in a similar mode also on 20 November 2015, concurrent with LAP1, coincident with the period of observations presented above. Figure 3 shows observations of the prominent waves between 05:55-06:02 by four different instruments. Panel a shows the probe currents of LAP1 and 2 (black and blue lines, to be read off the left-hand y axis, green segments indicate data gap interpolation, see section A1), density measurements by MIP and spacecraft potential estimates by ICA (red and yellow lines, to be read off the correspondingly colored right-hand y axes). The latter has been identified here as the first (lowest) energy bin that has at least 2 counts and for which the following bin also has at least two counts. The cut-off energies thereby obtained are also shown in the ICA ion energy spectra in Panel b (white line). The resulting spacecraft potential estimate is subject to an unknown additive offset of a few volts Stenberg Wieser et al., 2017) but should accurately capture the variations of V S/C , albeit with a time resolution of 4 s and with somewhat limited precision due to the rather coarse energy resolution (cf. section 2.4). The density oscillations are clearly observed in all the instruments, including LAP2 most of the time. However, their signature is much weaker in LAP2 than in LAP1; the LAP2 current has been multiplied by a factor of 3 here to yield comparable values to LAP1. We note, however, that the phase of the waves in LAP1 and LAP2 are not always the same (e.g., around 06:00). The reason for this is not well understood at this time; possible explanations include effects of wave electric fields or variations in the spacecraft potential modulating the ion flux reaching LAP2; we may note that the S/C potential in Panel b also is oscillating. The phases of the other measurements are on the other hand in good agreement with each other, thus we conclude that the phase of the oscillations in the LAP1 current can be trusted, but that the phase (and magnitude) of the oscillations in the LAP2 current is not to be taken at face value. For times when suitable LAP1 data is not available, we may then use LAP2 as an indicator of the existence of the waves, but not to quantify their properties.  during this time. As suggested by the above discussion, and reinforced by the MIP density measurements in Figure 4, the weaker amplitudes are not physical; in fact MIP indicates remarkably large-amplitude oscillations during the first, most prominent, wave-train. Both the latter events also exhibit magnetic field signatures, albeit quite weak (Panel f). In at least the last case there is some elevated coherence in Panel g, whereas for the preceding event it is hard to distinguish any elevation in coherence over the broadband transient signatures due to variations in the background magnetic field. We note also that the angle between B 0 and the comet radial direction is no longer completely perpendicular: here the angle is close to 60 • .

Minimum Variance Analysis
We have again selected three ∼2 min intervals for minimum variance analysis (indicated by dashed-dotted rectangles in Figures 4c and 4d); the results are shown in Figures 4h-4p. As before, the principal component variances are shown below Figures 4j, 4m, and 4p, respectively, for the three subintervals. They all have ratios of maximum to minimum variance of about 7-8. The first interval corresponds to the first, most prominent event (Panels h-j). We observe wave vectors perpendicular to the background field; however, unlike before, the intermediate instead of the maximum variance direction appears to be (close to) aligned with B 0 . The maximum variance component here appears to be perpendicular to B 0 instead. However, there appears to be a change of direction of the background field during this time, with two components of the magnetic field crossing each other in the middle and two components approaching each other near the end of the interval. Thus, the maximum variance component obtained here is probably not a very good estimate of the principal axis of polarization of the wave. The second event (Panels k-m) deviates even more from expectations: now we even appear to no longer have perpendicular wave vectors. However, the weakness of the magnetic field signature, and the dubious coherence with the density fluctuations, leads us to believe that this event is not representative of the magnetic properties of the density waves. The third event on the other hand (Panels n-p) appears more reliable, with stronger magnetic signature and clearer coherence with the density fluctuations. Here things look more familiar: perpendicular wave vectors and B 1 (close to) aligned with B 0 . Figure 5 shows the results of polarization analysis for the time interval shown in Figures 4c-4g (with the events selected for minimum variance analysis indicated by rectangles as in Figure 2). In the first rectangle in Panel d, we see why establishing the direction of B 1 w.r.t. B 0 from the minimum variance analysis of the first interval in Figure 4 was difficult: There appears to be multiple components present here, with different maximum variance directions. We note in Panel g that the coherence is strongest in the middle left portion of the rectangle, which coincides with high-cosine values in Panel d. This suggests that the part of the magnetic field oscillations most closely related to the density fluctuations has a principal axis of polarization close to aligned with B 0 . The magnitude of the ellipticity is again ∼0.5, with ambiguous sense of polarization. Somewhat surprisingly, the degree and planarity of polarization (Panels e and f) appear to be quite low for this event, in spite of the substantial power in the field and coherence with the LAP2 current fluctuations. Perhaps this can also be explained by the presence of multiple modes at similar frequencies; at least for the planarity the largest values (∼0.6) seem to coincide with the region of the rectangle exhibiting the strongest coherence with the density fluctuations. The direction of the wave vector (Panel c) is clearly perpendicular to B 0 . The phase (Panel h) is w.r.t. the LAP2 current and is therefore not reliable for determining the actual phase difference between magnetic field and density.

Concluding Remarks
The kind of density oscillations presented here are commonly observed in the region surrounding the diamagnetic cavity, although not always as clearly and persistently as in Figure 1; this represents the most prominent case we have come across in the data so far. They generally appear to be associated with the asymmetric magnetized structures, occurring predominantly on their descending slopes. Clear magnetic field signatures are not a general feature of these oscillations; more often than not, magnetic field signatures are very weak or absent. We have not found any clear distinguishing circumstances between occurrences of the oscillations with and without clear magnetic signatures, although this should be further investigated in future works.

On the Nature of the Fluctuations
An obvious question is the physical nature of the observed narrowband fluctuations. To start with, we consider the highly correlated detection of the signal in three very different instruments (LAP, MAG, and MIP, to which we may also add the spacecraft potential signature from ICA) to make any interpretation in terms of spurious oscillations caused by spacecraft-plasma interactions highly unlikely. The preferred location of the waves, on the falling flank of the sawtooth-like pulses in plasma density, magnetic field and also ion flux (Odelstad et al., 2018;Stenberg Wieser et al., 2017) bordering the diamagnetic cavity is also a strong sign of a real physical phenomenon, so we will disregard any artificial spacecraft related mechanisms.

Characterizing the Plasma Environment
The plasma surrounding the diamagnetic cavity is characterized by total densities in the range 500-2,000 cm −3 , an electron temperature T e ∼ 5 eV (Odelstad et al., , 2018, at least intermittently interspersed with a cold electron population with T e ≲ 0.1 eV (Engelhardt et al., 2018;Eriksson et al., 2017), and magnetic field strengths |B 0 | ∼ 10-20 nT. The relative abundance of cold electrons, when present, is not known. We will here use the lower end of the aforementioned range of densities (500 cm −3 ) as a lower bound for the warm electrons, giving an electron plasma e = n e k B T e ∕(B 2 ∕2 0 ) ≳ 1. Attributing total densities in excess of 500 cm −3 to cold electrons gives relative abundances of the latter in the range 0-75%, roughly in line with preliminary estimates by (some) previous authors (Engelhardt et al., 2018;Gilet et al., 2017Gilet et al., , 2020Wattieaux et al., 2020). The electron gyro-radius is r Le ∼ 400-800 m (for T e ∼ 5 eV).
In terms of number density, the ions are dominated by H 2 O + (and H 3 O + ), predominantly of low energy (≲|eV S/C |) produced locally and not yet picked up by the solar wind electric field, but there are also pick-up ions at much higher energies likely ionized farther away and deflected back to the comet (Masunaga et al., 2019;Stenberg Wieser et al., 2017). (The total ion flux is actually dominated by this latter population.) Observations by ICA have shown that solar wind H + and He 2+ were entirely deflected away from the inner coma during the high-activity phase of the mission; such particles were consistently absent in ICA measurements 10.1029/2020JA028592 between 13 May and 11 December 2015  (except for very few brief events, in at least one case related to particular solar wind transient events Edberg et al., 2016). Assuming then an ion mass of 18 amu, the Alfvén velocity is v A ∼ 1-5 km/s, the ion plasma frequency pi ≳ 1 kHz, the LH frequency f LH ∼ 10-20 Hz and ion inertial length d i ∼ 20-40 km.
The ion velocity has been the subject of multiple studies (Odelstad et al., 2018;, indicating velocities likely ∼5 km/s (corresponding to energies ∼2 eV) in the region outside the diamagnetic cavity, although it is unclear to what extent this represents a bulk drift or thermal motion of the ions. Inside the diamagnetic cavity, the ions were inferred to flow radially outward from the nucleus with supersonic speed, at least w.r.t. the temperature in the direction perpendicular to the flow (Odelstad et al., 2018). Just outside the cavity the magnetic field is perpendicular to the radial direction, so one possibility is that the ion velocity, whether bulk or thermal, is largely in the direction perpendicular to the background magnetic field. This gives an ion gyro-radius r Li ∼ 50 km. If the flow of ions is diverted at the cavity boundary in such a way that a significant fraction of their energy goes into field-aligned motion, r Li could be somewhat lower, for example, ∼10 km at 0.1 eV (∼1 km/s). We note that Odelstad et al. (2018) did not constrain the ion temperature in the direction along the flow (i.e., the radial direction) inside the cavity. Qualitatively, the fact that the ions were found to be accelerated beyond the bulk flow speed of the neutral gas suggests that they were not completely collisionally coupled to the neutrals, but also subject to an accelerating electric field. Therefore, it does not appear likely that the ions would remain cold. Model calculations by  indeed suggest significant velocity spread in the radial direction, corresponding to an ion temperature T i ≳ 1 eV.
The (perpendicular) ion acoustic velocity c s = √ 2k B (T i + T e )∕m i is in the range 5-10 km/s, if T e is taken to be the temperature of the warm electron population. In the presence of an additional cold electron population this should be replaced by an "effective" temperature (Jones et al., 1975) T e,eff = n e,cold + n e,warm n e,warm ∕T e,warm + n e,cold ∕T e,cold . (4) Since we expect T e,warm ∕T e,cold ≳ 50, this requires n e,cold ∕n e,warm ≲ 1∕50 for the above ion acoustic velocity estimate to be accurate. The relative abundance of cold electrons is likely often much larger than that Gilet et al., 2017Gilet et al., , 2020, thus the presence of cold electrons may reduce the effective ion acoustic velocity, down to ∼1 km/s in the extreme scenario that the ions are also cold, but more moderately down to ∼3 km/s if the ions maintain thermal energies ∼1 eV.
The cometary plasma environment in general, and the region surrounding the diamagnetic cavity in particular, contains significant temporal and spatial inhomogeneities that may affect the existence and properties of many different plasma wave modes. The lack of multi-point measurements in the coma precludes unambiguously distinguishing between spatial and temporal variations. However, we may note from Panels c and d in Figures 1 and 4 that the wave observations do not generally coincide with the most rapid changes in background magnetic field or density in the spacecraft frame. In fact, during the precise time intervals corresponding to individual observed wave trains, their rates of change are much lower. Looking for example at the wave train between 06:05 and 06:07 on 20 November 2015, which has the steepest slope of the background magnetic field of the wave trains in Figure 1, ΔB B 0 ≈ 20−10 nT 15 nT ≈ 2∕3. This wave train contains roughly 12 wave periods, giving an apparent gradient length of roughly 18 wavelengths in the spacecraft frame. The apparent gradient length of the plasma density is typically even longer during the observed wave trains. While plasma drifts and a non-vanishing wave phase velocity may significantly complicate this picture, it is thus still of interest to first investigate if there are any linear homogeneous-plasma wave modes fitting the observations. The neutral density n n observed by the neutral gas spectrometer ROSINA-COPS (Balsiger et al., 2007) was around 5⋅10 7 cm −3 at the time of the observations (Odelstad et al., 2018). The momentum transfer cross-section in for ion-neutral collisions, assuming dominance of charge transfer H 2 O + + H 2 O interactions, is approximately given by in = 68∕ √ E ⋅ 10 −16 cm 2 , where E is the relative kinetic energy of the ions w.r.t. the neutrals. This gives an ion mean free path i = 1∕n n in of about 40 km for E = 2 eV and 10 km for E = 0.1 eV, that is, very close to the respective ion gyro-radii for the considered energies. Ion-neutral collisions may therefore play some role by for example decreasing the effective 10.1029/2020JA028592 magnetization of the ions. If there are several collisions between involved particles and neutrals during one wave period, this can seriously change the theory developed for homogenous and collisionless plasmas. We limit the scope of the following discussion to collisionless wave theory, while acknowledging that waves with wavelengths of a few ion gyro-radii or more may not be well described in this framework.

Doppler Shift
The waves reported in this study have typical frequencies in the spacecraft frame in between the water and proton gyro-frequencies ( f c ,H 2 O + ∼ 0.01 Hz and f c ,H + ∼ 0.2-0.4 Hz, respectively). However, these may be Doppler-shifted if the plasma has a drift velocity component along the direction of propagation, by a frequency shift Δ = k ⋅ v Di . Again, v Di is poorly constrained, especially its direction. The most generous assumption we dare to make here is |v Di | ≈ 5 km/s (cf. section 4.1.1), giving a maximum Doppler shift, obtained for v Di (anti-)parallel to k, of where + and -correspond to drift parallel and anti-parallel to k, respectively.
(1) can be ruled out since the electron gyro-frequency is in the kHz range; the spectral broadening ("Doppler spread") associated with such a large Doppler shift, orders of magnitude larger than the observed frequency in the S/C frame, would be enormous and unable to produce such narrow-band, almost sinusoidal signals as we observe in the ion cyclotron frequency range. Also, in (1) the response of ions to the waves is entirely negligible because of their large inertia, which is inconsistent with the large-amplitude fluctuations in ion density as observed, for example, by LAP. Thus, wave modes in this range are not further considered.
For waves in the LH range (2), with ≫ ci , the Doppler shift Δf would also have to be very large, many times larger than the the observed frequency in the S/C frame, to reduce the frequency down to the observed value. Furthermore, for a LH frequency ≳10 Hz, Equation 5 gives ≲ 500 m, which is close to the estimated electron gyro-radius (cf. section 4.1.1). LH waves are expected at wavenumbers kr Le ≲ 1 Graham et al., 2017Graham et al., , 2019Norgren, 2016;Norgren et al., 2012), corresponding to wavelengths ≳ 2 r Le ≈ 6r Le ≳ 2 km. Thus, they would appear to be out of Doppler shift range of the observed waves. It is of course possible that the herein assumed values of |v Di | and r Le are not entirely accurate; however, neither is likely off by more than a factor of 2, which would produce a maximum Doppler shift that falls just short of what would be required. In short, LH waves do not appear to be a likely candidate for the observed waves.
For waves in the MHD range (4), with ≪ ci in the plasma frame, the Doppler shift Δf would have to account for virtually the entire observed frequency of ∼0.1 Hz. From Equation 5 this requires a wavelength ≲ 50 km, which is close to the estimated ion gyro-radius. But the applicability of MHD is limited to ≫ r Li , thus the observed waves cannot fall into the MHD range.
It remains then, to look for possible wave modes in the ion cyclotron frequency range (3), that is, waves at the ion-kinetic scale. The waves in this range are typically classified into four different modes : (1 ′ ) the whistler mode (sometimes referred to as the ion-scale whistler mode to distinguish it from its namesake close to the electron gyro-frequency), which has typically been regarded as a kinetic extension of the fast magnetosonic mode (Gary, 1986), (2 ′ ) the kinetic Alfvén wave (also sometimes called the ion cyclotron wave for nearly parallel propagation), which is a kinetic extension of the MHD shear-Alfvén mode, (3 ′ ) the kinetic slow mode, similarly a kinetic extension of the slow magnetosonic mode, and (4 ′ ) ion Bernstein waves (IBW, or ion cyclotron harmonic waves (ICH) (André, 1985)), which appear as breakups of the whistler mode near the harmonics of the ion gyro-frequency Howes, 2009).
The kinetic slow mode (3 ′ ) is subject to strong ion Landau damping unless the ions are very cold (T i ≪ T e ) and represents an unphysical wave that does not exist in a weakly collisional or collisionless, finite ion temperature plasma (Howes, 2009;Krauss-Varban et al., 1994;Sahraoui et al., 2012). Thus, we will not consider it further here.

Whistler/Fast Magnetosonic Mode
The breakup of the whistler mode (1 ′ ) into ion Bernstein waves (4 ′ ) at large propagation angles to the magnetic field is due to finite Larmor orbit (FLR) effects (Swanson, 2012). Therefore, we can qualitatively expect this to occur when the wavelength approaches the order of the ion gyro-radius. Hence, the existence of (1 ′ ) as a continuous "whistler" mode for frequencies above the ion gyro-frequency at quasi-perpendicular propagation is contingent on it reaching these frequencies at sufficiently low wavenumbers and large wavelengths, that is, before FLR effects develop. Quantitatively, André (1985), Li et al. (2001), Howes (2009), andSahraoui et al. (2012) all indicate that this happens for k ⟂ r Li ≳ 1, which corresponds to ≲ 2 r Li . We can find rough constraints on the plasma parameters for (1 ′ ) to reach > ci before FLR effects kick in using the long-wavelength dispersion relation for this mode (Gary, 1986) ≈ k where is the angle between the wave vector k and the background magnetic field. The (perpendicular) ion acoustic velocity may be expressed as c 2 are the ion and electron (perpendicular) thermal velocities, respectively. Requiring > ci and ≫ 2 r Li then gives (in the quasi-perpendicular limit If r Li is taken to be the root mean square gyro-radius of the ion population, the right-hand side of Equation 7 is simply v th,i , giving where i is the ion plasma , that is, the ratio of ion thermal to magnetic pressure. Thus, (1 ′ ) exists as a continuous "whistler" mode for frequencies above the ion gyro-frequency (at quasi-perpendicular propagation) only in the limits of low i and/or small electron to ion temperature ratio. For the plasma at hand i ≳ 1, so we would require T e /T i ≫ 1 for Equation 8 to be fulfilled. As discussed in section 4.1.1, we expect T i ≳ 1 eV and T e ∼ 5 eV, giving T e ∕ T i ≲ 5. However, in the presence of both a warm and a cold electron population, T e in Equation 8 refers to the effective electron temperature T e,eff of Equation 4, which further decreases the temperature ratio, giving T e ∕T i ≲ 1 for n e,cold ∕n e,warm ≳ 0.1. Considering that the MIP detection threshold for cold electrons is n e,cold ∕n e,warm ≳ 1.5 (Gilet et al., 2020) (and likely even higher for LAP) and that cold electrons are observed so frequently (≳50% of the time) by these instruments in the region surrounding the diamagnetic cavity on the days of the wave observations Odelstad et al., 2018), it is tempting to infer the presence of some fraction (≳0.1) of cold electrons also during instances when they are not directly observed. This would then preclude the existence of a continuous "whistler" mode for frequencies above the ion gyro-frequency at quasi-perpendicular propagation. We therefore provisionally disregard (1 ′ ) as a possible wave mode here, in favor of the ion Bernstein waves (4 ′ ) that we expect to take its place in this plasma regime. Obviously, conclusively ruling out ion-scale whistler waves would require a more detailed investigation including numerical solutions of the dispersion relation (e.g., using WHAMP Rönnmark, 1982or PDRK Xie & Xiao, 2016, and perhaps also a more detailed analysis of the electron temperature data obtained by LAP and MIP during the specific events where the waves are observed. However, this is beyond the scope of this paper.

Kinetic Alfvén Wave
The kinetic Alfvén wave (2 ′ ) was investigated by Sahraoui et al. (2012) for large propagation angles to the background magnetic field in high-plasmas with warm ions, that is, fairly similar to what was observed outside the diamagnetic cavity of 67P. They found that for very oblique propagation, ≳ 89.9 • , (2 ′ ) extends all the way down to the electron gyro-scale with < ci and relatively weak damping. (For somewhat less oblique angles, it becomes dispersive at kr Li ≳ 1 and obtains frequencies larger than ci , but the damping also becomes more important.) Thus, we cannot rule out (2 ′ ) based on (non-)existence in relevant parts of the frequency-wavelength domain alone. However, we note that while (2 ′ ) does develop some degree of magnetic compression at smaller scales in this plasma regime, the compressive component doesn't become dominant at any scales, in contrast to what has been found for the waves we study here. It should also be noted that Sahraoui et al. found (2 ′ ) to be elliptically right-hand polarized, which is also what our observations suggest, but we still disregard this wave mode from further consideration here on account of its lacking magnetic compression.

Ion Bernstein Waves
Ion Bernstein waves (4 ′ ) typically refer to electrostatic waves propagating at (or near) right angles to the background magnetic field at (or near) harmonics of the ion cyclotron frequency. This places them squarely in the frequency range of the waves observed outside the diamagnetic cavity of 67P. In the absence of significant Doppler shift, the observed waves appear around the 8th to 10th harmonics of the cyclotron frequency. However, for v Di ∼ 5 km/s as conjectured above, the first harmonic could be Doppler-shifted up to the observed frequencies for wavelengths ≲ 50 km ∼r Li .
The dispersive properties of ion Bernstein waves change depending on the precise angle of propagation. If the parallel phase velocity of the wave is much smaller than the thermal velocity of electrons along B 0 , /k || ≪ v th,e,|| , electrons can flow rapidly enough along the magnetic field lines to cancel charge separations and "neutralize' the wave, effectively being in Boltzmann equilibrium with the wave potential (Chen, 1984;Schmitt, 1973). In this regime the waves are referred to as neutralized ion Bernstein waves (NIBW). The opposite case, /k || ≫ v th,e,|| , is called pure ion Bernstein waves (PIBW). NIBW and PIBW are separated by an intermediate regime, /k || ∼ v th,e,|| , where electron Landau damping is important (Schmitt, 1973). In terms of angle of propagation between k and B 0 , for PIBW this corresponds to cos ≫ ∕kv th,e,|| , which for waves at ∼ ci and ∼ r Li gives ≫ arccos Conversely, NIBW requires ≪ 89.5 • . We make no attempt to quantify the uncertainty in the angle of propagation of the observed waves as derived from the minimum variance and polarization analyses in previous sections, but we presume that it is not nearly accurate enough to use for distinguishing between NIBW and PIBW.
The electrostatic nature of the waves only prevails for sufficiently low i and/or short wavelengths, being subject to the necessary condition m i m e i (kr Li ) 2 ≪ 1 (Callen & Guest, 1973). The nature of the lowest-order electromagnetic (EM) modifications of the waves is different for PIBW and NIBW. For PIBW, the electromagnetic modification is due to the current driven by the difference in E × B drift of the electrons and ions, caused by FLR effects, in the wave electric and background magnetic fields. Since the wave electric field E is longitudinal (i.e., parallel to k) in the electrostatic limit, this current will flow in the direction of k × B 0 and, through Ampère's law, give rise to a magnetic field perturbation B || parallel to the background magnetic field (Norgren et al., 2012). For NIBW, the electrons are free to move along the magnetic field lines and the resulting current gives rise to a magnetic field perturbation in the direction perpendicular to both k and B 0 . Thus, we note that the compressive nature of the magnetic field fluctuations associated with the observed waves would suggest PIBW as the more likely wave mode.
However, these are only the lowest-order EM modifications, still requiring |k × E| ≪ k ⋅ E. Modes having appreciable k × E (≳10 −2 k ⋅ E) have been dubbed generalized ion Bernstein waves by previous authors (Dougherty, 1975;Fredricks, 1968;Puri et al., 1973). Here, numerical solutions of the full kinetic EM dispersion relation is required. Referring again to the work of Sahraoui et al. (2012), which provides such solutions for the plasma parameters most comparable to ours that we have found in the surveyed literature, in a high-plasma the (generalized) IBWs exhibit dominance of parallel over perpendicular power for propagation angles ≳ 80 • , that is, well into the expected angular range of NIBW. Indeed, appreciable ellipticities ∼0.5, as generally observed in the polarization analyses of previous sections, are only obtained for propagation angles well into the NIBW range. Hojo et al. (1993) used the linearized Vlasov equation to obtain an expression for density fluctuations associated with electromagnetic waves in the ion cyclotron range of frequencies, that is applicable to IBW. They found that in the limit /k || ≫ v th,e,|| , that is, corresponding to PIBW, the ratio of density to magnetic field fluctuations | n/n|/ | B || /B 0 | approaches unity, and that the phase difference Δ = n − B approaches 0. In the opposite limit, that is, corresponding to NIBW, Δ → −90 • but the magnitude of the ratio tends to 0, that is, the incompressible limit. Maximum compressibility was attained for /k || v th,e,|| ≈ 1.5, giving | n/n|/ | B || /B 0 | ≈ 1.3 and Δ ≈ −20 • . Thus, the large values | n∕n|∕| B || ∕B 0 | ≳ 10 of the observed waves were not attainable in this model, and the ∼ −90 • phase difference was only achieved in the incompressible limit. Most likely, inhomogeneities, non-linearities, and/or unstable, non-Maxwellian velocity space distribution functions must be considered in order to properly explain the observed wave properties; however, such aspects are largely beyond the scope of this paper. We will, however, briefly mention a few instabilities which can lead to growth of IBWs and that could perhaps be operational near the diamagnetic cavity of 67P.
IBWs can be generated by the drift-cyclotron instability. The highest excited harmonic of the ion cyclotron frequency is l m ≈ 3r Li /L n , where L n is the natural scale length of the plasma (Treumann & Baumjohann, 1997). The wave frequencies here being in the range 4-12 times the H 2 O + cyclotron frequency, we would require L n ≲ r Li ∕4, if there is no appreciable Doppler shift. In light of the large r Li and the highly inhomogeneous nature of the plasma surrounding the diamagnetic cavity, this does not seem out of reach. However, it should be noted that the observed waves do not appear on the steepest gradients of the asymmetric plasma and magnetic field enhancements that dominate the region near the diamagnetic cavity, but usually some way down their descending slopes, where the rate of change in the spacecraft frame is lower. Unless there is substantial variation in plasma bulk velocity across these enhancements, so that their flatter "tails" are convected past the spacecraft faster than their steeper "heads," it would seem unlikely that the wave growth is triggered by spatial inhomogeneities.
IBWs can also be generated by various velocity space anisotropies, for example, when the distribution of perpendicular energy has a "hump" in it (Hall et al., 1965;Rosenbluth & Post, 1965). This includes ring, ring-beam and spherical shell distributions, which have been extensively investigated in existing literature (Crawford, 1968(Crawford, , 1965Harris, 1959;Kumar & Tripathi, 2012;Noreen et al., 2019) and have been invoked to explain ion cyclotron harmonic wave generation in many space plasmas (Chen, 2002;Joyce et al., 2012;McClements & Dendy, 1993;McClements et al., 1994, and references therein). Broadly speaking, it has been found that positive growth rate from such distributions requires the ion plasma frequency to gyro-frequency ratio pi / ci to be sufficiently large. Quantitive estimates of this threshold are scarce for IBWs, but similar distributions of electrons have been shown to be unstable to electron Bernstein waves for 2 pe ∕ 2 ce ≳ 6 (Tataronis & Crawford, 1970a& Crawford, , 1970b. This can likely be extended, at least qualitatively, to positive ions as well (Crawford, 1965). pi / ci being ∼10 5 in the plasma we consider here, this threshold should easily be met. (More detailed analytical and/or numerical calculations for the specific plasma conditions encountered by Rosetta outside the diamagnetic cavity would clearly be of interest here but is deferred to future studies.) These kinds of distributions are of particular interest here due to the proximity of the diamagnetic cavity boundary. Inside the cavity, the ions were inferred to flow radially outward from the nucleus with supersonic speed, at least w.r.t. the temperature in the direction perpendicular to the flow (Odelstad et al., 2018). The magnetic field is perpendicular to the radial direction outside the cavity, thus it is possible that these radially outflowing cometary ions form such distributions in the magnetic field outside the cavity. This might explain why the waves are not observed on the steepest gradients; if these are surmised to occur primarily very close to the cavity boundary, well within an ion gyro-radius of the boundary, the outflowing ions may not have formed ring-or shell-like distributions yet, since this would occur on spatial scales comparable to or greater than the ion gyro-radius. We remind the reader that similar distributions have previously been inferred at 67P by other authors (Volwerk et al., 2016). This is currently our preferred suggestion as to the generation mechanism of the waves. However, since the estimated ion mean free path here is on the order of the ion gyro-radius (cf. section 4.1.1), the possibility of ion-neutral collisions impeding the formation of such distributions may need to be considered, but will not be further addressed here.
Finally, we note that IBWs can produce heating of both electrons, through electron Landau damping, and ions, through ion cyclotron damping. Electrons can also be accelerated along the magnetic field lines by the non-vanishing parallel component of the wave electric field that develops when k is not exactly perpendicular to B 0 , and/or electromagnetic effects become non-negligible (e.g., at high plasma ). If generated by a velocity space anisotropy, the general effect of the waves will be to reduce the free energy associated with that anisotropy, so that these waves may contribute to pitch angle scattering of ring or ring-beam distributions, or thermalization of shell-like distributions. They may thus play an important role in redistributing energy between different particle populations and in the reconfiguration of the plasma environment that occurs as a result of the transition from unmagnetized plasma inside the diamagnetic cavity and (at least partially) magnetized plasma in the magnetic pile-up region outside.

Water-Proton Ion-Ion Hybrid Waves
One last possible wave mode we address is water-proton ion-ion hybrid waves, the frequency of which (the Buchsbaum frequency; André, 1985;Buchsbaum, 1960) lies between the cyclotron frequencies of water ions and protons and can be tuned to fit the observations by invoking a suitable relative abundance of protons in the plasma. The Buchsbaum frequency is given by 10.1029/2020JA028592 where f cp , m p and p are the proton cyclotron frequency, mass and relative abundance, respectively, and likewise for the water ions. For p = 0 or 1, this reduces to the proton or water gyro-frequencies, respectively. Solving for p yields giving p ≳ 25% for , which is about where we observe peak power, if there is no appreciable Doppler shift. The absence of ions of solar wind origin means that any protons present at the time of the wave observations in this study would have to be locally produced in the inner coma. ICA data shows no sign of such low-energy protons. However, such low-energy protons would be hard to observe for the ICA and ROSINA-DFMS (Balsiger et al., 2007) instruments, and would not be possible to separate from the low-energy water group ions in RPC-IES . The Alice far-ultraviolet spectrograph on Rosetta observed Lyman-alpha emissions attributed to electron impact dissociation of H 2 O by photoelectrons resulting from photoionization of H 2 O by solar EUV radiation (Feldman et al., 2015). In this process, the impact leads to an excited H 2 O molecule that dissociates, releasing a neutral hydrogen atom that may still be excited and then decays to the ground state. This indicates there are at least hydrogen atoms in the coma, which can be converted to H + by photoionization. However, based on models by Vigren and Galand (2013), Vigren et al. (2015), and Heritier, Altwegg, et al. (2017), ≳ 25% appears out of reach. In absence of chemical loss, ≲ 10% might be possible, but models including further chemical reactions predict < 1% and close to the nucleus < 0.1%. It is possible to form H + by collisions between H 2 O + (or H 3 O + ) with H 2 O if the ion has an energy of ≳ 13 eV, but cross-sections for such process are low (Lishawa et al., 1990). As no observations or models yet indicate significant abundance of protons generated in this way, we do not consider this mechanism further.
Finally, we note that the inhomogeneous plasma encountered outside the diamagnetic cavity may support additional or modified wave modes not covered by homogeneous plasma theory. This issue should be more thoroughly addressed in future works.

On a Possible Relationship With the "Singing Comet Waves"
The "singing comet waves" (hereafter SCWs) observed by Richter et al. (2015) have frequencies typically in the range 20-50 mHz, which is only slightly lower than those observed here (∼0.1 Hz). Thus, the question arises whether the herein observed waves are the same, or a similar, phenomenon. Richter et al. (2016) found that the SCWs had wave vectors roughly perpendicular to the background magnetic field, in agreement with our observations. However, the SWCs had large magnetic field amplitudes, B/B ∼ 1, about 3-6 nT in absolute terms, whereas the magnetic field signatures we observe rarely even reach 1 nT in amplitude, corresponding to B∕B ≲ 0.1 in the much stronger background magnetic field during our observations. The SCW wave activity was almost continuous , whereas the oscillations we observe are intermittent, only occurring specifically on the descending slopes of the plasma and magnetic field enhancements in the region surrounding the diamagnetic cavity, and the related magnetic field oscillations are weak and sporadic. Also, the suggested generation mechanism for the SCWs (Meier et al., 2016) does not work in this region due to the absence of solar wind protons. We note that the physical environments differ quite a lot between the SCW observations and the ones presented in this paper, on account of the different comet activity levels. Therefore, one should not precipitate toward invoking the same explanation in both cases just because the frequencies are similar. Thus, we conclude that the herein presented wave observations constitute the detection of a new type of plasma waves at the comet.

Summary and Conclusions
We have investigated strong low-frequency density fluctuations observed in the region surrounding the diamagnetic cavity of comet 67P. We find that they are large-amplitude ( n/n ∼ 1) quasi-harmonic oscillations at frequencies on the order of 0.1 Hz, about 10 times the water ion cyclotron frequency and less than half the proton cyclotron frequency. They occur predominantly on the descending slopes of asymmetric plasma 10.1029/2020JA028592 and magnetic field enhancements that are characteristic of the region surrounding the diamagnetic cavity (Goetz, Koenders, Hansen, et al. 2016;Hajra et al., 2018;Henri et al., 2017). We have detected weak ( B∕B ≲ 0.1) signatures of them in the magnetic field as well. The magnetic fluctuations are generally elliptically polarized and with wave vectors perpendicular to the background magnetic field. Their principal axis of polarization is closely aligned with the background field and lags the density fluctuations by about 90 • . We have considered several possible wave modes, limiting ourselves here to linear homogeneous collisionless plasma wave theory, but discarded most of them for various reasons. We finally land on ion Bernstein waves (IBWs) as the most likely candidate. This constitutes a new type of wave not previously observed at comets (previous observations of IBWs in space plasmas are largely limited to the magnetic reconnection outflow region Narita et al., 2016 and fluctuations in the solar wind Perschke et al., 2013Perschke et al., , 2014. At 67P, they are possibly generated by unstable ring, ring-beam or spherical shell distributions of the cometary ions in the plasma near the diamagnetic cavity. Such waves may heat both electrons and ions, accelerate electrons along the magnetic field lines, reduce anisotropies and gradients in the plasma and reshape the plasma environment of the comet.

A1. Resampling LAP Fixed-Bias Current to MAG Sample Times
For the wave analysis in this study, we make use of LAP fixed-bias current measurements sampled at 57.8 Hz. LAP operational modes are organized in 32 s long sequences with a brief gap (∼1 s) at the end of each sequence. Also, every fifth 32-s sequence starts with a Langmuir probe bias voltage sweep, increasing the length of the gap to about 6 s. Thus, the fixed-bias current data obtained is not uniformly sampled over time periods longer than about half a minute. To prepare the data for spectral (wavelet) analysis, we downsample these measurements to the MAG 20 Hz sample times (forward-backward filtering with a third-order elliptic anti-aliasing filter with cut-off frequency at 8 Hz, then using linear interpolation), then filling the LAP data gaps by interpolation from forward and reverse autoregressive fits of the surrounding samples. The maximum length of prediction sequences is set to 300 points (at 20 Hz sampling rate) and an autoregressive model order is selected that minimizes the Akaike information criterion (fillgaps MATLAB, 2018). Examples of the results will be shown in Figure 3.

A2. Wavelet Analysis
Wavelet scalograms of power spectral density of LAP current and magnetic field presented in section 3 were computed using the processed version of the LAP fixed-bias current described in section A1 and the original 20 Hz magnetic field data, using a filter-bank implementation of the continuous wavelet transform (cwt MATLAB, 2018). Here, we have used the standard analytic Morlet wavelet (Eriksson, 1998;Morlet et al., 1982) with main central frequency 0 = 6 rad/s, 10 voices per octave and  2 normalization. However, the subsequent cross-wavelet analysis yielding coherence and polarization parameters is based on wavelet cross-spectra that have been computed using a Fourier transform based algorithm (wcoherence MATLAB, 2018), on LAP current and magnetic field data that have been further downsampled to 6.25 Hz, to reduce computational cost, using a polyphase antialiasing filter (resample MATLAB, 2018). Here, we have used 12 voices per octave (eight octaves between ∼0.01 and ∼3 Hz) and smoothed over eight wave periods (with a gaussian window) and 6 scales (rectangular window) (Grinsted et al., 2004).

A3. Obtaining the Background Magnetic Field
For the background magnetic field, we use the high-quality 1 Hz data produced by Goetz, Koenders, Hansen, et al. (2016). For obtaining relevant directions in the wavelet-based polarization analysis of section 3, this has been lowpass-filtered at 0.02 Hz, using a third-order elliptic filter applied in the forward and backward directions, and then upsampled to the aforementioned 6.25 Hz sample rate by linear interpolation.

Data Availability Statement
The data used in this paper have been delivered to the ESA Planetary Science Archive (http://archives.esac. esa.int/psa) where it will be made available shortly (MAG data are already available there). In the meantime, the LAP, MIP, and ICA data used in this paper can be found in the Zenodo repository (Odelstad et al., 2020).