# Analysis of a JWST NIRSpec Lab Time Series: Characterizing Systematics,
Recovering Exoplanet Transit Spectroscopy, and Constraining a Noise Floor

Zafar Rustamkulov<EMAIL_ADDRESS>Department of Earth and Planetary Science,
Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA
David K. Sing Department of Earth and Planetary Science, Johns Hopkins
University, 3400 N. Charles Street, Baltimore, MD 21218, USA  Department of
Physics and Astronomy, Johns Hopkins University, 3400 N. Charles Street,
Baltimore, MD 21218, USA  Rongrong Liu Department of Physics and Astronomy,
Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA
Ashley Wang Department of Physics and Astronomy, Johns Hopkins University,
3400 N. Charles Street, Baltimore, MD 21218, USA

###### Abstract

The James Webb Space Telescope’s NIRSpec instrument will unveil the nature of
exoplanet atmospheres across the wealth of planet types, from temperate
terrestrial worlds to ultrahot Jupiters. In particular, the 0.6-5.3 $\mu$m
PRISM mode is especially well-suited for efficient spectroscopic exoplanet
observations spanning a number of important spectral features. We analyze a
lab-measured NIRSpec PRISM mode Bright Object Time Series (BOTS) observation
from the perspective of a JWST user to understand the instrument performance
and detector properties. We create two realistic transiting exoplanet time
series observations by performing injection-recovery tests on the lab-measured
data to quantify the effects of real instrument jitter, drift, intrapixel
sensitivity variations, and 1/$f$ noise on measured transmission spectra. By
fitting the time series systematics simultaneously with the injected transit,
we can obtain more realistic transit depth uncertainties that take into
account noise sources that are currently not modeled by traditional exposure
time calculators. We find that sources of systematic noise related to
intrapixel sensitivity variations and PSF motions are apparent in the data at
the level of a few hundred ppm, but can be effectively detrended using a low-
order polynomial with detector position. We recover the injected spectral
features of GJ 436 b and TRAPPIST-1 d, and place a 3-$\sigma$ upper limit on
the detector noise floor of 14 ppm. We find that the noise floor is consistent
with $<$10 ppm at the 1.7-$\sigma$ level, which bodes well for future
observations of challenging targets with faint atmospheric signatures.

planets and satellites: atmospheres — techniques: spectroscopy

††facilities: JWST(NIRSpec)††software: Matplotlib (Hunter, 2007), Numpy (van
der Walt et al., 2011), Scipy (Virtanen et al., 2020), astropy (Astropy
Collaboration et al., 2013), dynesty (Speagle, 2020), lmfit (Newville et al.,
2014)

##  1 Introduction

The James Webb Space Telescope (JWST), now in its L2 orbit, will soon shed new
light on transiting exoplanet atmospheres. It will probe at high signal-to-
noise ratio (SNR) planets spanning the phase space of mass, radius,
temperature, and age, yielding our most comprehensive view yet. Recent
discoveries by the Transiting Exoplanet Survey Satellite (TESS) (Ricker et
al., 2015) are beginning to add to an ever-growing catalog of nearby
exoplanets that are excellent targets for efficient spectral observations with
JWST (e.g., Winters et al., 2019; Jenkins et al., 2020; Hobson et al., 2021;
Vanderspek et al., 2019). A variety of spectral modes are available with
NIRSpec, each with unique wavelength coverage and spectral resolution.
NIRSpec’s versatility made it an extremely popular mode for Cycle 1 GO
proposers, with roughly two-thirds of awarded exoplanet time series programs
making use of one or more of its modes111https://www.stsci.edu/jwst/science-
planning/user-committees/jwst-users-committee. Several targets on the
Guaranteed Time Observations (GTO) and Early Release Science (ERS) lists also
widely make use of NIRSpec.

NIRSpec’s low-resolution ($\lambda/\Delta\lambda$ $\sim$ 20-100) PRISM mode
has a uniquely broad spectral range (0.6 - 5.3 $\mu$m), offering a doubling in
throughput relative to other JWST mode combinations, which require a minimum
of two separate transits to achieve the same wavelength coverage. This
property makes PRISM the ideal instrument to carry out efficient, ‘one-shot’
measurements of a number of carbon, oxygen, and even nitrogen-bearing
molecules in exoplanet atmospheres, in all three of the main observational
geometries (i.e. Batalha & Line, 2017a; Morley et al., 2017; Lustig-Yaeger et
al., 2019). Measurements of these information-rich chemical species will allow
for exquisite comparative studies between planets, and will provide important
insights into the mechanisms that sculpt planetary atmospheres and
populations.

In this work we used a 3-hour Bright Object Time Series (BOTS) dataset
measured from the NIRSpec instrument in the lab to look at the data from the
perspective of a JWST observer to quantify how the instrument’s systematics
and noise performance affect the final scientific product. The Pandexo
(Batalha & Line, 2017b) and Pandeia (Pontoppidan et al., 2016) packages have
served as excellent community tools for estimating exposure times and building
synthetic JWST observations. We build off the recent work by Birkmann et al.
(2022) and Jakobsen et al. (2022), who characterize many important
engineering-side detector diagnostics, and show how these nuances will affect
a real exoplanet observation. Similar characterization work with NIRCam has
shown that the detector’s noise floor from known systematics sources is 9 ppm,
and that its systematics are similar in nature and amplitude to the ones
outlined in our work (Schlawin et al., 2020, 2021).

As is often seen in time series data from Spitzer IRAC (e.g., Deming et al.,
2015; May & Stevenson, 2020), TESS, and Kepler/K2 (Luger et al., 2016, 2018),
undersampled PSFs often introduce a non-negligible source of systematic noise
due to the PSF’s motion across pixels with nonuniform inter- and intrapixel
sensitivity variations. Thankfully, with knowledge of the PSF’s position over
time, generalized pixel-level decorrelation methods (PLD) are able to largely
remove these signatures in photometry. With exoplanet time series
spectroscopy, analogous jitter decorrelation techniques (e.g., Sing et al.,
2019) have shown promise in removing systematic trends and greatly improving
light curve fits. In our work, we leverage a $\sim$0.3-pixel, $\sim$30-mas
drift of the light source on the detector to measure the amplitude of shift-
induced systematics in the light curve and the extent to which they can be
decorrelated. Further, we fit the observed systematics jointly with the
transit light curve at each wavelength channel of our injected transit
spectrophotometry, and perform a uniform model intercomparison to predict
which systematics vectors are important in the detrending of real JWST time
series data. Lastly, we select a wavelength channel that displays weak
systematic noise and bin it down through time to place an empirical upper
limit on NIRSpec’s noise floor.

##  2 Overview of the CV3 PRISM Time Series

We analyze an experimental dataset collected in January of 2016 during the
cryo-vacuum testing campaign (CV3) of the JWST Integrated Science Instrument
Module (ISIM) at NASA’s Goddard Space Flight Center. The single BOTS exposure
was taken with NIRSpec held at a temperature of $\sim$40 K (its L2 operating
temperature), and consists of 12,000 integrations with three groups each,
corresponding to an effective integration time of 0.45 seconds. The exposure
used the SUB512 subarray, with size 512$\times$32 pixels, which is wide enough
to encapsulate the PSF of the light source and the unilluminated region of the
detector. The NRSRAPID readout pattern and S1600A1 (1.6”$\times$1.6”) slit
were used. The data was processed by the NIRSpec ramp-to-slopes pipeline to
fit the up-the-ramp count rate, the statistical variance, and quality flag for
each pixel (see Birkmann et al., 2022, Section 7). The CV3 BOTS prism data is
publicly available through the ESA
website222https://www.cosmos.esa.int/web/jwst-nirspec/test-data. Full details
of the NIRSpec CV3 observation can be found in Birkmann et al. (2022), which
we briefly summarize here.

Figure 1: Overview of the lamp spectrum. Top: a single 0.45s integration of
the lamp spectrum summed along the cross-dispersion axis. Prominent frost
absorption features in NIRSpec’s optical path sculpt the spectrum’s shape.
Bottom: a single log-stretched raw NIRSpec PRISM frame. The aspect ratio of
the detector is vertically exaggerated by a factor of $\sim$2 for clarity.

The NIRSpec instrument was illuminated with a tungsten filament source with a
similar PSF width as expected from JWST (see Fig. 1). The PSF displays a
spatial drift of $\sim$30 mas ($\sim$0.3 px) over the course of the exposure,
which $\sim$4$\times$ larger than the 1-$\sigma$ RMS stability requirement for
JWST. The light source displays a long term drift in the flux level at the
percent-level over several hours, as well as chromatic flickers at the
$\sim$0.2$\%$-level, active on the timescale tens of minutes. As is apparent
in Fig. 2c, these flickers have a stronger wavelength-dependence at the blue-
end of the detector, suggesting that they are related to sudden temperature
changes of the filament. The spectrum further shows deep and shallow
absorption features presumably due to condensed ices along the optical path.
Some of the absorption features change somewhat in depth throughout the
exposure, with a notable concave-down profile at 3.55-3.7 $\mu$m that imparts
an erroneous excess in the extracted transit depths in this channel. The first
hour of the time series shows a strong, chromatic increase in flux due to the
lamp warming up, and the chromaticity of the flickers is stronger here. We
therefore use the last 213 = 8192 integrations, or $\sim$2.1 hours, of data to
mitigate the strong spectral systematics and to better distill the detector’s
performance. The light source, whose spectrum is far removed from a typical
stellar spectrum in shape, has an equivalent J-band magnitude of roughly
J$\sim$11, which is approaching PRISM’s saturation limit.

##  3 Data Extraction and Reduction

In order to prepare our data for injection-recovery testing, we first extract
the spectrophotometry from the time series of the raw frames. The overall
approach to reduce and calibrate NIRSpec data has been described by Rawle et
al. (2016) and Alves de Oliveira et al. (2018) and is also summarized in
Birkmann et al. (2022) and Jakobsen et al. (2022). For each integration, we
perform the following steps: we first apply the flat field correction and save
the wavelength solution. We next subtract the 1/$f$ noise, as measured by the
vertical median of the unilluminated regions of the detector at each pixel
column. This process is also referred to as “destriping.” The 1/$f$ noise
component, whose contribution can be seen in Fig. 2a, imparts a near-doubling
to the white noise level to the spectrophotometry if not correctly removed.
Next, we measure the vertical and horizontal pixel shift of the trace in each
integration relative to the first by performing a 1-D cross correlation
between the vertically- and horizontally-summed counts (Fig. 2b). We then
remove these shifts and align spectra by interpolating each frame with a 3rd-
order spline and re-gridding the image data onto the shifted pixel array. By
aligning the spectra, we remove a common time series systematic effect where
the flux of a pixel is observed to change with time as spectral features drift
across a pixel.

Next, we fit a Gaussian profile to the PSF at each pixel column to measure the
trace position and width. We smooth the trace profile with a median filter,
fit it with a 6th-order polynomial, and perform an sub-pixel extraction of the
counts enclosed within 5$\times$FWHM of the trace center. Lastly, we convert
from count rate to counts by multiplying by the total integration time of each
frame (0.45 s). At this stage, the flat-fielded, 1/$f$ corrected, shift-
stabilized, and spatially summed spectrophotometry resembles the true,
distilled behavior of the temperamental light source with most of the
confounding artifacts removed. The extracted spectrophotometry can be seen in
Fig. 2c, where light source flickers and flux variations imparted by x-shifts
across sharp spectral features show up strongly in the relative flux at the
$\pm$1$\%$-level.

Before injecting transit spectra into the time series, we first apply a
common-mode correction to each wavelength pixel in order to suppress the
effects of the lamp flicker. This common-mode technique to remove systematic
errors is an often-used method to remove time series trends which are the same
or similar at each wavelength (e.g., Berta et al., 2012; Deming et al., 2013;
Sing et al., 2013). While simply dividing the time series of each pixel by the
white light curve can remove most of the observed lamp flickers and long-term
flux drifts, the amplitude of these variations are found to also have a color
dependence, which we remove by applying the common-mode correction on a pixel-
by-pixel basis. This is done by performing a linear regression fit of the
white light curve (see Fig. 2a) to the spectral light curves at each
wavelength pixel, and dividing it out. This correction removes the majority of
the light source flicker while preserving the detector’s noise properties,
intrapixel systematics, and most of the light source’s chromaticity. The self-
correction effect imposed by division of the white light curve can only be
seen in the post-corrected white light curve, which by definition is a flat
line at 1. The common-mode-detrended spectrophotometry is shown in Fig. 2d,
where modest residual chromatic light source systematics are seen. At this
stage, the spectrophotometry is largely free from confounding trends and
flickers from the light source.

##  4 Transit Spectrum Injection

Injection-recovery tests are often repeated with several draws from a large
dataset or noise model, but our study relies on a single real observation, and
therefore a single noise instance. We perform our tests on two different
planet cases, as two separate manifestations of the same noise instance, to
show the consistency of our injection-recovery process. We inject the transit
spectrophotometry of two planets: the warm Neptune GJ 436 b (T${}_{eq}\sim$700
K, J = 6.9) and the temperate terrestrial planet TRAPPIST-1 d
(T${}_{eq}\sim$300 K, J = 11.35). We choose these planets due to their short
$\leq$1 hr transit duration and their broad interest to the exoplanet
atmosphere community.

GJ 436 b is slated for transmission spectroscopy with NIRSpec G395H (Stevenson
and Lustig-Yeager et al., GO 1981), and TRAPPIST-1 d will be observed in
transit with the NIRSpec PRISM (Lafreniere et al., GTO 1201). With
TRAPPIST-1’s J-band magnitude of 11.4, the count rates from the CV3 light
source are similar in our dataset to what will be seen on-sky with JWST,
making it a compelling comparison. GJ 436A, however, is 4 magnitudes brighter
(prohibitively bright for PRISM) and the expected precision of its measured
transmission spectrum using other modes is significantly higher than what we
present. For GJ 436 b, we use a transmission spectrum from the Goyal Grid of
generic models (Goyal et al., 2020) with parameters most closely matching
those of GJ 436 b. The TRAPPIST-1 d spectrum is taken from the work of
Lincowski et al. (2018) and Lustig-Yaeger et al. (2019) where a CO2-dominated
10-bar atmosphere is assumed, with and without high-altitude H2SO4 clouds.

For each planet, we bin our model transmission spectra to the wavelength grid
of each PRISM pixel, and compute 4-parameter limb darkening coefficients
between the pixel edges using values from the 3D STAGGER-grid (Magic et al.,
2015). These LD coefficients are weighted by NIRSpec’s response function and
the respective stellar spectrum. Given the narrow wavelength bins used here,
we find that the errors in the LD coefficients introduced by using a stellar
spectral weighting rather than the lamp’s spectrum are very small. We then
compute the transit light curve time series at each wavelength pixel using
batman (Kreidberg, 2015).

Finally, we inject the spectral transit model while simultaneously dividing
out the common-mode correction into the raw data, and reinsert the original
time series’ 1/$f$ noise and jitter. At this stage, our injected files mimic a
JWST observation as closely as possible, with the jitter exaggerated by a
factor of 4.

Figure 2: Overview of the NIRSpec PRISM CV3 time series observation. (a) The
mean-normalized white light curve with and without 1/$f$ correction. (b) The
x- and y-shifts of the PSF relative to the middle integration. The grey region
shows the 1-$\sigma$ RMS jitter requirement of JWST over 104 seconds. (c) The
raw relative spectrophotometry after shift-stabilization, 1/$f$ noise
subtraction, and extraction. Note the strong light source flickers and trends.
(d) same as (c), but after common-mode correction that mostly removes the
systematics imparted by the unstable light source.

##  5 Transit Spectrum Recovery

We reduce our transit-injected data files following the procedure described in
the first paragraph of Section 3. We bin the transit spectrophotometry into
variable-width wavelength channels with roughly equal counts. We use 51 and 19
wavelength channels for GJ 436 and TRAPPIST-1, respectively.

We model the systematic errors with a parameterized deterministic model, where
the photometric trends are found to correlate with a number of external
parameters (or optical state parameters, x). These parameters describe changes
in the instrument or other external factors as a function of time during the
observations, and are fit with coefficients, $p_{n}$, for each optical state
parameter to detrend the light curves. We then fit each spectral light curve
jointly with a transit model and systematics model $S(\mathbf{x})$ of the
form:

$S(\mathbf{x})=p_{1}X+p_{2}Y+p_{3}X^{2}+p_{4}Y^{2}+p_{5}XY+p_{6}\Delta_{\mathrm{CM}}$
(1)

where $X$ and $Y$ are the mean-subtracted, standard deviation-normalized x-
and y-shifts, respectively, and $\Delta_{\mathrm{CM}}$ is the pre-injection
normalized common-mode light curve, normalized in the same way. The full light
curve flux model, $f(t)$, is then assembled with the equation

$f(t)=f_{0}\times T(t,\boldsymbol{\theta})\times S(\mathbf{x})$ (2)

where $f_{0}$ is the baseline flux and $T(t,\theta)$ is the theoretical
transit model which depends on the transit parameters $\theta$. This model
captures the behavior of jitter-induced intrapixel variations, as well as
leftover light source variability. The detrending vectors are normalized so
that their fit coefficients may be directly compared. We perform joint fits of
the systematics model and the planet’s transit depth (while keeping all other
transit parameters fixed) using the Python-based Levenburg-Marquardt least
squares algorithm lmfit (Newville et al., 2014). Using Python’s
multiprocessing package to parallelize the fitting routine, we fit all 51 GJ
436 b channels in 30 seconds using 16 cores. We initially assume photon noise
uncertainties for the light curve, re-fit with the residual standard deviation
as the uncertainty, and then inflate the final lmfit parameter uncertainties
on the transit depths by a factor $\beta(\lambda)$. $\beta(\lambda)$ is
computed from the fit-residuals of each spectral light curve taking into
account both white and red noise sources. Its median value is 1.35, and it
peaks where there are residual systematics due to chromatic lamp flickers and
changing cryogenic ice features. $\beta(\lambda)$ is near 1 in wavelength
regions where these two effects are small. If $\beta(\lambda)$ were not
multiplied into our transit depth uncertainties, the residual chromatic
flicker effects would not be taken into account, giving unrealistic errorbars.
Further, we avoid removing these residuals using stochastic models such as
Gaussian Processes (GPs) since they could remove systematics due to detector’s
noise properties which we aim to characterize.

Table 1: Systematics Fitting Results Summary

Scenario | Description | $N_{\textrm{params}}$ |  $\sigma_{\mathrm{TS\ resid}}$ (ppm) | $\chi^{2}_{\nu}$ |  $\Delta$BIC |  $\Delta$AIC
---|---|---|---|---|---|---
1 | Fiducial | 8 | 148.70 | 1.32 | – | –
2 | No $X^{2}$ | 7 | 148.74 | 1.25 | -6.15 | -2.48
3 | No $XY$ | 7 | 149.30 | 1.31 | -5.01 | -1.35
4 | No $Y^{2}$ | 7 | 151.21 | 1.28 | -4.87 | -1.20
5 | No $X^{2},Y^{2}$ | 6 | 152.38 | 1.28 | -7.37 | -0.08
6 | No $Y$ | 7 | 171.28 | 2.32 | 39.09 | 42.75
7 | No $X$ | 7 | 184.86 | 2.25 | 36.15 | 39.85
8 | Fiducial with $X^{3},Y^{3}$ | 10 | 243.70 | 3.36 | 88.69 | 81.27
9 | No $X^{2},Y^{2},XY$ | 5 | 267.5 | 4.33 | 130.17 | 141.02
10 | No $X^{2},Y^{2},XY,X$ | 4 | 285.93 | 4.84 | 154.76 | 169.10
11 | No $X^{2},Y^{2},XY,Y$ | 4 | 290.39 | 6.05 | 211.69 | 226.03
12 | No $X^{2},Y^{2},XY,X,Y$ | 3 | 293.37 | 5.66 | 195.09 | 212.79
13 | Fiducial, no flat field | 8 | 390.98 | 7.12 | 248.97 | 248.97
14 | Fiducial, no 1/$f$ subtraction | 8 | N/A | 555.69 | 23837.78 | 223837.78

Results from a suite of tests for goodness-of-fit and systematics model
quality sorted in decreasing order of residual scatter. Each row represents an
independent transmission spectrum extraction from our injected transit
spectrophotometry of GJ 436 b, fitted with some linear combination of 1st,
2nd, and 3rd order $X$ and $Y$ jitter vectors, and the cross-term $XY$.
$\sigma_{\mathrm{resid}}$ is the standard deviation of the residual between
the injected and recovered transmission spectrum, and $\chi^{2}_{\nu}$ is the
reduced $\chi^{2}$ statistic of the fit. The scatter in Scenario 14 is not
reported since the noise level was too overwhelming for the fit to prefer a
non-zero transit depth.

We also use the nested sampling package dynesty to fit the light curves
(Speagle, 2020), and find equivalent results in both the fitted depth and the
uncertainties. Moreover, we find that there are no complicated degeneracies
between the transit depth and the optical state vectors, indicating that the
posterior is well-described by a multivariate Gaussian distribution. This
means that a Levenburg-Marquardt least squares minimizer is sufficient to
quickly estimate the posterior distribution with accurate uncertainties.

This fit to the transmission spectra of both planets, which we term the
“fiducial” fit hereafter, match the injected spectra well, with a reduced chi-
squared $\chi^{2}_{\nu}$ of 1.32 for GJ 436b, and 1.15 for TRAPPIST-1 d (see
Fig. 5). In the latter case, the extracted spectrum qualitatively shows
increases in transit depth at the three main absorption features. The
$\chi^{2}_{\nu}$ for GJ 436 b’s fiducial fit without taking into account
errorbar inflation is 3.13. The residuals are largely Gaussian, with the
notable exception of the peaky 3.6$\mu$m feature, which is imparted by a
downward concavity in the ice deposition profile there. This same feature can
be seen in Fig. A.2, panel f, of Birkmann et al. (2022). The average amplitude
of the $X$ and $Y$ shift vectors are found to be 455 and 445 ppm,
respectively, across both spectra. This is in excellent agreement with the
jitter systematic noise modeled by Birkmann et al. (2022). Rescaling these
amplitudes linearly by a factor proportional to the expected jitter of JWST,
we estimate the on-sky jitter-induced photometric signatures to be at the
$\sim$100 ppm level.

Figure 3: The recovered transmission spectra of GJ 436 b (top), and TRAPPIST-1
d (bottom) recovered from the injected data using our novel pipeline. We
include a cloudy TRAPPIST-1 d atmosphere model for an illustrative comparison.

##  6 Systematics Model Intercomparison

We explore the relevance of each term in the systematics model, as well as the
flat field and 1/$f$-subtraction, by performing a model intercomparison on the
GJ 436 b data. We re-fit the transit spectrophotometry using the methods
described in Section 5, but with key elements added or removed to see how the
goodness-of-fit of the final transmission spectrum changes. We recommend this
technique as an informative method for systematics model selection, and to
explore the response of the final spectrum to different procedures in the data
reduction. At each scenario, we compute the standard deviation of the fit
residual, the reduced Chi-squared statistic $\chi^{2}_{\nu}$, and the change
in Bayesian and Akaike Inference Criteria ($\Delta$BIC and $\Delta$AIC,
respectively) relative to the fiducial model. We note that in each instance,
we keep the uncertainties on the transit depths fixed to those of the fiducial
fit in order to avoid erroneously inflating the errorbars due to poor model
selection, and to preserve the uniformity of the comparison. The results of
this intercomparison are in Table 1.

Decorrelation against the linear jitter vectors $X$, $Y$, and $XY$ all
substantially improve goodness-of-fit over a model without them. The second-
order terms $X^{2}$ and $Y^{2}$ are marginally, if at all, important, as
evidenced by the negative changes to their Bayesian and Akaike Information
Criteria relative to the fiducial model ($\Delta$BIC and $\Delta$BIC,
respectively). This indicates that the intrapixel sensitivity profile of the
typical NIRSpec pixel is well-approximated by a low-order polynomial. Adding
the third-order terms $X^{3}$ and $Y^{3}$ to the systematics model results in
a much poorer goodness-of-fit, likely due to the even nature of the detector’s
intrapixel sensitivity function.

Figure 4: (a) The injected transit light curve photometry of GJ 436 b in the
relatively systematics-free 3.77-4.17 $\mu$m channel (tomato points), and the
best-fitting joint model (grey line). b Residual scatter as a function of bin
size for the 3.77-4.17 $\mu$m channel (black points) with the best-fitting
noise components (solid lines) and the photon shot noise limit (dashed line).

Flat fielding has a major impact on the quality of the resulting fit, nearly
halving the resulting residual scatter when the flat field is applied. From
this large effect, we determine that significant efforts should be made to
obtain several, accurate flat fields for the subarrays used in JWST BOTS data.
Sub-pixel flat fields could help further reduce the impact of jitter on the
photometry. Perhaps most important, though, is the 1/$f$ correction–this added
noise component overwhelms the transit depths, resulting in non-detections of
the transit feature in all spectroscopic channels when not removed.

A number of systematic sources that are not explored by our work, or
potentially fully captured by this dataset, are likely to be present in on-sky
JWST data. Birkmann et al. (2022) discusses many potential sources of
systematic sources in more depth, including charge trapping-induced
persistence in the detector, which results in exponential “ramps” in measured
flux as is regularly seen in HST/WFC3 data (e.g Deming et al. 2013; Wakeford
et al. 2013; Kreidberg et al. 2014; Zhou et al. 2017). This effect, which
cannot be measured in the CV3 data due to the instability of the lamp, is
expected to be greatly reduced for the relatively new JWST detectors (Rauscher
et al., 2014; Birkmann et al., 2022). High-gain antenna moves, increased
cosmic rays, and thermal drift of the telescope may also introduce systematics
and are not captured in the CV3 data. Birkmann et al. (2022) estimates that
slit losses and intrapixel sensitivity variations will both be on the order of
several hundred ppm with bright object time series data, pending JWST’s real
on-sky PSF profile and pointing stability. The CV3 data displays these trends
at the expected levels given the $X$ and $Y$ position dependence of the
spectrophotometry, which are well-decorrelated by our jitter detrending
technique.

##  7 Searching for the Noise Floor

To search for the detector’s noise floor, we first select a wide wavelength
region that is largely free of systematics and light source spectral features.
A bin sufficiently narrow to minimize the enclosed chromatic systematics, and
wide enough to encapsulate maximal counts, is ideal. The relatively clean
3.77-4.17 $\mu$m region is chosen for these reasons. We then fit the light
curve of GJ 436 b at this channel as described previously, and study its
residual by binning it through time (see Fig. 4). The residual RMS bins down
closely with $1/\sqrt{N}$ (where N is the number of points per bin). We fit
this RMS curve with a quadrature-sum of the white noise component, a noise
floor, and a logistic decay law to capture the leftover chromatic flickers. We
find that the fitted white noise at this bin closely matches the expected
photon shot limit, and that the residual flicker noise contributes 35 ppm at
the $\sim$2-minute timescale, which is consistent with the flicker timescale
seen in the white light curve. The RMS curve bins to $\sim$15 ppm. Most
importantly, we recover a 3-$\sigma$ upper bound of 14 ppm for the noise
floor, with the fit preferring a $<$10 ppm floor at 1.7-$\sigma$. Given the
joint probability distribution with our systematics model, the transit depth
uncertainty at this channel is 1.35$\times$ larger than what is expected due
to photon noise alone, or 74$\%$ of the photon noise limit. We expect on-sky
JWST, with its expected weaker systematics, to reach closer to the photon
limit, making our value somewhat pessimistic.

##  8 Summary and Conclusion

In this work we analyzed NIRSpec lab time series data, injecting and
recovering transmission spectra including real detector noise, making it the
closest analog available before science data is collected. The dataset,
complete with the detector’s actual noise performance and intrapixel
sensitivity variations, informed our optimistic expectations as future JWST
users. We demonstrate that subtraction of the 1/$f$ noise is a crucial step in
the data reduction, and caution that this step may be complicated in the
narrower SUB512S subarray which is only 16 pixels tall. Further, we show that
neglecting the flat field increases the scatter of the final spectrum by a
factor of 2.6 (Table 1). Given JWST’s 7 mas (0.07 pixel) on-sky stability
requirement we expect the amplitude of jitter-induced systematics to be at the
$\sim$100 ppm-level. We demonstrate that such jitter systematics can be
readily detrended using a low-order polynomial with detector position, binning
the data down close to the photon limit. While we do not confidently detect a
noise floor, we show that it is $<$14 ppm, and likely below 10 ppm. Using our
fast reduction and analysis pipeline, we successfully recover our injected
transmission spectra, and show that Levenburg-Marquardt least-squares fitters
are sufficient to quickly estimate the posterior distribution. We also find
that, when simulating observations, transit depth uncertainties should be
scaled conservatively by a factor of 1.35 relative to the photon noise limit
expectation. All of these facts considered, we believe NIRSpec PRISM is well-
suited to characterize challenging exoplanet atmosphere targets with $\sim$10
ppm spectral amplitudes. We expect that the precision of such spectra will be
limited by the photon limit and confounding astrophysical signatures such as
stellar activity rather than instrumental systematics.

We thank Giovanna Giardino, Stephan Birkmann, and Pierre Ferruit for their
useful discussions and the NIRSpec team for their efforts in collecting and
sharing the CV3 dataset.

## References

  * Alves de Oliveira et al. (2018) Alves de Oliveira, C., Birkmann, S. M., Böker, T., et al. 2018, 10704, 107040Q
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Batalha & Line (2017a) Batalha, N. E., & Line, M. R. 2017a, The Astronomical Journal, 153, 151
  * Batalha & Line (2017b) —. 2017b, The Astronomical Journal, 153, 151
  * Berta et al. (2012) Berta, Z. K., Charbonneau, D., Désert, J.-M., et al. 2012, The Astrophysical Journal, 747, 35
  * Birkmann et al. (2022) Birkmann, S. M., Ferruit, P., Giardino, G., et al. 2022, arXiv:2202.03309 [astro-ph], arXiv:2202.03309, comment: Accepted for publication in A&A. 15 pages, 15 figures
  * Deming et al. (2013) Deming, D., Wilkins, A., McCullough, P., et al. 2013, The Astrophysical Journal, 774, 95
  * Deming et al. (2015) Deming, D., Knutson, H., Kammer, J., et al. 2015, The Astrophysical Journal, 805, 132
  * Goyal et al. (2020) Goyal, J. M., Mayne, N., Drummond, B., et al. 2020, Monthly Notices of the Royal Astronomical Society, 498, 4680
  * Hobson et al. (2021) Hobson, M. J., Brahm, R., Jordán, A., et al. 2021, The Astronomical Journal, 161, 235
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Jakobsen et al. (2022) Jakobsen, P., Ferruit, P., de Oliveira, C. A., et al. 2022, arXiv:2202.03305 [astro-ph], arXiv:2202.03305, comment: Accepted for publication in A&A. 22 pages, 19 figures
  * Jenkins et al. (2020) Jenkins, J. S., Díaz, M. R., Kurtovic, N. T., et al. 2020, Nature Astronomy, 1
  * Kreidberg (2015) Kreidberg, L. 2015, Publications of the Astronomical Society of the Pacific, 127, 1161
  * Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Nature, 505, 69
  * Lincowski et al. (2018) Lincowski, A. P., Meadows, V. S., Crisp, D., et al. 2018, The Astrophysical Journal, 867, 76
  * Luger et al. (2016) Luger, R., Agol, E., Kruse, E., et al. 2016, The Astronomical Journal, 152, 100
  * Luger et al. (2018) Luger, R., Kruse, E., Foreman-Mackey, D., Agol, E., & Saunders, N. 2018, The Astronomical Journal, 156, 99
  * Lustig-Yaeger et al. (2019) Lustig-Yaeger, J., Meadows, V. S., & Lincowski, A. P. 2019, The Astronomical Journal, 158, 27
  * Magic et al. (2015) Magic, Z., Chiavassa, A., Collet, R., & Asplund, M. 2015, Astronomy & Astrophysics, Volume 573, id.A90, $<$$>$10$<$$>$ pp., 573, A90
  * May & Stevenson (2020) May, E. M., & Stevenson, K. B. 2020, The Astronomical Journal, 160, 140
  * Morley et al. (2017) Morley, C. V., Kreidberg, L., Rustamkulov, Z., Robinson, T., & Fortney, J. J. 2017, The Astrophysical Journal, 850, 121
  * Newville et al. (2014) Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python (0.8.0), Zenodo
  * Pontoppidan et al. (2016) Pontoppidan, K. M., Pickering, T. E., Laidler, V. G., et al. 2016, 9910, 991016
  * Rauscher et al. (2014) Rauscher, B. J., Boehm, N., Cagiano, S., et al. 2014, Publications of the Astronomical Society of the Pacific, 126, 739
  * Rawle et al. (2016) Rawle, T. D., Alves de Oliveira, C., Birkmann, S. M., et al. 2016, 9904, 990446
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Schlawin et al. (2020) Schlawin, E., Leisenring, J., Misselt, K., et al. 2020, The Astronomical Journal, 160, 231
  * Schlawin et al. (2021) Schlawin, E., Leisenring, J., McElwain, M. W., et al. 2021, The Astronomical Journal, 161, 115
  * Sing et al. (2013) Sing, D. K., Lecavelier des Etangs, A., Fortney, J. J., et al. 2013, Monthly Notices of the Royal Astronomical Society, 436, 2956
  * Sing et al. (2019) Sing, D. K., Lavvas, P., Ballester, G. E., et al. 2019, The Astronomical Journal, 158, 91
  * Speagle (2020) Speagle, J. S. 2020, Monthly Notices of the Royal Astronomical Society, 493, 3132
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22
  * Vanderspek et al. (2019) Vanderspek, R., Huang, C. X., Vanderburg, A., et al. 2019, The Astrophysical Journal, 871, L24
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, doi:https://doi.org/10.1038/s41592-019-0686-2
  * Wakeford et al. (2013) Wakeford, H. R., Sing, D. K., Deming, D., et al. 2013, MNRAS, 435, 3481
  * Winters et al. (2019) Winters, J. G., Medina, A. A., Irwin, J. M., et al. 2019, The Astronomical Journal, 158, 152
  * Zhou et al. (2017) Zhou, G., Bakos, G. Á., Hartman, J. D., et al. 2017, The Astronomical Journal, 153, 211


# Evolution of switchbacks in the inner Heliosphere

Anna Tenerani  Department of Physics, University of Texas at Austin, TX 78712
<EMAIL_ADDRESS>Nikos Sioulas  Department of Earth, Planetary,
and Space Sciences, UCLA, Los Angeles, CA, 90095  Lorenzo Matteini  Imperial
College, London  Olga Panasenco  Advanced Heliophysics  Chen Shi  Department
of Earth, Planetary, and Space Sciences, UCLA, Los Angeles, CA, 90095  Marco
Velli  Department of Earth, Planetary, and Space Sciences, UCLA, Los Angeles,
CA, 90095

###### Abstract

We analyze magnetic field data from the first six encounters of PSP, three
Helios fast streams and two Ulysses south polar passes covering heliocentric
distances $0.1\lesssim R\lesssim 3$ au. We use this data set to statistically
determine the evolution of switchbacks of different periods and amplitudes
with distance from the Sun. We compare the radial evolution of magnetic field
variances with that of the mean square amplitudes of switchbacks, and quantify
the radial evolution of the cumulative counts of switchbacks per km. We find
that the amplitudes of switchbacks decrease faster than the overall turbulent
fluctuations, in a way consistent with the radial decrease of the mean
magnetic field. This could be the result of a saturation of amplitudes and may
be a signature of decay processes of large amplitude Alfvénic fluctuations in
the solar wind. We find that the evolution of switchback occurrence in the
solar wind is scale-dependent: the fraction of longer duration switchbacks
increases with radial distance whereas it decreases for shorter switchbacks.
This implies that switchback dynamics is a complex process involving both
decay and in-situ generation in the inner heliosphere. We confirm that
switchbacks can be generated by the expansion although other type of
switchbacks generated closer to the sun cannot be ruled out.

##  1\. Introduction

We have analyzed data from Parker Solar Probe (PSP), Helios and Ulysses to
investigate the evolution of switchbacks with heliocentric distance ranging
from about 0.1 to 3 au. Switchbacks are embedded within a continuous flux of
turbulent, Alfvénic fluctuations that permeates the solar wind and that
dominates the frequency range $f\simeq[10^{-6}-10^{-1}]$ Hz of the magnetic
field energy spectrum. Switchbacks make up the subset of largely Alfvénic
fluctuations whose amplitude is large enough that the magnetic field kinks
backwards on itself, leading to a local field polarity reversal and to a
corresponding radial velocity jet (Matteini et al., 2014). Although
switchbacks have been observed in the past, both by Helios at distances
$R\simeq 0.3-1$ au from the sun (Horbury et al., 2018), and by Ulysses beyond
1 au, one of the major findings of PSP is that this peculiar type of
fluctuation/structure appears to be ubiquitous and a prominent feature of the
solar wind closer to the sun, at distances $R\lesssim 0.2$ au. Understanding
what is the origin of switchbacks remains an important and open issue because
it may exclude or point to mechanisms taking place in the corona related to
coronal heating and solar wind acceleration. They may also contribute to
feeding the turbulent cascade via their nonlinear evolution.

Many mechanisms have been proposed to explain switchbacks. It has been
suggested that they are remnants of processes in the corona such as
interchange reconnection (Drake et al., 2020; Zank et al., 2020). In a high
beta plasma, they may also originate from instabilities like the firehose
(Tenerani & Velli, 2018), which is known to lead to highly kinked field lines.
Switchbacks may also form dynamically as turbulent fluctuations propagate into
the inner heliosphere, induced by shear flows or by the solar wind expansion
itself (Landi et al., 2006; Ruffolo et al., 2020; Squire et al., 2020; Mallet
et al., 2021; Schwadron & McComas, 2021).

Switchbacks have now been characterized extensively in their local properties
(Horbury et al., 2020; Woolley et al., 2020; Woodham et al., 2021; Laker et
al., 2021). They are three dimensional field-aligned structures displaying the
typical velocity-magnetic field correlation that characterizes Alfvén waves
propagating away from the sun. Just like the overall turbulent, Alfvénic
fluctuations, switchbacks are characterized by a nearly constant magnetic
field magnitude, a condition corresponding to spherical polarization
(Tsurutani et a., 2018; Matteini et al., 2015). A small level of
compressibility has been observed by PSP in some switchbacks (Farrell et al.,
2020; Larosa et al., 2021), especially at the boundaries, a property that may
be related to their evolutionary stage but that still remains to be
investigated.

Interestingly, large amplitude, Alfvénic fluctuations (in the sense of
velocity-magnetic field correlation) and switchbacks, with constant total
magnetic field magnitude, provide an exact, nonlinear dynamical state at the
magnetohydrodynamic scales (Goldstein et al., 1974; Barnes & Hollweg, 1974;
Hollweg, 1974). However, such a nonlinear state is thought to be unstable,
since large amplitude Alfvénic fluctuations tend to decay due to resonant
couplings with compressible modes via parametric and modulational
instabilities, proven to persist in various conditions, including expansion
effects and realistic temperature anisotropies. Theoretical investigations and
numerical simulations support the idea that Alfvénic fluctuations of coronal
origins should already begin to decay close to the sun, say, within tens of
solar radii (Tenerani and Velli, 2013; Réville et al., 2018; Shoda et al.,
2018). Recently, we have shown via numerical simulations that a localized
switchback is more stable than an Alfvénic wave of non constant total field
magnitude (like the one considered by, e.g., Landi et al. (2005)), or of a
periodic circularly or arc polarized shear Alfvén wave (Del Zanna, 2001), and
that it can propagate for distances up to a few tens of solar radii (Tenerani
et al., 2020). However, parametric decay eventually sets in, destroying the
switchback well before it can travel distances of the order of 1 au. This
suggests that in the absence of ongoing dynamical forcing – such as solar wind
expansion, for example – capable of continually replenishing switchbacks in-
situ, their occurrence rate should decrease with radial distance. Vice-versa,
the effect of a continuous forcing would lead to an increase in the occurrence
rate (or to a steady state) as we move further away from the sun.
Understanding how fluctuations and switchbacks evolve with radial distance is
relevant to understand not just how switchbacks can form out of the overall
turbulent fluctuations, but also to understand how they may contribute to the
evolution of the turbulent cascade in the solar wind.

In this paper, we study the radial evolution of switchbacks in the inner
heliosphere out to 3 au by combining data from PSP, Helios and Ulysses. We
track the radial evolution of the magnetic energy density carried by
switchbacks as compared to the overall turbulent fluctuations and we
investigate the occurrence rate of switchbacks. We finally discuss how our
results compare and fit into existing theories and observations of switchback
generation.

##  2\. Data sets and methods

In this work, we use magnetic field data from Helios 1 and 2 at 6 second
resolution (in SSE coordinates; Schwenn et al. (1975)) and consider a subset
of the intervals analyzed by Perrone et al. (2019), who classify data into
three fast Alfvénic streams crossed by Helios at various radial distances
during years 1975-1976. In this paper we will adopt their same notation to
denote intervals pertaining to stream A, stream B and stream C, respectively.

For Ulysses, magnetic field data from the VHM instrument are resampled at 6
seconds resolution (RTN; Balogh et al. (1992)). We consider the first and
third south polar passes, that occurred in 1994 and 2006-2007 denoted PP94 and
PP06, respectively, to compare data from different solar cycles.

Finally, we use high cadence magnetic field data from the FIELDS instrument
(in RTN coordinates; Bale et al. (2016)) onboard Parker Solar Probe (PSP)
resampled at 6 seconds during the first six orbits (E1 through E6), from year
2018 to year 2020. The same analysis was performed with 1-second resolution
and results were qualitatively similar. The PSP intervals have been identified
by visual inspection to avoid current sheet crossings (leading to north-south
main polarity reversals) and to select only those intervals with small values
of magnetic compressibility $\delta|{\bf B}|/|{\bf B}|<0.2$, on 5-minute
averages, although the majority of intervals had even lower values,
$\delta|{\bf B}|/|{\bf B}|\lesssim 0.1$. For the perihelia of E1, E2, E4 and
E6 we have also determined what intervals correspond to the same source region
via PFSS extrapolation (Panasenco et al., 2020), in order to compare
switchbacks counts within the same or nearby streams. Plasma data from Helios
(Schwenn et al., 1975), SWOOPS (Bame et al., 1992), and SPC (Case et al.,
2020) are also used to infer the mean radial speed of the wind for each time
interval.

Mean magnetic field, variances and mean amplitude of switchbacks, are obtained
via temporal averages over a sliding time window of length $\tau$ with
$\tau\in[30$ min $,12$ h], hence we consider time scales of fluctuations that
fall well into the inertial range and energy-containing range of Alfvénic
turbulence. For a fixed value of the time $\tau$ we calculate the mean
magnetic field $<{\bf B}>_{\tau}$ and the variances of the magnetic field
components $<\delta{B_{i}}^{2}>_{\tau}$. Switchbacks are identified as
follows. For each time scale $\tau$ we search for those intervals in which the
angle $\theta$ between the total magnetic field ${\bf B}$ and $<{\bf
B}>_{\tau}$ is larger than $\theta=90^{\circ}$. In this way we only select
fluctuations that correspond to a backwards kink of the magnetic field. By
using the direction of the total mean magnetic field rather than the radial
direction to identify switchbacks allows us to avoid including small radial
fluctuations in a non-dominantly radial mean magnetic field. The amplitude of
the switchbacks in the identified interval is defined as $(\delta
B_{i})_{\tau}^{2}=E[(B_{i}-<B_{i}>_{\tau})^{2}]$, where $E[\cdot]$ stands for
the mean, that in the case of switchbacks is evaluated inside the field
reversal. This procedure is also used to determine the cumulative counts of
switchbacks, essentially the number of switchbacks of duration less than
$\tau$. We adopt the same method to determine the duration of the field
reversal $\delta t$, and we use this information to investigate the evolution
of the probability distribution (pdf) of switchback duration with radial
distance. Finally, for each data set of $<\delta{B_{i}}^{2}>_{\tau}$ and
$(\delta B_{i})_{\tau}^{2}$ we then take the mean to obtain the radial trends
of such quantities. Count rates of switchbacks at given $\tau$ are obtained by
normalizing the number of switchbacks events in a given interval of duration
$\Delta t$ by the distance $\Delta R$ spanned in that interval by the
switchbacks, which is $\Delta R=\Delta t|V_{sw}-V_{sp}|$ with $V_{sp}$ the
average speed of the spacecraft, neglecting the contribution from the Alfvén
speed (we verified on PSP data that its inclusion does not change our
results).

##  3\. Results

###  3.1. Radial evolution of the energy density of fluctuations

We have analyzed the radial evolution of fluctuations amplitude from the
Helios data within each of the three fast streams, during the polar passes of
Ulysses, and by combining PSP data. Data from some Helios streams were also
used in the past by Villante & Vellante (1982) and Bavassano et al. (1982) to
infer the radial evolution of turbulence properties. In particular, they
showed that the rms amplitude of magnetic field fluctuations obey the WKB
prediction $\delta b_{rms}\simeq R^{-3/2}$ to a good approximation at the
lowest frequencies ($f\simeq 10^{-4}-10^{-3}$ Hz), whereas at higher
frequencies they decay faster with radial distance due to the ongoing
nonlinear cascade. Our results are consistent with such findings, and similar
trends are also found at PSP and Ulysses.

Figure 1 shows an example of the radial evolution of the variances of the
magnetic field for each component $<\delta B_{i}^{2}>$ and of the mean square
amplitude $(\delta B_{i})^{2}$ of the identified switchbacks for stream A of
Helios (in black color), PSP (in red) and Ulysses PP94 (in blue) for $\tau=8$
h. PSP data points from PSP (43 in total) have been binned in intervals of
width $\Delta R=0.05$ au, and the dashed lines represent reference slopes of
$\sim R^{-3}$ for the variances, and $\sim R^{-4}$ for the mean square
amplitude of switchbacks.The three panels from left to right show the radial
and two perpendicular components of the field.

By assuming a log-log scale linear relation of amplitudes with radial distance
$\log\delta B^{2}=A+\alpha\log R$, where $R$ is in units of au, we have
performed a best fit analysis of the variances of the magnetic field and of
the mean square amplitude of switchbacks at various values of $\tau$. Results
obtained from PSP E1-E6, Ulysses PP94 and from the average of the three Helios
streams A, B and C are displayed in Fig. 2 and Fig. 3. Fig. 2 shows that
variances follow the radial trend predicted by the WKB theory to a good
approximation, with $\alpha\lesssim-3$ for the tangential and normal
components (and $y,z$ in SSE coordinates) and $\alpha\gtrsim-3$ for the radial
component. While the coefficients $\alpha$ for the non-radial components are
consistent among the three different spacecraft, we notice that the radial
component of the fluctuations decays at different rates in Ulysses with
respect to PSP and Helios. In particular the difference between the obtained
angular coefficient $\alpha$ from Ulysses and Helios (or PSP) is larger than
the estimated error, pointing to the fact that at 1 au and beyond, other
mechanisms, such as nonlinearities and possibly interactions with large scale
shear flows, affect the evolution of fluctuations also at the larger scales.
In general, however, we recover the trends found by (Villante & Vellante,
1982; Bavassano et al., 1982), with the fluctuation energy at shorter
timescales decaying with radial distance faster than the larger timescales.
The coefficient $\alpha$ displays a monotonic dependence with $\tau$ and a
change of slope at around $\tau_{*}\simeq 10^{2}$ min can be seen in all the
three components, and at the three spacecraft. The timescale $\tau_{*}\simeq
10^{2}$ min roughly corresponds to the spectral break that separates the
inertial range from the $f^{-1}$ range of the turbulent spectrum. This
particular time scale depends on radial distance, however, we had to consider
a wide range of heliocentric distances in order to find the best fit for the
angular coefficient $\alpha$. As a consequence, we cannot capture clearly the
shift of $\tau_{*}$ with $R$. Fig. 3 shows the same quantities but for the
mean square amplitude of switchbacks only. Contrary to the overall variances,
the coefficient $\alpha$ does not depend strongly on $\tau$. The mean values
of $\alpha$ for the radial, tangential ($y$ for SSE) and normal ($z$ for SSE)
are listed in Table 1. For comparison, in the table we report also the mean
$\alpha$ values for the variances obtained for values $\tau<100$ min (left
values) and $\tau>100$ min (right values). As can be seen, the radial
component of switchbacks decays much faster than the overall turbulent
spectrum with a coefficient of about $\alpha\simeq-4$, while the non-radial
components have a decay rate $-4<\alpha<-3$ which is consistent with the decay
rate of transverse turbulent fluctuations. This suggests that switchbacks
behave as a separate population within the turbulent flux of Alfvénic
fluctuations in the sense that their energy density overall decays faster,
pointing to an ongoing continuous dynamical evolution of these structures.

Table 1Mean value of the power-law coefficient $\alpha$ for switchbacks (first three columns) and of the variances (last three columns). The mean value of $\alpha$ for the variances has been evaluated for $\tau<100$ min (left values) and $\tau>100$ min (right values). | $<\alpha_{r}>_{sb}$ | $<\alpha_{\bot,1}>_{sb}$ | $<\alpha_{\bot,2}>_{sb}$ | $<\alpha_{r}>$ | $<\alpha_{\bot,1}>$ | $<\alpha_{\bot,2}>$
---|---|---|---|---|---|---
PSP E1-E6 | $-3.6$ | $-3.1$ | $-3.1$ | $-3.2,-3.0$ | $-3.7,-3.4$ | $-3.5,-3.2$
Helios | $-3.8$ | $-3.4$ | $-3.4$ | $-3.2,-2.7$ | $-3.6,-3.1$ | $-3.6,-3.2$
Uly PP94 | $-3.9$ | $-3.6$ | $-3.0$ | $-3.5,-3.3$ | $-3.5,-3.2$ | $-3.5,-3.1$
Figure 1.— Variance (dots) and mean square amplitude of switchbacks
(triangles) as a function of radial distance for the radial (left panel), and
perpendicular (middle and right panels) components of the magnetic field
obtained for $\tau=8$ h. The dashed line shows the scaling $\sim R^{-3}$ and
the dot-dashed $R^{-4}$. Figure 2.— Angular coefficient $\alpha$ at different
$\tau$ obtained from the best fit of magnetic field variances as a function of
radial distance in log-log scale. The dashed horizontal lines mark the value
$\alpha=-3$. Figure 3.— Same as Fig. 2 but for the mean square amplitude of
switchbacks only. The dashed line in the top panel marks the value
$\alpha=-4$.

###  3.2. Cumulative counts of switchbacks/km and pdf

Figure 4 shows the cumulative counts of switchbacks per km as a function of
$R$ at different values of $\tau$. We have chosen three representative values,
$\tau=30$ min, 3 h and 8 h, to show how trends change as longer switchbacks
are considered. Data points from PSP (48 in total) have been binned inside
intervals of width $\Delta R=0.05$ au, and the error bars represent the range
of points within each bin. Count rates from stream A, B, C, and PP94, PP06 are
much less scattered than PSP and on average they lie systematically below PSP
rates as $\tau$ increases above 1 h. Count rates for $\tau=30$ min in both PSP
and Helios show a decreasing rate with radial distance, although PSP data
points are quite scattered. As timescale increases, and thus switchbacks of
longer duration are included, Helios displays a net radial increase in the
count rates while PSP rates have a more oscillatory behavior. At Ulysses,
count rates appear to somewhat decrease with radial distance regardless of the
scale considered. This suggests that switchbacks undergo a complex evolution,
where expansion or other in-situ effects are efficient in inducing magnetic
field folds at larger scales while others decay, especially the shorter ones.
We notice that such trends are clear particularly in Helios, where switchbacks
can be followed at several radial distances within a given stream. We argue
that the combined effect of the presence of decaying and forming switchbacks
and the mixing of different source regions and streams can explain the highly
scattered data points found with PSP. For this reason we have also calculated
the cumulative counts per km for E1, E2, E4 and E6 by considering only those
intervals that correspond to streams originating from nearby regions on the
sun (not shown here). However, due to the few points satisfying such
conditions, it is difficult to infer solid trends for the count rates in PSP.
We find indeed that some streams display an increasing rate of switchbacks
with radial distance while other a decreasing trend (not shown). Further away
from the sun, within the relatively steady fast streams observed by Ulysses,
effects of the magnetic spiral may also affect the count rates because the
non-radial components of the mean magnetic field decrease slower than
fluctuations amplitude. At Ulysses, the tangential mean magnetic field can
reach up to half the mean radial magnetic field. On average, the mean
tangential field forms an angle of about $20^{\circ}$ with the radial
direction. This deviation from a purely radial mean magnetic field may be
sufficient to affect the strength of the fluctuation amplitude relative to the
mean magnetic field, inhibiting the formation of switchbacks by expansion
effects (see discussion).

In order to gain insight on the evolution of switchbacks with radial distance,
in Fig. 5 we show the probability distribution of switchbacks as a function of
their duration $\delta t$. In this case we have chosen a relatively long value
of $\tau$, $\tau=8$ h, to include a wide range of time scales. Pdfs for larger
values of $\tau$ (not shown here) display the same qualitative features. Since
count rates for PSP are very scattered, a finer radial resolution of the pdf
provides a better representation of its radial evolution. For PSP data we
therefore show the pdfs for 6 mean radial distances. For Helios and Ulysses,
instead, we show the pdfs at the closest and furthest heliocentric distance
for stream A and PP94, respectively (other Helios streams and PP06 are
similar). Inspection of the pdfs shows that shorter duration switchbacks are
more common than the longer ones, in agreement with previous work (Dudok de
Wit et al., 2020). However, our analysis also shows that the fraction of
switchbacks of longer duration increases with radial distance, while the
opposite trend is typically seen at the shorter scales, for $\delta t\approx
10-100$ s.

Figure 4.— Cumulative counts of switchbacks per $km$ obtained for three
different values of $\tau$.

Figure 5.— Probability distribution function of switchback duration $\delta t$
for $\tau=8$ h.

##  4\. Discussion

A well known property of shear Alfvén waves in the expanding solar wind is
that, for frequencies higher than the expansion rate, they evolve according to
the conservation law of wave action (Bretherton & Garrett, 1969; Jacques,
1977). In the radially expanding solar wind this leads to fluctuation
amplitudes scaling as $\delta B_{i}^{2}\propto R^{-3}$. Such a scaling holds
true for any Alfvénic structure at arbitrary amplitude, involving also radial
fluctuations, provided the magnitude of the total magnetic field remains
constant, so that nonlinearities and coupling with compressible modes are
quenched. On the other hand, this condition implies that the tip of the
magnetic field is bound to rotate on the surface of a sphere of radius $B$,
where $B$ is the magnitude of the total magnetic field, thus imposing a
constraint on the maximum possible excursion of the magnetic field (Matteini
et al., 2018).

Our data analysis shows that at the larger scales, those roughly lying in the
$f^{-1}$ range of frequencies, the overall turbulent fluctuations follow to a
good approximation the WKB trend in all three components of the magnetic
field, thus, including the radial component. However, switchbacks behave as a
separate population in that the radial component of their fluctuations
decreases with radial distance as the mean magnetic field, $\delta
B_{r}^{2}\sim R^{-4}$. In other words, while the expansion can naturally drive
relative large amplitudes for a radial mean magnetic field ($\delta
B_{r}/<B_{r}>\propto R^{1/2}$) out from the turbulent bath of fluctuations,
some other physical processes must come into play to prevent relative
amplitudes from increasing without bound and maintaining the observed nearly
constant-$B$ constraint. Compressible effects leading to parametric
instabilities and wave steepening/collapse become dynamically important at
large amplitudes of about $10\%$ the mean magnetic field, or more. Coupling
with compressible modes may provide a possible way to control fluctuations
amplitudes whereas kinetic effects have been shown to play an important role
in controlling the level of plasma and magnetic compressibility (Tenerani and
Velli, 2013; González et al., 2021). Our finding of a faster decay of the
radial amplitudes may be an indication of such type of processes, although a
focused numerical study on this problem needs to be done. Recent simulations
in the expanding solar wind have indeed shown that switchbacks can form in-
situ (Squire et al., 2020; Shoda et al., 2021), however a study of how
switchbacks evolve once they are generated still needs to be carried out and
it will be of interest to compare numerical results with our observations.

Mozer et al. (2021) have performed a study of the occurrence rate of
switchbacks using data from encounters 3 through 7 of PSP. The authors adopt a
Poisson regression model to infer the dependence of the occurrence rate
(counts per hour) of switchbacks on solar wind speed and radial distance, and
conclude that the occurrence rate depends on the wind speed (higher rates for
higher wind speed) and does not depend on radial distance, thus, they exclude
in-situ generation mechanisms. However, our analysis, which includes the wind
speed in the switchbacks count normalization, shows that switchbacks counts
with PSP are extremely scattered and it is difficult to infer a dominant trend
with radial distance. The highly scattered occurrence of switchbacks is likely
due to the mixing of different streams, as also mentioned by Mozer et al.
(2021), but also to the presence of switchbacks that decay and reform in the
wind. Our analysis indeed also shows that the occurrence of switchbacks is
scale-dependent, a trend particularly clear in the cumulative counts in the
Helios streams, and also shown explicitly in the pdfs from PSP, Helios and
Ulysses data. Specifically, we have found that the fraction of switchbacks of
duration of a few tens of seconds and longer increases with radial distance
and that the fraction of those of duration below a few tens of seconds instead
decreases. Switchbacks in the solar wind can decay and reform in the expanding
solar wind, with in-situ generation being more efficient at the larger scales.

##  5\. Summary

We have analyzed magnetic field data from PSP, Helios and Ulysses covering
heliocentric distances $0.15\lesssim R\lesssim 3$ au at 6 second resolution.
We have determined the radial evolution of the amplitudes of turbulent
fluctuations by comparing the magnetic field variances with the mean square
amplitudes of switchbacks, and quantified the radial evolution of the
cumulative counts of switchbacks per km. Our results can be summarized as
follows: (1) radial amplitudes of switchbacks ($\delta B_{r}$) decrease faster
than that of the overall turbulent fluctuations. Specifically, their amplitude
square scales approximately as $\delta B_{r}^{2}\sim R^{-4}$ at all scales as
opposed to the bath of turbulent fluctuations scaling between $R^{-3}$ at the
larger scales (in agreement with the WKB prediction), and faster at the
shorter scales. Thus, switchbacks amplitudes saturate with respect to the
radial mean field; (2) the evolution of Alfvénic fluctuations at larger scales
(obeying WKB trends) in the expanding wind is a likely driver for in-situ
switchback generation; (3) the occurrence of switchbacks in the solar wind
depends on the duration of the switchbacks. Both in-situ formation and decay
are at play in the solar wind; (4) in-situ switchback generation is more
efficient at larger scales, however, we observe a net decrease of the
cumulative counts per km of switchbacks beyond 1 au at Ulysses, regardless of
their duration. We conclude that it is possible that switchbacks of two or
more type can coexist: those generated close to the sun that decay as they
propagate away, and those that are reformed or maintained alive in the inner
heliosphere by a combination of expansion effects (or other in-situ drivers)
and decay processes. At increasing distances $R\gtrsim 1$ au, effects of the
Parker spiral and possibly nonlinearities and large scale shears may be the
cause of the observed net decrease of switchbacks cumulative counts per km.

This research was supported by NASA grant #80NSSC18K1211 and the HERMES NASA
DRIVE Science Center Grant #80NSSC20K0604.

## References

  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., et al., Space Science Reviews, 204, 49 (2016)
  * Balogh et al. (1992) Balogh A., Beek T. J., Forsyth R. J., Hedgecock P. C., Marquedant R. J., Smith E. J., Southwood D. J., Tsurutani B. T., 1992 A&A, 92, 2.
  * Bame et al. (1992) Bame, S. J., McComas, D. J., Barraclough, B. L., Phillips, J. L., Sofaly, K. J., Chavez, J. C., 1992 A&A Supplement Series 92, 2.
  * Bavassano et al. (1982) Bavassano, B., Dobrowolny, M., Mariani, F., and Ness, N. F. 1982, Journal of Geophys. Research, 87, 3617.
  * Barnes & Hollweg (1974) Barnes, A. & Hollweg, J. V. J. Geophys. Res., 79, 2302–2318 (1974).
  * Bretherton & Garrett (1969) Bretherton, F. P. & Garrett, J. R. 1969, Proc. Roy. Soc. 302, 529.
  * Case et al. (2020) W. Case et al 2020 ApJS 246 43
  * Del Zanna (2001) Del Zanna L., Geophys. Res. Letters, 28, 13 (2001).
  * Drake et al. (2020) Drake, J. F., O. Agapitov, M. Swisdak, S. T. Badman, S. D. Bale, T. S. Horbury et al. 2020, A & A 650 A2.
  * Dudok de Wit et al. (2020) Dudok de Wit, T. et al. 2020 ApJS 246 39
  * Farrell et al. (2020) W. M. Farrell et al 2020 ApJS 249 28
  * Goldstein et al. (1974) Goldstein, M. L., Klimas, A. J., Barish, F. D., Solar wind three; Proceedings of the Third Conference, Pacific Grove, Calif., March 25-29, 1974. (A75-28001 12-92) Los Angeles, University of California, 1974, p. 385-387.
  * González et al. (2021) González, C. Tenerani, A., Matteini, L., Hellinger, P. and Velli, M., 2021 ApJL 914 L36
  * Hollweg (1974) Hollweg, J. V., JGR, 79, 10 (1974)
  * Horbury et al. (2018) Horbury, T.S., Matteini, L. and Stansby, 2018 D. Monthly Notices of the Royal Astronomical Society, 478, 1980–1986
  * Horbury et al. (2019) Horbury, T. et al. 2019, The Astrophysical Journal Supplement Series, 246:45, 2020
  * Horbury et al. (2020) Timothy S. Horbury et al 2020 ApJS 246 45
  * Jacques (1977) Jacques, S. A. 1977, The Astrophysical Journal, 215, 942.
  * Landi et al. (2005) Landi, S., Hellinger, P. and Velli, M., Proceedings of the Solar Wind 11 / SOHO 16, “Connecting Sun and Heliosphere” Conference (ESA SP-592). 12 - 17 June 2005 Whistler, Canada. Editors: B. Fleck, T.H. Zurbuchen, H. Lacoste, p. 785.
  * Landi et al. (2006) Landi, S., Hellinger, P. and Velli, M., 2006 Geophys. Res. Lett. 3, L14101
  * Larosa et al. (2021) Larosa, A et al, 2021 A&A 650, A3
  * Laker et al. (2021) Laker, R et al. 2021, A&A 650 A1
  * Landi et al. (2005) Landi, S., Hellinger, P. and Velli, M., 2005 Proc. Solar Wind 11
  * Malara and Velli (1996) Malara, F., and Velli, M., 1996 Phys. Plasmas, 3 4427
  * Mallet et al. (2021) Mallet et al. 2021 arXiv:2104.08321 [physics.space-ph]
  * Matteini et al. (2014) Matteini, L., Horbury, T. S., Neugebauer, M., and Goldstein, B. E., 2014 Geophys. Res. Lett., 41, 259.
  * Matteini et al. (2015) Matteini, L. et al. 2015, ApJ 802, 1
  * Matteini et al. (2018) L. Matteini et al 2018 ApJL 869 L32.
  * Mozer et al. (2021) Mozer F., et al. arXiv:2105.07601
  * Primavera et al. (2019) Primavera L., Malara, F., Servidio, S, and Nigro, G., 2019 ApJ 880 156
  * Réville et al. (2018) Réville V., Tenerani A., and Velli M., 2018, ApJ 866, 38
  * Panasenco et al. (2020) ,Panasenco, O., Velli, M., D’Amicis, R., Shi, C., . Réville, V. , Bale, S. D. , Badman, S. T., Kasper, J., Korreck, K., Bonnell, J. W., Dudok de Wit, T.,Goetz, K., Harvey, P. R., MacDowall, R. J., Malaspina, D. M., Pulupa, M., Case, A. W., Larson, D.,Livi, R.,Stevens, M., Whittlesey, P., 2020, APJ SS 246, 54 doi = 10.3847/1538-4365/ab61f4,
  * Perrone et al. (2019) Perrone, D., Stansby, D., Horbury, T. S. and Matteini, L. MNRAS 483, 3730–3737 (2019)
  * Ruffolo et al. (2020) Ruffolo, D., Matthaeus, W. H, Chhibber R. et al. The Astrophysical Journal, 902:94, 2020
  * Scearce et al. (1975) Scearce, C., S. Cantarano, N. F. Ness, F. Mariani, R. Terenzi, and L. F. Burlaga, 1975 Rome-GSFC magnetic field experiment for HELIOS A and B, NASA/GSFC Tech. Rep..X-692-75-112, Goddard Space Flight Center, Greenbelt, Md.,
  * Schwadron & McComas (2021) N. A. Schwadron and D. J. McComas 2021 ApJ 909 95
  * Schwenn et al. (1975) Schwenn, R., Rosenbauer, H., Miggenrieder, H.: 1975, The plasma experiment on board Helios E1. Raumfahrtforschung 19, 226.
  * Shoda et al. (2018) Shoda, M. ,Yokoyama, T. and Suzuki, K., 2018, ApJ 860 17
  * Shoda et al. (2021) Shoda, M. et al, 2021, arXiv:2101.09529
  * Squire et al. (2020) Squire, J., Chandran, B. D. G., and Meyrand, R., ApJL 891 L2 2020
  * Tenerani and Velli (2013) Tenerani, A. and Velli, M., 2013, J. Geophys. Res. 118, 7507-7516
  * Tenerani & Velli (2018) Tenerani A. & Velli, M., 2018 ApJ 867 (2), L26
  * Tenerani et al. (2020) Tenerani, A. et al. 2020, ApJS 246 32
  * Tsurutani et a. (2018) Tsurutani, B T. et al., 2018, JGR, 123 4.
  * Villante & Vellante (1982) Villante & Vellante, M 1982, Solar Physics, 81, 2.
  * Woodham et al. (2021) Woodham, L. D., et al. 2021, A&A 650, L1
  * Woolley et al. (2020) Woolley, T. et al. 2020, MNRAS, 498, 4
  * Zank et al. (2020) G. P. Zank et al. 2020 ApJ 903 1