* [72] M. J. Wainwright, High-dimensional statistics: A non-asymptotic viewpoint, vol. 48.  Cambridge University Press, 2019.


# Bayesian Solar Wind Modeling with Pulsar Timing Arrays

Jeffrey S. Hazboun NANOGrav Physics Frontiers Center Postdoctoral Fellow
Physical Sciences Division, University of Washington Bothell, 18115 Campus Way
NE, Bothell, WA 98011, USA  Joseph Simon Department of Astrophysical and
Planetary Sciences, University of Colorado, Boulder, CO 80309, USA  Dustin R.
Madison Department of Physics, University of the Pacific, 3601 Pacific Avenue,
Stockton, CA 95211, USA  Zaven Arzoumanian  X-Ray Astrophysics Laboratory,
NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA  Kathryn
Crowter Department of Physics and Astronomy, University of British Columbia,
6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada  Megan E. DeCesar George
Mason University, resident at the Naval Research Laboratory, Washington, DC
20375, USA  Paul B. Demorest National Radio Astronomy Observatory, 1003
Lopezville Rd., Socorro, NM 87801, USA  Timothy Dolch Department of Physics,
Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA  Eureka
Scientific, Inc. 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017  Justin
A. Ellis  Infinia ML, 202 Rigsbee Avenue, Durham NC, 27701  Robert D. Ferdman
School of Chemistry, University of East Anglia, Norwich, NR4 7TJ, United
Kingdom  Elizabeth C. Ferrara Department of Astronomy, University of Maryland,
College Park, MD, 20742, USA  Center for Exploration and Space Studies
(CRESST), NASA/GSFC, Greenbelt, MD 20771, USA  NASA Goddard Space Flight
Center, Greenbelt, MD 20771, USA  Emmanuel Fonseca Department of Physics and
Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA
Center for Gravitational Waves and Cosmology, West Virginia University,
Chestnut Ridge Research Building, Morgantown, WV 26505, USA  Peter A. Gentile
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315,
Morgantown, WV 26506, USA  Center for Gravitational Waves and Cosmology, West
Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505,
USA  Glenn Jones  Rigetti Computing, Inc., 775 Heinz Ave. Berkeley, CA 94710,
USA  Columbia Astrophysics Laboratory, Columbia University, NY 10027, USA
Megan L. Jones Center for Gravitation, Cosmology and Astrophysics, Department
of Physics, University of Wisconsin-Milwaukee,
P.O. Box 413, Milwaukee, WI 53201, USA  Michael T. Lam School of Physics and
Astronomy, Rochester Institute of Technology, Rochester, NY 14623, USA
Laboratory for Multiwavelength Astrophysics, Rochester Institute of
Technology, Rochester, NY 14623, USA  Lina Levin Jodrell Bank Centre for
Astrophysics, School of Physics and Astronomy, The University of Manchester,
Manchester M13 9PL, UK  Duncan R. Lorimer Department of Physics and Astronomy,
West Virginia University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for
Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge
Research Building, Morgantown, WV 26505, USA  Ryan S. Lynch Green Bank
Observatory, P.O. Box 2, Green Bank, WV 24944, USA  Maura A. McLaughlin
Department of Physics and Astronomy, West Virginia University, P.O. Box 6315,
Morgantown, WV 26506, USA  Center for Gravitational Waves and Cosmology, West
Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505,
USA  Cherry Ng Dunlap Institute for Astronomy and Astrophysics, University of
Toronto, 50 St. George St., Toronto, ON M5S 3H4, Canada  David J. Nice
Department of Physics, Lafayette College, Easton, PA 18042, USA  Timothy T.
Pennucci Institute of Physics, Eötvös Loránd University, Pázmány P. s. 1/A,
1117 Budapest, Hungary  Scott M. Ransom National Radio Astronomy Observatory,
520 Edgemont Road, Charlottesville, VA 22903, USA  Paul S. Ray U.S. Naval
Research Laboratory, Washington, DC 20375, USA  Renée Spiewak Jodrell Bank
Centre for Astrophysics, School of Physics and Astronomy, The University of
Manchester, Manchester M13 9PL, UK  Ingrid H. Stairs Department of Physics and
Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver,
BC V6T 1Z1, Canada  Kevin Stovall Department of Physics and Astronomy,
University of New Mexico, 210 Yale Blvd NE, Albuquerque, NM 87106, USA  Joseph
K. Swiggum NANOGrav Physics Frontiers Center Postdoctoral Fellow  Department
of Physics, Lafayette College, Easton, PA 18042, USA  Weiwei Zhu National
Astronomical Observatories, Chinese Academy of Science, 20A Datun Road,
Chaoyang District, Beijing 100012, China

###### Abstract

Using Bayesian analyses we study the solar electron density with the NANOGrav
11-year pulsar timing array (PTA) dataset. Our model of the solar wind is
incorporated into a global fit starting from pulse times-of-arrival. We
introduce new tools developed for this global fit, including analytic
expressions for solar electron column densities and open source models for the
solar wind that port into existing PTA software. We perform an ab initio
recovery of various solar wind model parameters. We then demonstrate the
richness of information about the solar electron density, $n_{E}$, that can be
gleaned from PTA data, including higher order corrections to the simple
$1/r^{2}$ model associated with a free-streaming wind (which are informative
probes of coronal acceleration physics), quarterly binned measurements of
$n_{E}$ and a continuous time-varying model for $n_{E}$ spanning approximately
one solar cycle period. Finally, we discuss the importance of our model for
chromatic noise mitigation in gravitational-wave analyses of pulsar timing
data and the potential of developing synergies between sophisticated PTA solar
electron density models and those developed by the solar physics community.

solar wind, pulsar timing, pulsar timing arrays, dispersion, solar electron
density

††journal: ApJ††facilities: AO, GBO††software: Numpy (Harris et al., 2020),
Scipy (Virtanen et al., 2020), MatPlotLib (Hunter, 2007), Pint (Luo et al.,
2021), Enterprise (Ellis et al., 2019), enterprise_extensions (Taylor et al.,
2021), Tempo2 (Hobbs & Edwards, 2012), PTMCMCSampler (Ellis & van Haasteren,
2017)

##  1 Introduction

Radio observations of distant astrophysical sources have long been used to
study the content and characteristics of the solar wind. For example, Hewish &
Dennison (1967) used observations of scintillation from quasars as a way to
probe the structure of the solar wind. The dispersion of pulsed radio emission
from pulsars, due to the diffuse ionized medium along the observational line
of sight (LOS), has been known from their first radio observations (Hewish et
al., 1968). The utility of these radio pulses for investigations of electron
density was used in Counselman et al. (1970), only two years after the
discovery of pulsars, to measure the solar electron density ($10\;{\rm
cm}^{-3}$ at $1{\rm au}$) based on a spherically symmetric model for the solar
wind.

Since these early investigations, pulsar astronomers have often included a
model for the solar electron density as a part of pulsar ephemerides—see, for
example, the three main pulsar timing software packages Tempo, Tempo2, and
Pint (Nice et al., 2015; Hobbs & Edwards, 2012; Luo et al., 2021). In Lommen
et al. (2006) and Splaver et al. (2005), the solar wind signal was shown to be
highly covariant with astrometric components of the timing models for
individual pulsars, in particular the parallax and sky position, because these
parameters contain strong Fourier components at $1/{\rm yr}$ and higher
harmonics. Beyond the importance of a solar wind model for accuracy in pulsar
astronomy, You et al. (2007) and You et al. (2012) showed the sensitivity of
pulsar data sets to more complex (than the spherically symmetric $1/r^{2}$
wind) features in the solar electron density111We present more detail on these
models in §2. More recent work (Madison et al., 2019; Tiburzi et al., 2019,
2021) has shown the potential of pulsar timing arrays (PTAs) as independent
probes of the solar wind and its behavior as a function of time.

In this paper we introduce new methods, expanding upon those used in Madison
et al. (2019) and Tiburzi et al. (2021), to obtain information about the solar
electron density from PTA data. Using a fully Bayesian framework, we show that
much more information about the solar wind can be obtained from the same set
of pulsar data used in Madison et al. (2019), the NANOGrav 11-year Data Set
Arzoumanian et al. (2018a, henceforth referred to as NG11). As in these recent
studies, it is important that we are using an array of pulsars. This allows us
to separate the variations in electron density of the ionized interstellar
medium (ISM) along the kiloparsec distances to pulsars from the local
fluctuations in electron density due to variations in the solar wind.

The NANOGrav 11-year Data Set consists of high-precision time-of-arrival (TOA)
measurements 45 millisecond pulsars spanning up to 11.4 years. Each pulsar was
observed at approximately a monthly cadence, over widely separated radio
frequencies in order to measure both TOAs and dispersion of the pulsar signal
at each observing epoch. Observations were made using the Arecibo Observatory
and the Green Bank Observatory. Typical TOA measurement precision was in the
range $0.1$ to 1 $\mu$s. Further details of the data set are in NG11.

###  1.1 Monitoring of Solar Electron Density

Since the initial use of astrophysical radio sources to measure solar wind
density via scintillation (Hewish & Dennison, 1967), a program for monitoring
solar wind densities, especially at higher solar latitudes, has continued
(Coles, 1978; Manoharan, 2010; Tokumaru, 2013). In addition to ground-based
monitoring, a considerable number of resources have been used to study the
solar wind from space. Many of these spacecraft have included instruments for
electron density measurements, including Ulysses’s Solar Wind Observations
Over the Poles of the Sun (SWOOPS; Bame et al., 1992), the Orbiter Retarding
Potential Analyzer (ORPA) instrument mounted on the NASA Pioneer Venus Orbiter
spacecraft (Knudsen et al., 1980) and the Advanced Composition Explorer (ACE)
with its Solar Wind Electron Proton Alpha Monitor (McComas et al., 1998).
Similar to pulsar timing measurements, the Viking mission used dual frequency
delays to make early measurements of the latitudinal dependence of the solar
electron density (Muhleman & Anderson, 1981). The Parker Solar Probe (Bale et
al., 2016) with its Solar Wind Electrons Alphas and Protons (SWEAP) instrument
(Kasper et al., 2016) is revolutionizing our understanding of the solar wind
in the inner solar system (Bale et al., 2019).

Apart from the Viking data, these space missions take in situ measurements of
the solar electron density. These allow for one variety of long-term
monitoring of the solar wind and Ulysses data have been used to study its
structure and behavior at high to mid latitudes (Issautier et al., 2001).
Obviously, this monitoring requires the use of costly spacecraft, and while
the data from pulsar timing does not provide the same type of fine-grained
local spatial information, the many lines of sight to pulsars allow for
omnidirectional integral monitoring of the solar electron density from ground-
based facilities and the scrutinization of various models for its structure
and evolution. PTAs provide a probe of the solar wind that is distinct from
and complementary to space missions.

###  1.2 Noise Mitigation in PTAs

The main goal of PTA experiments is to use precise long-term measurements of
millisecond pulsars to observe gravitational-waves (GWs) in the nanohertz
regime. Many pulsars are used since the unique noise properties of individual
pulsars necessitates the corroboration of common222“Common” is used here in
the sense that the signal is present, in whole or in part, in observations of
all pulsars in the timing array. astrophysical signals across multiple sources
and because lower-amplitude common signals can be drawn out of the noise of
more pulsars. Additionally, PTAs allow us to search for other common signals,
for instance the motion of the solar system barycenter (Vallisneri et al.,
2020) or errors in terrestrial time standards (Hobbs et al., 2020).

One motivation for this work is to develop a pan-PTA solar wind model as one
component of the next generation of PTA noise models. Dispersion measure (DM)
variations are one of the largest noise sources in PTA data, (Lam et al.,
2016; Cordes et al., 2016; Cordes & Shannon, 2010; Jones et al., 2017), and
modeling these variations is an important part of PTA noise mitigation
strategies—it is, in fact, a major driver of PTA observing strategies (Ransom
et al., 2019; Lam et al., 2018a). DM values in NG11 vary from $\sim
3-300\;{\rm pc}/{\rm cm}^{3}$, while the variations are usually on the order
of $10^{-3}\;{\rm pc}/{\rm cm}^{3}$. While a piecewise binned model of DM
variations (DMX, see Arzoumanian et al., 2015) has served NANOGrav well, it
removes up to a third of the power through the timing model transmission
function (Hazboun et al., 2019; Arzoumanian et al., 2015; Jones et al., 2017),
a problem endemic to any DM variation model that is part of the timing model
fit. Additionally it has been shown that asynchronous multi-band observations
can lead to misestimation of DM variations (Lam et al., 2015) when binned
together, either in DMX or an interpolation basis. Niu et al. (2017)
demonstrated that observations only a day apart can be insufficient to remove
sharp solar wind effects, or “cusps”. In Hazboun et al. (2020) and Lam et al.
(2018c) the effect of sharp unmodeled cusps in DM variations was demonstrated.
Mismodeling of solar wind cusps could then also present as broadband white
noise in pulsar datasets. In particular, this could adversely affect GW
searches for single sources which are especially dependent on the high
frequency noise floor of PTAs (Lam, 2018; Lam et al., 2018b; Lam & Hazboun,
2021).

Perhaps most importantly for PTAs, in addition to the noise introduced by
mismodeling DM variations, it was shown in Tiburzi et al. (2016) that the
solar wind can actually manifest spatial correlations amongst the pulsars, a
potential source of confusion when searching for GWs, especially a GW
background which presents in the data as a low-frequency spatially-correlated
process. In this paper, we demonstrate how a fully Bayesian framework, built
upon the extensive analysis infrastructure developed by PTA collaborations to
search for GWs, can isolate the signal from the solar wind in PTA data. This
allows for precise noise mitigation as well as monitoring of the overall solar
electron density both as a function of time and observational LOS.

The paper is organized as follows. In §2 the DM parameter will be derived
along with various analytical models for a spherically symmetric and static
solar electron density. In §3 we present the full Bayesian model for pulsar
timing data and how we incorporate our solar wind models. In §4 three models
for the behavior and structure of the solar electron density will be presented
along with the results of using these models to analyze NG11. In §5 we discuss
the incorporation of these models into PTA data analyses as well as how these
models might be used in the future to bolster efforts to study solar physics.

##  2 Solar Electron Density Models

The simplest and most common model for the solar electron density used in
pulsar timing software packages like Tempo, Tempo2 and Pint (Nice et al.,
2015; Hobbs & Edwards, 2012; Luo et al., 2021) is the spherically-symmetric,
time-independent expression $n_{e}(\vec{r})=n_{E}\left(1\,{\rm
au}/r\right)^{2}$, where $n_{E}$ is the electron density at $1\,{\rm au}$,
herein measured using units of $1/{\rm cm}^{3}$ to match pulsar timing
software.

Figure 1: Solar Wind Geometry: The figure defines the variables used for the
solar wind electron density line integral. Here we have adopted some of the
nomenclature of scattering calculations, as in Aksim et al. (2019), using $b$
as the impact parameter and defining $\theta_{i}$ as the “impact angle”, i.e
the solar elongation. The vector $\vec{p}$ points from the Sun to the pulsar,
while $\vec{r}_{\oplus}$ points from the Sun to the observatory on Earth. The
distance to the pulsar is effectively infinite when compared to the distance
to the Earth, so $z_{p}$ and $\vec{p}$ are effectively parallel.

The DM parameter used in pulsar timing is the column density, i.e. the
integral of $n_{e}(\vec{r})$ along the line-of-sight to the radio source

$\displaystyle{\rm DM}_{\odot}$ $\displaystyle=$
$\displaystyle\int_{r_{\oplus}}^{r_{p}}n_{e}\left(\vec{r}\right)ds$ (1)
$\displaystyle=$
$\displaystyle\int_{r_{\oplus}}^{r_{p}}n_{E}\left(\frac{1\,{\rm
au}}{r}\right)^{2}ds$ (2) $\displaystyle=$ $\displaystyle n_{E}\left(1\,{\rm
au}\right)^{2}\frac{1}{b}\left[\frac{\pi}{2}+\tan^{-1}\left(\frac{z_{\oplus}}{b}\right)\right]$
(3) $\displaystyle=$ $\displaystyle n_{E}\left(1\,{\rm
au}\right)^{2}\frac{\pi-\theta_{i}}{r_{\oplus}\sin\theta_{i}},$ (4)

where the first line is the generic expression, and the last equality uses the
relationship
$\frac{z_{\oplus}}{b}=\cot(\theta_{i})=\tan\left(\frac{\pi}{2}-\theta_{i}\right)$
and $b=r_{\oplus}\sin\theta_{i}$. See Figure 1 for the definitions of these
variables333Note that the $\rho$ used in Edwards et al. (2006, (and so
forth)); Tiburzi et al. (2019, (and so forth)) and $\theta_{i}$ are
supplementary angles, i.e., $\rho=\pi-\theta_{i}$.. The time delay is then
dependent on the radio frequency and the DM by

$\Delta t=\frac{e^{2}}{2\pi m_{e}c^{2}}\frac{\rm DM}{\nu^{2}}\,.$ (5)

DM is usually measured in units of ${\rm pc}/{\rm cm}^{3}$, and the constants
of nature in front combine to $4.15\times 10^{3}\,{\rm MHz}^{2}\,{\rm
pc}^{-1}\,{\rm cm}^{3}\,{\rm s}$.

In You et al. (2007) the authors use a two phase model assuming that the LOS
to the pulsar crosses the fast and slow solar wind. Higher order terms, in
addition to the simple $1/r^{2}$ model are used to ameliorate the effect of
realistic solar wind time delays. The slow solar wind model444Here we have
converted the expressions in You et al. (2007) to units more familiar to
pulsar timing astronomers.,

$\displaystyle n_{\odot}\times{\rm cm}^{3}$ $\displaystyle=$ $\displaystyle
8.867\left(\frac{1\,{\rm au}}{r}\right)^{2}+1.186\left(\frac{1\,{\rm
au}}{r}\right)^{2.7}$ (6) $\displaystyle+1.517\times
10^{-6}\left(\frac{1\,{\rm au}}{r}\right)^{6}$ $\displaystyle+1.431\times
10^{-29}\left(\frac{1\,{\rm au}}{r}\right)^{16}$

was developed in You et al. (2007) by combining the near and far distance
models in Allen (1947) and Muhleman & Anderson (1981). A separate model was
used for the fast solar wind

$\displaystyle n_{\odot}\times{\rm cm}^{3}$ $\displaystyle=$ $\displaystyle
2.498\left(\frac{1\,{\rm au}}{r}\right)^{2}$ (7) $\displaystyle+1.860\times
10^{-4}\left(\frac{1\,{\rm au}}{r}\right)^{4.39}$ $\displaystyle+4.067\times
10^{-30}\left(\frac{1\,{\rm au}}{r}\right)^{16.25}$

developed in Guhathakurta & Fisher (1998, 1995). These models were then
combined with measurements of the portion of the line-of-sight in the fast and
slow wind to better estimate the cusps to the pulsar timing data. In You et
al. (2012) the fraction parameter for the two phases of the model is fit for
as part of the analysis. Additionally a time dependent component is
investigated. The model in You et al. (2007, 2012) has potential as a useful
tool for removing the solar wind signal from pulsar timing data, especially
for single pulsars, where disentangling the DM variations from the ISM is
difficult.

Recently, however, Tiburzi et al. (2019, 2021) showed that the simple
$1/r^{2}$ model for $n_{E}$ commonly used in pulsar timing software packages
and the You et al. (2012) model are both insufficient for accurately modeling
dispersion delays from the solar wind. Even with allowances for time
variability, these models still fall short of the mitigation needed for the
low frequency radio data taken by the LOw-Frequency ARray (LOFAR) (Stappers et
al., 2011; van Haarlem et al., 2013) used in the analyses.

In Aksim et al. (2019) the authors compare their measurements of DM using Very
Long Baseline Interferometry (VLBI) observations of quasars to both a
spherically symmetric model and a numerical integral of the Alfvén Wave Solar
Model (AWSoM) (van der Holst et al., 2014). As in other studies the AWSoM
matches well with the astrophysical measurements, however this model is
currently too cumbersome to be used effectively by PTAs. Aksim et al. (2019)
also fit for a variable power spherically symmetric model and find reasonable
agreement with their time delay measurements.

###  2.1 General Spherically Symmetric Solar Wind Density

Here we derive generic expressions for the dispersion measure for any
spherically symmetric electron density fall off. While these expressions will
assume the source of the electrons is the sun, they are general enough to be
used for any spherically symmetric source of streaming electrons along the LOS
to a pulsar, a stellar wind from a star orbiting a pulsar, for instance.
Relaxing the assumption that the electron density around the sun drops off as
$1/r^{2}$, we model the drop off with a more generic power law dependence.
These calculations are similar to those in Aksim et al. (2019), where the
analytical expression for the time delay was calculated directly.

Assuming that $p$ is the exponent in the power law dependence we go through a
similar calculation as in Equation 4,

$\displaystyle\text{DM}_{\odot}^{p}$ $\displaystyle=$
$\displaystyle\int_{\text{Earth}}^{\text{Pulsar}}n\left(\vec{r}\right)ds=\int_{r_{\oplus}}^{r_{p}}n^{(p)}_{E}\left(\frac{1\,{\rm
au}}{r}\right)^{p}ds$ (8) $\displaystyle=$
$\displaystyle\int_{-z_{\oplus}}^{z_{p}}n^{(p)}_{E}\left(\frac{1\,{\rm
au}}{\sqrt{b^{2}+z^{2}}}\right)^{p}dz$ (9) $\displaystyle=$ $\displaystyle
n^{(p)}_{E}\left(\frac{1\,{\rm
au}}{b}\right)^{p}\left[z\times\,_{2}\mathcal{F}_{1}\\!\left(\frac{1}{2},\frac{p}{2},\frac{3}{2},\text{-}\frac{z^{2}}{b^{2}}\right)\right]_{\text{-}z_{\oplus}}^{z_{p}}$
(10)

where ${}_{2}\mathcal{F}_{1}$ is a hypergeometric function. This expression
can be further simplified by assuming the pulsar is sufficiently distant that
we can take the limit $z_{p}\rightarrow\infty$. This gives

$\text{DM}_{\odot}^{p}=n^{(p)}_{E}\left(\frac{1\,{\rm
au}}{b}\right)^{p}\left(z_{\oplus}\times\,{}_{2}\mathcal{F}_{1}\\!\left(\frac{1}{2},\frac{p}{2},\frac{3}{2},\text{-}\frac{z_{\oplus}^{2}}{b^{2}}\right)+\frac{b\sqrt{\pi}}{2}\frac{\Gamma\left(\frac{p}{2}-\frac{1}{2}\right)}{\Gamma\left(\frac{p}{2}\right)}\right)\;,$
(11)

where $\Gamma\left(x\right)$ is the Gamma function. Compare this expression to
the time delay equation given in Aksim et al. (2019). Equation 11 simplifies
to Equation 4 when $p=2$. See Appendix A for expressions with higher integer
values of $p$. The relationships below Equation 4 can be used to write this in
terms of $\theta_{i}$. One can choose to build a solar electron density model
using multiple summed terms, as in You et al. (2007, 2012); Allen (1947);
Muhleman & Anderson (1981); Guhathakurta & Fisher (1998, 1995), or allow the
index to vary, as in Aksim et al. (2019).

##  3 Bayesian Methods

The solar wind modeling presented here relies on the analysis infrastructure
developed by PTAs for GW searches (van Haasteren & Levin, 2013; Lentati et
al., 2016; Taylor et al., 2013; van Haasteren & Vallisneri, 2014; van
Haasteren et al., 2009). While a gravitational wave background only causes
residuals of tens of nanoseconds, the solar wind signal is much stronger
($\sim 1.7\mu{\rm s}$ at 1 GHz with $\theta_{i}=10^{\circ}$). This allows us
to undertake an in depth study of the solar wind by modeling the solar
electron density across all of the pulsars simultaneously.

To decouple the variations within the ISM along the LOS to each of the pulsars
from the solar electron density that is varying locally, we use a
deterministic solar wind model based on the expressions in §2.1. The solar
wind electron density is fit as a global parameter across all pulsars, while
the geometry of the solar wind model is dependent on the details of the
individual Sun-Earth-Pulsar angles, $\theta_{i}$, shown in Figure 1. In
addition to this model each pulsar is also fit with a quadratic polynomial in
time over the data set for the DM, encoded with the pulsar timing parameters,
DM1 and DM2, as well as an additional Gaussian process model that emulates the
variations of the ISM, which will be discussed in §3.2.

The PTA analysis framework includes a timing model marginalization, based on a
linearized timing model (van Haasteren & Levin, 2013) and a full noise
treatment that includes a generalized covariance matrix, parametrized by TOA
errors, three additional white noise parameters and a power law red noise
model (Lentati et al., 2016; Taylor et al., 2013). The PTA likelihood model is
implemented using the Enterprise software suite (Ellis et al., 2019) and the
various models and extensions compiled in enterprise_extensions (Taylor et
al., 2021). We use the Parallel Tempering Markov-Chain Monte Carlo Sampler,
PTMCMCSampler, for numerically integrating our likelihood (Ellis & van
Haasteren, 2017).

It should be noted that while this formalism is focused on high precision
millisecond pulsars, the code infrastructure can be used on data from any type
of pulsar and its timing data. Observations of the more numerous canonical
pulsars may be useful for these types of solar investigations in the future.

###  3.1 ACE SWEPAM Prior

As mentioned in §1, a number of past and ongoing space missions have collected
solar electron density data over the past few decades. In order to take
advantage of this large base of knowledge about the solar electron density we
used the solar electron density data from the ACE Solar Wind Electron Proton
Alpha Monitor (McComas et al., 1998) to build an informative prior for
$n_{E}$, the solar electron density at $1\,{\rm au}$. We binned the data from
the same time span as NG11, corrected the density from the L1 Lagrange point
using a $1/r^{2}$ density model555$n_{e}^{\rm ACE}=n_{E}\left(\frac{1\;{\rm
au}}{0.997\,{\rm au}}\right)^{2}$ and used an empirical distribution to build
a random variable in SciPy. The distribution and median are shown in Figure 2.

Figure 2: ACE SWEPAM Prior. We use in situ measurements of the solar electron
density at Earth’s orbit to construct a prior for our analysis. Note that the
prior extends all the way up to $\sim 50{\rm\,cm}^{-3}$.

The ACE SWEPAM in situ data are a useful prior for limiting the extent of
parameter space at which we expect to be working, roughly
$0.01-50{\rm\,cm}^{-3}$. However, the electron flow in the ecliptic plane is
known to be different than the fast flow at other solar latitudes (Muhleman &
Anderson, 1981; Bame et al., 1992; Guhathakurta & Fisher, 1998). We do not
expect the ACE measurements to match our results perfectly, but it is a good
comparison in these analyses as a sanity check. In practice the prior was only
informative when compared to the data in some individual pulsar runs — usually
pulsars without strong solar wind signals in their TOAs. In tests, the full
PTA runs retrieved the same posteriors from uninformative uniform priors, but
we use the ACE prior throughout the work presented here to decrease
convergence times in our Markov Chain Monte Carlo runs.

Uniform priors were used for the spectral index in all higher order searches
and log-uniform priors were used for $n^{(p)}_{E}$ priors when $p\neq 2$.

###  3.2 Dispersion Measure Variation Model

Gaussian processes (GPs) have been used extensively to model time-correlated
noise in pulsar datasets (Arzoumanian et al., 2020, 2018b; Goncharov et al.,
2021). Commonly, they are used to model the GW background, achromatic red
noise, such as that caused by intrinsic spin noise, and chromatic noise due to
the bulk movement of the ISM (Goncharov et al., 2021; Lentati et al., 2016).

The standard prior function of GPs in these datasets is a red (negative
spectral-indexed) power law model with a Fourier-basis. The solar wind has
strong spectral characteristics at Fourier frequencies of $1/{\rm yr}$ and
higher modes (Madison et al., 2019). In early testing these modes showed
themselves to be highly covariant with some frequencies in a Fourier-basis GP,
so we have instead adopted for our DM variations model another common ansatz
for GPs: a square exponential kernel,

$\displaystyle k_{\mathrm{SE}}(t_{1},t_{2})=$
$\displaystyle\sigma^{2}\exp{\left(-\frac{|t_{1}-t_{2}|^{2}}{2\ell^{2}}\right)}+\left(\frac{\sigma}{500}\right)^{2},$
(12)

where $\sigma$ is an overall variance, $\ell$ is a time scale of variation,
and a constant variance is also added as a regularizing term for stable
inversions and to ensure a minimal value of the variance. This kernel
represents the realization average of covariance matrices described by this
type of variation. They are similar to the structure functions used in the
pulsar timing literature for characterizing noise due to variations in the
ISM—for instance, see Equation (7.1) in Foster & Cordes (1990) or Section 5 of
Lam et al. (2016). Such structure functions were also used as a check in
Tiburzi et al. (2019) to ensure the iterative process used therein full
disentangled the solar wind signal form the DM variations in the ISM. Our
Bayesian framework allows for simultaneous fitting across multiple pulsars and
frees us from the iterative process used in this latter publication. Tiburzi
et al. (2021) also uses a Bayesian analysis to separate these signals, but
only within single pulsar datasets in order to keep track of solar wind
dependencies on ecliptic latitude. The GP model above (along with our solar
wind model) replaces the use of a short-time-span DMX model normally used in
NANOGrav datasets (Arzoumanian et al., 2018a; Alam et al., 2020a, b), since
DMX cannot distinguish between the ISM and solar wind contributions to the
measured DM value at each epoch. Here a standard $15$ day linear interpolation
basis in the time domain is used for constructing realizations of the DM
variations.

###  3.3 Noise Analyses

Since we are implementing entirely new DM variation models from the original
NG11 analysis, we redid single pulsar noise analyses using our solar wind and
DM GP models. This allowed us to recover new values for the standard white
noise parameters, EFAC, EQUAD and ECORR, (see e.g., Arzoumanian et al., 2015;
Lam et al., 2017) and power law achromatic red noise parameters for each
pulsar. These single pulsar noise runs then consist of a set of three white
noise parameters for each backend/receiver combination used for observations,
two DM GP parameters, $\sigma$ and $\ell$, a solar electron density, $n_{E}$,
and two power law red noise parameters, $A_{\rm RN}$ and $\gamma$.

Using only one pulsar’s dataset it is impossible to fully disentangle the
solar wind electron density from the ionized ISM. Any given observation’s DM
variation can be attributed to either changes in the interplanetary or
interstellar media. The key to these analyses is a global fit where the DM GP
and red noise parameters are allowed to vary individually for each pulsar,
while the solar electron density is set as a global parameter across all of
the pulsars. Since pulsars have upwards of $20+$ white noise parameters it is
cumbersome to vary all of these parameters across the entire PTA, so we adopt
the usual practice in PTA analyses (Arzoumanian et al., 2020) of setting their
values constant in the full PTA analysis following single-pulsar noise
modeling.

##  4 Bayesian Model Implementation and Analysis

We use the expressions in §2.1 to build models for a spherically symmetric
solar electron density and analyze NG11. We explore three different models for
the behavior and structure of the solar electron density: higher order
spherically symmetric solar wind models, time-binned solar wind density
models, and continuous time-dependent solar wind perturbations. Each model
demonstrates how our full Bayesian framework provides access to the unique
probes of the solar wind available in PTA datasets.

###  4.1 Higher Order Models

Any number of higher order spherically symmetric models can be constructed
using the tools developed here. We present only a few to demonstrate the
feasibility of detecting these higher order signals using PTA datasets. In
general we follow the methods of Aksim et al. (2019) and fit for a varying
power dependence, $\sim 1/r^{p}$, of the solar electron density, though we
still investigate models of the type in You et al. (2007, 2012) with multiple
summed fixed-power terms. The fit for the proportion of the two-phases (slow
and fast solar wind) in You et al. (2012) is dependent on the individual
lines-of-sight for the individual pulsars which is out of the scope of the
current tests for these modeling tools.

Table 1: Table of higher order spherically symmetric solar wind density terms. Parameters without sub(super)scripts are set constant. The median of free parameters is reported along with the $68\%$ credible interval. Parameters with a superscript 95% represent $95\%$ upper limits of the parameter. | First Term | Second Term | Third Term
---|---|---|---
Model |  +$p$ | $n^{(p)}_{E}\;[{\rm cm}^{-3}]$ |  +$p$ | $\log_{10}n^{(p)}_{E}$ |  +$p$ | $\log_{10}n^{(p)}_{E}$
1 |  +$2$ | $6.9^{+0.13}_{-0.13}$ | — | —
2 |  +$2.29^{+0.01}_{-0.02}$ | $4.7^{+0.20}_{-0.20}$ | — | —
3 |  +$2$ | $2.5^{+0.56}_{-0.63}$ |  +$2.4^{+0.04}_{-0.05}$ |  +$0.38^{+0.10}_{-0.07}$ | —
4 |  +$2.30^{+0.02}_{-0.01}$ | $4.5^{+0.18}_{-0.22}$ |  +$4.3^{+0.58}_{-0.73}$ |  +$-4.1^{95\%}$ | —
5 |  +$2$ | $6.8^{+0.18}_{-0.19}$ |  +$4.39$ |  +$-2.87^{+0.03}_{-0.02}$ |  +$16.25$ |  +$-24.8^{95\%}$

In Table 1 we summarize results from a number of similar analyses. The aim
here was not to exhaustively test the various powers from Equation 6 and
Equation 7, as the real model is expected to be a mixture of the two (You et
al., 2007), but rather to test the sensitivity of NG11 to some of these higher
order terms. Table 1 shows a progression from the overly simplistic $1/r^{2}$
model, where we find $n^{(2)}_{E}=6.9^{+0.13}_{-0.13}\frac{1}{{\rm cm^{3}}}$,
to models with more terms, and more free parameters. Model 2, where the
exponent in the density relation is allowed to vary, shows that there is
support for a higher order model for the solar electron density,
$p=2.29^{+0.01}_{-0.02}$. This power is in agreement with the best fit power
from Aksim et al. (2019) of $p=2.3$. Their value of $n^{2.3}_{E}=2.59\pm
0.13\frac{1}{{\rm cm^{3}}}$ is smaller than our recovered value, but the Aksim
et al. (2019) value is from only one observation, while our value is an $11.4$
year average. We will see in §4.2 that time varying values can differ by up to
a factor of 4 or more.

Models 3 and 4 show that there is broad support for a second term in the
density model when the power of the first term is kept constant, but the
analysis only retrieves upper limits for a third term when the first two are
kept constant. These analyses broadly support what has been previously shown
in You et al. (2012, 2007) that pulsar timing data, especially for an array as
we are treating here, is sensitive to a much more complicated model of the
solar electron density than is routinely assumed. Model 5 shows the limits of
the data set’s sensitivity, as we do not seem to detect a third term with the
large index predicted in You et al. (2007).

###  4.2 Time-Binned Solar Wind Density

As in Madison et al. (2019) and Tiburzi et al. (2021), we search for separate
values of $n_{E}$ with $p=2$ during different periods of time by fitting
values of $n_{E}$ for discrete time bins of the dataset. In the Madison et al.
(2019) analysis the bins are set to be a year long and a fit is done to the
DMX times series of the pulsars. DMX is a piece-wise time-binned analysis of
the DM variations, (see, e.g. Arzoumanian et al., 2015). We present our
results and those of Madison et al. (2019) in Figure 3, along with binned ACE
SWEPAM data for the same time period. The ACE data are in situ measurements of
the local electron density, and so can be considerably noisier than the large
effective average taken by the column density measurements to which the pulsar
timing data are sensitive. In addition, since ACE orbits in the ecliptic
plane, the SWEPAM instrument is really only sampling the slow solar wind, and
would not track changes in the higher altitude fast wind. This is supported by
scintillation studies that track high altitude solar wind density, see for
instance Porowski et al. (2021) for a model based on scintillation data. Their
model for solar electron density shows strong positive correlations on the
high latitude wind density with the solar cycle. Lastly, the ACE data are
taken by the same instrument across the 11+ year time span shown here, while
the NANOGrav PTA continually increased the number of pulsar in the array
(hence LOSs) and sensitivity, including a “high cadence” observing campaign in
the latter portion of the dataset (Arzoumanian et al., 2018a).

The largest differences between our values and those of Madison et al. (2019)
occur in the first part of the dataset where the observation cadence was
significantly less.

Figure 3: Year binned solar electron density. The blue plusses show the values
of our binned solar electron density search. The horizontal error bars show
the width of the bins used, mostly one year here, except the last bin which is
1.4 years long. The orange lines with vertical error bars show the results
from Madison et al. (2019) that used this same dataset, (NG11). The green
shaded blocks show median and 68% confidence intervals from the in situ
measurements of the ACE satellite’s SWEPAM instrument. The pink and purple
shaded vertical regions show the minimum and maximum, respectively, months of
solar cycle 24.

The fit in Madison et al. (2019) was done using a $\chi^{2}$-minimization, and
did not include any priors for the solar electron density. Hence the
unphysical negative value of the second bin is a result of the
$\chi^{2}$-minimization. The large value in the fifth bin occurs during a time
when there is a gap in the PTA data. This will be explored more in the finer
binning used next. The Madison et al. (2019) analysis is not as robust as the
analysis here, since rather than fitting a secondary data product (DMX) we are
fitting the solar wind as part of a full PTA-wide analysis, which includes a
noise model fit and timing model marginalization.

In Figure 4 we show Bayesian binned results for bin sizes of three months
again with $p=2$, showing that NG11 has enough information to fit $n_{E}$ with
finer resolution than is done in Madison et al. (2019). The $68\%$ credible
intervals for the binned values are small compared to the same intervals in
the ACE data. In addition to the binned fit for $n_{E}$, we also run an
analysis with a second density, $n^{(4.39)}_{E}$, in the solar wind density
profile model with $p=4.39$. We chose this index from the fast solar wind
model of Equation 7 in order to investigate how a higher order term might
effect the binned values of $n_{E}$ specifically from the few pulsars with
very small $\theta_{i}$. We use a time-constant value of $n^{(4.39)}_{E}$ in
this application for ease of analysis.

This second analysis, with an $n^{(4.39)}_{E}$ parameter highlights the
importance of these higher order terms in the density model when there are
observations with lines of sight close to the Sun. The geometric factor for
the ${\rm DM}_{\odot}$ delay, $(\pi-\theta_{i})/\sin\theta_{i}$, of all TOAs
is shown in the bottom panel of Figure 4. There are two occasions in the
11-year dataset where $\theta_{i}$ is very small for PSR J0030+0451 and PSR
J1614$-$2230\. These can be seen as the points where the geometric factor is
larger than $\sim 100$. The two bins, marked by vertical dots, where these
occur correspond to the two bins where the recovered values for $n_{E}$ differ
by many standard deviations from the usual $1/r^{2}$ model for the density,
over estimating $n_{E}$. While the $p=4.39$ model seems to ameliorate the
effect of these small $\theta_{i}$ observations, the direct cause of the
anomalously high measurements of $n_{E}$ could also be due to the highly
inhomogeneous streaming of the solar wind at these small distances or from
individual events, i.e. coronal mass ejections. See Appendix B for more
discussion.

Two other bins warrant attention. The first bin has only a few TOAs with
narrow frequency coverage, and hence parameters recovered have large error
bars. The bin near MJD 54250 has no TOAs due to concurrent down time at both
the Arecibo Observatory and the Green Bank Observatory. In this bin, the ACE
prior is returned in the posterior. Hence no information beyond the prior is
gained.

Figure 4: Three month binned solar electron density. The horizontal error bars
show the width of the bins used, while the vertical error bars show the inner
$68\%$ credible region of the binned $n_{E}$ posteriors. The pink and purple
shaded vertical regions show the minimum and maximum, respectively, months of
solar cycle 24. The vertical dotted lines show the boundaries of a few
important bins. The first shows the end of the first bin which contains the
fewest number of TOAs. The bin near MJD $54250$ has no TOAs due to down time
at both the Arecibo Observatory and the Green Bank Observatory. The dotted
vertical lines near MJD $55800$ and MJD $56400$ show bins where the values for
$n_{E}$ in the two models differ by a substantial amount. The bottom panel
shows the geometric part of the solar DM obtained from a $1/r^{2}$ model for
the solar wind density calculated for all TOAs in the dataset. Note the log
scale on the y-axis. Only two pulsars, J0030+0451 and J1614$-$2230, have
geometric factors larger than $100$ and the factors for these have been
highlighted with colored x’s and +’s, respectively.

###  4.3 Continuous Time-Dependent Solar-Wind Perturbations

In addition to the piecewise constant models for $n_{E}$ discussed above, we
also implemented a continuous Fourier-basis model for $n_{E}(t)$ as a
perturbation to a mean $n_{E}$. The model is reminiscent of the free spectral
models used by PTAs to describe achromatic red noise (Arzoumanian et al.,
2020). It is constructed using the TOAs for each pulsar, parametrized with a
separate value for each component of the basis,

$n^{F}_{e}(t)=\sum^{N_{f}}_{j=0}a_{j}\sin(2\pi f_{j}t)+b_{j}\cos(2\pi
f_{j}t)\;.$ (13)

Rather than trying to reproduce the fine structure recoverable by the You et
al. (2012) model, here we are concerned with long to mid scale variations of
$n_{E}$, therefore we use a set of frequencies based on the time span,
$T=11.4$ yr, of NG11 and include 30 linearly spaced frequencies ranging from
$[1/T,30/T]$, using 2 parameters per frequency, $a_{j}$ and $b_{j}$.

In various test analyses we recorded the same type of behavior for these
continuous models as was seen in Figure 4, i.e., the fit returned large values
of $n_{E}$ during periods of time when J0030+0451 and J1614$-$2230
observations had small values of $\theta_{i}$. Here we only report on the
continuous model for $n_{E}$ where we also fitted for the same static, higher
order $n^{(4.39)}_{E}$ term discussed in §4.2. Figure 5 shows the resulting
recovery of $n_{E}(t)=n^{mean}_{E}+n^{F}_{E}(t)$ from an analysis that uses
the model in Equation 13.

Figure 5: Continuous solar wind density model. The black dots show the same
binned values for $n_{E}$ shown in orange in the top panel of Figure 4. The
blue trace and shading show the median and inner $68\%$ credible interval of
1000 realizations of the continuous solar electron density perturbation model.

The use of the Fourier basis gives us immediate access to an analogue of the
power spectral density for the perturbations to the solar electron density.
Figure 6 shows violin plots for the power in the various frequencies
calculated using the posteriors of the Fourier decomposition coefficients,
$P=a_{i}^{2}+b_{i}^{2}$. Notice that a majority of the frequencies show small
amounts of power, while the lowest frequency, $1/(11.4{\rm\,yr})$, recovers
the most power — very close to an inverse average solar cycle length of
$11.75\,{\rm yr}$. It will be interesting to see which frequencies have the
most power in longer datasets.

Figure 6: Power in the continuous solar wind density model. The violin plots
show a simple representation of a power posterior calculated as
$a^{2}_{i}+b^{2}_{i}$ for each set of Fourier coefficients. The horizontal
bars show the median values. The orange dotted line shows a frequency of
$1/{\rm yr}$ while the blue dashed vertical line shows $1/(11.75\,{\rm yr})$,
an inverse average solar cycle.

A continuous signal like this would be a good candidate for a Gaussian
process, however current code bases for PTAs construct different realizations
of a process for each pulsar in an array. Since this work reveals a great deal
of potential for gleaning information about the solar wind from PTAs we plan
on developing per LOS solar wind perturbations for these pulsars in future
work by using a simple deterministic inhomogeneous model.

##  5 Discussion

This work demonstrates a new set of Bayesian tools that allow for the use of
PTA datasets, or any collection of pulsar timing data, as probes of the solar
wind. The main innovations include a completely general Bayesian fit that
includes timing model marginalization, generic spherically symmetric models
for the solar electron density, and time dependent models. These models can be
assembled in myriad ways in order to mitigate the solar wind as a noise source
in pulsar data, or to use these sensitive datasets to study the Sun’s coronal
behavior across the solar system. The NANOGrav data set is not aimed at
monitoring the solar wind, nonetheless it is sensitive to interesting solar
phenomena. The tools we have developed, if applied to a PTA data set built for
the purposes of solar science, have great potential.

One important aspect of this work not yet discussed is how these models can
enable more in depth generic ISM studies. These solar wind models allow the
removal of the DM variation signal from the local interplanetary medium in a
principled way. This would better enable the use of these pulsar timing
datasets for studying the ISM unencumbered by local variations in electron
density.

The solar wind signal has been known to cause issues (Lommen et al., 2006;
Splaver et al., 2005) when trying to measure some timing effects in pulsar
data sets—for instance, timing parallax, critical for constraining pulsar
distances. The covariances with measuring parallax imply that improvements in
the solar wind model could ultimately improve single source GW searches, as
the distance to pulsars is an important aspect of pulsar term searches for GW
from supermassive black binaries (Aggarwal et al., 2019). The solar wind has
also been implicated as one of the systematics in the “triple system”
(Archibald et al., 2018), impeding better limits of Einstein’s equivalence
principle. Better modeling of the solar wind might allow for more accurate
measurements of these limits, in addition to allowing for better noise
mitigation when these data are used for GW searches. Improvements to the
handling of the solar wind in pulsar timing investigations stand to improve
precision pulsar science in a wide variety of ways.

PTA studies of the solar wind have the potential to add a completely
independent set of measurements to the myriad space missions studying solar
weather. PTA observations survey a much wider swath of the Sun’s environment
than any spacecraft would be able to, taking column density information in 70+
lines of sight every month, probing beyond the outermost reaches of the solar
environment. The observational campaigns of lower frequency observatories like
the Canadian Hydrogen Intensity Mapping Experiment (CHIME/Pulsar Collaboration
et al., 2021), LOFAR (Stappers et al., 2011; van Haarlem et al., 2013) and the
Long Wavelength Array (LWA, Stovall et al., 2015) are continually adding data
that is extremely useful for DM variability studies, since these lower
frequencies allow for more accurate measurement of the variations. These data
can add important tests to the well developed solar weather modeling efforts
undertaken by the solar physics community (van der Holst et al., 2014). These
large scale measurements would allow for long timescale monitoring of the fast
and slow solar wind, giving access to continuous measurements of the solar
wind at all solar latitudes.

Author contributions. This paper is the result of the work of many people and
uses data from over a decade of pulsar timing observations. We list specific
contributions below. JSH developed and ran the solar wind analyses and led the
paper writing. JS, and DRM contributed substantially to paper writing,
discussion and interpretation of results. ZA, KC, MED, PBD, TD, JAE, RDF, ECF,
EF, PAG, GJ, MLJ, MTL, LL, DRL, RSL, MAM, CN, DJN, TTP, SMR, PSR, RS, IHS, KS,
JKS, and WZ ran observations and developed the 11-year data set.
Acknowledgments. The NANOGrav project receives support from National Science
Foundation (NSF) as a Physics Frontier Center award number 1430284. We thank
Caterina Tiburzi, Aditya Parthasarathy, Joseph Lazio and Patrick O’Neill for
comments on an early draft of this manuscript. TD and MTL are supported by an
NSF Astronomy and Astrophysics Grant (AAG) award number 2009468. ECF is
supported by NASA under award number 80GSFC17M0002. The Dunlap Institute is
funded by an endowment established by the David Dunlap family and the
University of Toronto. TTP acknowledges support from the MTA-ELTE
Extragalactic Astrophysics Research Group, funded by the Hungarian Academy of
Sciences (Magyar Tudományos Akadémia), which was used during the development
of this research. SMR is a CIFAR Fellow. Portions of this work performed at
NRL were supported by ONR 6.1 basic research funding. Pulsar research at UBC
is supported by an NSERC Discovery Grant and by CIFAR. WWZ is supported by the
CAS Pioneer Hundred Talents Program and the Strategic Priority Research
Program of the Chinese Academy of Sciences, grant No. XDB23000000.

##  Appendix A Generic Spherically Symmetric DM for Integer Powers

The expression in Equation 11 simplifies considerably for integer values of
$p>1$. Here we show the first few expressions for powers greater than $p=2$
for reference.

$\displaystyle\text{DM}_{\odot}^{3}$ $\displaystyle=$ $\displaystyle
n_{E}^{\left(3\right)}\left(1\,\text{au}\right)^{3}\frac{1+\cos\theta_{i}}{R_{E}^{2}\sin^{2}\theta_{i}}$
(A1) $\displaystyle\text{DM}_{\odot}^{4}$ $\displaystyle=$ $\displaystyle
n_{E}^{\left(4\right)}\left(1\,\text{au}\right)^{4}\frac{2\pi-2\theta_{i}+\sin
2\theta_{i}}{4R_{E}^{3}\sin^{3}\theta_{i}}$ (A2)
$\displaystyle\text{DM}_{\odot}^{5}$ $\displaystyle=$ $\displaystyle
n_{E}^{\left(5\right)}\left(1\,\text{au}\right)^{5}\frac{2-\cos\theta_{i}}{12R_{E}^{4}\sin^{4}\frac{\theta_{i}}{2}}$
(A3)

##  Appendix B Coronagraph Images

As discussed in §4.2, three observation epochs from NG11 have fairly small
solar impact angles. In Figure 7 we show images from the Large Angle
Spectroscopic Coronagraph (LASCO, Brueckner et al., 1995) on board the Solar
and Heliospheric Observatory (SOHO) from the same time as the pulsar
observations, with the approximate pulsar positions marked. These white light
images from the C2 camera can be related to the solar electron density
(Quémerais & Lamy, 2002), however such calculations are out of the scope of
the current work.

Figure 7: Coronagraph images from the LASCO/C2 instrument onboard the SOHO
spacecraft. The inset image is from the Extreme-Ultraviolet Imaging Telescope
(EIT, Delaboudinière et al., 1995). These coincide with the three pulsar
timing observations discussed in §4.2. The approximate pulsar positions are
marked in each image with a white crosshairs. Images a) and b) are for the PSR
J1614$-$2230 observations and are separated by one day. Image c) is for the
PSR J0030+0451 observation. The dates and time in UT are given in the images.

These images are provided to highlight how close to the Sun these observations
were taken ($<6R_{\odot}$) and how inhomogeneous the inner solar system can
be. Testing various estimates of 2 and 3-dimensional electron density from
these white light images would be an interesting use of dedicated pulsar
observations close to the Sun.

## References

  * Aggarwal et al. (2019) Aggarwal, K., Arzoumanian, Z., Baker, P. T., et al. 2019, ApJ, 880, 116, doi: 10.3847/1538-4357/ab2236
  * Aksim et al. (2019) Aksim, D., Melnikov, A., Pavlov, D., & Kurdubov, S. 2019, ApJ, 885, 159, doi: 10.3847/1538-4357/ab499a
  * Alam et al. (2020a) Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2020a, arXiv e-prints, arXiv:2005.06490.  https://arxiv.org/abs/2005.06490
  * Alam et al. (2020b) —. 2020b, arXiv e-prints, arXiv:2005.06495.  https://arxiv.org/abs/2005.06495
  * Allen (1947) Allen, C. W. 1947, MNRAS, 107, 426, doi: 10.1093/mnras/107.5-6.426
  * Archibald et al. (2018) Archibald, A. M., Gusinskaia, N. V., Hessels, J. W. T., et al. 2018, Nature, 559, 73, doi: 10.1038/s41586-018-0265-1
  * Arzoumanian et al. (2015) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2015, ApJ, 813, 65
  * Arzoumanian et al. (2018a) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018a, ApJS, 235, 37, doi: 10.3847/1538-4365/aab5b0
  * Arzoumanian et al. (2018b) Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2018b, ApJ, 859, 47, doi: 10.3847/1538-4357/aabd3b
  * Arzoumanian et al. (2020) Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, ApJ, 905, L34, doi: 10.3847/2041-8213/abd401
  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., et al. 2016, Space Sci. Rev., 204, 49, doi: 10.1007/s11214-016-0244-5
  * Bale et al. (2019) Bale, S. D., Badman, S. T., Bonnell, J. W., et al. 2019, Nature, 576, 237, doi: 10.1038/s41586-019-1818-7
  * Bame et al. (1992) Bame, S. J., McComas, D. J., Barraclough, B. L., et al. 1992, A&AS, 92, 237
  * Brueckner et al. (1995) Brueckner, G. E., Howard, R. A., Koomen, M. J., et al. 1995, Sol. Phys., 162, 357, doi: 10.1007/BF00733434
  * CHIME/Pulsar Collaboration et al. (2021) CHIME/Pulsar Collaboration, Amiri, M., Bandura, K. M., et al. 2021, ApJS, 255, 5, doi: 10.3847/1538-4365/abfdcb
  * Coles (1978) Coles, W. A. 1978, Space Sci. Rev., 21, 411, doi: 10.1007/BF00173067
  * Cordes & Shannon (2010) Cordes, J. M., & Shannon, R. M. 2010, arXiv e-prints, arXiv:1010.3785.  https://arxiv.org/abs/1010.3785
  * Cordes et al. (2016) Cordes, J. M., Shannon, R. M., & Stinebring, D. R. 2016, ApJ, 817, 16, doi: 10.3847/0004-637X/817/1/16
  * Counselman et al. (1970) Counselman, C. C., Rankin, J. M., & Richards, D. W. 1970, in Bulletin of the American Astronomical Society, Vol. 2, 189
  * Delaboudinière et al. (1995) Delaboudinière, J. P., Artzner, G. E., Brunaud, J., et al. 1995, Sol. Phys., 162, 291, doi: 10.1007/BF00733432
  * Edwards et al. (2006) Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, MNRAS, 372, 1549, doi: 10.1111/j.1365-2966.2006.10870.x
  * Ellis & van Haasteren (2017) Ellis, J., & van Haasteren, R. 2017, jellis18/PTMCMCSampler: Official Release, doi: 10.5281/zenodo.1037579
  * Ellis et al. (2019) Ellis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T. 2019, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE.  http://ascl.net/1912.015
  * Foster & Cordes (1990) Foster, R. S., & Cordes, J. M. 1990, ApJ, 364, 123, doi: 10.1086/169393
  * Goncharov et al. (2021) Goncharov, B., Reardon, D. J., Shannon, R. M., et al. 2021, MNRAS, 502, 478, doi: 10.1093/mnras/staa3411
  * Guhathakurta & Fisher (1998) Guhathakurta, M., & Fisher, R. 1998, ApJ, 499, L215, doi: 10.1086/311371
  * Guhathakurta & Fisher (1995) Guhathakurta, M., & Fisher, R. R. 1995, Geophys. Res. Lett., 22, 1841, doi: 10.1029/95GL01603
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hazboun et al. (2019) Hazboun, J. S., Romano, J. D., & Smith, T. L. 2019, Phys. Rev. D, 100, 104028, doi: 10.1103/PhysRevD.100.104028
  * Hazboun et al. (2020) Hazboun, J. S., Simon, J., Taylor, S. R., et al. 2020, ApJ, 890, 108, doi: 10.3847/1538-4357/ab68db
  * Hewish et al. (1968) Hewish, A., Bell, S. J., Pilkington, J. D. H., Scott, P. F., & Collins, R. A. 1968, Nature, 217, 709, doi: 10.1038/217709a0
  * Hewish & Dennison (1967) Hewish, A., & Dennison, P. A. 1967, J. Geophys. Res., 72, 1977, doi: 10.1029/JZ072i007p01977
  * Hobbs & Edwards (2012) Hobbs, G., & Edwards, R. 2012, Tempo2: Pulsar Timing Package.  http://ascl.net/1210.015
  * Hobbs et al. (2020) Hobbs, G., Guo, L., Caballero, R. N., et al. 2020, MNRAS, 491, 5951, doi: 10.1093/mnras/stz3071
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Issautier et al. (2001) Issautier, K., Hoang, S., Moncuquet, M., & Meyer-Vernet, N. 2001, Space Sci. Rev., 97, 105, doi: 10.1023/A:1011878228168
  * Jones et al. (2017) Jones, M. L., McLaughlin, M. A., Lam, M. T., et al. 2017, ApJ, 841, 125, doi: 10.3847/1538-4357/aa73df
  * Kasper et al. (2016) Kasper, J. C., Abiad, R., Austin, G., et al. 2016, Space Sci. Rev., 204, 131, doi: 10.1007/s11214-015-0206-3
  * Knudsen et al. (1980) Knudsen, W. C., Spenner, K., Bakke, J., & Novak, V. 1980, IEEE Transactions on Geoscience and Remote Sensing, 18, 54, doi: 10.1109/TGRS.1980.350261
  * Lam (2018) Lam, M. T. 2018, ApJ, 868, 33, doi: 10.3847/1538-4357/aae533
  * Lam et al. (2015) Lam, M. T., Cordes, J. M., Chatterjee, S., & Dolch, T. 2015, ApJ, 801, 130, doi: 10.1088/0004-637X/801/2/130
  * Lam et al. (2016) Lam, M. T., Cordes, J. M., Chatterjee, S., et al. 2016, ApJ, 821, 66, doi: 10.3847/0004-637X/821/1/66
  * Lam & Hazboun (2021) Lam, M. T., & Hazboun, J. S. 2021, ApJ, 911, 137, doi: 10.3847/1538-4357/abeb64
  * Lam et al. (2018a) Lam, M. T., McLaughlin, M. A., Cordes, J. M., Chatterjee, S., & Lazio, T. J. W. 2018a, ApJ, 861, 12, doi: 10.3847/1538-4357/aac48d
  * Lam et al. (2018b) —. 2018b, ApJ, 861, 12, doi: 10.3847/1538-4357/aac48d
  * Lam et al. (2017) Lam, M. T., Cordes, J. M., Chatterjee, S., et al. 2017, ApJ, 834, 35, doi: 10.3847/1538-4357/834/1/35
  * Lam et al. (2018c) Lam, M. T., Ellis, J. A., Grillo, G., et al. 2018c, ApJ, 861, 132, doi: 10.3847/1538-4357/aac770
  * Lentati et al. (2016) Lentati, L., et al. 2016, Mon. Not. Roy. Astron. Soc., 458, 2161, doi: 10.1093/mnras/stw395
  * Lommen et al. (2006) Lommen, A. N., Kipphorn, R. A., Nice, D. J., et al. 2006, ApJ, 642, 1012, doi: 10.1086/501067
  * Luo et al. (2021) Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45, doi: 10.3847/1538-4357/abe62f
  * Madison et al. (2019) Madison, D. R., Cordes, J. M., Arzoumanian, Z., et al. 2019, ApJ, 872, 150, doi: 10.3847/1538-4357/ab01fd
  * Manoharan (2010) Manoharan, P. K. 2010, Sol. Phys., 265, 137, doi: 10.1007/s11207-010-9593-5
  * McComas et al. (1998) McComas, D. J., Bame, S. J., Barker, P., et al. 1998, Space Sci. Rev., 86, 563, doi: 10.1023/A:1005040232597
  * Muhleman & Anderson (1981) Muhleman, D. O., & Anderson, J. D. 1981, ApJ, 247, 1093, doi: 10.1086/159119
  * Nice et al. (2015) Nice, D., Demorest, P., Stairs, I., et al. 2015, Tempo: Pulsar timing data analysis, Astrophysics Source Code Library.  http://ascl.net/1509.002
  * Niu et al. (2017) Niu, Z.-X., Hobbs, G., Wang, J.-B., & Dai, S. 2017, Research in Astronomy and Astrophysics, 17, 103, doi: 10.1088/1674-4527/17/10/103
  * Porowski et al. (2021) Porowski, C., Bzowski, M., & Tokumaru, M. 2021, arXiv e-prints, arXiv:2110.15847.  https://arxiv.org/abs/2110.15847
  * Quémerais & Lamy (2002) Quémerais, E., & Lamy, P. 2002, A&A, 393, 295, doi: 10.1051/0004-6361:20021019
  * Ransom et al. (2019) Ransom, S., Brazier, A., Chatterjee, S., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 195.  https://arxiv.org/abs/1908.05356
  * Splaver et al. (2005) Splaver, E. M., Nice, D. J., Stairs, I. H., Lommen, A. N., & Backer, D. C. 2005, ApJ, 620, 405, doi: 10.1086/426804
  * Stappers et al. (2011) Stappers, B. W., Hessels, J. W. T., Alexov, A., et al. 2011, A&A, 530, A80, doi: 10.1051/0004-6361/201116681
  * Stovall et al. (2015) Stovall, K., Ray, P. S., Blythe, J., et al. 2015, ApJ, 808, 156, doi: 10.1088/0004-637X/808/2/156
  * Taylor et al. (2021) Taylor, S. R., Baker, P. T., Hazboun, J. S., Simon, J., & Vigeland, S. J. 2021, enterprise$\\_$extensions.  https://github.com/nanograv/enterprise_extensions
  * Taylor et al. (2013) Taylor, S. R., Gair, J. R., & Lentati, L. 2013, Phys. Rev. D, 87, 044035, doi: 10.1103/PhysRevD.87.044035
  * Tiburzi et al. (2016) Tiburzi, C., Hobbs, G., Kerr, M., et al. 2016, MNRAS, 455, 4339, doi: 10.1093/mnras/stv2143
  * Tiburzi et al. (2019) Tiburzi, C., Verbiest, J. P. W., Shaifullah, G. M., et al. 2019, MNRAS, 487, 394, doi: 10.1093/mnras/stz1278
  * Tiburzi et al. (2021) Tiburzi, C., Shaifullah, G. M., Bassa, C. G., et al. 2021, A&A, 647, A84, doi: 10.1051/0004-6361/202039846
  * Tokumaru (2013) Tokumaru, M. 2013, Proceedings of the Japan Academy, Series B, 89, 67, doi: 10.2183/pjab.89.67
  * Vallisneri et al. (2020) Vallisneri, M., Taylor, S. R., Simon, J., et al. 2020, ApJ, 893, 112, doi: 10.3847/1538-4357/ab7b67
  * van der Holst et al. (2014) van der Holst, B., Sokolov, I. V., Meng, X., et al. 2014, ApJ, 782, 81, doi: 10.1088/0004-637X/782/2/81
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2, doi: 10.1051/0004-6361/201220873
  * van Haasteren & Levin (2013) van Haasteren, R., & Levin, Y. 2013, Mon. Not. Roy. Astron. Soc., 428, 1147, doi: 10.1093/mnras/sts097
  * van Haasteren et al. (2009) van Haasteren, R., Levin, Y., McDonald, P., & Lu, T. 2009, MNRAS, 395, 1005, doi: 10.1111/j.1365-2966.2009.14590.x
  * van Haasteren & Vallisneri (2014) van Haasteren, R., & Vallisneri, M. 2014, Phys. Rev. D, 90, 104012, doi: 10.1103/PhysRevD.90.104012
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * You et al. (2012) You, X. P., Coles, W. A., Hobbs, G. B., & Manchester, R. N. 2012, MNRAS, 422, 1160, doi: 10.1111/j.1365-2966.2012.20688.x
  * You et al. (2007) You, X. P., Hobbs, G. B., Coles, W. A., Manchester, R. N., & Han, J. L. 2007, ApJ, 671, 907, doi: 10.1086/522227


# Measurement of the temperature dependence of dwell time and spin relaxation
probability of Rb atoms on paraffin surfaces using a beam-scattering method

Kanta Asakawa1<EMAIL_ADDRESS>Yutaro Tanaka1 Kenta Uemura1 Norihiro
Matsuzaka2 Kunihiro Nishikawa1 Yuki Oguma1 Hiroaki Usui2 Atsushi Hatakeyama1
<EMAIL_ADDRESS>1Department of Applied Physics, Tokyo University of
Agriculture and Technology, Koganei, Tokyo 184-8588, Japan  2Department of
Organic and Polymer Materials Chemistry, Tokyo University of Agriculture and
Technology, Koganei, Tokyo 184-8588, Japan

###### Abstract

The scattering of Rb atoms on an anti-relaxation coating was studied. No
significant change in the spin relaxation probability of Rb atoms by single
scattering from a tetracontane surface was observed by cooling the film from
305 to 123 K. The mean surface dwell time was estimated using a time-resolved
method. Delay-time spectra, from which mean surface dwell times can be
estimated, were measured at 305, 153, and 123 K, with a time window of
$9.3\times 10^{-5}$ s. The increase in mean surface dwell time with cooling
from 305 to 123 K was smaller than $4.4\times 10^{-6}$ s, which is
significantly smaller than the value expected from the mean dwell time at room
temperature measured using the Larmor frequency shift. These results can be
explained by assuming a small number of scattering components, with a mean
surface dwell time at least three orders of magnitude longer than the majority
component.

Spin

###### pacs:

73.20.At,68.47.Gh

##  I Introduction

Anti-relaxation coatings are used to reduce the spin relaxation of alkali
metal atoms resulting from wall collisions in the alkali-metal vapor cells of
atomic clocks Risley  _et al._ (1980); Robinson and Johnson (1982); Frueholz
_et al._ (1983) and atomic magnetometers Budker  _et al._ (1998); Balabas  _et
al._ (2006); Wasilewski  _et al._ (2010). ParaffinRobinson  _et al._ (1958);
Bouchiat and Brossel (1966), octadecyltrichlorosilane (OTS)Seltzer  _et al._
(2007); Zhao  _et al._ (2008), and polydimethylsiloxane (PDMS)Atutov  _et al._
(2016) are representative materials for anti-relaxation coatings. It has been
reported that paraffin-coated surfaces can support $10^{4}$ spin-preserving
collisions for Rb atoms Bouchiat and Brossel (1966). The performance of an
anti-relaxation coating depends on the surface dwell time, as well as the
strength of the interaction between alkali metal spins and the surface. The
mean dwell time $\tau_{\textrm{s}}$ can be described by the Arrhenius formula:

$\tau_{\textrm{s}}=\tau_{0}\exp\left(\frac{E_{\textrm{des}}}{k_{\textrm{B}}T_{\textrm{s}}}\right),$
(1)

where $\tau_{0}$ is the pre-exponential factor, $E_{\textrm{des}}$ is the
desorption energy, $k_{\textrm{B}}$ is Boltzmann constant, and
$T_{\textrm{s}}$ is the temperature of the surface. In the case of Rb atoms on
tetracontane (C40H82), which is a representative type of paraffin, the
experimentally obtained desorption energy is 0.06 eVRahman and Robinson
(1987); Budker  _et al._ (2005). By adopting a commonly used assumption, i.e.,
that the pre-exponential factor is $1\times 10^{-12}$ s, which is the typical
period of thermal vibration of atoms, we can roughly estimate
$\tau_{\textrm{s}}$ as $1\times 10^{-11}$ s at 300 K. However, using the
Larmor frequency shift, Ulanski  et al. reported a mean dwell time of $(1.8\pm
0.2)\times 10^{-6}$ s for Rb atoms on paraffin coatings at 345 K Ulanski and
Wu (2011). The reason for the large difference in the mean dwell time
calculated from the desorption energy compared to that measured by experiments
is still unclear Atutov and Plekhanov (2015). One possibility is that the
assumption $\tau_{0}\simeq 1\times 10^{-12}$ s is incorrect. By substituting
$E_{\textrm{des}}=0.06$ eV, $T_{\textrm{s}}=345$ K, and
$\tau_{\textrm{s}}=1.8\times 10^{-6}$ s into equation (1) and regarding
$\tau_{0}$ as a variable, we obtain $\tau_{0}=2.4\times 10^{-7}$ s. However,
this is two orders of magnitude larger than the pre-exponential factor
$\tau_{0}=2.2\times 10^{-9}$ s of 87Rb atoms on Pyrex glass surfaces coated
with OTS, estimated from the temperature dependence of the mean dwell time
Zhao  _et al._ (2009) and is five orders of magnitude larger than the typical
period of thermal vibration. Therefore, this issue requires further
investigation.

Equation (1) shows that the mean dwell time increases with cooling, which
makes it easier to measure dwell times using time-resolved methods. If we
assume that $\tau_{0}=2.4\times 10^{-7}$ s (which is obtained by substituting
$E_{\textrm{des}}=0.06$ eV Rahman and Robinson (1987); Budker  _et al._
(2005), $\tau_{\textrm{s}}=1.8\times 10^{-6}$ s, and $T_{\textrm{s}}=345$ K
Ulanski and Wu (2011)) is correct, then it can be seen from Eq. (1) that
$\tau_{\textrm{s}}$ will increase by $6.7\times 10^{-5}$ s with cooling of a
sample from 305 to 123 K, which is sufficient time to detect using time-
resolved methods. The spin relaxation probability with surface scattering is
also expected to increase at low temperatures due to an increased dwell time.

In this study, we investigated the temperature dependence of the spin
relaxation probability and dwell time of Rb atoms on tetracontane coatings. A
beam-scattering method and X-ray photoelectron spectroscopy (XPS) were
employed to analyze the surface dwell time and surface chemical composition.
Using an atomic beam and optical hyperfine pumping, the dwell time can be
measured more directly compared to the methods used in earlier studies Ulanski
and Wu (2011); Zhao  _et al._ (2009). The results show that the increase in
mean dwell time (averaged over the majority of scattered atoms) with cooling
from 305 to 123 K was shorter than $4.4\times 10^{-6}$ s, which is
significantly shorter than the value of $6.7\times 10^{-5}$ s expected from
the previously reported desorption energy $E_{\textrm{des}}=0.06$ eV Rahman
and Robinson (1987); Budker  _et al._ (2005) and the mean dwell time
$\tau_{\textrm{s}}=1.8\times 10^{-6}$ s at 345 K Ulanski and Wu (2011). This
indicates the existence of minor scattering components with dwell times at
least three orders of magnitude larger than that of the major component.

##  II Experimental

Figure 1: (a) Schematic illustration of the experimental setup and (b) the
energy-level diagram of 85Rb. Figure 2: (a) Schematic diagram of the optical
system for generating pump and probe light. AOM, EOM, BS, PBS, NPBS, RM, and
QWP denote the acousto-optic modulator, electro-optic modulator, beam
splitter, polarizing beam splitter, non-polarizing beam splitter, removable
mirror, and quarter-wave plate, respectively. (b) Schematic diagram of the
signal-processing system used for the delay-time measurements.

Figure 1(a) shows the experimental setup. A tetracontane-coated quartz
substrate was mounted in an ultra-high vacuum (UHV) chamber with a base
pressure lower than $3\times 10^{-7}$ Pa. The tetracontane film was deposited
on the substrate in another high-vacuum chamber, the base pressure of which
was $1.0\times 10^{-5}$ Pa, by evaporating tetracontane at 513 K for 10 min.
The thickness and average roughness, Ra, of the film were measured to be
$0.93\pm 0.17$ µm and 50 nm, respectively, using atomic force microscopy. The
Rb beam was generated using a multi-channel effusive atomic beam source. The
full width at half maximum of the atomic beam at the position of the sample
was estimated to be 8.2 mm from a fluorescence image taken with a charge-
coupled device (CCD) camera.

The two pump light beams were directed perpendicular to the Rb beam. The pump
and probe frequencies are shown in Fig. 1(b). The optical system used to
generate the pump and probe light is illustrated in Fig. 2(a). The frequency
of the first pump light was tuned to the $F=3\rightarrow F^{\prime}=3$
transition frequency of the 85Rb $D_{1}$ transition line using polarization
spectroscopy Wieman and Hänsch (1976); Harris  _et al._ (2006), where $F$ and
$F^{\prime}$ are the total angular momentum of atoms in the $5^{2}S_{1/2}$ and
$5^{2}P_{1/2}$ states, respectively. The second pump light, the frequency of
which was tuned to the $F=2\rightarrow F^{\prime}=3$ transition frequency of
the $D_{1}$ line, was generated by blue-detuning the first pump light by 3,036
MHz Arimondo  _et al._ (1977); Schultz  _et al._ (2008) using an electro-optic
modulator (EOM) and an acousto-optic modulator (AOM), and was pulsed to
$5\times 10^{-6}$ s by the AOM.

The probe light was used to selectively excite Rb atoms in the $F=3$ state to
the $F^{\prime\prime}=4$ state through the $D_{2}$ transition, resulting in
fluorescence. Here, $F^{\prime\prime}$ is the total angular momentum of the
atoms in the $5^{2}P_{3/2}$ state. The fluorescence was detected by a CCD
camera and a photomultiplier tube (PMT), which were equipped with interference
filters that were designed to transmit only the probe light and fluorescence.
The frequency of the probe light could be tuned, which enabled selective
excitation of the incident or scattered atoms. For probing scattered atoms,
the probe light was blue-detuned from the $F=3\rightarrow F^{\prime\prime}=4$
transition frequency of the $D_{2}$ transition line by 160 MHz, such that it
did not excite incident atoms. With this frequency, the probe light excites Rb
atoms whose velocity component along the probe light is 125$\pm$5 m/s; these
are abundant among scattered atoms Sekiguchi  _et al._ (2018) but negligible
among incident atoms. Here, the natural line width of the Rb $D_{2}$
transition line (6.06 MHz) Volz and Schmoranzer (1996) was used to calculate
the uncertainty. When probing the incident atoms, the direction of the probe
light beam was the same as that used to probe the scattered atoms; however,
the frequency was red-detuned by 60.3 MHz from the $F=3\rightarrow
F^{\prime\prime}=4$ transition frequency. With this frequency, the probe light
excites atoms with a velocity component along the probe light of $-47\pm 5$
m/s. Because atoms with this velocity component are found among both incident
and scattered atoms, the fluorescence intensity of the incident atoms was
estimated by subtracting the contribution of scattered atoms from the measured
fluorescence intensity, as discussed below.

Mean dwell-time estimates were based on time-of-flight (TOF) measurements
obtained using the pump and probe light. The incident Rb atoms were first
pumped to the $F=2$ state by the first pump light and subsequently irradiated
with the second pump light, which was pulsed; atoms that were irradiated with
the second pump light were momentarily pumped to the $F=3$ state. The incident
atoms pumped to the $F=3$ state by the second pump light enhanced the
fluorescence induced by the probe light when they reached it. The delay time
in the fluorescence enhancement induced by irradiation by the second pump
light is the sum of the TOF of the Rb atoms (from the second pump light to the
probe light via the film surface) and the surface dwell time. The probe-light-
induced fluorescence was detected by the PMT. The signal from the PMT was
processed using the system shown in Fig. 2(b). Delay-time spectra were
acquired by accumulating the time intervals between the irradiation of the
second probe light and the detection of fluorescence by PMT. The delay-time
distribution can be treated as the distribution of the sum of the TOF and the
dwell time only when the hyperfine relaxation by a single collision is
negligibly small. When the hyperfine relaxation by a single collision is
significant, hyperfine polarization of incident atoms is lost over the surface
dwell time. In this case, a large percentage of the scattered atoms experience
spin relaxation while on the surface; as such, they do not contribute to the
delay-time spectra. Therefore, the probability of spin relaxation resulting
from surface scattering must be estimated prior to measurement of the dwell
time.

To evaluate the spin relaxation resulting from a single collision, we used the
first pump light and the probe light. The first pump light polarizes incident
Rb atoms in the beam to the $F=2$ state. Given that the probe-light-induced
fluorescence of the incident and scattered atoms reflects the number of atoms
in the $F=3$ state, the fluorescence intensity decreases when the first pump
light is introduced. The population fraction $f_{\textrm{i}}$ of the $F=2$
state of the incident atoms pumped by the first pump light can be written as

$\displaystyle f_{\textrm{i}}$ $\displaystyle=$
$\displaystyle\frac{N_{2}}{N_{2}+N_{3}}$ (2) $\displaystyle=$ $\displaystyle
1-\frac{N_{3}}{N},$ (3)

where $N_{2}$ and $N_{3}$ are the numbers of atoms in the $F=2$ and $F=3$
states in the incident atoms, respectively, and $N=N_{2}+N_{3}$. Because the
fluorescence intensity is proportional to the number of atoms in the $F=3$
state,

$\displaystyle N_{3}$ $\displaystyle=$ $\displaystyle CI_{\textrm{i,p}},$ (4)

where $C$ is a constant and $I_{\textrm{i,p}}$ is the intensity of the
fluorescence of incident atoms induced by the first pump light.
$I_{\textrm{i,p}}$ can be obtained by

$\displaystyle I_{\textrm{i,p}}$ $\displaystyle=$ $\displaystyle
i_{-60.3\textrm{ MHz,p}}$ (5) $\displaystyle-i_{160\textrm{
MHz,p}}\times\frac{M(-47\textrm{ m/s},T_{\textrm{s}})}{M(125\textrm{
m/s},T_{\textrm{s}})}$

where $i_{\delta,\textrm{p}}$ is the fluorescence intensity measured with the
first pump light introduced with the probe light blue-detuned by $\delta$ from
the $F=3\rightarrow F^{\prime\prime}=4$ transition frequency, and
$M(v,T_{\textrm{s}})$ is the Maxwell distribution given by

$\displaystyle M(v,T_{\textrm{s}})=\sqrt{\frac{2m}{\pi
k_{\textrm{B}}T_{\textrm{s}}}}\exp\left(-\frac{mv^{2}}{2k_{\textrm{B}}T_{\textrm{s}}}\right),$
(6)

where $m$ is the mass of an 85Rb atom. The first term in Eq. (5) includes
contributions from both incident and scattered atoms. To subtract the latter,
the second term is introduced. The second term is the fluorescence intensity
of the scattered atoms excited by the probe light red-detuned by 60.3 MHz
estimated from the fluorescence intensity measured with the probe light blue-
detuned by 160 MHz (based on the fact that the velocity distribution of the
scattered atoms can be expressed as a Maxwell distribution Sekiguchi  _et al._
(2018)).

Because $F=2$ and $F=3$ states have five- and seven-fold degeneracy, $N_{3}$
and $N$ can be written as

$\displaystyle N$ $\displaystyle=$
$\displaystyle\frac{12}{7}CI_{\textrm{i,np}}.$ (7)

Here, $I_{\textrm{i,np}}$ is the intensity of the fluorescence of incident
atoms without the first pump light, which can be obtained by

$\displaystyle I_{\textrm{i,np}}$ $\displaystyle=$ $\displaystyle
i_{-60.3\textrm{ MHz,np}}$ (8) $\displaystyle-i_{160\textrm{
MHz,np}}\times\frac{M(47\textrm{ m/s},T_{\textrm{s}})}{M(125\textrm{
m/s},T_{\textrm{s}})},$

where $i_{\delta,np}$ is the intensity of fluorescence measured with the probe
light blue-detuned by $\delta$ with the first pump light blocked. Therefore,
we can experimentally determine $f_{\textrm{i}}$ as

$\displaystyle
f_{\textrm{i}}=1-\frac{7}{12}\frac{I_{\textrm{i,p}}}{I_{\textrm{i,np}}}.$ (9)

Similarly, the population fraction $f_{\textrm{s}}$ of the $F=2$ state of the
scattered atoms is written as

$\displaystyle
f_{\textrm{s}}=1-\frac{7}{12}\cdot\frac{I_{\textrm{{s,p}}}}{I_{\textrm{{s,np}}}},$
(10)

where $I_{\textrm{s,p}}$ and $I_{\textrm{s,np}}$ are the intensities of the
fluorescence of the scattered atoms measured with the probe light blue-detuned
by 160 MHz with and without the pump light, respectively. The fluorescence of
the incident atoms is negligible when the probe laser is blue-detuned by 160
MHz, which was confirmed by the fluorescence intensity being lower than the
detection limit when the sample was removed from the atomic beam position. By
comparing $f_{\textrm{i}}$ and $f_{\textrm{s}}$, the proportion of atoms whose
total angular momentum is changed by scattering at the surface can be
estimated.

The delay-time and spin-relaxation measurements were conducted using different
samples prepared by the same procedure, to minimize aggregation of Rb atoms,
which may contaminate the surface. The aggregation of Rb on the surface
resulting from atomic Rb beam irradiation was checked by XPS using Al
k$\alpha$ radiation with a photon energy of 1486.6 eV.

##  III Results and discussion

###  III.1 X-ray photoelectron spectroscopy

Aggregation of Rb atoms on the surface was investigated by XPS. The
temperature at which aggregation was noticeable was approximately 123 K; this
varied slightly among samples. Figure 3 shows the XPS spectra of the as-
prepared sample and the spectra taken after exposing the sample to an atomic
Rb beam with a flux of $10^{11}$–$10^{12}$ s-1 at 153 K for 2 h and 123 K for
2 h. The spectrum of the as-prepared sample displayed a strong peak at
$E-E_{\textrm{F}}=-285$ eV, which was assigned to the C $1s$ state. After
exposure to the Rb beam at low temperature, new peaks appeared at
$E-E_{\textrm{F}}=-531$, $-246$, $-238$, and $-109$ eV; these were assigned to
the O $2p$, Rb $3p_{1/2}$, Rb $3p_{3/2}$, and Rb $3d$ states, respectively.
The existence of the O $2p$ peak indicates that some of the adsorbed Rb atoms
had become oxidized by the residual O2 or H2O in the UHV chamber.

Figure 3: X-ray photoelectron spectroscopy spectra of the clean sample and the
sample exposed to an Rb beam at 153 K for 2 h and 123 K for 2 h.  Figure 4:
(a)Temperature dependence of the $F=2$ population fraction of the incident and
scattered atoms, $f_{\textrm{i}}$ and $f_{\textrm{s}}$, respectively and (b)
the temperature dependence of the ratio $f_{\textrm{s}}/f_{\textrm{i}}$ of the
incident and scattered atoms.

###  III.2 Spin relaxation resulting from surface scattering

The spin relaxation caused by surface scattering was evaluated at 305, 153,
and 123 K. The temperature of the Rb oven of the Rb beam source was set to 393
K. Under these conditions, the flux intensity was estimated to be $10^{11}$ –
$10^{12}$ atoms per second based on the designed value and the fluorescence
induced by the probe light. Below 123 K, the number of scattered atoms was
significantly smaller than at above 123 K, indicating the initiation of Rb
atom adsorption at around 123 K, which is consistent with the XPS results.
Figure 4(a) shows the population fractions of the $F=2$ state for the incident
and scattered atoms. Uncertainties in the population fractions were estimated
by repeating the measurements three to five times. Because the spin relaxation
induced by a single scattering process is negligibly small at room temperature
Bouchiat and Brossel (1966), the difference between $f_{\textrm{s}}$ and
$f_{\textrm{i}}$ at 305 K is attributed to the incident atoms away from the
beam center, as opposed to relaxation due to scattering. When measuring the
population fraction of the $F=2$ state of incident atoms, atoms that pass the
edge of the atomic beam are difficult to pump or detect due to the large
deviation in velocity direction with respect to the major component of the
incident atoms, which the pump and probe light frequencies are tuned to
excite. However, when the probe light is blue-detuned to detect only the
scattered atoms, atoms that were not pumped can be detected, as surface
scattering changes the direction of the translational movement of atoms. The
temperature dependence of the ratio of the population fraction of the incident
and scattered atoms $f_{\textrm{s}}/f_{\textrm{i}}$ is shown in Fig. 4(b); no
increase in spin relaxation probability induced by cooling was observed above
123 K within the experimental error. Given that Rb atoms experience $10^{4}$
collisions before their spins relax Bouchiat and Brossel (1966) in paraffin-
coated cells at room temperature, the low spin relaxation probability at 123 K
means that the mean dwell time at 123 K is smaller than $10^{4}$ times the
mean dwell time at 300 K. Thus, from Eq. (1),
$\tau_{0}\exp\left(\frac{E_{\textrm{des}}}{k_{\textrm{B}}\cdot 123\textrm{
K}}\right)\leq
10^{4}\cdot\tau_{0}\exp\left(\frac{E_{\textrm{des}}}{k_{\textrm{B}}\cdot
305\textrm{ K}}\right)$. By solving this, we obtain $E_{\textrm{des}}\leq
0.16$ eV, which agrees with $E_{\textrm{des}}=0.06$ eV from previous reports
Budker  _et al._ (2005); Rahman and Robinson (1987).

###  III.3 Delay-time spectra

Figure 5: Delay-time spectra (a) at 303 K taken with a high-intensity Rb beam
and long repetition period (b) at 305, 153, and 123 K. Dots with error bars
represent the experimental data and solid lines represent the simulation
results.

Figure 5(a) shows the delay-time spectrum at $T_{\textrm{s}}=303$ K. The
intervals of the second pump light pulses were $2.00\times 10^{-4}$ s. During
the measurement, spectra with and without the second pump light were acquired
by switching the second pump light repeatedly using a shutter. The delay-time
spectrum was obtained by subtracting the latter from the former. For this
measurement, the temperature of the Rb oven of the beam source was set to 453
K, which was 60 K higher than that used for the low-temperature measurement,
to achieve higher signal intensity. As a result, the signal intensity was
enhanced by a factor of $\sim 20$. The feature at 0–$1.0\times 10^{-5}$ s is
attributed to the second pump light, which partially penetrated the
interference filter, and the fluorescence of the incident atoms; this is
excluded from the intensity integration discussed below. The time origin was
defined by the rising edge of the second pump light-derived feature. The
feature peaking at around $2\times 10^{-5}$ s was attributed to the enhanced
fluorescence of scattered atoms caused by the second pump light. The peak area
obtained by integrating the intensity in the region $1.0\times
10^{-5}$–$1.89\times 10^{-4}$ s and subtracting the average level of the
region $1.89\times 10^{-4}$–$1.94\times 10^{-4}$ s as the base level, which
corresponds to the sum of the area of the regions represented by $\alpha$ and
$\beta$ in Fig. 5 (a), accounted for $91\pm 5\%$ of the total signal
intensity. Here, the total signal intensity is the sum of regions $\alpha$,
$\beta$, and $\gamma$. If we adopt the average level of the region $9.3\times
10^{-5}$ –$9.8\times 10^{-5}$ s as the base level, the integrated intensity in
the region $1.0\times 10^{-5}$–$9.3\times 10^{-5}$ s, which corresponds to the
area of region $\alpha$, accounts for $79\pm 2\%$ of the total intensity. This
implies that $79\pm 2\%$ of the scattered atoms will contribute to the peak
intensity in the delay-time spectra, if we regard the sum of the region
$t<9.3\times 10^{-5}$ s as the peak intensity and adopt the average level of
the region $9.3\times 10^{-5}$ –$9.8\times 10^{-5}$ s as the base level. For
the temperature-dependence measurement, we used $1.00\times 10^{-4}$ s as the
second pump light interval and subtracted the average of the region near the
back edge of the time window from the whole spectra, instead of subtracting
the background spectra taken without the second pump light. This dramatically
reduced the measurement time, which was essential to prevent Rb aggregation
during measurements at low temperatures.

Figure 5 (b) shows the delay-time spectra taken at 305, 153, and 123 K. The
temperature of the Rb oven of the Rb beam source was set to 393 K. By cooling
the tetracontane film from 305 to 123 K, the delay-time spectrum shifted to
the longer side and the mean delay time $\tau_{\textrm{M}}$ increased by
$(7.0\pm 3.2)\times 10^{-6}$ s from $(2.54\pm 0.21)\times 10^{-5}$ to
$(3.24\pm 0.24)\times 10^{-5}$ s. Here, $\tau_{\textrm{M}}$ is defined by

$\displaystyle\tau_{\textrm{M}}=\frac{\sum_{i}t_{i}I_{i}}{\sum_{i}I_{i}},$
(11)

where $t_{i}$ and $I_{i}$ are the delay time and the intensity at the $i$th
point, respectively. The uncertainty in $\tau_{\textrm{M}}$ originates from
the uncertainty in $I_{i}$ at each point, which is $\sqrt{I_{i}}$. Data points
in the region $6\times 10^{-6}\leq t_{i}<9.3\times 10^{-5}$ s were included in
the summation. The average of the region $9.5\times 10^{-5}\textrm{ s}\leq
t_{i}<9.9\times 10^{-5}\textrm{ s}$ was adopted as the base level. Intensities
between 0 and $6\times 10^{-6}$ s, which include the peak originating from the
second pump light, were not included in the summation. Because the velocity
distribution of the scattered beam also depends on the film temperature
Sekiguchi  _et al._ (2018), we cannot simply attribute the increase in
$\tau_{\textrm{M}}$ to the increase in mean dwell time. To evaluate the
increase in TOF due to the change in velocity distribution, we simulated TOF
spectra without taking the dwell time into account.

The simulation considered the TOF from the second pump light to the surface,
and from the surface to the probe light. The velocity distribution
$d_{\textrm{b}}(v)$ of the incident atoms, and the TOF distribution from the
second pump light to the surface $D_{\textrm{b}}(t,L_{1})$, were calculated
using equations Sekiguchi  _et al._ (2018)

$\displaystyle d_{\textrm{b}}(v)$ $\displaystyle=$
$\displaystyle\frac{m^{2}}{2k_{\textrm{B}}^{2}T_{\textrm{b}}^{2}}v^{3}\exp\left(-\frac{mv^{2}}{2k_{\textrm{B}}T_{\textrm{b}}}\right),$
(12) $\displaystyle D_{\textrm{b}}(t,L_{1})$ $\displaystyle=$ $\displaystyle
d_{\textrm{b}}\left(\frac{L_{1}}{t}\right)\frac{\textrm{d}}{\textrm{dt}}\left(\frac{L_{1}}{t}\right)$
(13) $\displaystyle=$
$\displaystyle\frac{m^{2}L_{1}^{4}}{2k_{\textrm{B}}^{2}T_{\textrm{b}}^{2}t^{5}}\exp\left(-\frac{m}{2k_{\textrm{B}}T_{\textrm{b}}}\left(\frac{L_{1}}{t}\right)^{2}\right),$
(14)

where $T_{\textrm{b}}$ is the temperature of the incident beam determined by
the temperature of the capillary of the beam source, $v$ is the velocity of
atoms, $L_{1}$ is the distance between the second pump light and the surface
along the atomic beam direction, and $t$ is the time. $T_{\textrm{b}}$ was 453
K, and $L_{1}$ was roughly estimated to be $1.8\times 10^{-3}$ m. The
distribution $d_{\textrm{s}}(v_{\perp s})$ of the velocity component of the
scattered atoms perpendicular to the surface, which are in thermal equilibrium
with the film, and the TOF distribution $D_{\textrm{s}}(t,L_{2})$ are given by

$\displaystyle d_{\textrm{s}}(v_{\perp\textrm{s}})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2m}{\pi
k_{\textrm{B}}T_{\textrm{s}}}}\exp\left(-\frac{mv_{\perp\textrm{s}}^{2}}{2k_{\textrm{B}}T_{\textrm{s}}}\right),$
(15) $\displaystyle D_{\textrm{s}}(t,L_{2})$ $\displaystyle=$ $\displaystyle
d_{\textrm{s}}\left(\frac{L_{2}}{t}\right)\frac{\textrm{d}}{\textrm{dt}}\left(\frac{L_{2}}{t}\right)$
(16) $\displaystyle=$ $\displaystyle\sqrt{\frac{2m}{\pi
k_{\textrm{B}}T_{\textrm{s}}}}\frac{L_{2}}{t^{2}}\exp\left(-\frac{m}{2k_{\textrm{B}}T_{\textrm{s}}}\left(\frac{L_{2}}{t}\right)^{2}\right),$
(17)

where $v_{\perp\textrm{s}}$ is the velocity component perpendicular to the
film surface, and $L_{2}$, which was estimated to be $1.42\times 10^{-3}$ m,
is the height of the probe light from the film surface Sekiguchi  _et al._
(2018). The total TOF spectrum $S(t)$ is given by

$\displaystyle S(t)$ $\displaystyle=$
$\displaystyle\frac{1}{w}\int_{t-w}^{t}dr\int_{0}^{r}ds\textrm{
}u_{\textrm{b}}(s,L_{1})$ (18) $\displaystyle\times
u_{\textrm{s}}(r-s,L_{\textrm{s}}),$

where

$\displaystyle u_{\textrm{b}}(t,L_{1})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{1}{\sqrt{2\pi\sigma_{1}^{2}}}\exp\left(-\frac{(x-L_{1})^{2}}{2\sigma_{1}^{2}}\right)$
(19) $\displaystyle\times$ $\displaystyle D_{\textrm{b}}(t,x)dx,$
$\displaystyle u_{\textrm{s}}(t,L_{2})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{1}{\sqrt{2\pi\sigma_{2}^{2}}}\exp\left(-\frac{(x-L_{2})^{2}}{2\sigma_{2}^{2}}\right)$
(20) $\displaystyle\times$ $\displaystyle g_{\textrm{s}}(t,x)dx,$

where $w$ is the duration of the second pump light and $\sigma_{1}$ and
$\sigma_{2}$ are the $1/\sqrt{e}$ half width of the second pump light and
probe light, respectively. $\sigma_{1}$ and $\sigma_{2}$ were $2.4\times
10^{-4}$ and $1.7\times 10^{-4}$ m, respectively.

The simulation results are indicated by solid lines in Fig. 5. The intensities
of the simulation results were adjusted to fit the experimental results. The
simulated TOF spectra are in good agreement with the experimental results.
According to the simulation results at 305 K, the peak area calculated by
integrating the spectral intensity between $1.0\times 10^{-5}$ and $9.3\times
10^{-5}$ s and subtracting the average of the region $9.3\times
10^{-5}$–$9.8\times 10^{-5}$ s accounts for 81% of the total intensity, which
is in good agreement with the experimental results shown in Fig. 5(a). The
simulated mean TOF calculated from the region 0–9.3 $\times 10^{-5}$ s
increased by $5.8\times 10^{-6}$ s from $2.48\times 10^{-5}$ to $3.06\times
10^{-5}$ s with cooling from 305 to 123 K. Here, the average of the region
$9.5\times 10^{-5}\textrm{ s}\leq t_{i}<9.9\times 10^{-5}\textrm{ s}$ was
adopted as the base level. The experimentally observed shift of
$\tau_{\textrm{M}}$, which is $(7.0\pm 3.2)\times 10^{-6}$ s, is the sum of
the increase in mean dwell time and mean TOF. Therefore, the increase in mean
dwell time induced by cooling can be obtained by subtracting the increase in
the simulated mean TOF from the experimentally obtained increase in the mean
delay time. Therefore, we can see from the experimental and simulation results
that

$\displaystyle\tau_{\textrm{s,$t<93$ \textmu s}}(123\textrm{
K})-\tau_{\textrm{s,$t<93$ \textmu s}}(305\textrm{ K})$ $\displaystyle=(1.2\pm
3.2)\times 10^{-6}\textrm{ s},$ (21)

which means

$\displaystyle\tau_{\textrm{s,$t<93$ \textmu s}}(123\textrm{
K})-\tau_{\textrm{s,$t<93$ \textmu s}}(305\textrm{ K})$ $\displaystyle\leq
4.4\times 10^{-6}\textrm{ s.}$ (22)

Here, $\tau_{\textrm{s},t<t_{\textrm{max}}}(T_{\textrm{s}})$ is the mean dwell
time of the scattering component, with a delay time of less than
$t_{\textrm{max}}$ at a surface temperature of $T_{\textrm{s}}$.

###  III.4 Discussion

The discrepancy between our result ($\tau_{\textrm{s}}(123\textrm{
K})-\tau_{\textrm{s}}(305\textrm{ K})\leq 4.4\times 10^{-6}$ s) and the value
($\tau_{\textrm{s}}(123\textrm{ K})-\tau_{\textrm{s}}(305\textrm{
K})=6.7\times 10^{-5}$ s) obtained by substituting the desorption energy
$E_{\textrm{des}}=0.06$ eV Budker  _et al._ (2005); Rahman and Robinson (1987)
and mean dwell time $\tau_{\textrm{s}}=1.8\times 10^{-6}$ s at $T=345$ K
Ulanski and Wu (2011) into Eq. (1) may be explained by assuming multiple
scattering components with different mean dwell times. Using the method
described in this study, the scattering components with dwell times larger
than $9.3\times 10^{-5}$ s is not detected. In Ref. Ulanski and Wu (2011), on
the other hand, the mean dwell time was estimated based on the Larmor
frequency shift caused by the interaction with the surface and evanescent pump
light.

The reason for the difference in mean dwell times between the scattering
components can be attributed to differences in pre-exponential factors. It has
been reported that a certain proportion of the incident atoms penetrate the
PDMS film, diffuse into the bulk, and desorb from the surface Atutov and
Plekhanov (2015), which makes the mean dwell time about a million times larger
than that calculated from the desorption energy and film temperature. If the
diffusion barrier in the bulk is significantly smaller than the desorption
energy, the temperature dependence of the diffusion time can be neglected so
that the temperature dependence of the mean dwell time is almost entirely
determined by the desorption energy. By assuming that two scattering
components with different mean dwell times exist, we can approximate the
temperature dependence of the mean dwell time as

$\displaystyle\tau_{\textrm{s}}=(1-p)\tau_{1}\exp\left(\frac{E_{\textrm{des}}}{k_{\textrm{B}}T_{\textrm{s}}}\right)+p\tau_{2}\exp\left(\frac{E_{\textrm{des}}}{k_{\textrm{B}}T_{\textrm{s}}}\right),$
(23)

where $p$ is the proportion of scattering events with longer mean dwell times,
and $\tau_{1}$ and $\tau_{2}$ are the pre-exponential factors for the
scattering events with shorter and longer mean dwell times, respectively. We
suppose that
$\tau_{2}\exp\left(\frac{E_{\textrm{des}}}{k_{\textrm{B}}T_{\textrm{s}}}\right)$
is significantly larger than the time window of $9.3\times 10^{-5}$ s and only
the component with a shorter mean dwell time, which corresponds to the first
term in Eq. (23), contributes to the delay-time spectra. By substituting Eq.
(1) into our results ($\tau_{\textrm{s}}(123\textrm{
K})-\tau_{\textrm{s}}(305\textrm{ K})\leq 4.4\times 10^{-6}$ s), we obtain
$0<\tau_{1}\leq 1.6\times 10^{-8}$ s. From
$\tau_{2}\exp\left(\frac{E_{\textrm{des}}}{k_{\textrm{B}}T_{\textrm{s}}}\right)\gg
9.3\times 10^{-5}$ s at $T_{\textrm{s}}\leq 305$ K, we get $\tau_{2}\gg
9.47\times 10^{-6}$ s. By substituting $\tau_{\textrm{s}}=1.8\times 10^{-6}$ s
at 345 K Ulanski and Wu (2011) and $E_{\textrm{des}}=0.06$ eV Rahman and
Robinson (1987); Budker  _et al._ (2005) into Eq. (23), $p=\frac{2.4\times
10^{-7}\textrm{s}-\tau_{1}}{\tau_{2}-\tau_{1}}$. From $0<\tau_{1}$ and
$\tau_{2}\gg 9.47\times 10^{-6}$ s, $p<0.025$, which means that the component
with a shorter mean dwell time is the major component. This is consistent with
the fact that the observed $79\pm 2$ % fraction within $9.3\times 10^{-5}$ s
is nearly the same as the 81 % fraction obtained from the simulation without
dwell times.

##  IV Conclusions

Scattering of Rb atoms on tetracontane surfaces was investigated. No
significant spin relaxation was observed with a single scattering process down
to 123 K. The temperature evolution of delay time showed that the increase in
mean surface dwell time induced by cooling from 305 to 123 K was less than
$4.4\times 10^{-6}$ s. Taken together, the results indicate the existence of
multiple scattering sites. The pre-exponential factor $\tau_{0}$ of the minor
components is at least three orders of magnitude larger than that of the major
component, which means that the mean dwell time of the minor scattering
components is at least three orders of magnitude larger than that of the major
component.

##  V Acknowledgments

This work was supported by JSPS KAKENHI Grant Number JP17H02933.

## References

  * Risley  _et al._ (1980) A. Risley, S. Jarvis Jr, and J. Vanier, Journal of Applied Physics 51, 4571 (1980).
  * Robinson and Johnson (1982) H. Robinson and C. Johnson, Applied Physics Letters 40, 771 (1982).
  * Frueholz  _et al._ (1983) R. Frueholz, C. Volk, and J. Camparo, Journal of Applied Physics 54, 5613 (1983).
  * Budker  _et al._ (1998) D. Budker, V. Yashchuk, and M. Zolotorev, Physical Review Letters 81, 5788 (1998).
  * Balabas  _et al._ (2006) M. Balabas, D. Budker, J. Kitching, P. Schwindt, and J. Stalnaker, Journal of the Optical Society of America B 23, 1001 (2006).
  * Wasilewski  _et al._ (2010) W. Wasilewski, K. Jensen, H. Krauter, J. J. Renema, M. Balabas, and E. S. Polzik, Physical Review Letters 104, 133601 (2010).
  * Robinson  _et al._ (1958) H. Robinson, E. Ensberg, and H. Dehmelt, Bulletin of the American Physical Society 3 (1958).
  * Bouchiat and Brossel (1966) M. A. Bouchiat and J. Brossel, Physical Reveview 147, 41 (1966).
  * Seltzer  _et al._ (2007) S. Seltzer, P. Meares, and M. Romalis, Physical Review A 75, 051407 (2007).
  * Zhao  _et al._ (2008) K. Zhao, M. Schaden, and Z. Wu, Physical Review A 78, 034901 (2008).
  * Atutov  _et al._ (2016) S. Atutov, F. Benimetskiy, A. Plekhanov, and V. Sorokin, The European Physical Journal D 70, 1 (2016).
  * Rahman and Robinson (1987) C. Rahman and H. Robinson, IEEE journal of quantum electronics 23, 452 (1987).
  * Budker  _et al._ (2005) D. Budker, L. Hollberg, D. F. Kimball, J. Kitching, S. Pustelny, and V. V. Yashchuk, Physical Review A 71, 012903 (2005).
  * Ulanski and Wu (2011) E. Ulanski and Z. Wu, Applied Physics Letters 98, 201115 (2011).
  * Atutov and Plekhanov (2015) S. Atutov and A. Plekhanov, Journal of Experimental and Theoretical Physics 120, 1 (2015).
  * Zhao  _et al._ (2009) K. Zhao, M. Schaden, and Z. Wu, Physical Review Letters 103, 073201 (2009).
  * Wieman and Hänsch (1976) C. Wieman and T. W. Hänsch, Physical Review Letters 36, 1170 (1976).
  * Harris  _et al._ (2006) M. Harris, C. Adams, S. Cornish, I. McLeod, E. Tarleton, and I. Hughes, Physical Review A 73, 062509 (2006).
  * Arimondo  _et al._ (1977) E. Arimondo, M. Inguscio, and P. Violino, Reviews of Modern Physics 49, 31 (1977).
  * Schultz  _et al._ (2008) B. Schultz, H. Ming, G. Noble, and W. Van Wijngaarden, The European Physical Journal D 48, 171 (2008).
  * Sekiguchi  _et al._ (2018) N. Sekiguchi, A. Hatakeyama, K. Okuma, and H. Usui, Physical Review A 98, 042709 (2018).
  * Volz and Schmoranzer (1996) U. Volz and H. Schmoranzer, Physica Scripta 1996, 48 (1996).