Modification of Velocity Power Spectra by Thermal Plasma Instrumentation

The upcoming Solar Probe Plus mission (Launch 2018) will launch with the newest and fastest space plasma instrumentation to date. The Solar Wind Electrons, Alphas, and Protons (SWEAP) instrument suite, which measures thermal plasma, will make measurements faster than the local gyro-frequency and proton plasma frequency. By developing an end-to-end computer model of a SWEAP instrument, this work explores the specific instrumental effects of thermal space plasma measurement, particularly in the reproduction of velocity power spectra, or Power Spectral Densities (PSDs). This model reproduces the slowest measurement cadence of the Solar Probe Cup (SPC), a Faraday cup (FC) style instrument on that will measure thermal plasma density, velocity, and temperature on SPP. By using the calibrated model to model measurement of fully determined and synthetic turbulent time series data, a consistent underestimation of the velocity power spectral indices has been quantified, with possible implications for previous missions flying similar instrumentation.


Introduction
The spectral index is important in space physics because it describes the distribution of energy as a function of frequency. Typically magnetic field data are used to generate PSDs because magnetic field measurements are made at higher cadence, but in order to determine how Alfvénic a plasma is, comparison between the velocity PSDs and the magnetic PSDs is necessary. The amplitude of the PSD yields the overall power as a function of frequency, and the spectral index describes the relative distribution of power.
Power Spectra are often generated from magnetic field measurements, where the instrumentation often permits higher cadence measurements (Wüest et al., 2007). Power Spectra from velocity measurements often reveal a lower power spectral index of P (f ) ∝ f −3/2 , whereas magnetic field measurements usually reveal a P (f ) ∝ f −5/3 (Alexandrova et al., 2013). This disparity between spectral indices indicates an unequal distribution of energy between magnetic field and velocity fluctuations δV and δB, violating usual assumptions of Alfvénicity and equipartition.
This work explores the influence that instrumental effects may have upon disagreement in spectral indices between velocity power spectra and magnetic power spectra. The upcoming Solar Probe Plus mission has instrumentation that will produce the fastest thermal plasma measurements to date. By using the calibrated model to model measurement of fully determined and synthetic turbulent time series data, a consistent underestimation of the velocity power spectral indices has been quantified, with possible implications for previous missions flying similar instrumentation. This paper is split into six parts: Section one is this introduction. Section two is a primer on the traditional spectral analysis technique used in this paper. Section three describes the method of creating synthetic turbulent time series data. Section four includes details of the Solar Probe Cup and the computer model. Section five includes the results from feeding the SPC model the turbulent time series, and section six is a discussion of the results.

Spectral Analysis Background
The Power Spectrum, or more appropriately, the Power Spectral Density (PSD), is defined as the Fourier transform of the autocorrelation function R uu , where the autocorrelation function of a general function u(t) is, which is the time-averaged autocorrelation value (Zank, 2014, Blackman andTukey, 1958). Although u(t) in Equation 2 is a continuous function in nature, it is discretized by instrumentation to a finite time series of point measurements. Thus, Equations 1 and 2 must be discretized. The autocorrelation function is rewritten for the case of discrete data with time Δt between measurements, where n is the number of velocity measurements, h is some integer greater than zero and Δτ = hΔt, r = 0, 1, 2, ...m, where mh < n and describes the longest lag mΔt between measurements (Blackman and Tukey, 1958). The PSD is now estimated via the real-valued Fourier transform of R uu . The domain of the PSD in frequency space is determined by two values: the total time duration of the data and the time between measurements Δt. The Nyquist, or folding, frequency determines the minimum detectable frequency, is directly related to the time-resolution of the data record, where Δt is the sampling rate. For example, to detect a signal with a 2000Hz frequency, the sampling rate must be at least 1000Hz. Nyquist frequencies for different levels of SPC data are shown in Table 2 below.  Level 0 is the least sophisticated data level, and requires Fourier processing before it can be interpreted physically. Level 1 data contains data for distinct, discrete sections of velocity space, and is in instrument units of current rather than density, velocity, or temperature. Level 1 data can be used to determine plasma flow angles only. Level 2 data is in physical units of density, velocity, and temperature, and is the product of fitting Level 1 data to the Faraday Cup instrument response function. See Section 4 for details of instrument operation.

Generation of Turbulent Velocity Fluctuations
To examine turbulent fluctuations in SPC measurement, synthetic time-series velocity measurements are generated. Here the estimated PSD, or autocorrelation tensor, is defined in terms of velocity fluctuation values, Several simplifying assumptions are made in deriving an expression for δV as a function of solar distance. First, assume fast solar wind conditions. The idealized threepart turbulent spectrum consisting of the energy-containing, inertial, and dissipation ranges is most often and consistently associated with fast solar wind (Bruno and Carbone, 2005). Second, assume equipartition between magnetic fluctuations δB and velocity fluctuations δV . Assuming equipartition implies a purely Alfvénic plasma, a condition which is expected to become more frequent with increasing solar proximity as seen from Helios data Carbone, 2005, Bruno et al., 1985). Also assume that the magnetic field vector is well described by the Parker spiral configuration, thus the background magnetic field has a radial component pointing anti-sunward. Simplify coordinates by assuming the velocity fluctuations δV are antiparallel to the magnetic field fluctuations δB, as is often seen in spacecraft observations (Horbury et al., 2005). Finally, assume that the spacecraft is in the ecliptic plane.
The absolute direction of the velocity and magnetic field fluctuation vectors depends on whether the turbulence is slab or quasi-2D. In the case of slab turbulence, the magnetic field variation vector δB is parallel or antiparallel to the background magnetic field B 0 , and the sense of δV is antiparallel to δB. For two-dimensional turbulence, the vector δB is in the plane perpendicular to the background magnetic field B 0 , with the sense of δV again antiparallel to δB.
First assume a general expression for the spectrum energy contained in velocity fluctuations of form, Now examine the form of the estimated PSD. The estimated PSD can be described in three regions: the energy-containing range, the inertial range, and the dissipation range ( Figure 8). The transition between the energy-containing range is called the 'bendover scale' or 'bendover length', whereas the transition between the inertial range and the dissipation range is called the 'breakpoint frequency'. In this chapter's work the energycontaining range is neglected because the primary interest is in reproducing the inertial and dissipation ranges. Thus a range of frequencies significantly below the bendover scale are chosen. Typically the breakpoint frequency is on the order of either the ion cyclotron frequency, where m p is the proton mass, or the ion-inertial length, where n p is the proton mass density and ω p,i is the plasma ion frequency. The quantities B and n p can be written as functions of radius. Here, where the value n p,0 = 5.6 × 10 7 m −3 is the number density at 0.29AU, given in Adhikari et al. (2015) as calculated from Chandran and Hollweg (2009). Likewise the radial and azimuthal magnetic field magnitude are expressed as functions of radius,  Containing Range usually has k −1 form, and the Dissipation range is thought to be of ∼ k −3 form. The Inertial Range above is shown with Kolmogorov k −5/3 form.
where B R,0 is the magnetic field at r = r 0 , and θ is the colatitude angle (Owens and Forsyth, 2013). Making the simplifying assumption that all measurements occur in the ecliptic plane reduces sin θ = 1, For completeness, assume a purely radial flow. Thus the θ component of the magnetic field, With the solar wind density and magnetic field defined as functions of radius, the remaining quantity V r is left as a free parameter as long its value is within the range of acceptable solar wind radial speeds.
From Equation 7 the ion-cyclotron frequency is shown to be inversely proportional to the square of radial distance, and from Equations 8 and Equation 9 that the ion inertial length is proportional to radial distance, Taking the inverse of the ion inertial length and multiplying by the phase velocity yields the equivalent frequency for use in the following analysis, Thus, both the ion cyclotron frequency and the equivalent frequency for the ion-inertial length are expected to increase with growing solar proximity. More significantly, the separation between the two will also increase compared to their values at 1AU. At 1AU the frequencies are similar enough that there remains ambiguity as to which value is the true break-point frequency. One of the science goals of Solar Probe Plus is to distinguish which value is the actual breakpoint frequency, thereby determining the roles that turbulent heating and ion-cyclotron heating have in coronal temperatures.
To generate a synthetic time-series, a PSD is defined that describes the desired δV fluctuations. The PSD equation will differ depending on the turbulence model and dissipation mechanism. Once the PSD is defined by a continuous function (see Eq 6), the PSD's domain is discretized according to, where n time is the number of measurements in the desired time series result, and n P SD is the number of discretized values from the estimated PSD (Blackman and Tukey, 1958). Once the number of discrete PSD values is known, the discrete Fourier transform (DFT) solves for frequencies at, where T 0 is the total time length of the desired time series. These frequencies designate the values of the mean squared amplitudes A(k) 2 in the PSD. The square root of the A(k) 2 yields the average amplitude as a function of wave number or frequency. Now construct a wave function for δV in wave number space, where φ(k j ) is a randomly-chosen phase angle on [0, 2π) for each wavenumber or frequency j. The PSD by definition has no phase information, so a randomized phase must be substituted in the wave function. Thus, taking the inverse real DFT of the series in Equation 19 yields the synthetic time series δV (t) fluctuations, which are then summed with the bulk velocity and used as input to the SPC model. The upper and lower boundaries of the PSD are determined by 1) the length of the observation desired for analysis, and 2) the Nyquist frequency. The total length of the observation determines the slowest detectable frequency. The inverse of the length of the observation yields the slowest detectable frequency, since a signal of that frequency can be contained only once over the length of the observation (Blackman and Tukey, 1958).   The rate is the only free parameter; magnetic field, density, temperature, and bulk velocity are constant. The proton gyrofrequency Ω p is indicated by a vertical dashed line, and is also constant.
Here, 0.0005Hz is chosen as the slowest frequency; thus data is generated for a 2000 second period (33.3 minutes). Thus, the time series includes 75000000 measurements at 37500Hz. The chosen 0.0005 Hz frequency is sufficiently below the turbulence bendover scale separating the energy-containing range from the inertial range (Borovsky, 2012, Leamon et al., 1998.
To vary the PSD as a function of radial distance from the sun, first construct the PSD function using an isotropic Kolmogorov turbulent plasma as an example. The kinetic energy PSD is expressed as, where C K is an experimentally derived constant as described in (Vasquez et al., 2007), and is the turbulent heating or dissipation rate. Determining the radial evolution of the heating rate allows us to study the radial evolution of the corresponding velocity fluctuations as measured by SPC. To determine this dissipation or heating rate as a function of radius, the example set byVasquez et al. (2007) and developed in Verma et al. (1995) is followed, in which the solar wind expansion experiences adiabatic cooling, in-situ heating, and all of the energy dissipation manifests as proton temperature T pr . Equation 21 implies a radial proton temperature dependence of T pr = T pr,0 r 0 r ξ where ξ has hitherto been determined by spacecraft measurement, and varies depending on the heliospheric location. For a system experiencing adiabatic cooling alone, ξ = 4/3. For radial distances  Figure 4: An example of the PSDs that are used to generate synthetic turbulent time series. In this example, the inertial range spectral index is −5/3, the dissipation range spectral index is −3, and the breakpoint frequency is the local cyclotron frequency. The bulk velocity for the data above is 600km/s. between 0.3 AU and 1.0 AU however, ξ = 0.9, in support of in-situ heating (Totten et al., 1995, Vasquez et al., 2007. Solving for the heating rate, and substitute ξ = 0.9 for r = (0.3, 1.0) AU to arrive at the radial expression for the heating rate, where T pr,0 is the proton temperature at r 0 . An illustration of the PSD's variance with respect to radius is shown in Figure 3. The synthetic turbulent time-series are generated using all combinations of parameters shown in Table 3. Although the density and magnetic field were expressed as functions of radius previously to solve for radial expressions of the cyclotron frequency Ω p and ion-plasma frequency ω p,i , we choose to keep them constant here. The density

Faraday Cup Model
The Faraday cup's enduring success may be attributed to a simple design, robust against deterioration over time. An advantage that the Faraday cup has over other thermal plasma instrumentation is that it can make angular measurements simultaneously to energy measurements. Likewise, the instrument response can be modeled with an analytic response function. Successful science missions that have flown Faraday cup instruments include Explorer 10, IMP-8, Voyager 2, Wind, DSCOVR (formerly Triana), and others (Bonetti et al., 1963, Bellomo and Mavretic, 1978, Bridge et al., 1977, Ogilvie and Chornay, 1995, Valero et al., 1999. A summary of the differences between well known Faraday cups is shown in Table 4. To produce a single ion spectrum, a Faraday cup scans incrementally through a preset number of energy windows, each window bounded by specific voltage boundaries (V i , V i + ΔV i ) induced on the modulator grid 7. A predetermined number of windows span SPC's full energy range without overlap. Usually, the window width increases for each successive window as SPC scans upward through the energy windows. Window spacing is defined by choosing a minimum measurable energy and a desired window width as a percentage of the bottom of each energy window, Hence, the width of the windows is a constant percentage ΔE/E of each E bottom . SPC uses narrow windows according to Eq 24 at the bottom of its dynamic energy range, changing to fewer windows with broader energy widths at the upper end of its energy range.  Thus, the Faraday cup instrument response function is dependent on the size of, and rate of scanning through, successive energy windows. For a fully detailed derivation of the Faraday cup response function for a square wave modulated Faraday cup, see the dissertations of Kasper (2003) and Maruca (2012). In their derivation they begin by assuming the plasma is well described as a Maxwellian, where the total plasma velocity v has been rewritten in terms of thermal speed w = 2k B T m , and the bulk speed U , so that v − U = w. Defining a segment of phase space in terms of a current Equation (26) describes the instrument response in units of instrument current within one window of width ΔV between V o and V o + ΔV during square wave modulation, assuming the plasma sufficiently resembles a Maxwellian distributed plasma. However, to measure kinetic scale phenomena, which may have frequencies nearing 60Hz at SPP perihelion, SPC uses a sine wave modulation technique (Kasper et al., 2015). Thus, the SPC computer model analytically solves Equation 26 for each point measured at 37,500 Hz. Because the range of each energy window is controlled only by the voltage bias on the modulator grid, the energy windows are highly configurable during the mission operation phase. For SPC, three main operational modes have been identified:  (i) Full Scan (FS) mode: SPC scans through all the energy windows incrementally from lowest to highest, and then back down to lowest before repeating the scan cycle. (ii) Peak Tracking (PT) mode: SPC executes one Full Scan, identifies the energy window with the largest current response, and scans incrementally through a predefined number of windows centered about the peak of the distribution. (iii) Flux or Flow Angle (FA) mode: SPC executes one Full or Peak Tracking scan, identifies the energy window with the largest current response, and samples the current in that window alone.
The time required for each of the scans outlined above is a function of the modulation frequency, the number of windows used, and the number of modulation cycles in each window. Typically each window consists of 8 integration cycles, the first two of which are discarded, leaving 6 cycles' worth of samples for each window. Therefore, the number of seconds for each ion spectrum to be taken is, SPC's sine wave modulation is sampled discretely, 32 times per cycle, yielding a total of 256 samples per 8-cycle energy window. Only the last 192 samples are used to create the measured ion spectrum. For each energy window, these 192 measured points of Level 0 data are processed using Fourier analysis to create Level 1 data (Kasper et al., 2015).
Of the three operational modes, the Flux Angle (FA) mode contains the least information, followed by Peak Tracking (PT), while the Full Scan (FS) mode obtains the full electron and ion spectrum. Table 5 details the expected data output from each of these scans. In a worst case operational scenario where the PT mode does not function as expected, the FS mode is the default operational mode. For this work we consider only the FS mode.  Shown above are characteristics of each SPC operational mode. Electrons are only measured during the Full Scan mode. Bulk velocity is measured in every mode, but thermal velocity is not measured during Flux Angle mode. These values are subject to change before launch.

Results
Using the method described in section 3, and the model in section 4, all possible combinations of the values in Table 3 are used to generate synthetic turbulent spectra. The velocity fluctuations generated from Table 3 were used as input for the SPC model. The resultant Level 2 data ( Figure 2) is used to generate PSDs, an example of which is shown in Figure 10. Tables 7 and 8 show the results when the first third and last third of the data are fit according to a power law. The spectral index is indicated by k. The lower and upper thirds of the PSD frequencies are fit separately to prevent data points close to the breakpoint frequency contributing to the fit parameters. Spectral indices corresponding to the last third of the data roughly correspond to the dissipation range spectral indices; spectral indices corresponding to the first third of the data roughly correspond to the inertial range spectral indices.  Table 6 lists fit values for simulated SPC data produced by exposure to a plasma where the spectral index of −5/3 is extended to the entire frequency range, e.g. without dissipation. From Table 6 it is shown that the SPC does not produce accurate fits for 2D turbulence, as expected since SPC is insensitive to fluctuations parallel to the instrument normal vector. Thus, 2D turbulence is consciously neglected in this chapter's remaining analysis. However, Table 6 shows that the SPC model produces data that underestimates the B parameter of the inertial range of the slab Kolmogorov spectrum, where B = −1.66 or −5/3 is expected. Reasons for this under-estimation of the PSD power are explored further in this section.  Two breakpoint frequencies are used for the remaining analysis of slab turbulence in this chapter: the ion-inertial length, and the ion-cyclotron frequency. Table 7 contains the results from Kolmogorov turbulence in the inertial range, with a −3 spectral index imposed in the dissipation range. Again, the simulated SPC data consistently underestimates the magnitude of the spectral index in both the inertial and dissipation ranges. Note that the uncertainty is significantly lower in the inertial range than those for a straight Kolmogorov spectrum ( Table 6), but that the uncertainty for the dissipation range is much higher. The results for a both breakpoint frequencies are similar (Table 7). Introducing Iroshnikov-Kraichnan (IK) turbulence (spectral index = −3/2) into the initialization plasma produces similar results as seen for Kolmogorov turbulence, as seen in Table 8. The ideal spectral index B = −1.5, but instead the spectral index is again undervalued by the analysis. Likewise, the dissipation range spectral index magnitudes are also underestimated. See Figure 10 15th AIAC: "The Science of  Likewise, for IK turbulence with an ion inertial length breakpoint frequency, the spectral index fits for the inertial range are lower in magnitude than expected (Table 8).
To isolate the source of the PSD underestimation within the inertial range, the simulated SPC Level 0 data are compared to the input synthetic turbulent plasma. The Level 0 data contains samples taken at 37,500 Hz. The SPC energy windows in the Level 0 data with the highest signal for each input plasma are selected and five windows on either side of the maximum are used for a total of 11 energy windows analyzed. To keep the cadences identical between Level 0 and Level 2 data, the first Level 0 measurement in each window are used. These values are compared to the idealized plasma input at the same measurement time as the maximum modulator value. The results are shown in Table 9.
As the results in Table 9 show, the simulated Level 0 data does not accurately reproduce Kolmogorov or IK parameters. The inaccuracy here is attributed to the Level 0 measurement units, which are in instrument units of current and not in physical units as the synthetic turbulent time series are. When a single point from each L0 data energy window is used to to generate a PSD, the resultant power is grossly in error, 25.6% and 31.2% from the ideal Kolmogorov or IK inertial range value. However, when the corresponding samples from the initial synthetic turbulent time series are analyzed in an identical manner, the resultant power is only 2% in error, which is in agreement with the model's established accuracy (Whittlesey, 2016). Furthermore, the fit values from the window-sampled initialization plasma are remarkably similar to the fit values from the Level 2 data analysis (Table 10).
The resultant spectral index in the inertial range for all values are on the order of 2%. Considering the limitations of the model are within a 2% margin of error (Whittlesey,    2016), these values are considered a good metric for the accuracy of both the model and instrument's positive ability to resolve the spectral index in both the inertial and dissipation ranges. However, the source of the error in spectral index is yet unfound. To explore the source of error, a set of synthetic turbulent data have been generated for a range of spectral indices ranging between -1.0 and -4.0 for R = 0.2AU , with the same density, bulk speed, thermal speed, and as previously used. The synthetic data are sampled at the same rate as in the comparison to Level 0 SPC data. The results are plotted in Figures 11 and 12. It is apparent that undersampling the synthetic turbulent time series at the SPC Level 2 data rate results in a ∼2% error for inertial range spectral indices of magnitude greater than 1.5, but that the error is greater for indices of lesser magnitude.

Conclusion: Velocity PSD Spectral Indices are Always Undermeasured
This paper has established that the SPC model and the spectral analysis detailed therein reproduce inertial range spectral indices to a high degree of precision, but a low degree of absolute accuracy, for slab and turbulence. Except for circumstances where the bulk flow is at highly oblique angles to the instrument normal axis, 2D turbulence is not distinguishable from noise in SPC. Examining simulated Level 2 slab turbulence data indicates that the magnitude of the inertial range spectral index is consistently underestimated by 2% of the input value. This underestimation is consistent for both spectral indices chosen. The dissipation range accuracy and precision are considerably less accurate, and with higher variance, likely due to the higher noise levels observed in the dissipation range of the PSD. Significantly, the spectral index is consistently underestimated in both the inertial and the dissipation range of the inertial range when discretely sampled in a style consistent with thermal plasma instrumentation. This underestimation has been quantified at 2% for SPC's measurement cadence for spectral indices of magnitude greater than 1.5, with possible extension to other thermal plasma instruments of differing cadence.
Finally, note that traditional turbulent spectral analysis typically uses magnetic field measurements instead of velocity measurements (e.g. Leamon et al. (1998), Wicks et al. (2010), Bruno and Carbone (2005) and others) due to faster measurement cadences prevalent in magnetic field instrumentation. Thus, it is possible that the higher measurement cadence in measurement by these instruments prevents such a spectral under estimation.