# Testing Gravity with Frequency-Dependent Overlap Reduction Function in
Pulsar Timing Array

Qiuyue Liang<EMAIL_ADDRESS>Ippei Obata<EMAIL_ADDRESS>Misao
Sasaki<EMAIL_ADDRESS>Kavli Institute for the Physics and Mathematics
of the Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan

###### Abstract

The positive evidence of a nano-hertz gravitational wave background recently
found by several pulsar timing array (PTA) collaborations opened up a window
to test modified gravity theories in a unique frequency band in parallel to
other gravitational wave detection experiments. In particular, the overlap
reduction function (ORF) in PTA observation is sensitive to the phase velocity
of gravitational waves. In this work, we provide analytical expressions for
the coefficients of the multipole moments in the ORF, and utilize these
analytical results to study constraints on the phase velocity from the
frequency dependent overlap reduction function obtained from the Chinese PTA
(CPTA) data. While the data contain large error bars yet, interesting
constraints are found in the frequency-dependent ORF in the case of subluminal
phase velocity. This makes us expect that the nano-hertz band gravitational
wave background will become one of the important arenas for exploring modified
gravity theories.

††preprint: YITP-24-55

##  I Introduction

Pulsars are compact neutron stars emitting pulsed lights with very stable
rotation periods Hewish  _et al._ (1968). Thanks to this feature, their pulse
arrival time measurements can test the fluctuations in the gravitational field
at frequencies on the order of 1-100nHz Detweiler (1979) with long temporal
baseline data. Moreover, we can construct the correlation function among the
pulsar pairs in the system to study the underlying gravitational wave
background (GWB). For an isotropic GWB, we can split the signal into power
spectrum and the overlap reduction function, which is a characteristic angular
correlation function between pulsar pairs, dubbed the Hellings-Downs (HD)
curve Hellings and Downs (1983). Recently, an evidence of a nano-hertz GWB in
the pulsar timing array signals, over 3 $\sigma$ statistical significance, has
been reported by NANOGrav, PPTA, EPTA, and CPTA Agazie  _et al._ (2023); Afzal
_et al._ (2023); Antoniadis  _et al._ (2023); Reardon  _et al._ (2023); Xu
_et al._ (2023). This unique probe opened up the possibility to test
fundamental theories, including modifications to general relativity.

Constraining modified gravity theories with extra non-Einstein polarization
modes and/or modified dispersion relations by using PTA has been well
discussed in the literature Lee  _et al._ (2010); Gair  _et al._ (2015); Qin
_et al._ (2021); Liang and Trodden (2021); Liang  _et al._ (2023); Bernardo
and Ng (2023a); Bi  _et al._ (2024); Bernardo and Ng (2023b); Wu  _et al._
(2024, 2023). Constraining the graviton mass is, among others, of particular
research interest Finn and Sutton (2002); Abbott  _et al._ (2016); Liang and
Trodden (2021); Abbott  _et al._ (2021); Bernardo and Ng (2023c); Wang and
Zhao (2023); Wu  _et al._ (2024); Bernardo and Ng (2023b). A modified
dispersion relation of gravitational waves will give rise to a deviation in
the overlap reduction function from the HD curve. As pointed out in Liang  _et
al._ (2023), the phase velocity other than the group velocity
$v_{p}=\omega(k)/k$ is the parameter entering the observable, where
$\omega(k)$ and $k$ are the angular frequency (or energy) and the wavenumber
(or momentum) of the gravitational wave, respectively. While we have $v_{p}>1$
in canonical massive gravity, theories with $v_{p}<1$, which we refer to as
subluminal gravity, can naturally arise from various effective field theories
of gravity Carrillo Gonzalez  _et al._ (2022); Ezquiaga  _et al._ (2021); de
Rham and Tolley (2020). In this paper, we derive analytical expressions for
the coefficients of the multipole moments of ORF for both $v_{p}>1$ and
$v_{p}<1$ cases.

It is of great interest to apply these analytical expressions to the current
data sets and see if we can obtain constraints on the phase velocity. In this
respect, as one generically expects that the dispersion relation is frequency
dependent, CPTA offers a unique opportunity in that they report an analysis of
the overlap reduction function based on three different frequency bins, in
contrast to those by the other experimental groups that provide only
frequency-averaged data analyses. In this paper, as an application of our
analytical formula, we discuss the frequency-dependent angular correlation
obtained from CPTA. We mention, however, CPTA has been operating only for 3.4
years, and only one frequency bin, $f=1.5/(3.4\,{\rm yrs})$, shows significant
statistical evidence of a HD signal. Therefore, we will focus on this specific
frequency, while hoping for the possibility that future PTA observations will
eventually offer frequency dependent data at different frequency bins with
high statistical significance that will enable us to quantitatively test
modified dispersion relations.

The paper is organized as follows. In Sec. II we briefly review the formalism
of the decomposition of ORF into Legendre polynomials in the PTA system. In
Sec. III, we give a full analytical expression for the coefficients of these
Legendre polynomials. It is based on the assumption that distances to the
pulsars are much longer than the wavelength. Then we also study the finite
distance effect and find that its effect is non-perturbative in the subluminal
case. In Sec. IV, we study the constraints on the phase velocity from CPTA
data. Throughout the paper we use the metric signature $(-,+,+,+)$, and set
$\hbar=c=1$.

##  II Introduction on the pulsar timing array system

In this section, we review how to obtain a gravitational wave background
signal in the PTA system. As mentioned above, each pulsar can be treated as a
standard clock, and their pulse frequency variation induced by the underlying
metric fluctuations will accumulate in the residual of arrival time, which is
the observable for each pulsar. It is defined as

$R(t,\hat{p})\equiv\int_{0}^{t}dt^{\prime}\,\left(\frac{\nu_{0}(\hat{p})-\nu(t^{\prime},\hat{p})}{\nu_{0}(\hat{p})}\right)=\int_{0}^{t}dt^{\prime}z(t^{\prime},\hat{p})\
,$ (1)

where $\nu$ is the pulse frequency, and the integrand is the difference in the
frequencies propagating in the direction of the unit vector $\hat{p}$ due to
the presence of gravitational waves, which we will refer to as the redshift,
$z(t,\hat{p})$. By measuring the arrival time residual of each individual
pulsar, one can construct the two point correlation function of the pulsar
pairs.

We evaluate the redshift by solving the geodesic equation for pulsar light. In
the presence of gravitational waves $h_{ij}$, it is given by

$\dfrac{d\nu}{d\lambda}\simeq\frac{\nu^{2}}{2}\dfrac{dh_{ij}}{dt}\hat{p}^{i}\hat{p}^{j}\
,$ (2)

where $\lambda$ is the affine parameter defined such that $dt/d\lambda$ is
equal to the frequency, $\nu$. To solve this, we assume that the gravitational
wave signal is approximated as a plane wave solution
$h_{ij}(t,\vec{x})=h_{ij}(t-\hat{\Omega}\cdot\vec{x}/v_{p}(k))$, where
$\hat{\Omega}$ is a unit vector from which gravitational waves propagate and
$v_{p}(k)=\omega(k)/k$ is the phase velocity with respect to the angular
frequency of gravitational waves, $\omega(k)$. Then, in terms of the following
relationship of time and space variables,

$t=\nu_{0}(\lambda-\lambda_{e})+t_{e}\ ,\ \
\hat{\Omega}\cdot\vec{x}=\hat{\Omega}\cdot\vec{x}_{e}-\nu_{0}(\lambda-\lambda_{e})\hat{\Omega}\cdot\hat{p}\,,$
(3)

between the pulsar and earth (labeled with $e$), (2) is solved as

$\log\left(\dfrac{\nu}{\nu_{0}}\right)=\dfrac{v_{p}~{}\hat{p}_{i}\hat{p}_{j}}{2\left(v_{p}+\hat{\Omega}\cdot\hat{p}\right)}\Delta
h_{ij}\rightarrow~{}z(t,\hat{p})\simeq-\dfrac{v_{p}~{}\hat{p}_{i}\hat{p}_{j}}{2\left(v_{p}+\hat{\Omega}\cdot\hat{p}\right)}\Delta
h_{ij}\ ,$ (4)

where

$\Delta h_{ij}\equiv
h_{ij}\left(t_{p}-\dfrac{\hat{\Omega}\cdot\vec{x}_{p}}{v_{p}}\right)-h_{ij}\left(t_{e}-\dfrac{\hat{\Omega}\cdot\vec{x_{e}}}{v_{p}}\right)\
,$ (5)

is the difference of the gravitational wave at pulsar and at earth.

For a gravitational wave background, $z(t,\hat{p})$ needs to be averaged over
the direction from which gravitational waves propagate. It reads

$\displaystyle z(t,\hat{p})$ $\displaystyle=$
$\displaystyle\int\frac{dk}{2\pi}d^{2}\hat{\Omega}~{}\frac{p^{i}p^{j}h_{ij}(k,t_{0})}{2}e^{i\omega(k)(t-t_{0})}\frac{v_{p}(k)}{\left(v_{p}(k)+\hat{\Omega}\cdot\hat{p}\right)}\left(1-e^{-ikL\left(v_{p}(k)+\hat{\Omega}\cdot\hat{p}\right)}\right)\
.$ (6)

Notice that we apply the Fourier transformation in momentum $k$ instead of
frequency, which can be related back through $2\pi f=\omega(k)$. The time
residual is given by

$\displaystyle R(t,\hat{p})=\int_{0}^{t}z(t^{\prime})dt^{\prime}$ (7)
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{dk}{2\pi}\int_{S^{2}}d^{2}\hat{\Omega}~{}\frac{p^{i}p^{j}h_{ij}(k,t_{0})}{-2i\omega(k)}\left(e^{i\omega(k)t}-1\right)\frac{v_{p}(k)}{\left(v_{p}(k)+\hat{\Omega}\cdot\hat{p}\right)}\left(1-e^{-ikL\left(v_{p}(k)+\hat{\Omega}\cdot\hat{p}\right)}\right)\
,$

where we have set the emission time of the gravitational wave, $t_{0}$, to be
zero. We can read off the Fourier transformation, $\tilde{R}(k,\hat{p})$, as

$\displaystyle\tilde{R}(k,\hat{p})$ $\displaystyle=$
$\displaystyle\frac{\hat{p}^{i}\hat{p}^{j}}{-2i\omega(k)}\int_{S^{2}}d^{2}\hat{\Omega}~{}\frac{v_{p}(k)}{\left(v_{p}(k)+\hat{\Omega}\cdot\hat{p}\right)}\left(1-e^{-ikL\left(v_{p}+\hat{\Omega}\cdot\hat{p}\right)}\right)h_{ij}\left(k\right)$
(8) $\displaystyle=$
$\displaystyle\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}\left(R^{E}_{(\ell
m)}(k,\hat{p})a_{(\ell m)}^{E}(k)+R^{B}_{(\ell m)}(k,\hat{p})a_{(\ell
m)}^{B}(k)\right)$ $\displaystyle=$
$\displaystyle\int_{S^{2}}d^{2}\hat{\Omega}\left(R^{+}(k,\hat{p},\hat{\Omega})h^{+}(k,\hat{\Omega})+R^{\times}(k,\hat{p},\hat{\Omega})h^{\times}(k,\hat{\Omega})\right)$

where $R^{E,B}_{(\ell m)}(k,\hat{p})$ and
$R^{+,\times}(k,\hat{p},\hat{\Omega})$ are the response functions.
$a_{(lm)}^{E,B}$ and $h^{+,\times}(k,\hat{\Omega})$ are the amplitudes of
gravitational waves in spherical harmonics and polarization basis,
respectively. We defer details of the basis and coefficients transformation
Gair  _et al._ (2014); Liang  _et al._ (2023) to Appendix (A).

The two point correlation functions of the arrival time $R(k,\hat{p})$ can be
obtained through

$\displaystyle\braket{R(t,\hat{p}_{1})R(t^{\prime},\hat{p}_{2})}=\int_{-\infty}^{\infty}\frac{dk}{2\pi}e^{i\omega(k)(t-t^{\prime})}H(k)\Gamma(k,\xi)\
,$ (9)

where $H(k)$ is the power spectrum of the GWB defined by

$\displaystyle\braket{h^{+}(k,\hat{\Omega})h^{+*}(k^{\prime},\hat{\Omega}^{\prime})}=\braket{h^{\times}(k,\hat{\Omega})h^{\times*}(k^{\prime},\hat{\Omega}^{\prime})}=\pi
H(k)\delta(k-k^{\prime})\delta^{2}(\hat{\Omega}-\hat{\Omega}^{\prime})\,,$
(10)

or

$\displaystyle\left\langle a_{(\ell
m)}^{E}(k)a_{\left(\ell^{\prime}m^{\prime}\right)}^{E*}\left(k^{\prime}\right)\right\rangle=\left\langle
a_{(\ell
m)}^{B}(k)a_{\left(\ell^{\prime}m^{\prime}\right)}^{B*}\left(k^{\prime}\right)\right\rangle=2\pi\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}H(k)\delta(k-k^{\prime})\,,$
(11)

and $\Gamma(k,\xi)$ is the momentum-dependent overlap reduction function (ORF)
of the angle between pairs of pulsars
$\xi\equiv\cos^{-1}(\hat{p}_{1}\cdot\hat{p}_{2})$, which contains information
about the angular correlation.

In polarization basis, the overlap reduction function takes the form,

$\displaystyle\Gamma(k,\xi)=\beta\int_{S^{2}}d^{2}\Omega\sum_{A=+,\times}R^{A}(k,\hat{p}_{1},\hat{\Omega})R^{A*}(k,\hat{p}_{2},\hat{\Omega})\
,$ (12)

where $\beta$ is a normalization factor, and in harmonic basis, it can be
expressed as

$\displaystyle\Gamma(k,\xi)$ $\displaystyle=$
$\displaystyle\mathcal{C}\sum_{\ell=2}^{\infty}\Gamma_{12,\ell}(k,\xi)=\mathcal{C}\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}\sum_{P=E,B}R_{(\ell
m)}^{P}(k,\hat{p}_{1})R_{(\ell m)}^{P*}(k,\hat{p}_{2})$ (13) $\displaystyle=$
$\displaystyle\mathcal{C}\sum_{\ell=2}^{\infty}(2\ell+1)\frac{2(\ell-2)!}{(\ell+2)!}|c_{\ell}(k)|^{2}P_{\ell}(\cos\xi)\equiv\mathcal{C}\sum_{\ell=2}^{\infty}a_{\ell}P_{\ell}(\cos\xi)\
,$

where $\mathcal{C}$ is the overall normalization factor, and the coefficients
can be computed as

$\displaystyle c_{\ell}(k)$ $\displaystyle=$
$\displaystyle\int_{-1}^{1}dx\frac{v_{p}(k)}{v_{p}(k)+x}\left[1-e^{-ikL\left(v_{p}+x\right)}\right](1-x^{2})^{2}\frac{d^{2}}{dx^{2}}P_{\ell}(x)\
.$ (14)

In both GR and massive gravity (superluminal phase velocity) cases, the
exponential factor in Eq. (14) can be dropped to a good approximation in the
PTA system, where the typical distances between pulsar and earth are kpc, and
the frequency band is 1-100 nHz. The coefficients without the exponential term
give $c_{\ell}=4(-1)^{\ell}$ in GR, which will give the well-known HD curve
when substituted to Eq. (13). However, as pointed out in Liang  _et al._
(2023), the subluminal phase velocity case has to be treated differently. If
we neglect the exponential term, the integrand is ill-defined as there appears
a singularity at $x=-v_{p}$, corresponding to a “surfing” angle
$\theta=-\cos^{-1}(-v_{p})$. Integrating this exponential factor is
computationally expensive, but necessary to compare with data to extract
useful constraints. In the next section, we present an analytical expression
for these coefficients.

##  III Computing the Overlap Reduction Function Beyond GR: Examples

The main result of this paper is to give an analytical expression for the
coefficients $c_{\ell}$ in Eq. (14). After some manipulations, we obtain

$\displaystyle
c_{\ell}=\left\\{\begin{array}[]{@{}ll@{}}-2v_{p}(1+\ell)\left(((2+\ell)v_{p}^{2}-\ell)Q_{\ell}(-v_{p})+2v_{p}Q_{\ell+1}(-v_{p})\right)+\mathcal{O}\left(\frac{1}{kL}\right);{}&v_{p}(k)>1\,,\\\
4(-1)^{\ell}+\mathcal{O}\left(\frac{1}{kL}\right);{}&v_{p}(k)=1\,,\\\
-2v_{p}(1+\ell)\left(((2+\ell)v_{p}^{2}-\ell)\tilde{Q}_{\ell}(-v_{p})+2v_{p}\tilde{Q}_{\ell+1}(-v_{p})\right)\\\
~{}+\pi
iv_{p}(1+\ell)\left(((2+\ell)v_{p}^{2}-\ell)P_{\ell}(-v_{p})+2v_{p}P_{\ell+1}(-v_{p})\right)+\mathcal{O}\left(\frac{1}{kL}\right);{}&v_{p}(k)<1\,,\end{array}\right.$
(19)

where $P_{\ell}(x)$ is the Legendre P polynomial, and $Q_{\ell}$ and
$\tilde{Q}_{\ell}$ are Legendre Q functions defined as

$\displaystyle
Q_{\ell}(z)=-\frac{1}{2}P_{\ell}(z)\log\frac{z+1}{z-1}+\frac{1}{\ell!2^{\ell}}\frac{d^{\ell}}{dz^{\ell}}\left((z^{2}-1)^{\ell}\log\frac{z+1}{z-1}\right)\quad\text{for}\
z\notin[-1,1]\ ,$
$\displaystyle\tilde{Q}_{\ell}(x)=-\frac{1}{2}P_{\ell}(x)\log\frac{1+x}{1-x}+\frac{1}{\ell!2^{\ell}}\frac{d^{\ell}}{dx^{\ell}}\left((x^{2}-1)^{\ell}\log\frac{1+x}{1-x}\right)\quad\text{for}\
x\in[-1,1]\ .$ (21)

The difference between the two Legendre Q functions is the choice of the
branch cut. In our definition, $Q_{\ell}(z)$ has branch cut over $[-1,1]$,
while $\tilde{Q}_{\ell}(x)$ has branch cut on $(-\infty,-1)\cup(1,\infty)$.
Note that the choice of the Legendre Q functions in (19) is consistent with
the fact that the real part of $c_{\ell}$ must be real. We save the detail of
the algebraic derivation to Appendix (B).

In Fig. 1, a comparison of the coefficients $a_{\ell}$ defined in Eq. (13) is
made between the numerical and analytical results for both massive gravity and
subluminal cases. To show the finite distance effect in the plot, we
exaggeratedly chose the parameter $kL=6$, which is usually the order of $1000$
in a real PTA system. One can see that the two match well. The numerical
result starts to deviate from the analytical expression for $\ell>kL$. Thus,
we can claim that the analytical expression for the coefficients of the
overlap reduction function in harmonic decomposition is a good approximation
to the response function for PTA measurements.

Figure 1: A comparison of the coefficients of the first $8$ Legendre P
polynomials $a_{\ell}$ defined in Eq. (13) obtained from numerical integration
of Eq. (14) and from the analytical expression given in Eq. (19). The solid
lines represent the analytical results where the infinite distance limit
$kL\to\infty$ is taken. The red, blue and green colors represent the cases of
subluminal phase velocity, general relativity, and massive gravity,
respectively. The dashed lines are the numerical results, assuming $kL=6$ to
exaggerate the finite distance effect. In a real system, $kL$ is usually an
order of thousands, therefore, the deviation from our analytical expression to
the numerical integration can be neglected.

It is worthwhile pointing out that the imaginary part of $c_{\ell}$ is non-
vanishing when $v_{p}<1$. In this case, we can express the real and imaginary
parts of the coefficients at the large $\ell$ limit as

$\displaystyle\lim_{\ell\to\infty}\Re\,c_{\ell}$ $\displaystyle~{}=-\pi
v_{p}\sqrt{\frac{\theta}{\sin\theta}}(\ell+1)\left((-\ell+(2+\ell)v_{p}^{2})Y_{0}\left(\left(\ell+\frac{1}{2}\right)\theta\right)+2v_{p}Y_{0}\left(\left(\ell+\frac{3}{2}\right)\theta\right)\right)\,,$
(22) $\displaystyle\lim_{\ell\to\infty}\Im\,c_{\ell}$ $\displaystyle~{}=\pi
v_{p}\sqrt{\frac{\theta}{\sin\theta}}(\ell+1)\left((-\ell+(2+\ell)v_{p}^{2})J_{0}\left(\left(\ell+\frac{1}{2}\right)\theta\right)+2v_{p}J_{0}\left(\left(\ell+\frac{3}{2}\right)\theta\right)\right)\,,$
(23)

where $\cos\theta=-v_{p}$, and $J_{0}$ and $Y_{0}$ are the Bessel functions of
the 0-th order. Here, we have used the asymptotic behavior of the Legendre
functions for large arguments. Notice that for
$c_{\ell}=\Re\,c_{\ell}+i\,\Im\,c_{\ell}$, we have

$\displaystyle\lim_{\ell\to\infty}c_{\ell}$ $\displaystyle~{}=i\pi
v_{p}\sqrt{\frac{\theta}{\sin\theta}}(\ell+1)\left((-\ell+(2+\ell)v_{p}^{2})H^{(1)}_{0}\left(\left(\ell+\frac{1}{2}\right)\theta\right)+2v_{p}H^{(1)}_{0}\left(\left(\ell+\frac{3}{2}\right)\theta\right)\right),$
(24)

where $H^{(1)}_{0}$ is the 0-th order Hankel function of the first kind.

Let us focus on the large $\ell$ behavior in the subluminal case. Using the
asymptotic behavior of the Hankel function for large arguments, we obtain

$\displaystyle\lim_{\ell\to\infty}|c_{\ell}|^{2}=2\pi
v_{p}^{2}\frac{1}{\sin\theta}(\ell+1)^{2}\Bigg{|}\frac{(-\ell+(2+\ell)v_{p}^{2})}{\sqrt{\ell+\frac{1}{2}}}+\frac{2v_{p}e^{i\theta}}{\sqrt{\ell+\frac{3}{2}}}\Bigg{|}^{2}\,.$
(25)

Then, in the limit $kL\to\infty$, the coefficients $a_{\ell}$ approaches a
constant at large $\ell$,

$\displaystyle\lim_{\ell\to\infty}a_{\ell}=\lim_{\ell\to\infty}(2\ell+1)\frac{2(\ell-2)!}{(\ell+2)!}\left|c_{\ell}(k)\right|^{2}=8\pi
v_{p}^{2}(1-v_{p}^{2})^{3/2}+\mathcal{O}(1/\ell)\ ,$ (26)

as found in Liang  _et al._ (2023). One may wonder if this indicates the
divergence at $\xi\to 0$. It turns out that this is a fake divergence due to
the assumption of infinite distance. When one takes account of the finite
distance effect, one can see from Fig. 1 that the numerical calculated
coefficients starts to deviate from the analytical prediction at $\ell\sim
kL$. It is of interest (and perhaps necessary) to study the finite distance
effect in the PTA system. We explicitly calculate the perturbative expansion
of Eq. (14) in $1/(kL)$ in Appendix (B), and show that the corrections to the
coefficients $a_{\ell}$ are exactly zero at $\mathcal{O}(1/kL)$ and
$\mathcal{O}(1/(kL)^{2})$. This indicates that the finite distance effect
might be a non-perturbative effect, which gives an exponential decay like
$e^{-\ell/kL}$. This exponential decay behavior means a truncation in the
summation of the Legendre P polynomials, and eventually leads to a finite
result at $\xi\to 0$.

##  IV Constraints from current CPTA data

In this section, as an application of our analytical formula, we consider
observational constraints on the phase velocity. We use the data set given by
CPTA collaboration since they provide the frequency dependent angular
correlation function in their 3.4 years result. Although they provide data
sets at three different frequencies, we use the $1.5/3.4yr$ data since it has
the best statistical significance. 111The use of the lowest frequency data is
mostly affected by systematic errors due to modelling the pulsar timing. On
the other hand, the highest frequency also has a low significance due to a
suppressed signal in the power law modelling. We expect a future data release
may give us freqeuncy-dependent ORFs at various frequencies.

We use the significance function $\mathcal{S}$ discussed in CPTA collaboration
Xu  _et al._ (2023); Jenet  _et al._ (2005) to study the statistical
significance of the modified Hellings-Downs curve in modified gravity against
a constant correlation,

$\displaystyle\mathcal{S}\equiv\sqrt{\frac{N(N-1)}{2}}\frac{\sum_{i<j}\left(c_{ij}-\overline{c}_{ij}\right)\left(\Gamma_{ij}-\overline{\Gamma}_{ij}\right)}{\sqrt{\sum_{i<j}\left(c_{ij}-\overline{c}_{ij}\right)^{2}\sum_{i<j}\left(\Gamma_{ij}-\overline{\Gamma}_{ij}\right)^{2}}}\
,$ (27)

where $i,j$ sum from $1$ to $N$ through all pulsars in the system.
$\Gamma_{ij}$ is the theoretical ORF in a modified gravity theory with a
pulsar pair $\xi_{ij}=\hat{p}_{i}\cdot\hat{p}_{j}$ obtained throught Eq. (13),
while $c_{ij}$ is the observed angular correlation of this pulsar pair.
$\bar{c}_{ij}$ and $\bar{\Gamma}_{ij}$ represent the average
$\bar{c}_{ij}=\sum_{i<j}c_{ij}/N_{p}$, $N_{p}=N(N-1)/2$ is the number of
pulsar pairs in the system.

We show in Fig. 2 the significance function as a function of the phase
velocity. It is interesting that we can provide quite a good constraint for
the subluminal phase velocity case. When $v_{p}<0.88$, the significance drops
below $4\sigma$, and when the phase velocity approaches $v_{p}\lesssim 0.4$,
the significance drops below $3\sigma$. While in the $v_{p}>1$ case, shown in
the green curve of Fig. 2, the significance does not decrease rapidly, leaving
room for the possibility of massive gravity.

It is also interesting to notice that there are two peaks in the significance
function Eq.(27), in both $v_{p}<1$ and $v_{p}>1$ regions: $v_{p}=0.98$ and
$v_{p}=1.04$, although they do not have any statistical significance. As an
exercise, however, let us see if these peaks can be explained by modified
gravity. Let us assume the modified dispersion relation,

$\displaystyle\omega^{2}=c_{s}^{2}k^{2}+m^{2}\
,v_{p}=\frac{\omega}{k}=\sqrt{c_{s}^{2}+\frac{m^{2}}{k^{2}}}\,.$ (28)

Apparently a massless gravity theory with $c_{s}=0.98$ is a candidate for the
peak $v_{p}=0.98$, though constraints from other observational data would
probably exclude such a case. If we focus on the peak at $v_{p}=1.04$, this
corresponds to a graviton mass $m=1.02\times 10^{-22}$eV if $c_{s}=1$. It
seems this massive gravity theory is still viable Finn and Sutton (2002);
Abbott  _et al._ (2016); Liang and Trodden (2021); Abbott  _et al._ (2021);
Bernardo and Ng (2023c); Wang and Zhao (2023); Wu  _et al._ (2024); Bernardo
and Ng (2023b); Wu  _et al._ (2024, 2023).

Figure 2: We show the statistical significance of CPTA data as a function of
phase velocity. One can see a significant drop down for $v_{p}<1$.

##  V Conclusions and Discussions

In this paper we gave an analytical expression for the coefficients of the
multipole moments of the overlap reduction function for the angular
correlation function in the PTA system. Then as an application, adopting the
recent data from CPTA, which provides a frequency dependent data set, we
studied constraints on the phase velocity. Although the current data has large
uncertainty and could not provide good constraints to the phase velocity, we
find it disfavors the subluminal phase velocity case since the significance
function drops down drastically.

We conclude that in the future data release, the joint constraints from
different frequency bands should be available, and the PTA collaborations
should provide the frequency-dependent ORF as part of the results.

###### Acknowledgements.

We thank Kejia Lee for providing CPTA data for us. We thank Siyuan Chen, Meng-
Xiang Lin, Mark Trodden, and Wayne Hu for useful discussions. This work is
supported by the World Premier International Research Center Initiative (WPI),
MEXT, Japan, and also supported in part by JSPS KAKENHI Nos. JP19K14702 (IO),
JP20H05859 (IO), JP20H05853 (MS), and JP24K00624 (MS).

##  Appendix A Spherical harmonics decomposition

In this section, we follow Gair  _et al._ (2014); Liang  _et al._ (2023) and
list the basis transformation between the polarization basis and the spherical
harmonic function basis.

The harmonic basis $Y^{E,B}_{(lm)ij}(\hat{\Omega})$ and the polarization basis
$e^{+,\times}_{ij}(\hat{\Omega})$ are related through the spin-weighted
spherical harmonics $~{}_{\pm 2}Y_{lm}(\hat{\Omega})$,

$\displaystyle
Y^{E}_{(lm)ij}(\hat{\Omega})=\frac{1}{2}\left(~{}_{2}Y_{lm}(e_{ij}^{+}(\hat{\Omega})-ie_{ij}^{\times}(\hat{\Omega}))+~{}_{-2}Y_{lm}(e_{ij}^{+}(\hat{\Omega})+ie_{ij}^{\times}(\hat{\Omega}))\right)\
,$ $\displaystyle
Y^{B}_{(lm)ij}(\hat{\Omega})=\frac{1}{2}\left(~{}_{2}Y_{lm}(e_{ij}^{+}(\hat{\Omega})+ie_{ij}^{\times}(\hat{\Omega}))+~{}_{-2}Y_{lm}(e_{ij}^{+}(\hat{\Omega})-ie_{ij}^{\times}(\hat{\Omega}))\right)\
.$ (29)

The coefficients of gravitational wave $h_{ij}$ in these two basis are related
as

$\displaystyle
h^{+}(k,\hat{\Omega})=\frac{1}{2}\sum_{lm}\left(~{}_{2}Y_{lm}(a_{lm}^{E}(\hat{\Omega})+ia_{lm}^{B}(\hat{\Omega}))+~{}_{-2}Y_{lm}(a_{lm}^{E}(\hat{\Omega})-ia_{lm}^{B}(\hat{\Omega}))\right)\
,$ $\displaystyle
h^{\times}(k,\hat{\Omega})=\frac{1}{2}\sum_{lm}\left(~{}_{2}Y_{lm}(a_{lm}^{E}(\hat{\Omega})-ia_{lm}^{B}(\hat{\Omega}))+~{}_{-2}Y_{lm}(a_{lm}^{E}(\hat{\Omega})+ia_{lm}^{B}(\hat{\Omega}))\right)\
.$ (30)

##  Appendix B Legendre Q functions

We first consider the change of the integration variable from $x$ to $y$,

$y\equiv kL(x+v),\ \
c_{\ell}=v\int_{kL(v-1)}^{kL(v+1)}dy\frac{\left(1-e^{-iy}\right)}{y}\left[(1-x^{2})^{2}\frac{d^{2}}{dx^{2}}P_{\ell}(x)\right]_{x=-v+\frac{y}{kL}}\
,$ (31)

where we note that

$(1-x^{2})^{2}\frac{d^{2}}{dx^{2}}P_{\ell}(x)=(\ell+1)\left((-\ell+(2+\ell)x^{2})P_{\ell}(x)-2xP_{\ell+1}(x)\right)\,.$
(32)

The exponential factor in Eq. (31) can be split into two parts,

$\frac{1-e^{-iy}}{y}=\frac{1-\cos y}{y}+i\frac{\sin y}{y}\ .$ (33)

The first term is odd in $y$ while the second term is even. As $kL\to\infty$,
the high oscillation behavior for the even term $\sin y/y$ will pick up the
phase at $y=0$, provided that $v<1$. Therefore, the imaginary part of
$c_{\ell}$ is

$\Im
c_{\ell}=v\pi(\ell+1)\left((-\ell+(2+\ell)v^{2})P_{\ell}(-v)+2vP_{\ell+1}(-v)\right)\
.$ (34)

Numerically this matches the result well.

One can also verify that the imaginary part vanishes when $v\geq 1$. This is
because, when $v>1$, the phase in the exponent oscillates infinitely many
times in the limit $kL\to\infty$, due to the fact that there exists no
gravitational wave whose phase can remain constant along the null geodesic
from a pulsar to the earth. In the case of $v=1$, i.e., the case of general
relativity, the only gravitational wave whose phase can remain constant is the
one that propagates along the same null geodesic from the pulsar to the earth,
which is an isolated point in the solid angle
$\cos^{-1}(\hat{\Omega}\cdot\hat{p})$ and hence does not contribute to the
integral. On the other hand, if $v<1$, there exists a circle in the solid
angle in the momentum space toward which the gravitational wave phase remains
constant. This makes the integral along the circle finite. In other words, the
imaginary term is due to the coherent sum of the gravitational waves whose
propagation directions are toward this circle on which the phase remains
constant.

We now focus on the real part of $c_{\ell}$. One may naively guess that since
$(1-\cos y)/y$ is odd, as $kL\to\infty$, this would vanish under the
integration. However, expanding the Legendre polynomial in series of
associated Legendre polynomial,

$P_{\ell}\left(-v+\frac{y}{kL}\right)=\sum_{m=0}^{\ell}\frac{(-1)^{m}P^{m}_{\ell}(-v)}{(1-v^{2})^{m/2}m!}\left(\frac{y}{kL}\right)^{m}\
,$ (35)

we see that it contains both even and odd terms. For even $m$, since $(1-\cos
y)/y$ is odd, they vanish after integration. However, the odd $m$ terms are
non-vanishing and they give the real part of the coefficients. In principle,
one can integrate each term in the series and resum to obtain the result
analytically. Unfortunately, the expression is too complicated to resum for
general $\ell$ since each term contributes at the same order.

Nevertheless, in the limit, $kL\to\infty$, one can resum the leading order
terms to obtain

$\displaystyle\int_{kL(v-1)}^{kL(v+1)}dy\frac{1-\cos
y}{y}P_{\ell}\left(-v+\frac{y}{kL}\right)$ (36) $\displaystyle=$
$\displaystyle\int_{kL(v-1)}^{kL(v+1)}dy\frac{1-\cos
y}{y}P_{\ell}(-v)+\int_{kL(v-1)}^{kL(v+1)}\frac{dy}{kL}\sum_{m=1}^{\ell}P^{m}_{\ell}(-v)(-1)^{m}\frac{1-\cos
y}{m!(1-v^{2})^{m/2}}\left(\frac{y}{kL}\right)^{m-1}$ $\displaystyle=$
$\displaystyle\log\left(\frac{1+v}{1-v}\right)P_{\ell}(-v)+\sum_{m=1}^{l}P^{m}_{\ell}(-v)(-1)^{m}\frac{(1+v)^{m}-(v-1)^{m}}{m!m(1-v^{2})^{m/2}}+\mathcal{O}(1/kL)$
$\displaystyle=$ $\displaystyle-2\tilde{Q}_{\ell}(-v)+\mathcal{O}(1/kL)\ ,$

where we have used the definitions of the Legendre Q functions in Eq. (III).
We note that, in the limit $kL\to\infty$, there is an interesting similarity
of the above to Neumann’s integral representation of the Legendre Q functions
Bateman and Project (1953),

$\displaystyle
Q_{\ell}(z)=\frac{1}{2}\int_{-1}^{1}\frac{P_{\ell}(v)}{z-v}dv\,;~{}z\notin[-1,1]\,.$
(37)

The difference is that this $Q_{\ell}(z)$ is defined with a different choice
of branch cut (see (III) in the main context).

We can also evaluate the finite distance corrections in Eq. (36) through the
series expansion,

$\displaystyle\int_{kL(v-1)}^{kL(v+1)}dy\frac{1-e^{-iy}}{y}P_{\ell}\left(-v+\frac{y}{kL}\right)$
(38) $\displaystyle=$ $\displaystyle-2\tilde{Q}_{\ell}(-v)+\pi
iP_{\ell}(-v)+\int_{kL(v-1)}^{kL(v+1)}\frac{dy}{kL}\sum_{m=0}^{\ell}\frac{-e^{iy}P^{m}_{\ell}(-v)(-1)^{m}}{m!(1-v^{2})^{m/2}}\left(\frac{y}{kL}\right)^{m-1}$
$\displaystyle=$ $\displaystyle-2\tilde{Q}_{\ell}(-v)+\pi iP_{\ell}(-v)$
$\displaystyle\quad+\sum_{m=0}^{\ell}\frac{P^{m}_{\ell}(-v)(-1)^{m}}{m!(1-v^{2})^{m/2}}\Bigl{(}\gamma\bigl{(}m,ikL(v+1)\bigr{)}-\gamma\bigl{(}m,ikL(v-1)\bigr{)}\Bigr{)}\
,$

where $\gamma(m,z)$ is the incomplete gamma function. It has the asymptotic
expansion for large $z$,

$\displaystyle\lim_{z\to\infty}\gamma(m,z)=\Gamma(m)-e^{-z}z^{m-1}\left(1+\frac{m-1}{z}+\frac{(m-1)(m-2)}{z^{2}}+\mathcal{O}(1/z^{3})\right)\
.$ (39)

With the above expansion, we could obtain the corrections due to the finite
distance effect. We find

$\displaystyle\int_{kL(v-1)}^{kL(v+1)}dy\frac{\left(1-e^{-iy}\right)}{y}P_{\ell}\left(-v+\frac{y}{kL}\right)$
$\displaystyle=-2\tilde{Q}_{\ell}(-v)+\pi
iP_{\ell}(-v)+\frac{i}{kL}\left(\frac{e^{-ikL(v-1)}(-1)^{\ell}}{v-1}-\frac{e^{-ikL(v+1)}}{v+1}\right)+\mathcal{O}(1/k^{2}L^{2}).$
(40)

It is easy to see that the $\mathcal{O}(1/kL)$ correction in Eq.(32) vanishes.
While the next order $\mathcal{O}(1/(kL)^{2})$ is less trivial to see, one can
nevertheless check that their contribution to $c_{\ell}$ still vanishes. This
indicates that the finite distance contribution is perturbatively zero, and
therefore, it is a non-perturbative effect. Deriving this potentially non-
perturbative behavior for large $\ell$ is left for future studies.

## References

  * Hewish  _et al._ (1968) A. Hewish, S. J. Bell, J. D. H Pilkington, P. F. Scott, and R. A. Collins, “Observation of a rapidly pulsating radio source,” Nature 217, 709–713 (1968).
  * Detweiler (1979) Steven L. Detweiler, “Pulsar timing measurements and the search for gravitational waves,” Astrophys. J. 234, 1100–1104 (1979).
  * Hellings and Downs (1983) R. w. Hellings and G. s. Downs, “UPPER LIMITS ON THE ISOTROPIC GRAVITATIONAL RADIATION BACKGROUND FROM PULSAR TIMING ANALYSIS,” Astrophys. J. Lett. 265, L39–L42 (1983).
  * Agazie  _et al._ (2023) Gabriella Agazie _et al._ (NANOGrav), “The NANOGrav 15-year Data Set: Evidence for a Gravitational-Wave Background,” Astrophys. J. Lett. 951 (2023), 10.3847/2041-8213/acdac6, arXiv:2306.16213 [astro-ph.HE] .
  * Afzal  _et al._ (2023) Adeela Afzal _et al._ (NANOGrav), “The NANOGrav 15-year Data Set: Search for Signals from New Physics,” Astrophys. J. Lett. 951 (2023), 10.3847/2041-8213/acdc91, arXiv:2306.16219 [astro-ph.HE] .
  * Antoniadis  _et al._ (2023) J. Antoniadis _et al._ (EPTA, InPTA:), “The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals,” Astron. Astrophys. 678, A50 (2023), arXiv:2306.16214 [astro-ph.HE] .
  * Reardon  _et al._ (2023) Daniel J. Reardon _et al._ , “Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array,” Astrophys. J. Lett. 951, L6 (2023), arXiv:2306.16215 [astro-ph.HE] .
  * Xu  _et al._ (2023) Heng Xu _et al._ , “Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I,” Res. Astron. Astrophys. 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .
  * Lee  _et al._ (2010) Kejia Lee, Fredrick A. Jenet, Richard H. Price, Norbert Wex, and Michael Kramer, “Detecting massive gravitons using pulsar timing arrays,” Astrophys. J. 722, 1589–1597 (2010), arXiv:1008.2561 [astro-ph.HE] .
  * Gair  _et al._ (2015) Jonathan R. Gair, Joseph D. Romano, and Stephen R. Taylor, “Mapping gravitational-wave backgrounds of arbitrary polarisation using pulsar timing arrays,” Phys. Rev. D 92, 102003 (2015), arXiv:1506.08668 [gr-qc] .
  * Qin  _et al._ (2021) Wenzer Qin, Kimberly K. Boddy, and Marc Kamionkowski, “Subluminal stochastic gravitational waves in pulsar-timing arrays and astrometry,” Phys. Rev. D 103, 024045 (2021), arXiv:2007.11009 [gr-qc] .
  * Liang and Trodden (2021) Qiuyue Liang and Mark Trodden, “Detecting the stochastic gravitational wave background from massive gravity with pulsar timing arrays,” Phys. Rev. D 104, 084052 (2021), arXiv:2108.05344 [astro-ph.CO] .
  * Liang  _et al._ (2023) Qiuyue Liang, Meng-Xiang Lin, and Mark Trodden, “A test of gravity with Pulsar Timing Arrays,” JCAP 11, 042 (2023), arXiv:2304.02640 [astro-ph.CO] .
  * Bernardo and Ng (2023a) Reginald Christian Bernardo and Kin-Wang Ng, “Stochastic gravitational wave background phenomenology in a pulsar timing array,” Phys. Rev. D 107, 044007 (2023a), arXiv:2208.12538 [gr-qc] .
  * Bi  _et al._ (2024) Yan-Chen Bi, Yu-Mei Wu, Zu-Cheng Chen, and Qing-Guo Huang, “Constraints on the velocity of gravitational waves from the NANOGrav 15-year data set,” Phys. Rev. D 109, L061101 (2024), arXiv:2310.08366 [astro-ph.CO] .
  * Bernardo and Ng (2023b) Reginald Christian Bernardo and Kin-Wang Ng, “Beyond the Hellings-Downs curve: Non-Einsteinian gravitational waves in pulsar timing array correlations,” (2023b), arXiv:2310.07537 [gr-qc] .
  * Wu  _et al._ (2024) Yu-Mei Wu, Zu-Cheng Chen, Yan-Chen Bi, and Qing-Guo Huang, “Constraining the graviton mass with the NANOGrav 15 year data set,” Class. Quant. Grav. 41, 075002 (2024), arXiv:2310.07469 [astro-ph.CO] .
  * Wu  _et al._ (2023) Yu-Mei Wu, Zu-Cheng Chen, and Qing-Guo Huang, “Search for stochastic gravitational-wave background from massive gravity in the NANOGrav 12.5-year dataset,” Phys. Rev. D 107, 042003 (2023), arXiv:2302.00229 [gr-qc] .
  * Finn and Sutton (2002) Lee Samuel Finn and Patrick J. Sutton, “Bounding the mass of the graviton using binary pulsar observations,” Phys. Rev. D 65, 044022 (2002), arXiv:gr-qc/0109049 .
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “Tests of general relativity with GW150914,” Phys. Rev. Lett. 116, 221101 (2016), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Abbott  _et al._ (2021) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), “Tests of General Relativity with GWTC-3,” (2021), arXiv:2112.06861 [gr-qc] .
  * Bernardo and Ng (2023c) Reginald Christian Bernardo and Kin-Wang Ng, “Constraining gravitational wave propagation using pulsar timing array correlations,” Phys. Rev. D 107, L101502 (2023c), arXiv:2302.11796 [gr-qc] .
  * Wang and Zhao (2023) Sai Wang and Zhi-Chao Zhao, “Unveiling the Graviton Mass Bounds through Analysis of 2023 Pulsar Timing Array Data Releases,” (2023), arXiv:2307.04680 [astro-ph.HE] .
  * Carrillo Gonzalez  _et al._ (2022) Mariana Carrillo Gonzalez, Claudia de Rham, Victor Pozsgay, and Andrew J. Tolley, “Causal effective field theories,” Phys. Rev. D 106, 105018 (2022), arXiv:2207.03491 [hep-th] .
  * Ezquiaga  _et al._ (2021) Jose Maria Ezquiaga, Wayne Hu, Macarena Lagos, and Meng-Xiang Lin, “Gravitational wave propagation beyond general relativity: waveform distortions and echoes,” JCAP 11, 048 (2021), arXiv:2108.10872 [astro-ph.CO] .
  * de Rham and Tolley (2020) Claudia de Rham and Andrew J. Tolley, “Speed of gravity,” Phys. Rev. D 101, 063518 (2020), arXiv:1909.00881 [hep-th] .
  * Gair  _et al._ (2014) Jonathan Gair, Joseph D. Romano, Stephen Taylor, and Chiara M. F. Mingarelli, “Mapping gravitational-wave backgrounds using methods from CMB analysis: Application to pulsar timing arrays,” Phys. Rev. D 90, 082001 (2014), arXiv:1406.4664 [gr-qc] .
  * Jenet  _et al._ (2005) Fredrick A. Jenet, George B. Hobbs, K. J. Lee, and Richard N. Manchester, “Detecting the stochastic gravitational wave background using pulsar timing,” Astrophys. J. Lett. 625, L123–L126 (2005), arXiv:astro-ph/0504458 .
  * Bateman and Project (1953) Harry Bateman and Bateman Manuscript Project,  _Higher Transcendental Functions [Volumes I-III]_ (McGraw-Hill Book Company, 1953) funding by Office of Naval Research (ONR).


# Characterizing and Mitigating Telluric Absorption in Precise Radial
Velocities II: A Study of an M2 Type Star

Natasha Latouf Department of Physics and Astronomy, George Mason University,
4400 University Drive MS 3F3, Fairfax, VA, 22030, USA  Sharon Xuesong Wang
Department of Astronomy, Tsinghua University, 30 Shuangqing Rd, Haidian
District, Beijing, China  Bryson Cale Department of Physics and Astronomy,
George Mason University, 4400 University Drive MS 3F3, Fairfax, VA, 22030, USA
Peter Plavchan Department of Physics and Astronomy, George Mason University,
4400 University Drive MS 3F3, Fairfax, VA, 22030, USA Natasha Latouf, Sharon
Xuesong Wang<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Telluric absorption lines impact measuring precise radial velocities (RVs)
from ground-based, high-resolution spectrographs. In this paper, we simulate
the dependence of this impact on stellar spectral type and extend the work of
the first paper in this series, which studied a G type star, to a synthetic M
dwarf star. We quantify the bias in precise RV measurements in the visible and
near-infrared (NIR) from the presence of tellurics in a simulated set of
observations. We find that M dwarf RVs are more impacted by tellurics compared
to G type stars. Specifically, for an M dwarf star, tellurics can induce RV
errors of up to 16 cm/s in the red-optical and in excess of 220 cm/s in the
NIR. For a G dwarf, comparable RV systematics are 3 cm/s in the red optical
and 240 cm/s in the NIR. We attribute this relative increase for M dwarfs
stars to the increased concordance in wavelength between telluric lines and
stellar Doppler information content. We compare the results of our simulation
to data collected on Barnard’s star from the iSHELL spectrograph at the NASA
Infrared Telescope Facility (IRTF). This study was conducted as a follow-up to
the NASA probe mission concept study EarthFinder.

methods: numerical; techniques: radial velocities, high angular resolution

##  1 Introduction

Doppler spectroscopy is currently the most prolific method for measuring the
masses and characterizing the orbits of exoplanets. Many factors contribute to
the measured radial velocity (RV) precision, including limitations from photon
noise, instrument precision, telluric absorption, and stellar activity
(Plavchan et al., 2015; Fischer et al., 2016). New precise radial velocity
(PRV) spectrometers such as ESPRESSO, NEID, MAROON-X and EXPRES now have the
requisite instrument stability to be sensitive to habitable zone (HZ) Earth-
mass planets (Jurgenson et al., 2016; Kanodia et al., 2020; Pepe et al., 2021;
Seifahrt et al., 2021). Additionally, these spectrometers have broadened their
spectral grasps into the red-optical out to $\sim$900 nm, and several near-
infrared (NIR) PRV spectrometers have recently demonstrated $<$3 m/s long-term
RV precision at wavelengths out to 2.4 $\mu$m, such as IRD, HPF, CARMENES,
SPIRou, and iSHELL (Tamura et al., 2012; Mahadevan et al., 2014; Quirrenbach &
Carmenes Consortium, 2016; Donati et al., 2018; Cale et al., 2019).

Figure 1: Top: Fluxes with arbitrary offsets of the Kurucz solar model, BT-
Settl M dwarf model (Allard et al. (2012a, b)), and normalized telluric
spectra gathered from TAPAS (Bertaux et al. (2014)), shown in blue, orange,
and green, respectively, plotted as a function of wavelength in nanometers.
This plot ranges from 350 - 850 nm. There is more spectral information content
in the M dwarf red-optical spectrum, where there are more telluric lines. By
contrast, the solar spectrum has more information content in the optical, with
less telluric lines. Middle: The same as the top panel, but from 850 - 1650
nm. Bottom: The same as the other two panels, but from 1650 - 2500 nm.

Stellar activity challenges our current sensitivity to HZ Earth-mass planets,
the latter of which produce $\sim$8-25 cm/s Doppler reflex semi-amplitudes for
G type to M dwarf stars. By comparison, spots, plages/faculae, the suppression
of convective blueshift, and other photospheric effects induce apparent RV
changes on the order of a few m/s (Dumusque et al., 2011; Haywood et al.,
2020; Cegla et al., 2019; Dumusque, 2018). Stellar activity is inherently
chromatic, and can be larger in amplitude in the visible wavelengths than in
the red-optical and NIR (Marchwinski et al., 2015; Robertson et al., 2020;
Plavchan et al., 2020; Tran et al., 2021; Cale et al., 2021). Investigating
the impact of stellar activity on PRVs as a function of wavelength may help to
isolate stellar activity from Keplerian (achromatic) signals. Increased
spectral grasp of PRV spectrometers have also enabled more RV searches for
exoplanets orbiting M dwarf stars, especially the cooler, fainter, late types
(Morales et al., 2019). Unlike G dwarfs, M dwarfs have the highest RV
information content in the red-optical and are brightest in the NIR (Reiners
et al., 2018, 2010). M dwarfs stars are also the most abundant main sequence
stars and exoplanet hosts (Henry et al., 2018; Benedict et al., 2016; Dressing
& Charbonneau, 2013; Raghavan et al., 2010).

Consequently, we are motivated by the chromaticity of stellar activity and RV
information content for G and M dwarfs alike to investigate the challenges to
RV precision across a broad spectral range. At almost all wavelengths,
tellurics are present in the spectra recorded with ground-based PRV
spectrometers. At visible wavelengths, which we define in this work to be
400-700 nm, telluric lines appear to be limited within certain oxygen bands,
but in reality, there are copious micro-telluric lines. These shallow lines,
typically ignored by telluric masks when extracting radial velocities (e.g.,
see Cunha et al. 2014), are dominated by water, a molecule with a strong
electric dipole moment, producing absorption features that are small in depth
($\sim$1%) but highly variable with time and airmass. At red-optical
wavelengths, which we define to be $\sim$700-1000 nm, both micro-telluric and
deeper macro-telluric lines are present from water and molecular oxygen (O2).
At near-infrared wavelengths, which we define to be $\sim$1-2.5 $\mu$m, more
prevalent and deeper telluric absorption lines are formed from water in
addition to methane, carbon dioxide, nitrous dioxide, and oxygen, producing
lines that can be saturated and remain variable with time and airmass. These
lines are challenging to model to an accuracy of 1-2% in flux, (e.g., Smette
et al., 2015; Gullikson et al., 2014; Seifahrt et al., 2010; Ulmer-Moll et
al., 2019). However, when analyzing real RV observations, micro-tellurics are
typically comparable to or below the noise level and are either masked or
ignored altogether.

In this work we specifically consider the impact tellurics have on PRV
measurements of M dwarfs. This project is a direct continuation from Wang et
al. (2022 in prep, hereafter W22). In W22, we quantified the error due to
tellurics for a solar twin star (G2V), using the Kurucz solar model (Kurucz,
2005). Using simulated spectra, we extracted the RVs and assessed the RV
precision as set by telluric contamination alone. By generating noiseless
simulated spectra, we isolated the telluric impact on PRV measurements from
other sources such as photon noise and stellar jitter. We extracted the RVs
with the cross-correlation function (CCF) and the Forward Modeling method in a
broad wavelength range (350 - 2500 nm).

In this project, we follow the same methodology as in W22, with the
substitution of a synthetic M dwarf model instead of the solar model. There is
significant wavelength overlap between the stellar RV information content and
the telluric absorption lines for M dwarfs, and thus we are motivated to
quantify the RV error budget term due to the telluric absorption. In $[\S]$ 2
we present the methodology of our experiment, also providing a brief summary
of W22. In $[\S]$ 3, we present the results of our experiment with the M dwarf
stellar model. In ${\S}$ 4, we analyze the differences between the G type star
results presented in W22 with the results of the M dwarf contained herein. We
present a comparison to our simulation of an M dwarf with an analysis of
iSHELL data of Barnard’s star in subsection 4.4. We present our conclusions
and future work in $[\S]$ 5.

##  2 Methodology

In this section, we describe the simulations used in our analyses. We follow a
similar approach to W22. Herein we present a summary of the main steps in our
analysis. Our simulation’s code is summarized in a flow chart in W22. The code
uses several inputs, including:

  * •

A model for the stellar M dwarf spectrum. We use a BT-Settl model,
rotationally broadened to a $v$sin$i=2$ km/s, which is representative of a
relatively quiescent M dwarf (Allard et al., 2012a, b).

    * –

We used the BT-Settl atmospheric models from Allard et al. (2012a, b); Baraffe
et al. (2015) to generate the synthetic spectrum, with version 15.5 of the
PHOENIX stellar atmosphere modeling code. The spectrum has been continuum
normalized by fitting a blackbody curve initially and then refining it by
using a power-law to continuum normalize certain sections that were not well
normalized in the first pass.

    * –

The stellar parameters used to generate the synthetic spectrum are as follows:
Effective temperature of 3563K, Metallicity of $[M/H]=-0.38$, $\log{g}=4.8$ to
best match the stellar spectrum of the RV standard star HD 95735 as used in
Wang et al. (2019) (see section 2.1 therein for more detail). This will enable
us to compare results in the V band in the optical directly with Wang et al.
(2019), which studied the effects of telluric lines in the iodine-calibrated
PRVs in the V band, at roughly 500-600 nm.

  * •

A suite of models for the Earth’s atmospheric absorption provided by TAPAS,
containing data that were computed based on the atmospheric conditions of a
random night per month over a year from a location on Kitt Peak, i.e. we
generated 12 spectra suites, generally from the 15th day of each month for
consistency, unless data was unavailable (Bertaux et al., 2014). The TAPAS
suite models used for the purposes of this project was identical to the suites
used in W22.

  * •

A single Gaussian model for the line spread function (LSF) equivalent to a
spectral resolving power of $R=120,000$, which is the approximate spectral
resolution commonly adopted by the next generation, ground-based RV
instruments, e.g., NEID, EXPRES, ESPRESSO. This is kept unchanged throughout
our simulations.

  * •

The airmass and optical depth ($\tau$) for 6 molecular species (water,
methane, carbon dioxide, nitrous oxide, oxygen, and ozone molecules), as well
as a random value drawn from a distribution of real site values for the amount
of precipitable water vapor (PWV). The barycentric velocity correction (BC)
was calculated as in W22:

$BC=sin(date/(2\pi\cdot 365.15)-\pi/2)\cdot 30km/s$ (1)

The telluric spectrum was shifted according to the calculated BC value to
mimic Earth’s motion. We shifted the telluric spectrum rather than the stellar
spectrum due to the higher fidelity of TAPAS relative to the fidelity of the
stellar spectrum. The barycentric velocity, airmass, and PWV distributions are
the same between W22 and this experiment to allow for a direct comparison.

A direct comparison of the Kurucz solar model, the BT-Settl M dwarf model, and
the telluric model from TAPAS is shown in Figure 1. As seen and described
further, the solar model has more RV information content in the optical
wavelengths, less plagued by tellurics. For the M dwarf, however, the peak RV
information content occurs in the red-optical, where there are deeper telluric
absorption lines.

We make several assumptions in the setup of our analysis. In particular, we
assume we are observing from a telescope on Kitt Peak. We assume that we
observe 365 nights per year, once for each night of the year, e.g. a
circumpolar best case scenario. As in W22, we compute our observed spectra in
the rest frame of the star. Thus we shift our model telluric spectrum by the
negative of the barycenter velocity. We approximate the barycentric velocity
correction as a sine wave with a one year period, phase, and semi-amplitude of
10.7 km/s that approximately matches the archival observations of Tau Ceti
with Keck/HIRES (Howard & Fulton, 2016; Tuomi et al., 2013), to match with the
assumption in W22 since, as can be seen later, the individual RV measurement
errors correlate strongly with the barycentric RVs.

For these simulations, we used the BT-Settl M dwarf spectral model. We assume
no stellar activity modifies the stellar spectra and line shapes. We then
multiplied the M dwarf spectra with the TAPAS atmosphere, which varies by
night for each of the 365 nights in our simulations according to the airmass,
the PWV, and the barycentric correction (BC) of each night. Lastly, we
convolved the spectra at a resolution R = 120,000, as expressed above. We did
not add photon noise or order blaze functions, in order to further isolate the
errors due to tellurics.111For more information, please see the discussion on
caveats in W22. Then, we extracted the RVs using two methods: CCF and Forward
Modeling. Each of these two methods was run with two sets of assumptions: with
perfectly known telluric line profiles and imperfectly known telluric line
profiles.

The first method used for extracting RVs was the CCF (Baranne et al., 1996):

  1. 1.

No Correction: In this method, RVs were extracted using the method similar to
Zechmeister et al. (2018) using a stellar template instead of a line mask. We
also assumed, unrealistically, that we have a perfectly known wavelength
solution, LSF, and stellar template. We do not apply any telluric line-masking
before extracting RVs, as is done with SERVAL and other RV analyses with the
CCF method at visible wavelengths.

  2. 2.

Correction by Division: In this method, we assume a perfectly known telluric
line that is the same as the telluric absorption that we put in for each
night. We correct for tellurics by dividing the telluric lines (broadened to a
resolution of 120,000) before extracting the RVs. This is the most commonly
used method of correction for instruments such as HARPS or CARMENES, which are
stabilized spectrographs (e.g., Cunha et al. 2014; Artigau et al. 2014). While
this method works fairly well, it is also technically mathematically
incorrect, involving division after, rather than before, convolution by the
LSF. For more information on the incorrect mathematical treatment of dividing
an observed spectrum by a telluric model, consult W22.

  3. 3.

Correction by Division with Wrong Telluric Line Profile: In this method, we
divide out a set of telluric lines generated with a line profile that is
mismatched to the telluric line profile we used when multiplying to form our
total simulated spectra. Here, the wrong telluric line profile refers to the
intrinsic line profiles of the telluric lines and not the spectrograph line
spread function. The mismatched telluric profile is based on the atmospheric
model of a different observing site than our assumed site of Kitt Peak. Same
as in W22, we used Maunakea as our mismatched profile location, hereafter
known as the K-profile. We selected Maunakea due to this difference in
elevation from our assumed observation site of Kitt Peak, which would mimic
the lack of precise knowledge of the telluric line profile and molecular
abundances in ground-based observations.

The second method for extracting RVs we used was Forward Modeling:

  1. 1.

Perfectly Known Line Profile: In this method, we expect to get a result of 0
m/s error with our RV extraction. This method is used as a control to reaffirm
our analysis yields the injected RVs. These results are as expected and not
presented herein.

  2. 2.

Imperfectly Known Line Profile: In this method, the Forward Modeling algorithm
is performed with imperfect knowledge of the existing telluric line profile,
the same as in the method of using CCF with correction by division with
K-profile as described above.

For the RVs of each order, we calculate the overall scatter over time (RMS),
which is used as the quantification of the errors induced by the tellurics
from this order. The spectrum is divided into 230 orders and this is done for
all orders. The spectrum division is conducted to mimic real observations,
where the RVs tend to be individually computed for each order independently
and combined afterwards. The orders are weighted to account for the Doppler
content, and the order sizes are set to linearly increase with increasing
wavelength to remain constant in $\delta\lambda/\lambda$. In addition, we take
a slice through time, getting RVs from all orders for any given night, and
statistically combine them according to the stellar Doppler content and its
overall scatter due to tellurics. In this way, the final nightly multi-order
RV is calculated.

Ground-based | No Correction | Division | Division | Modeling | Photon-limited
---|---|---|---|---|---
Instrument | (m/s) | (m/s) | K-Profile (m/s) | K-Profile (m/s) | Precision (m/s)
EFE Visible Arm |  |  |  |  |
(380-900 nm) | 0.167 | 0.125 | 0.120 | 0.140 | 0.102
EFE NIR Arm |  |  |  |  |
(900-2500 nm) | 2.158 | 0.636 | 1.622 | 1.636 | 0.160
ESPRESSO |  |  |  |  |
(380-788 nm) | 0.166 | 0.137 | 0.119 | 0.137 | 0.109
EXPRES |  |  |  |  |
(380-680 nm) | 0.149 | 0.042 | 0.119 | 0.104 | 0.138
NEID |  |  |  |  |
(380-930 nm) | 0.167 | 0.124 | 0.112 | 0.127 | 0.104
CARMENES |  |  |  |  |
Visible Arm |  |  |  |  |
(520-960 nm) | 0.302 | 0.170 | 0.171 | 0.198 | 0.159
CARMENES |  |  |  |  |
NIR Arm |  |  |  |  |
(960-1710 nm) | 2.122 | 0.685 | 1.149 | 1.156 | 0.293
Table 1: Summary of RV RMS of an M dwarf for different instruments including
an EarthFinder ground equivalent instrument using the different RV extraction
methods previously discussed in Section 2. Columns 2-5 remain the same as
Columns 3-6 in Table 3. Column 6: theoretically computed SNR = 100 photon-
limited precision (m/s) for our chosen M dwarf in each band.

Figure 2: Example spectral orders shown on the left with their corresponding
RVs vs. the barycentric correction on the right when specifically ignoring
tellurics during RV extraction, i.e. CCF + no correction. The left plots show
the M dwarf spectrum in blue and the telluric spectrum in orange for each
order, with the telluric absorption on order #28, shown in the top left panel
amplified by 10 in flux for better visibility.

##  3 Results

Herein we present our results for the simulation of the BT-Settl M dwarf
model, as compared to the prior results of the G star simulation. All plots
are replicated in similar forms to W22 for consistency, and comparative
material has been added where necessary.

First, we present example spectral orders with corresponding RV vs barycentric
correction in Figure 2. Looking at the two left panels (top and bottom), we
can see the telluric spectrum and M dwarf spectrum between 445 - 448.5 nm and
1099 - 1107 nm, respectively. The corresponding right plots show how varying
depths of telluric lines affect RVs to different degrees and that the RV
biases caused by tellurics in the temporal space are highly structured due to
Earth’s barycentric motion. This is caused by the changes in wavelength
separation between the telluric lines and the stellar lines due to the Earth’s
barycentric motion. This rastering is also discussed in W22 and other prior
work, e.g. Figure 3 in Cunha et al. (2014) and references therein. This shows
how tellurics add to the RV error budget, by adding scatter to the measured
RVs and adding coherent systematic signal. This is described further in W22,
presenting the peak pulling effect and the dependence of amplitude on the
relative Doppler content of the stellar and telluric lines.

The full suite of extracted RVs as a function of airmass, PWV, and BC is
presented in Figure 3. Each data point presented is a weighted average of
reported RVs from 350 - 2500 nm, with assumed SNR of 400 per pixel and
combined photon noise with telluric-induced RV RMS. The lack of correlation
between RVs and PWVs is expected, as the input PWV value would not match the
K-profile PWV value and the water lines are largely mitigated, although not
perfectly. This also shows that the scatter in RVs is weakly related to the
airmass, due to the final RVs being dominated by RVs in the optical region.
However, we can see that the BC is the most strongly correlated to the RV RMS,
due to the highly structured systematics.

Figure 3: Extracted RVs from 365 nights as functions of airmass, PWV, and BC for three methods: No correction (left), CCF$+$division with the K-profile for tellurics (middle), and Forward Modeling with the K-profile for tellurics (right). Each data point is a weighted average of reported RVs from all orders from 350 - 2500 nm, the entire wavelength range in our simulations. The weights assume SNR = 400 per pixel and combine photon noise with telluric-induced RV RMS, as described in W22, Section 2, and as used in the presented data in Table 1.  Photometric | No Correction | Division | Division | Modeling | Photon-limited
---|---|---|---|---|---
Band | (m/s) | (m/s) | K-Profile (m/s) | K-Profile (m/s) | Precision (m/s)
B (398–492 nm) | 0.066 | 0.000 | 0.050 | 0.058 | 0.451
V (507–595 nm) | 0.233 | 0.021 | 0.154 | 0.169 | 0.339
R (589–727 nm) | 0.556 | 0.259 | 0.523 | 0.428 | 0.255
I (732–881 nm) | 1.523 | 0.525 | 1.513 | 1.551 | 0.632
Z (0.8–1.0 $\mu$m) | 2.624 | 0.886 | 2.223 | 2.306 | 0.625
Y (0.96–1.1 $\mu$m) | 2.131 | 0.419 | 2.032 | 1.910 | 0.706
J (1.1-1.3 $\mu$m) | 17.769 | 2.859 | 11.275 | 9.872 | 1.341
H (1.5–1.8 $\mu$m) | 6.176 | 1.383 | 3.901 | 3.355 | 0.769
K (2.0–2.4 $\mu$m) | 33.987 | 1.757 | 22.424 | 6.941 | 1.056
Table 2: Summary of RV RMS of an M dwarf for several photometric bands.
Columns 2-6 remain the same as columns 2-6 in Table 1.

Figure 4: Telluric line induced RV errors as a function of wavelength for each
order in our simulation for the RV extraction methods listed in $[\S]$ 2. In
the top two plots, each point shows the RMS value of RVs extracted from one
spectral order each across one year (365 nights). The top left plot shows how
using an incorrect line profile introduces more error to the analysis. The top
right plot portrays how Forward Modeling generally has a lower RMS than
CCF$+$division, especially in the NIR. The bottom plot is the smoothed version
of the top plots for more clarity, with a window size of 11 orders and a
moving average.

Figure 5: Relationship between RV error and the stellar RV information content
for both the M dwarf and Kurucz solar model. The top left plot shows the M
dwarf RV error vs. the stellar RV information content and the bottom left plot
shows the solar RV error vs. the stellar RV information content. In the left
plot, the ‘no correction’ results are shown in gray, the K-profile division
results are shown in orange, and the K-profile Forward Model results are shown
in blue.

We study the comparison between an EarthFinder equivalent (EFE) instrument and
other currently existing telescopes in Table 1, broken down into visible and
NIR arms as necessary (Plavchan et al., 2020). We can see that the EFE
instrument has the widest wavelength range of any existing instrument, and the
lowest photon-limited precision in the visible wavelength bands, and the EFE
NIR arm outperforms or matches the currently available NIR instrument present
in CARMENES, while going out further into the NIR than CARMENES. The first two
rows presented in this table are reflected in Table 3 and the corresponding
Figure 4. These results show the effectiveness of EFE in the NIR for M dwarf
observations, which as discussed is also where the RVs are most affected by
tellurics. We present the relationship between RV error and the stellar RV
information content for both stellar models in Figure 5. We also present the
RV signal as a function of time for the M dwarf model only in Figure 6.
Further explanation for both can be found in Section 4.1.

We also present further analysis of RV RMS as a function of photometric band
in Table 2. As seen, the M dwarf RV RMS is significantly worse in the NIR than
for the G star, specifically in the JHK bands. Similarly to Table 3 presented
in this text, the results of Table 2 are also presented in Figure 7. This
figure contains all of the RV extraction methods discussed, and the horizontal
lines mark the photon-noise weighted precision in each band, representing the
last column in Table 2. Further explanation on Figure 7 can be found in
Section 4.2

##  4 Discussion

We present a discussion of comparison between M dwarf and G type,
effectiveness to mitigation methods, caveats to our simulation, and iSHELL
data comparisons.

###  4.1 M & G Dwarf Comparison

In Table 3 we present a comparative summary of the average nightly RV RMS
results for both the M dwarf and G type from W22 - with respect to wavelength
range retrieved with each method described in Section 2: No correction,
Division, Division with K-Profile, and Forward Modeling with the K-Profile.
Between an M dwarf and a G dwarf, RV RMS differences are present in all
wavelength bands present between 350 - 2500 nm. For an M dwarf (G star) in the
visible band, a simulated EFE instrument as shown in Table 3 had an RV error
RMS of 0.17 m/s (0.04 m/s) without correcting for tellurics, and 0.14 m/s
(0.03 m/s) using Forward Modeling with K-profile for M dwarfs. In the NIR for
the M dwarf (G star), there was an error of 2.2 m/s (2.4 m/s) without
correcting for tellurics, and 1.6 m/s (1.2 m/s) with the K-profile Forward
Modeling.

We graphically represent the results from Table 3 in Figure 4 for the M dwarf
specifically; the full results are contained within the top two panels, while
the bottom plot in Figure 4 depicts a smoothed version of the results for
clarity. As can be seen from this figure, the CCF division method with a
perfectly known telluric line profile yields the best results, but is an
unrealistic representation. However the forward modeling with K-profile, the
most realistic method of simulation, yields the next best results and is the
best realistic presented method between the other realistic scenarios of no
correction or CCF division with K-profile.

Based on the results of Table 3 and Figure 4, the M dwarf fares worse in the
optical (380-900 nm) wavelength range than the G type star. This is expected,
as the G type star has more stellar RV information content in the optical,
where there are naturally fewer and shallower telluric lines (Figure 1). For
the M dwarf, the peak in stellar RV information content appears in the red-
optical, at approximately 700 nm, where the telluric lines are deeper and more
frequent. This is explored in Figure 5. In the NIR (900-2500 nm), both models
return greater RV RMS as the telluric lines are deeply saturated in these
wavelengths.

For M dwarfs, the tellurics more significantly impact PRV observations at NIR
wavelengths. In the NIR, Forward Modeling with K-profile had slightly less RV
RMS than in the visible with no correction. This is expected since M dwarfs
are brightest in the NIR, with the highest RV information content in the red-
optical (Reiners & Zechmeister, 2020). Thus for M dwarfs, as expected, the
wavelengths with the most RV information are the most telluric-plagued. The
stellar information content was calculated using the following equation from
Butler et al. (1996), where $\sigma_{\bar{v}}$ is the error in mean velocity,
dI/dV is the slope obtained by the spectrometer and $\epsilon_{I}$ is the
fractional Poisson error in the number of photons that arrive into a pixel
$i$:

$\sigma_{\bar{v}}=\frac{1}{\sqrt{\Sigma(\frac{dI/dV}{\epsilon_{I}})^{2}}}$ (2)

For a larger amount of stellar RV information content (and thus lower in m/s),
the RV precision (or the RV limit set by stellar information content) gets
higher and vice versa, as illustrated in Figure 5. As can be seen from Figure
5, the G star sun-like model has very low telluric impact at the wavelengths
with the highest stellar information content. Conversely, the M dwarf has the
lowest telluric impact where there is low stellar information content.

Stellar | Wavelength | No Correction | Division | Division | Modeling
---|---|---|---|---|---
Type | Range (nm) | (m/s) | (m/s) | K-Profile (m/s) | K-Profile (m/s)
M dwarf | 380-900 | 0.167 | 0.125 | 0.120 | 0.140
| 900-2500 | 2.158 | 0.636 | 1.622 | 1.636
G dwarf | 380-900 | 0.035 | 0.021 | 0.021 | 0.032
| 900-2500 | 2.432 | 0.761 | 1.196 | 1.220
Table 3: Summary of the scatter in the nightly RVs retrieved via each method
for an M dwarf star and a G type star (calculated by W22). Column 1-2: Stellar
Type, Wavelength Range (nm). Column 3-6: RV RMS for each method – No
Correction, Division, Division K-Profile, Modeling K-Profile, in m/s. These
two wavelength ranges are for an EarthFinder Equivalent (EFE) instrument on
the ground.

Figure 6: RV signal as a function of time in units of days for the four RV
extraction methods described in Section 2 for the broad optical (left) and the
NIR (right) for an EFE spectrograph. Shown in purple is an RV analog of
Earth’s signal (semi-amplitude of 25 cm/s). We assume SNR = 100 per pixel for
R = 120,000.

In Figure 5, the relationship between RV error induced by tellurics and
stellar RV information content differs for the M dwarf and the G type star.
Whereas for the solar model there is very low RV error in wavelength ranges
with greater stellar RV information content, for the M dwarf the RV error is
still high in areas with high stellar RV information content due to the
tellurics present there as well. A larger stellar information content value
means lower RV precision, and vice versa. In the right plot, this phenomena
can also be seen in terms of telluric absorption line depth, with the M dwarf
in blue and the solar model in yellow. For the solar model, seen in the bottom
panel within the black box, there are orders in the bottom left corner with
little telluric absorption and high stellar RV information content. Notably,
these orders are absent from the lower left corner box in the top plot with
the M dwarf star. In the top plot, the lowest telluric RV error occurs where
there are far less stellar lines, and consequently less precise RVs. This is
also represented in terms of telluric absorption line depth in the right plot
of Figure 5. As seen in this plot, where there is virtually no telluric
absorption and a very high stellar RV information content and low stellar RV
RMS, there are an abundance of solar model orders that are absent for an M
dwarf. As expected, M dwarf RVs are more impacted by tellurics than G type
stars, and we have quantified by how much in Table 3.

We can also study the effectiveness of our extraction methods as a function of
time. We show in Figure 6, the RV signal as a function of time in units of
days for all four RV extraction methods discussed are presented for the broad
optical (350 - 900 nm) in the left panel, and the NIR (0.9 - 2.5 $\rm{\mu m}$)
in the right panel for an EFE spectrograph. Also included in this figure is an
RV analog of Earth’s signal with a semi-amplitude of 25 cm/s for an M dwarf.
This shows us precisely how the RVs in m/s vary as a function of time,
wavelength range, and mitigation method. As can be seen, in the optical, the
RV RMS is small compared to an Earth-like RV signal but when looking at the
NIR, there is significant RV RMS when compared to the Earth-like RV signal. As
expected and previously presented, this portrays the significant challenges
posed by tellurics in the NIR region (Reiners et al., 2018).

Figure 7: Telluric induced RV errors for each photometric band. The horizontal
lines mark the photon-limited precision for SNR=100 per pixel in each band.
The lower ends of the vertical lines mark the RV RMS shown in Table 2 and the
upper ends of the lines mark the range of RV systematics (95% percentile
range) from the simulation methods listed in Table 2.

###  4.2 Effectiveness of Mitigation Methods

In Table 1, we show the effectiveness of different mitigation methods for
different ground-based instruments and their wavelengths. Looking at the last
column, “Photon-limited Precision (m/s),” in comparison to the first column,
“No Correction (m/s)” with instruments operating, we can see that tellurics
did not considerably add to the error budget. Conversely, when we look at the
instruments that operate in the NIR band, tellurics greatly contribute to the
error budget in the optical band. Even division and Forward Modeling cannot
reach close to the photon-limited precision, with a scatter of about 2 m/s.
However, in general the three mitigation methods perform similarly. To explore
this further with greater wavelength granularity, we peer into the different
photometric bands to study the mitigation method effectiveness and photon-
limited precision. In Table 2, we show the simulated RV RMS for multiple
photometric bands. In the corresponding Figure 7, the horizontal lines show
the photometric bands, the range of RV systematics (95th percentile), and the
RV standard deviation of any band are the two ends to every vertical line.
This figure also illustrates that at SNR = 100 per pixel in the B and V bands,
we are photon-noise limited; however tellurics generate larger RV RMS than the
SNR = 100 photon noise would everywhere else. Conversely, this tells us the
SNR at which telluric errors start to matter at each wavelength. Our analysis
confirms a similar result for M dwarfs to that found in Wang et al. (2019), in
that in the V band (roughly 500 - 600 nm), tellurics do not matter as much;
they contribute errors on the order of $\sim$10 cm/s. Instead, the RV
uncertainties for such an M dwarf in the V band is dominated by photon noise,
as shown in Table 2 and Table 1 in Wang et al. (2019).

###  4.3 Additional Caveats to our Simulations

In our experimentation, we do not allow our code to fit for microtellurics
during the first round of Forward Modeling when we fit for the optical depth
for each species. Allowing fitting for microtellurics returns better results
than only turning on species with deep enough telluric lines. During our runs
in which we are fitting an incorrect line profile, fitting different orders
independently gives our code more flexibility to compensate for the incorrect
line profile for every molecular species. However, this approximation is
unrealistic. We assumed a noise-free simulated spectrum, which allows for
microtelluric line fitting down to incredibly small line depths, but in
reality this is not possible. For instance, the water line depth cannot be
recovered very well when they are obscured by photon noise, and the same logic
applies with most microtellurics that are below typical noise levels. For this
reason, we did not present the results of our trials with microtelluric
fitting as our fiducial result.

###  4.4 Empirically Fitting Tellurics with IRTF/iSHELL Data

Our K-profile approach is a good order-of-magnitude approximation to real
observations where tellurics cannot be modeled perfectly, although it
certainly does not represent the accurate and complete picture. In order to
compare the idealistic Forward Modeling simulated scenario to real
observations, we compare the simulated residuals in the Forward Modeling
method with those obtained from iSHELL on the 3.0 m NASA Infrared Telescope
Facility (IRTF) on Maunakea at K band wavelengths (2.2-2.5 $\mu$m). iSHELL is
an R $\sim$ 80,000 spectrograph and has shown 3-5 m/s precision on year long
timescales (Cale et al., 2019). iSHELL also has a methane isotopologue gas
cell (13CH4) for wavelength solution and LSF characterization. As in Cale et
al. (2019), we take high SNR ($>$ 100 per detector pixel near the continuum
and blaze wavelength) spectra of Barnard’s Star and Forward Model them using
an iterative approach to solve for the stellar template using the Python
package pychell. Relevant telluric species for the spectral range analyzed
(2.28-2.41$\mu$m) are $\mathrm{H}_{2}\mathrm{O}$ and $\mathrm{CH}_{4}$, with
water having the largest variability in opacity.

Figure 8: RMS of the flux residuals vs. telluric absorption depth from Forward
Modeling iSHELL spectra at K-band (2.8-2.41 $\micron$) plotted with our
simulations’ chunks that match these wavelengths. Each point represents a
residual flux RMS that corresponds to a small range of telluric absorption
($\approx\pm 1\%$).

For a given echelle order, the two telluric input templates are first raised
to their respective best fit optical depths from the Forward Model fit, and
then convolved with a Gaussian LSF of width 0.1 Å to take into account the
iSHELL line profile. Each modified telluric template is then linearly
interpolated onto the unique wavelength grid for each observation. The
modified water and methane templates are then separately binned using 50
equally spaced bins (in flux) with the minimum and maximum flux values as
endpoints. We match the corresponding residuals to each bin for water and
methane and compute an RMS. We directly sum residuals from (possibly) all
echelle ($m=216-232$) orders and observations since each observation (and
order) are at a similar S/N in the case of iSHELL, and are also identical at
infinite photon SNR for the simulated observations. We perform similar steps
on the simulated observations (across wavelength range 2284–2425 nm) using
chunks 219–227 which approximately match the relevant iSHELL spectral range.
We choose not to separate the methane and water analyses in this instance
since the simulated observations treat the two species identically.

For the simulated residuals, we find the residuals are largest around 50%
telluric absorption suggesting the wings of the telluric lines suffer more
from the “mismatch” of input line profiles, as shown in Figure 8. In the case
of iSHELL, there is no clear correlation with telluric absorption depth for
methane, but water possesses a higher RMS when absorption is $\gtrsim 70$%.
Further, the mismatches are on average larger in the case of the simulated
observations, indicating that our simulated induced line profile error is an
overestimation of the real differences between the TAPAS templates and
observations for Maunakea.

Figure 9: Histograms for the best fit telluric shift parameter for three
iSHELL orders and corresponding best-fit normal distributions to the
histograms (in units of $\mathrm{ms}^{-1}$). Order 5 ($m=216$) has
considerably more features from water than orders 9 ($m=220$) and 17
($m=228$).

Obtaining a robust stellar template to compare iSHELL observations against
pychell relies on isolating the intrinsic stellar spectrum present within the
data by co-adding residuals in the barycenter rest frame where they are nearly
coherent. This process is repeated up to 40 times for a given order before the
RVs tend to stabilize (Cale et al., 2019). Although other spectral features
are forward modeled and are further decoherent in the barycenter frame, we do
not always have the adequate barycenter sampling for a given target to fully
remove any potentially remaining errors. Therefore, we can not assume that
telluric features are fully excluded from the empirical stellar template.

In the case of iSHELL, we suspect the primary limiting factor in RV precision
is the overall complexity of the Forward Model rather than tellurics. The
solver iteratively optimizes 36 parameters, many of which are correlated in
non-intuitive ways. The variable iSHELL blaze function requires local
perturbations that we model with cubic splines, but these are commonly found
to be correlated with the LSF Hermite coefficients and telluric specific
parameters. This implies the blaze and LSF are working together to capture any
mismatches in the telluric lines. The telluric Doppler shift is also commonly
correlated with the stellar shift, although typically not until later
iterations. See Cale et al. (2019) for more details.

Other possibly limitations to this comparison between iSHELL and our
simulations pertain to the telluric Doppler shift for iSHELL and the overall
template “zero-points”. The best fit telluric shift parameters for each iSHELL
order follow a normal wavelength shift distribution; however the mean shifts
across orders are not consistent, as shown in Figure 9. Although we do not yet
fully understand why this is the case, we suspect the higher density of highly
variable water lines contribute to this effect. These offsets also hint that
the telluric model may benefit from allowing for species-specific Doppler
shifts.

##  5 Conclusions & Future Work

To summarize, tellurics add systematic error to our recovered RVs for every
method of RV extraction simulated herein for an M dwarf star, with the
exception of our control, perfect knowledge Forward Model. In the optical,
tellurics add at least 0.1 m/s RMS error, and in the NIR tellurics can add in
excess of 2 m/s, depending on the telluric mitigation method employed. By
comparison, the expected RV semi-amplitude for an Earth-mass planet orbiting
in the HZ of a 0.5 solar mass M dwarf star is $\sim$25 cm/s. This shows the
detrimental effect of tellurics on M dwarf observations, especially since the
most useful RV information content is found in red-optical and NIR
wavelengths. This problem would be mitigated with a space-based instrument as
shown in Plavchan et al. (2020). In the absence of such an instrument,
however, Forward Modeling partially mitigates the problem, bringing the
telluric induced error down to $\sim$ 0.16 m/s.

In future works, we will assess the dependence of our analysis on spectral
resolution, the assumed telluric line profile error, and the impact of time-
varying telluric absorption and airmass on exposures of finite duration. The
latter will become increasingly important for $<$10 cm/s photon noise
observations, as the exposure times will approach durations on the order of
the atmospheric dynamic time-scales. Higher spectral resolution may give us
the ability to mitigate telluric influence in a much better way by resolving
the telluric line profile and providing crucial information about the CCF
stellar line shape deformation from stellar activity. We will also extend the
specific work of this project across the full M star sequence from M0 to M9.

We thank the referee for their thoughtful and constructive review of this
manuscript. N. L. gratefully acknowledges support from an NSF GRFP. N. L. is
the lead author of the paper and performed the analyses and writing presented
herein. All iSHELL data were processed on ARGO, a research computing cluster
provided by the Office of Research Computing, and the exo computer cluster,
both at George Mason University, VA. PPP acknowledge support from NASA
(EarthFinder Probe Mission Concept Study 16-APROBES16-0020, Exoplanet Research
Program Award #80NSSC20K0251, JPL Exoplanet Exploration Program and Research
and Technology Development) and the NSF (Astronomy and Astrophysics Grants
#1716202 and 2006517). This work used the service provided by TAPAS, which is
being supported by CNES (Centre National des Etudes Spatiales) and CNRS
(Centre National de la Recherche Scientifique). TAPAS is a service maintained
by ETHER data center. ETHER acknowledges for TAPAS the use of HITRAN data base
and the LBLRTM radiative transfer code, the use of ECMWF data and the ETHER
data center. We acknowledge useful discussions with Iouli Gordon and Larry
Rothman.

## References

  * Allard et al. (2012a) Allard, F., Homeier, D., & Freytag, B. 2012a, Philosophical Transactions of the Royal Society of London Series A, 370, 2765, doi: 10.1098/rsta.2011.0269
  * Allard et al. (2012b) Allard, F., Homeier, D., Freytag, B., & Sharp, C. M. 2012b, in EAS Publications Series, Vol. 57, EAS Publications Series, ed. C. Reylé, C. Charbonnel, & M. Schultheis, 3–43
  * Artigau et al. (2014) Artigau, É., Astudillo-Defru, N., Delfosse, X., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9149, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 5
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42, doi: 10.1051/0004-6361/201425481
  * Baranne et al. (1996) Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373
  * Benedict et al. (2016) Benedict, G. F., Henry, T. J., Franz, O. G., et al. 2016, AJ, 152, 141, doi: 10.3847/0004-6256/152/5/141
  * Bertaux et al. (2014) Bertaux, J. L., Lallement, R., Ferron, S., Boonne, C., & Bodichon, R. 2014, A&A, 564, A46, doi: 10.1051/0004-6361/201322383
  * Butler et al. (1996) Butler, R. P., Marcy, G. W., Williams, E., et al. 1996, PASP, 108, 500, doi: 10.1086/133755
  * Cale et al. (2019) Cale, B., Plavchan, P., LeBrun, D., et al. 2019, AJ, 158, 170, doi: 10.3847/1538-3881/ab3b0f
  * Cale et al. (2021) Cale, B., Reefe, M., Plavchan, P., et al. 2021, Diving Beneath the Sea of Stellar Activity: Chromatic Radial Velocities of the Young AU Mic Planetary System.  https://arxiv.org/abs/2109.13996
  * Cegla et al. (2019) Cegla, H. M., Watson, C. A., Shelyag, S., Mathioudakis, M., & Moutari, S. 2019, ApJ, 879, 55, doi: 10.3847/1538-4357/ab16d3
  * Cunha et al. (2014) Cunha, D., Santos, N. C., Figueira, P., et al. 2014, A&A, 568, A35, doi: 10.1051/0004-6361/201423723
  * Donati et al. (2018) Donati, J.-F., Kouach, D., Lacombe, M., et al. 2018, SPIRou: A NIR Spectropolarimeter/High-Precision Velocimeter for the CFHT, ed. H. J. Deeg & J. A. Belmonte, 107
  * Dressing & Charbonneau (2013) Dressing, C. D., & Charbonneau, D. 2013, ApJ, 767, 95, doi: 10.1088/0004-637X/767/1/95
  * Dumusque (2018) Dumusque, X. 2018, A&A, 620, A47, doi: 10.1051/0004-6361/201833795
  * Dumusque et al. (2011) Dumusque, X., Udry, S., Lovis, C., Santos, N. C., & Monteiro, M. J. P. F. G. 2011, A&A, 525, A140, doi: 10.1051/0004-6361/201014097
  * Fischer et al. (2016) Fischer, D. A., Anglada-Escude, G., Arriagada, P., et al. 2016, PASP, 128, 066001, doi: 10.1088/1538-3873/128/964/066001
  * Gullikson et al. (2014) Gullikson, K., Dodson-Robinson, S., & Kraus, A. 2014, AJ, 148, 53, doi: 10.1088/0004-6256/148/3/53
  * Haywood et al. (2020) Haywood, R. D., Milbourne, T. W., Saar, S. H., et al. 2020, arXiv e-prints, arXiv:2005.13386.  https://arxiv.org/abs/2005.13386
  * Henry et al. (2018) Henry, T. J., Jao, W.-C., Winters, J. G., et al. 2018, AJ, 155, 265, doi: 10.3847/1538-3881/aac262
  * Howard & Fulton (2016) Howard, A. W., & Fulton, B. J. 2016, PASP, 128, 114401, doi: 10.1088/1538-3873/128/969/114401
  * Jurgenson et al. (2016) Jurgenson, C., Fischer, D., McCracken, T., et al. 2016, Ground-based and Airborne Instrumentation for Astronomy VI, doi: 10.1117/12.2233002
  * Kanodia et al. (2020) Kanodia, S., Ninan, J. P., Monson, A. J., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11447, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 1144740
  * Kurucz (2005) Kurucz, R. L. 2005, Memorie della Societa Astronomica Italiana Supplementi, 8, 14
  * Mahadevan et al. (2014) Mahadevan, S., Ramsey, L. W., Terrien, R., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, ed. S. K. Ramsay, I. S. McLean, & H. Takami, 91471G
  * Marchwinski et al. (2015) Marchwinski, R. C., Mahadevan, S., Robertson, P., Ramsey, L., & Harder, J. 2015, ApJ, 798, 63, doi: 10.1088/0004-637X/798/1/63
  * Morales et al. (2019) Morales, J. C., Mustill, A. J., Ribas, I., et al. 2019, Science, 365, 1441, doi: 10.1126/science.aax3198
  * Pepe et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96, doi: 10.1051/0004-6361/202038306
  * Plavchan et al. (2015) Plavchan, P., Latham, D., Gaudi, S., et al. 2015, Radial Velocity Prospects Current and Future: A White Paper Report prepared by the Study Analysis Group 8 for the Exoplanet Program Analysis Group (ExoPAG).  https://arxiv.org/abs/1503.01770
  * Plavchan et al. (2020) Plavchan, P., Vasisht, G., Beichman, C., et al. 2020, EarthFinder Probe Mission Concept Study: Characterizing nearby stellar exoplanet systems with Earth-mass analogs for future direct imaging.  https://arxiv.org/abs/2006.13428
  * Quirrenbach & Carmenes Consortium (2016) Quirrenbach, A., & Carmenes Consortium. 2016, in 41st COSPAR Scientific Assembly, Vol. 41, E1.21–11–16
  * Raghavan et al. (2010) Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1, doi: 10.1088/0067-0049/190/1/1
  * Reiners et al. (2010) Reiners, A., Bean, J. L., Huber, K. F., et al. 2010, ApJ, 710, 432, doi: 10.1088/0004-637X/710/1/432
  * Reiners & Zechmeister (2020) Reiners, A., & Zechmeister, M. 2020, ApJS, 247, 11, doi: 10.3847/1538-4365/ab609f
  * Reiners et al. (2018) Reiners, A., Zechmeister, M., Caballero, J. A., et al. 2018, A&A, 612, A49, doi: 10.1051/0004-6361/201732054
  * Robertson et al. (2020) Robertson, P., Stefansson, G., Mahadevan, S., et al. 2020, ApJ, 897, 125, doi: 10.3847/1538-4357/ab989f
  * Seifahrt et al. (2010) Seifahrt, A., Käufl, H. U., Zängl, G., et al. 2010, A&A, 524, A11, doi: 10.1051/0004-6361/200913782
  * Seifahrt et al. (2021) Seifahrt, A., Bean, J. L., Stürmer, J., et al. 2021, arXiv e-prints, arXiv:2106.02157.  https://arxiv.org/abs/2106.02157
  * Smette et al. (2015) Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77, doi: 10.1051/0004-6361/201423932
  * Tamura et al. (2012) Tamura, M., Suto, H., Nishikawa, J., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84461T
  * Tran et al. (2021) Tran, Q. H., Bowler, B. P., Cochran, W. D., et al. 2021, arXiv e-prints, arXiv:2101.11005.  https://arxiv.org/abs/2101.11005
  * Tuomi et al. (2013) Tuomi, M., Jones, H. R. A., Jenkins, J. S., et al. 2013, A&A, 551, A79, doi: 10.1051/0004-6361/201220509
  * Ulmer-Moll et al. (2019) Ulmer-Moll, S., Figueira, P., Neal, J. J., Santos, N. C., & Bonnefoy, M. 2019, A&A, 621, A79, doi: 10.1051/0004-6361/201833282
  * Wang et al. (2019) Wang, S., Wright, J. T., Bender, C., et al. 2019, AJ, 158, 216, doi: 10.3847/1538-3881/ab48e1
  * Wang et al. (2022 in prep) Wang, S. X., Latouf, N., Plavchan, P., Cale, B., & Cullen, B. 2022 in prep
  * Zechmeister et al. (2018) Zechmeister, M., Reiners, A., Amado, P. J., et al. 2018, A&A, 609, A12, doi: 10.1051/0004-6361/201731483