# Probing the solar interior with lensed gravitational waves from known
pulsars

Ryuichi Takahashi1, Soichiro Morisaki2 and Teruaki Suyama3 1 Faculty of
Science and Technology, Hirosaki University, 3 Bunkyo-cho, Hirosaki, Aomori
036-8561, Japan
2 Institute for Cosmic Ray Research (ICRR), KAGRA Observatory, The University
of Tokyo, Kashiwa City, Chiba 277-8582, Japan
3 Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama,
Meguro-ku, Tokyo 152-8551, Japan

###### Abstract

When gravitational waves (GWs) from a spinning neutron star arrive from behind
the Sun, they are subjected to gravitational lensing that imprints a
frequency-dependent modulation on the waveform. This modulation traces the
projected solar density and gravitational potential along the path as the Sun
passes in front of the neutron star. We calculate how accurately the solar
density profile can be extracted from the lensed GWs using a Fisher analysis.
For this purpose, we selected three promising candidates (the highly spinning
pulsars J1022+1001, J1730-2304, and J1745-23) from the pulsar catalog of the
Australia Telescope National Facility. The lensing signature can be measured
with $3\sigma$ confidence when the signal-to-noise ratio (SNR) of the GW
detection reaches $100\,(f/300{\rm Hz})^{-1}$ over a one-year observation
period (where $f$ is the GW frequency). The solar density profile can be
plotted as a function of radius when the SNR improves to $\gtrsim 10^{4}$.

###### Subject headings:

gravitational lensing: weak — gravitational waves — Sun: general

##  1\. Introduction

Since the direct detection of a gravitational-wave (GW) signal from a merging
black hole binary (GW150914; Abbott et al., 2016), GW astronomy has attracted
increasing interest. Currently, all GW sources detected by the ground-based
detectors of the Laser Interferometer Gravitational-Wave Observatory (LIGO),
the Virgo interferometer, and the Kamioka Gravitational Wave Detector (KAGRA)
are compact binary coalescences (CBCs) of stellar-mass black holes and neutron
stars. The LIGO-Virgo-KAGRA (LVK) collaboration reported $\sim 90$ candidates
of CBC events (The LIGO Scientific Collaboration et al., 2021). However, the
ground-based detectors are also expected to detect spinning neutron stars
(e.g., Abbott et al., 2022a). A neutron star that is non-axisymmetric around
its spin axis, resulting in a so-called mountainous profile, emits continuous
GWs. Non-axisymmetry may be caused by crustal deformation, magnetic fields,
and mass accretion from the star’s companion (e.g., see recent reviews by
Glampedakis & Gualtieri (2018) and Riles (2022)). Although the LVK
collaboration has been searching for continuous GW signals, no event has yet
been reported (Abbott et al., 2021, 2022a, 2022b; The LIGO Scientific
Collaboration et al., 2022).

If a GW signal encounters the Sun along its path, gravitational lensing
imprints a frequency-dependent modulation on the waveform. In geometrical
optics (i.e., the zero-wavelength limit of GWs), the solar density modulates
the amplitude with a magnification effect, while the gravitational potential
modulates the phase by imposing a potential (or Shapiro) time delay. These
modulations can be obtained along the transversal path of the Sun moving in
front of the source (the duration of this movement is approximately half a
day). Therefore, in principle, one can probe the solar interior using the
lensed signal. GW lensing by the Sun has been studied as a tool for amplifying
the strain amplitude of GWs from a distant source and probing the solar
structure (e.g., Cyranski & Lubkin, 1974; Sonnabend, 1979; Patla & Nemiroff,
2008; Marchant et al., 2020). Before the 1980s, these studies were based on
geometrical optics. Bontz & Haugan (1981) first demonstrated that the
diffraction effect caused by the finite wavelength of GWs suppresses
magnification near the focal point at lower frequencies ($f\lesssim
10^{4}\,{\rm Hz}$). Recently, Marchant et al. (2020) proposed that the solar
structure can be probed through the lensed GW signal from a pulsar behind the
Sun. However, they only roughly estimated the detectability of the lensing
signature, without calculating the measurement accuracy of the density
profile. Very recently, Jung & Kim (2022) reported that the Fresnel scale is
comparable to the solar radius within the frequency band of the ground-based
detectors; therefore, the density profile can (in principle) be probed with
the lensing modulation on a chirp signal from a CBC.

In this study, we investigate the accuracy of measuring the solar density
profile using the lensed signals of known pulsars. Because ground-based
detectors detect GW wavelengths longer than (or comparable to) the solar
Schwarzschild radius, gravitational lensing should employ wave optics (e.g.,
Ohanian, 1974; Bliokh & Minakov, 1975; Schneider et al., 1992; Nakamura &
Deguchi, 1999; Takahashi & Nakamura, 2003; Dai et al., 2018; Oguri, 2019;
Suvorov, 2022; Liao et al., 2022). We first calculate the lensed waveform
based on wave optics. We discuss the effects of frequency and impact parameter
on the waveform (Section 2). Then we extract the known pulsars crossing behind
the Sun from the Australia Telescope National Facility (ATNF) pulsar catalog
(Manchester et al., 2005). From the extracted list, we select suitable
candidates by calculating the lensing modulations of these pulsars (Section
3). Using these candidates, we calculate the detectability of the lensing
signature and the accuracy of measuring the solar density profile through a
Fisher analysis (Section 4). We roughly estimate the number of Galactic
millisecond pulsars (MSPs) behind the Sun, which are potentially detectable by
near-future radio surveys (Section 5). We discuss the lensing by a Galactic
star, which may be confused with the lensing by the Sun (Section 6). The study
findings are summarized in Section 7.

##  2\. Amplitude and phase modulations of solar lensing

This section introduces the amplitude and phase modulations of a lensed
waveform, based on wave optics.

###  2.1. Lensed waveform

Figure 1.— Configuration of the detector, Sun, and source. $D_{\rm L}$,
$D_{\rm LS}$, and $D_{\rm S}$ denote the distances between them. Solid line is
the GW path. Angular coordinates $\bm{\theta}$ and $\bm{\theta}_{\rm s}$
describe the incoming direction of the GWs and source position, respectively,
with respect to the solar center.

Figure 1 shows the configuration of the detector, Sun, and source. The solid
line describes a single GW path111 In wave optics, the lensed waveform is the
superposition of many GW paths from the source to the detector (Eq. (2)). The
solid line represents one of them.. We assume the thin-lens approximation, in
which the solar density and its gravitational potential are projected onto a
plane perpendicular to the line of sight. The thin-lens approximation is valid
within an error less than $0.5\,\%$ ($\simeq R_{\odot}/D_{\rm L}$) on the
lensed waveform (Appendix A; Suyama et al., 2005). Furthermore, we assume the
flat-sky approximation, in which $|\bm{\theta}|,|\bm{\theta}_{\rm s}|\ll 1$.
As the coordinate system is fixed at the solar position, the background source
(pulsar) moves along the ecliptic longitude, following the annual solar
motion. The distance to the Sun is $D_{\rm L}=1\,{\rm au}$ and the solar
angular radius is $\theta_{\odot}=\arctan[R_{\odot}/(1\,{\rm au})]\simeq 16$
arcmin. The annual modulation of $D_{\rm L}$ is ignored because it changes by
$<2\,\%$ as the Earth moves elliptically around the Sun.

The lensed waveform $\tilde{h}^{\rm L}(f)$ in the frequency domain is obtained
by multiplying the unlensed waveform $\tilde{h}(f)$ by a function
$F(f;\bm{\theta}_{\rm s})$ (e.g., Nakamura, 1998; Takahashi & Nakamura, 2003):

$\displaystyle\tilde{h}^{\rm L}(f;\bm{\theta}_{\rm s})$
$\displaystyle=F(f;\bm{\theta}_{\rm s})\tilde{h}(f),$
$\displaystyle=\left[1+A_{F}(f;\bm{\theta}_{\rm s})\right]{\rm e}^{{\rm
i}\varPhi_{F}(f;\bm{\theta}_{\rm s})}\tilde{h}(f).$ (1)

In the second line, $A_{F}\,(\equiv|F|-1)$ and
$\varPhi_{F}\,(\equiv\ln(F/|F|))$ represent the amplitude and phase
modulations, respectively. The GW polarization is ignored because the
polarization rotation due to lensing is negligibly small (e.g., Hou et al.,
2019; Ezquiaga et al., 2021; Dalang et al., 2022). The function $F$, called
the amplification factor or transmission factor, is given by the following
diffraction integral (e.g., Schneider et al., 1992):

$F(f;\bm{\theta}_{\rm s})=\frac{D_{\rm L}D_{\rm S}}{cD_{\rm LS}}\frac{f}{\rm
i}\int\\!{\rm d}^{2}\theta\exp\left[2\pi{\rm i}ft_{\rm
d}(\bm{\theta},\bm{\theta}_{\rm s})\right],$ (2)

where the time delay is given by

$t_{\rm d}(\bm{\theta},\bm{\theta}_{\rm s})=\frac{1}{c}\left[\frac{D_{\rm
L}D_{\rm S}}{2D_{\rm LS}}\left|\bm{\theta}-\bm{\theta}_{\rm
s}\right|^{2}-\frac{\hat{\psi}(\bm{\theta})}{c^{2}}\right].$ (3)

The first and second terms are the geometrical and potential (or Shapiro) time
delays, respectively. Because the distance to the Sun is much shorter than the
distance to the source, hereafter we assume $D_{\rm S}/D_{\rm LS}\simeq 1$.
The two-dimensional lens potential $\hat{\psi}(\bm{\theta})$ is determined
from the solar-projected density profile $\Sigma(\bm{\theta})$ using the
Poisson equation:

$\nabla^{2}_{\theta}\hat{\psi}(\bm{\theta})=8\pi GD^{2}_{\rm
L}\Sigma(\bm{\theta}).$ (4)

Outside the Sun, the potential is identical to that of a point mass:

$\hat{\psi}(\theta)=4GM_{\odot}\ln\left(\frac{\theta}{\theta_{\rm
E}}\right)~{}~{}{\rm for}~{}\theta>\theta_{\odot},$ (5)

where $\theta_{\rm E}$ is an arbitrary constant, here set to the angular
Einstein radius; i.e., $\theta_{\rm E}=[4GM_{\odot}/(c^{2}D_{\rm
L})]^{1/2}\simeq 0.043\,\theta_{\odot}$. The Sun is modeled using the BS05(OP)
spherical density model222The numerical table can be downloaded from
http://www.sns.ias.edu/~jnb/SNdata/sndata.html#bs2005. in Bahcall et al.
(2005). The model density agrees within $2\,\%$ of the helioseismological
results across the whole radius (their Fig. 1). It also agrees within $3\,\%$
and $4\,\%$ of the recent solar density models of B16-GS98 and B16-AGSS09met,
respectively, in Vinyoles et al. (2017). As the density profile is spherical,
the integration (2) can be performed over the azimuth and $F$ reduces to

$\displaystyle F(f;\theta_{\rm s})$ $\displaystyle=\frac{2\pi fD_{\rm L}}{{\rm
i}c}\int_{0}^{\infty}\\!\\!{\rm d}\theta\,\theta J_{0}\left(\frac{2\pi fD_{\rm
L}}{c}\theta\theta_{\rm s}\right)$ $\displaystyle\times\exp\left[2\pi{\rm
i}f\left\\{\frac{1}{2c}D_{\rm L}\left(\theta^{2}+\theta_{\rm
s}^{2}\right)-\frac{\hat{\psi}(\theta)}{c^{3}}\right\\}\right],$ (6)

where $J_{0}$ is the zero-order Bessel function. Equation (6) is numerically
integrated using integration by parts (Appendix A of Takahashi (2004); Guo &
Lu (2020)).

Let us comment on the validity of the spherical density model of the Sun. The
spin rotation makes the equatorial radius larger than the polar radius, but
the observed relative difference between these surface radii is very small
$\approx 10^{-5}$ (Chapter 7 of Stix, 2004). The quadrupole mass moment of the
Sun ($J_{2}$)333Modeling the Sun as a spheroidal body about its spin axis ($z$
axis), $J_{2}$ is simply $J_{2}=(I_{zz}-I_{xx})/(M_{\odot}R_{\odot}^{2})$ with
the moment of inertia $I_{ij}$. A positive (negative) $J_{2}$ means an oblate
(prolate) shape. generates the precession of Mercury’s orbital perihelion
(which is much smaller than the known general-relativistic effect). The
precession was precisely measured by the MESSENGER (MErcury Surface, Space
ENviroment, GEochemistry, and Ranging)
spacecraft444https://solarsystem.nasa.gov/missions/messenger/in-depth/, which
gave a constraint on $J_{2}\simeq 2\times 10^{-7}$ (Park et al., 2017). These
observations suggest that $F$ for the spherical model would be valid with an
error of $\lesssim 10^{-5}$.

###  2.2. High- and low-frequency limits

Figure 2.— Solar deflection angle, convergence, and potential (or Shapiro)
time delay as functions of angular radius $\theta$ (normalized by the solar
radius $\theta_{\odot}$). The convergence $\kappa$ in the middle panel is the
dimensionless projected density profile. The potential time delay in the right
panel is the projected gravitational potential divided by $c^{3}$. Near and
outside the surface ($\theta/\theta_{\odot}\gtrsim 1$), the result approaches
that of Eq. (5) for a point mass (the dashed line).

This subsection presents the high- and low-frequency limits of the
amplification factor. In the high-frequency limit (i.e., the geometrical-
optics limit), because the exponential term of $F$ in Eq. (2) oscillates
violently, a stationary point of $t_{\rm d}(\bm{\theta},\bm{\theta}_{\rm s})$
contributes to the integral (e.g., Schneider et al., 1992). This point is a
solution of $\nabla_{\theta}t_{\rm d}(\bm{\theta},\bm{\theta}_{\rm s})=0$ :

$\bm{\theta}_{\rm s}=\bm{\theta}-\hat{\bm{\alpha}}(\bm{\theta}),$ (7)

from Eq. (3), where $\hat{\bm{\alpha}}=\nabla_{\theta}\hat{\psi}/(c^{2}D_{\rm
L})$ is the deflection angle. Equation (7) is the lens equation. The image
position $\bm{\theta}_{\rm i}$ is obtained by solving Eq. (7). Note that only
a single image is formed irrespective of $\bm{\theta}_{\rm s}$. In the high-
frequency limit, $F$ reduces to the following simple form:

$\lim_{f\rightarrow\infty}F(f;\bm{\theta}_{\rm s})=\left|\mu(\bm{\theta}_{\rm
i},\bm{\theta}_{\rm s})\right|^{1/2}\exp\left[2\pi ift_{\rm
d}(\bm{\theta}_{\rm i},\bm{\theta}_{\rm s})\right],$ (8)

where $\mu=[\,\det(\partial\bm{\theta}_{\rm s}/\partial\bm{\theta}_{\rm
i})\,]^{-1}$ is the magnification of the image.

In the weak-gravitational-field limit, the magnification can be approximated
as $\mu\simeq 1+2\kappa$, where $\kappa=\Sigma/\Sigma_{\rm cr}$ is the
convergence and $\Sigma_{\rm cr}=c^{2}/(4\pi GD_{\rm L})$ is the critical
density (e.g., Bartelmann & Schneider, 2001). Accordingly, $A_{F}$ and
$\varPhi_{F}$ reduce to

$\displaystyle\lim_{f\rightarrow\infty}\\!A_{F}(f;\bm{\theta}_{\rm s})$
$\displaystyle\simeq\kappa(\bm{\theta}_{\rm s})=\frac{4\pi GD_{\rm
L}}{c^{2}}\Sigma(\bm{\theta}_{\rm s}),$
$\displaystyle\lim_{f\rightarrow\infty}\\!\varPhi_{F}(f;\bm{\theta}_{\rm s})$
$\displaystyle=2\pi ft_{\rm d}(\bm{\theta}_{\rm i},\bm{\theta}_{\rm
s})\simeq-2\pi f\frac{\hat{\psi}(\bm{\theta}_{\rm s})}{c^{3}}.$ (9)

Therefore, $A_{F}$ and $\varPhi_{F}$ trace the projected density profile and
gravitational potential, respectively.

Figure 2 plots the deflection angle, convergence, and gravitational potential
as functions of $\theta$ (as these variables are functions of $\theta$, we
plot $\theta$ rather than $\theta_{\rm s}$ along the x-axis). The deflection
angle is consistent with the known result
$\hat{\alpha}=4GM_{\odot}/(c^{2}R_{\odot})\simeq 1.75\,{\rm arcsec}$ at the
surface. The maximum deflection is $\hat{\alpha}\simeq 4.57\,{\rm arcsec}$ at
$\theta\simeq 0.23\,\theta_{\odot}$. Because $\hat{\alpha}$ is much smaller
than the solar radius ($\theta_{\odot}\simeq 16\,{\rm arcmin}$), we can safely
set $\bm{\theta}_{\rm i}\simeq\bm{\theta}_{\rm s}$ from Eq. (7). Therefore,
the second term of $t_{\rm d}$ in Eq. (3) exceeds the first term in the
geometrical-optics limit. Meanwhile, the convergence $\kappa$ is maximized at
$\sim 0.04$ near the center. The profile is flat at the core ($\theta\lesssim
0.03\,\theta_{\odot}$), but it drops steeply at $\theta\gtrsim
0.1\,\theta_{\odot}$. The typical gravitational time delay is $\approx
10^{-5}$ s. Near the center, the potential in Eq. (4) can be approximated as
$\hat{\psi}(\theta)$$={\rm const.}$$+(c^{2}D_{\rm
L}/2)\,\kappa(\theta\\!=\\!0)\,\theta^{2}$$+\mathcal{O}(\theta^{3})$. The
results of Fig. 2 are consistent with previous works on lensing by the
transparent Sun (Bontz & Haugan, 1981; Patla & Nemiroff, 2008; Marchant et
al., 2020; Jung & Kim, 2022).

In the large-angle limit $\theta_{\rm s}\gg\theta_{\odot}$, the geometrical-
optics result (9) is recovered because $2\pi ft_{\rm d}\gg 1$. As the
gravitational potential in Eq. (5) is valid only for $\theta_{\odot}<\theta\ll
1$, it must be replaced with (e.g., Backer & Hellings, 1986)

$\hat{\psi}(\theta)=2GM_{\odot}\ln\left\\{\frac{2(1-\cos\theta)}{\theta_{\rm
E}^{2}}\right\\},$ (10)

when $\theta\gtrsim 1$.

We now investigate the low-frequency limit. Changing the variables to
$\bm{\theta}^{\prime}=\sqrt{f}\bm{\theta}$ and $\bm{\theta}_{\rm
s}^{\prime}=\sqrt{f}\bm{\theta}_{\rm s}$ in Eq. (2), $F$ is rewritten as

$\displaystyle F(f;\bm{\theta}_{\rm s})$ $\displaystyle=\frac{D_{\rm L}}{{\rm
i}c}\int\\!{\rm d}^{2}\bm{\theta}^{\prime}\exp\left[\frac{2\pi{\rm
i}}{c}\left\\{\frac{D_{\rm
L}}{2}\left|\bm{\theta}^{\prime}-\bm{\theta}^{\prime}_{\rm
s}\right|^{2}\right.\right.$
$\displaystyle~{}~{}~{}~{}\left.\left.-\frac{f}{c^{2}}\hat{\psi}\left(\frac{\bm{\theta}^{\prime}}{\sqrt{f}}\right)\right\\}\right].$
(11)

In the low-frequency limit, $\hat{\psi}$ can be replaced with its value at
infinity:

$\lim_{f\rightarrow
0}\hat{\psi}\left(\frac{\theta^{\prime}}{\sqrt{f}}\right)=4GM_{\odot}\ln\left(\frac{\theta^{\prime}}{\sqrt{f}\theta_{\rm
E}}\right),$

from Eq. (5). Inserting this $\hat{\psi}$ into Eq. (11) and expanding $f$ as a
Taylor series, we have

$\displaystyle\lim_{f\rightarrow 0}A_{F}(f;\bm{\theta}_{\rm s})$
$\displaystyle=\frac{2\pi^{2}GM_{\odot}}{c^{3}}f,$
$\displaystyle\lim_{f\rightarrow 0}\varPhi_{F}(f;\bm{\theta}_{\rm s})$
$\displaystyle=\frac{4\pi GM_{\odot}}{c^{3}}f\left[\gamma+\ln\left(\frac{4\pi
GM_{\odot}}{c^{3}}f\right)\right],$ (12)

where $\gamma=0.5772\cdots$ is Euler’s constant.

###  2.3. Amplitude and phase modulations

Figure 3.— Amplitude and phase modulations as functions of frequency at
various source positions $\theta_{\rm s}/\theta_{\odot}=0$–$2$. Horizontal
dashed lines and dot-dashed orange line represent the analytical results in
the high- and low-frequency limits, respectively, given by Eqs. (9) and (12).
Up arrows indicate the rough boundary between the wave- and geometrical-optics
regions for each $\theta_{\rm s}$ (see the description after Eq. (13) in the
text). In the left panel, dotted curves denote negative values ($A_{F}<0$). In
the right panel, the phase modulation is divided by $2\pi f$, which
corresponds to an arrival time delay.

Figure 4.— Same as Fig.3, but plotted as functions of $\theta_{\rm
s}/\theta_{\odot}$ for various $f$. Red curve overlaps the purple curve
everywhere except near the center ($\theta_{\rm s}/\theta_{\odot}\lesssim
0.2$).

Figure 3 plots $A_{F}$ and $\varPhi_{F}$ as functions of $f$ for various
$\theta_{\rm s}$. In the high- and low-frequency limits, the plots are
consistent with the analytical results of Eqs. (9) and (12), respectively. As
the frequency of a chirp signal from an inspiral binary sweeps from low to
high, $A_{F}$ increases proportionally to $f$ in the wave-optics regime but
approaches a constant value of $\kappa$ in the geometrical-optics regime.
Similarly, $-\varPhi_{F}/(2\pi f)$ decreases logarithmically, but approaches a
constant value of the potential time delay ($\hat{\psi}/c^{3}$). These
modulations are observed even outside the Sun ($\theta_{\rm
s}=2\theta_{\odot}$). Therefore, within the frequency range of ground-based
detectors, the solar lensing imprints unique frequency-dependent modulations
on both the amplitude and phase. The solar structure can (in principle) be
extracted from the chirp signal (see the recent study by Jung & Kim, 2022).
However, in the high-frequency limit, solar structure extraction from the
signal (even from a chirp signal) is impossible at a single $\theta_{\rm s}$
because the constant $A_{F}$ is degenerate with an intrinsic amplitude of GWs.
Similarly, the constants $\varPhi_{F}$ and $\varPhi_{F}/(2\pi f)$ are
degenerate with an intrinsic phase and arrival time, respectively. Therefore,
in geometrical optics, the lensing modulations are indistinguishable from the
intrinsic source properties (unless lensed signals with various $\theta_{\rm
s}$ from a single source are available). Figure 3 is consistent with Fig. 2 of
Jung & Kim (2022).

To discuss the behavior of these modulations in more detail, let us introduce
the angular Fresnel scale (Macquart, 2004; Takahashi, 2006):

$\displaystyle\theta_{\rm F}$ $\displaystyle=\left(\frac{c}{2\pi
f}\frac{1}{D_{\rm L}}\right)^{1/2},$ $\displaystyle\simeq
0.38\,\theta_{\odot}\left(\frac{f}{100\,{\rm Hz}}\right)^{-1/2}.$ (13)

In the low-frequency limit, the lensed signal probes a circle of radius
$\theta_{\rm F}$ around the source position (Section 2 and Fig. 1 of Choi et
al., 2021). In other words, the Fresnel scale can be interpreted as the
effective source radius (Oguri & Takahashi, 2020). The effective source radius
shrinks with increasing $f$ in the chirp signal. In the high-frequency limit,
the lensing probes a small region around the image position. As the
modulations at different $f$ can probe different regions, one can in principle
probe the density profile (Jung & Kim, 2022). Here, $\theta_{\rm
F}=\theta_{\rm s}$ (unless $\theta_{\rm s}=0$) roughly represents the boundary
between the wave- and geometrical-optics regions. This boundary condition is
derived from $2\pi ft_{\rm d}=1$ with neglecting the potential term in Eq. (3)
(Section II.B of Choi et al. (2021) and Section 3 of Jung & Kim (2022))555This
condition can also be derived from the intersection of the low- and high-
frequency limits of the phase modulation. Setting $\lim_{f\rightarrow
0}\varPhi_{F}=\lim_{f\rightarrow\infty}\varPhi_{F}$ in Eqs. (9) and (12) with
the point-mass potential (5) (which is valid only for $\theta_{\rm s}\gtrsim
0.3\,\theta_{\odot}$ in the right panel of Fig. 2), one obtains $\theta_{\rm
F}\simeq 0.94\,\theta_{\rm s}$.. The boundary is plotted as a series of up
arrows in Figs. 3 and 4.

Figure 4 plots $A_{F}$ and $\varPhi_{F}$ as functions of $\theta_{\rm s}$ for
various fixed $f$. As the GWs from a pulsar are continuous, one can (in
principle) measure $A_{F}$ and $\Phi_{F}$ as functions of $\theta_{\rm s}$ for
a pulsar moving behind the Sun. For smaller (larger) $\theta_{\rm s}$, the
exponential term of the amplification factor ($=2\pi ft_{\rm d}$) is smaller
(larger) and the results approach the low-frequency (high-frequency) limit in
Eq. (12) (Eq. (9)). At $f=1000$ Hz, the result almost matches that of
geometrical optics, except in the near-center region ($\theta_{\rm
s}/\theta_{\odot}\lesssim 0.2$). As $\Phi_{F}$ reflects the gravitational
potential, it can be measured even outside the Sun. However, $A_{F}$, which
reflects the mass density, falls steeply near the surface. Figure 4 is
consistent with Figs. 6 and 7 of Marchant et al. (2020).

Figure 5.— Contour plot of the difference between the phase modulation and its
geometrical-optics limit ($\Phi_{F}-2\pi ft_{\rm d}$). Red, blue, and white
regions correspond to positive, negative, and zero differences, respectively.
Dotted line plots the rough boundary between the wave- and geometrical-optics
regions, calculated from $\theta_{\rm F}=\theta_{\rm s}$ (the description
after Eq. (13) in the text).

The LVK analysis accounts for the Shapiro time delay imposed by the Sun. In
searches of continuous GWs (Abbott et al., 2017, 2019, 2022a), the time delay
has been calculated using the TEMPO (Nice et al., 2015) and TEMPO2 (Edwards et
al., 2006; Hobbs et al., 2006) packages developed for timing analyses of
pulsar radio signals. However, as shown in Figs. 3 and 4, the Shapiro time
delay is less accurate at lower $f$ and smaller $\theta_{\rm s}$. To clarify
the accuracy of the Shapiro time delay, Fig. 5 plots a contour map of the
difference between the phase modulation and its geometrical-optics limit. The
difference ranges from $-0.015$ to $0.002$ rad and is most significant at
$f\simeq 100$–$1000\,{\rm Hz}$ near the center ($\theta_{\rm s}\lesssim
0.2\,\theta_{\odot}$). At lower frequencies ($f\lesssim 10\,{\rm Hz}$), the
difference reduces because both $\varPhi_{F}$ and $2\pi ft_{\rm d}$ are
proportional to $f$. As expected, the Shapiro time delay is recovered at high
$f$ or large $\theta_{\rm s}$. The dotted line, predicted by setting
$\theta_{\rm F}=\theta_{\rm s}$, is apparently consistent with the boundary.
Although the Shapiro time delay is less accurate inside the Sun, the error is
small ($0.015$ rad at most).

##  3\. Known pulsars

Table 1Pulsars with $f\\!>\\!10\,{\rm Hz}$ crossing behind the Sun (upper
four) and close to, but not behind, the Sun (lower eight)

name |  $f\,$(Hz) | $\theta_{\rm{\scriptscriptstyle LAT}}/\theta_{\odot}$ |  $D_{\rm s}\,$(kpc) | $h_{0}^{95\%}$ | $\epsilon^{95\%}$ | references
---|---|---|---|---|---|---
J1745$-$23 | $369$ | $-0.07$ | $7.94$ | —— | —— | C20
J1022+1001 | $122$ | $-0.24$ | $0.64$ | $7.7\\!\times\\!10^{-27}$ | $3.2\\!\times\\!10^{-7}$ | C96,R21
J1809$-$2332 | $14$ | $-0.45$ | $0.88$ | —— | —— | A09,R11
J1730$-$2304 | $246$ | $0.71$ | $0.47$ | $5.0\\!\times\\!10^{-27}$ | $3.7\\!\times\\!10^{-8}$ | L95,R21
J1858-2216 | $839$ | $1.80$ | $0.92$ | $7.8\\!\times\\!10^{-27}$ | $9.6\\!\times\\!10^{-9}$ | S16
J1142+0119 | $394$ | $-1.84$ | $2.17$ | $7.4\\!\times\\!10^{-27}$ | $9.8\\!\times\\!10^{-8}$ | S16
J2310-0555 | $766$ | $-1.99$ | $1.56$ | —— | —— | S16
J1756-2251 | $70$ | $2.14$ | $0.73$ | $6.1\\!\times\\!10^{-27}$ | $8.6\\!\times\\!10^{-7}$ | F04,F14
J1646-2142 | $342$ | $2.44$ | $0.97$ | —— | —— | R12,R13
J1811-2405 | $752$ | $-2.52$ | $1.83$ | $1.5\\!\times\\!10^{-26}$ | $4.5\\!\times\\!10^{-8}$ | K10,N20
J1836-2354B | $619$ | $-2.74$ | $3.20$ | —— | —— | L11
J1836-2354A | $596$ | $-2.78$ | $3.20$ | —— | —— | L11

Note. — Second column: GW frequency. Third column: ecliptic latitude
$\theta_{\rm{\scriptscriptstyle LAT}}$ (i.e., the minimum angular separation
to the solar center). Pulsars are listed in order of increasing
$|\theta_{\rm{\scriptscriptstyle LAT}}|$. Fourth column: distance $D_{\rm s}$,
measured by parallax for J1022+1001, J1730-2304, and J1756-2251 and estimated
from the dispersion measures with the Galactic free-electron distribution
model YMW16 (Yao et al., 2017) for the other pulsars. Fifth and sixth columns:
$95\%$ upper limits on the strain amplitude and ellipticity, respectively,
taken from the LVK results (Abbott et al., 2022a). Last column: references A09
(Abdo et al., 2009), C20 (Cameron et al., 2020), C96 (Camilo et al., 1996),
F04 (Faulkner et al., 2004), F14 (Ferdman et al., 2014), K10 (Keith et al.,
2010), L95 (Lorimer et al., 1995), L11 (Lynch et al., 2011), N20 (Ng et al.,
2020), R11 (Ray et al., 2011), R12 (Ray et al., 2012), R21 (Reardon et al.,
2021), R16 (Roy & Bhattacharyya, 2013), and S16 (Sanpa-Arsa, 2016).

This section describes some known pulsars moving behind the Sun. From the
results in Subsection 2.3, a pulsar with a higher frequency and a smaller
impact parameter to the Sun is more suitable for probing the solar structure.
We searched the ATNF pulsar
catalog666http://www.atnf.csiro.au/research/pulsar/psrcat/ for pulsars
satisfying the following criteria: i) $f\\!>\\!10\,$Hz and ii) closest
distance to the solar center is within $3\,\theta_{\odot}$. Twelve samples met
the criteria: four pulsars crossing behind the Sun and eight pulsars passing
near but not behind the solar surface. The information of these pulsars is
summarized in Table 1. The GW frequency $f$ is twice the spin frequency for
the GW emission from the mass-quadrupole moment. These frequencies are covered
by the sensitivity range of the current ground-based detectors. The LVK
collaboration obtained $95\%$ credible upper limits on the strain amplitudes
of some of these samples from the O2 and O3 runs (Abbott et al., 2022a). They
also presented the corresponding upper bounds on the ellipticity of the mass
distribution (estimated from Eq. (15)). The collaboration also performed an
all-sky search for continuous GWs coming from any (known and unknown) spinning
neutron stars in isolated and binary systems (Abbott et al., 2021; The LIGO
Scientific Collaboration et al., 2022). Their upper bound was
$h_{0}^{95\%}\sim 1.1\,(2.0)\times 10^{-25}$ at $f=100$–$200$
$(50$–$700)\,{\rm Hz}$ for isolated systems (and similarly for binary
systems).

Figure 6.— Trajectories of pulsars moving behind the Sun (represented by an
orange circle of radius $\theta_{\odot}$). $\theta_{\rm{\scriptscriptstyle
LAT}}$ is the ecliptic latitude and $\theta_{\rm{\scriptscriptstyle LON}}$ is
the ecliptic longitude with origin set at the solar center.

Figure 6 shows the trajectories of the upper four pulsars in Table 1. These
pulsars move along the horizontal axis with a velocity of
$\dot{\theta}_{\rm{\scriptscriptstyle LON}}\simeq-2\pi/(1{\rm yr})=-2.5\,{\rm
arcmin}/{\rm hour}=-0.15\,\theta_{\odot}/{\rm hour}$. The typical crossing
time is $2\theta_{\odot}/|\dot{\theta}_{\rm LON}|\simeq 13\,{\rm hours}$.
Marchant et al. (2020) included J1022+1001 and J1730-2304 as candidates, but
excluded the very recently discovered pulsar J1745-23 (Cameron et al., 2020).

Figure 7.— Amplitude and phase modulations of the four pulsars moving behind
the Sun. The lower $x$ axis is the pulsar position
$\theta_{\rm{\scriptscriptstyle LON}}/\theta_{\odot}$ along the ecliptic
longitude (Fig. 6), while the upper $x$ axis represents the corresponding
crossing time. Dotted curves in the left panel represent negative values.

Figure 7 plots $A_{F}$ and $\varPhi_{F}$ as functions of
$\theta_{\rm{\scriptscriptstyle LON}}$ for the four pulsars. J1745-23, with
the highest frequency and closest impact parameter, yields the maximum $A_{F}$
($\sim 2\,\%$). The highest $|\varPhi_{F}|$ is $\sim 0.1\,$rad for all pulsars
except J1809-2332. The anomalous result for J1809-2332 is attributable to the
low frequency ($14\,$Hz) of this pulsar. The typical $|\varPhi_{F}|$ is
estimated as $|\varPhi_{F}|\simeq 2\pi f|t_{\rm d}|$ $\simeq 0.1\,{\rm
rad}\,(f/200\,{\rm Hz})\,(|t_{\rm d}|/100\,\mu{\rm s})$. The right panel
clearly shows the potential well around $\theta_{\rm{\scriptscriptstyle
LON}}=0$. The effect is especially noticeable in the curves of J1745-23 and
J1022+1001. As evidenced in the figure, $|A_{F}|\ll|\varPhi_{F}|$.

Figure 8.— Phase modulations of the eight pulsars passing close to (but not
behind) the Sun.

Figure 8 plots $\varPhi_{F}$ as a function of $\theta_{\rm{\scriptscriptstyle
LON}}$ for the lower eight pulsars in Table 1. Because these pulsars (except
for J1756-2251) have relatively high frequencies, their $|\varPhi_{F}|$ are
relatively large ($>0.1$ rad). As these pulsars do not pass behind the Sun,
their $|A_{F}|$ are negligibly small.

##  4\. Parameter extraction from the lensed waveform

This section evaluates the accuracy to which the solar density profile can be
extracted from the lensed GW signal through a Fisher analysis.

###  4.1. Fisher matrix analysis

We consider a highly spinning neutron star with a small non-axisymmetry around
the spin axis as the GW source. Neglecting spin down, the source is assumed to
emit continuous monochromatic GWs during the observational period. A waveform
of frequency $f$ is then described as

$h(t;f)=h_{0}\cos\left(2\pi ft+\phi_{0}\right),$ (14)

where $\phi_{0}$ is a constant phase. The strain amplitude $h_{0}$ is
determined by the non-axisymmetry as (e.g., Riles, 2022)

$h_{0}=\frac{4\pi^{2}\epsilon GIf^{2}}{c^{4}D_{\rm s}},$ (15)

where $I$ is the moment of inertia about the spin axis and $\epsilon$ is the
ellipticity of the moment of inertia. A gravitationally lensed waveform
$h^{\rm L}(t;f)$ can be obtained from Eqs. (1) and (14) as777Here, we perform
a Fourier transform to obtain the unlensed waveform in the frequency domain,
then calculate the lensed waveform using Eq. (1), and finally repeat the
Fourier transform to obtain the lensed waveform (16) in the time domain.

$h^{\rm L}(t;f)=h_{0}\left[1+A_{F}(t;f)\right]\cos\left[2\pi
ft+\phi_{0}+\varPhi_{F}(t;f)\right].$ (16)

Here, $A_{F}$ and $\varPhi_{F}$ are numerically obtained from Eqs. (1) and
(6). When $\theta_{\rm s}$ is large ($\theta_{\rm s}>20\,\theta_{\odot}$), we
can apply the geometrical-optics results (Eqs. (9) and (10)), which are fully
valid at such large angles (Fig. 4).

Suppose that the pulsar is observed from time $t=-T/2$ to $T/2$, where $T$ is
the observational period and $t=0$ is the time at which the pulsar is closest
to the solar center. The signal-to-noise ratio (SNR) is calculated as

$\displaystyle{\rm SNR}^{2}$
$\displaystyle=\frac{1}{S_{n}(f)}\int_{-T/2}^{T/2}\\!dt\left[h^{\rm
L}(t;f)\right]^{2},$ $\displaystyle\simeq\frac{h_{0}^{2}}{2S_{n}(f)}\,T,$ (17)

where $S_{n}$ is the noise spectrum of the detector. Throughout this paper,
the measurement accuracies of the fitting parameters are normalized by the
SNR. Therefore, the results do not depend on a specific form of $S_{n}$. In
the second equality of Eq. (17), we neglect $A_{F}$ and use the approximation
$\cos^{2}(2\pi ft+\phi_{0}+\varPhi_{F})\simeq 1/2$ obtained by long-time
averaging. We comment on the validity of the approximation in Eq. (17). The
left panel of Fig. 7 shows $A_{F}\approx\mathcal{O}(0.01)$ for a few hours
just behind the Sun; therefore, its average over one year is $\langle
A_{F}\rangle_{\rm 1yr}\approx\mathcal{O}(0.01)\times({\rm
a~{}few~{}hours})/(1{\rm yr})\approx\mathcal{O}(10^{-6})$. Similarly, in the
phase of Eq. (16), the first term, $2\pi ft\approx 6\times 10^{10}(f/{\rm
300\,Hz})(t/1\,{\rm yr})$, is ten orders of magnitude larger than the third
term, $\varPhi_{F}\approx\mathcal{O}(1)$ rad from Figs. 7 and 8. This suggests
that $A_{F}$ and $\varPhi_{F}$ give negligible contributions to the SNR (but
note that $A_{F}$ does increase the SNR during the crossing). This also
suggests that $h_{L}(t;f)$ is rapidly oscillating during the observation and
thus the approximation $\cos^{2}(2\pi ft+\phi_{0}+\varPhi_{F})\simeq 1/2$ is
valid with a high precision.

The lensed waveform (16) depends on the solar density profile, the constant
amplitude $h_{0}$, and the constant phase $\phi_{0}$. Let the solar density
profile be characterized by $N$ parameters $A_{\rm i}$ (${\rm
i}=1,2,\cdots,N$). The total set of fitting parameters is $p_{\mu}=(A_{\rm
i},\ln h_{0},\phi_{0})$ with $\mu=1,2,\cdots,N+2$. Let us consider a GW signal
as a sum of the theoretical template (16) characterized by the fiducial
parameters $\bar{p}_{\mu}$ and a Gaussian noise. Here, our fiducial solar
model is the BS05(OP) in Bahcall et al. (2005). In a likelihood analysis, one
can obtain the best-fit parameters by fitting the theoretical template to the
signal. In the limit of high SNR, the posterior distribution of ${\bm{p}}$ is
multivariate Gaussian:
$P({\bm{p}})\propto\exp[-\sum_{\mu,\nu}(p_{\mu}-\bar{p}_{\mu})\Gamma_{\mu\nu}(p_{\nu}-\bar{p}_{\nu})/2\,]$
(for a detailed discussion, see Cutler & Flanagan, 1994). The measurement
accuracy of $p_{\mu}$ is given by the inverse of the Fisher matrix: $\Delta
p_{\mu}=[(\Gamma^{-1})_{\mu\mu}]^{1/2}$. Neglecting the degeneracy with the
other parameters, $\Delta p_{\mu}=(\Gamma_{\mu\mu})^{-1/2}$. The Fisher matrix
is given by (e.g., Cutler, 1998)

$\Gamma_{\mu\nu}=\frac{1}{S_{n}(f)}\int_{-T/2}^{T/2}\\!dt\,\frac{\partial
h^{\rm L}(t;f)}{\partial p_{\mu}}\frac{\partial h^{\rm L}(t;f)}{\partial
p_{\nu}}.$ (18)

Using Eqs. (16) and (17), all components of $\Gamma$ are given as follows:

$\displaystyle\Gamma_{A_{\rm i}A_{\rm j}}$
$\displaystyle=\int_{t}\,\left[\frac{\partial A_{F}(t;f)}{\partial A_{\rm
i}}\frac{\partial A_{F}(t;f)}{\partial A_{\rm j}}\right.$
$\displaystyle~{}~{}+\left\\{1+A_{F}(t;f)\right\\}^{2}\left.\frac{\partial\varPhi_{F}(t;f)}{\partial
A_{\rm i}}\frac{\partial\varPhi_{F}(t;f)}{\partial A_{\rm j}}\right],$
$\displaystyle\Gamma_{A_{\rm i}\ln\\!h_{0}}$
$\displaystyle=\int_{t}\,\left\\{1+A_{F}(t;f)\right\\}\frac{\partial
A_{F}(t;f)}{\partial A_{\rm i}},$ $\displaystyle\Gamma_{A_{\rm i}\phi_{0}}$
$\displaystyle=\int_{t}\,\left\\{1+A_{F}(t;f)\right\\}^{2}\frac{\partial\varPhi_{F}(t;f)}{\partial
A_{\rm i}},$ $\displaystyle\Gamma_{\ln\\!h_{0}\ln\\!h_{0}}$
$\displaystyle=\Gamma_{\phi_{0}\phi_{0}}=\int_{t}\,\left\\{1+A_{F}(t;f)\right\\}^{2},$
$\displaystyle\Gamma_{\ln\\!h_{0}\phi_{0}}$ $\displaystyle=0,$ (19)

with

$\int_{t}\,\equiv{\rm SNR}^{2}\,\frac{1}{T}\int_{-T/2}^{T/2}\\!dt.$

Similar to Eq. (17), we here applied the approximations $\sin^{2}{\rm
or}\cos^{2}(2\pi ft+\phi_{0}+\varPhi_{F})\simeq 1/2$ and $\sin(2\pi
ft+\phi_{0}+\varPhi_{F})\,\cos(2\pi ft+\phi_{0}+\varPhi_{F})\simeq 0$. From
Eq. (19), $\Delta p_{\mu}$ simply scales as $\Delta p_{\mu}\propto{\rm
SNR}^{-1}$.

The simple waveform (14) neglects the detector response which depends on other
source parameters such as the inclination of the spin axis and GW
polarization. Owing to Earth’s spin and orbital motion, the detector response
causes daily and yearly periodic modulations on both amplitude and phase.
Unlike the solar modulation (which is temporary but repeats annually), the
response modulations are trigonometric functions of time. We expect that the
solar modulation can be distinguishable from the response modulations and thus
the response does not significantly affect the measurement accuracies of the
solar model parameters.

###  4.2. Constraint on the potential amplitude

Figure 9.— Measurement accuracies of the solar-potential amplitude $A_{\psi}$
defined in Eq. (20), amplitude $h_{0}$, and phase $\phi_{0}$ of the pulsars
crossing behind the Sun. Dashed red curve plots $\Delta A_{\psi}$ neglecting
the parameter degeneracy. The horizontal axis $T$ denotes the observational
period. The signal-to-noise ratio (SNR) is set to $100\,(T/{\rm yr})^{1/2}$
and the results simply scale as $\Delta p\propto{\rm SNR}^{-1}$. Thick gray
line denotes ${\rm SNR}^{-1}$.

This subsection presents the measurement accuracy of the overall amplitude of
the solar gravitational potential. Replacing the potential as

$\hat{\psi}(\theta)\rightarrow A_{\psi}\hat{\psi}(\theta),$ (20)

we calculate the measurement error in $A_{\psi}$. The derivatives in the
Fisher matrix with respect to $A_{\psi}$ are numerically obtained by changing
$A_{\psi}$ by $\pm 1\,\%$.

Figure 9 plots the measurement accuracies of $A_{\psi}$, $\ln h_{0}$ and
$\phi_{0}$ for the pulsars moving behind the Sun at $f>100\,{\rm Hz}$. The
results are normalized with ${\rm SNR}=100$ over a one-year observation
period. Because $|A_{F}|\ll|\Phi_{F}|$, $A_{\psi}$ is mainly determined by the
phase modulation. The slope of $\Delta A_{\psi}$ changes slightly at the
surface (corresponding to $T\approx 10^{-3}\,{\rm yr}$); the slope is steeper
(shallower) inside (outside) the Sun. These results can be explained as
follows. In the Sun’s interior, $|\varPhi_{F}|$ increases with $\theta$ (in
contrast to the constant $\phi_{0}$), thus $A_{\psi}$ and $\phi_{0}$ can be
determined almost independently. Outside the Sun, $\varPhi_{F}$ is nearly
constant at large $\theta\,(\gg\theta_{\odot})$, thus $A_{\psi}$ and
$\phi_{0}$ are degenerate to some extent. The dashed red curve plots $\Delta
A_{\psi}$ without the parameter degeneracy (i.e., $\Delta
A_{\psi}=(\Gamma_{A_{\psi}A_{\psi}})^{-1/2}$), which is $\sim 10$ times better
than the solid red curve. Previously, Marchant et al. (2020) roughly estimated
the detectability of lensing signatures. As they ignored the parameter
degeneracy, they underestimated the measurement error. Inside the Sun, $h_{0}$
cannot be determined, because $h_{0}$ and $A_{F}$ are highly degenerate, but
outside the Sun, the degeneracy is broken and $\Delta\ln h_{0}\simeq{\rm
SNR}^{-1}$.

We caution that $\Delta\phi_{0}$ shown in Fig. 9 depends on an arbitrary
constant in the gravitational potential. If a constant term is added to the
potential, i.e.,
$\hat{\psi}(\theta)\rightarrow\hat{\psi}(\theta)+\hat{\psi}_{0}$, the phase
modulation changes as $\varPhi_{F}\rightarrow\varPhi_{F}-2\pi
f\hat{\psi}_{0}/c^{3}$ using Eqs. (1)–(3). Two components of the Fisher matrix
change
accordingly888$\Gamma_{A_{\psi}\phi_{0}}\rightarrow\Gamma_{A_{\psi}\phi_{0}}-(2\pi
f\hat{\psi}_{0}/c^{3})\Gamma_{\phi_{0}\phi_{0}}$ and
$\Gamma_{A_{\psi}A_{\psi}}\rightarrow\Gamma_{A_{\psi}A_{\psi}}-(4\pi
f\hat{\psi_{0}}/c^{3})\Gamma_{A_{\psi}\phi_{0}}+(2\pi
f\hat{\psi}_{0}/c^{3})^{2}\Gamma_{\phi_{0}\phi_{0}}$.. After some algebra, one
finds that $\Delta A_{\psi}$, $\Delta\ln h_{0}$, and
$(\Gamma^{-1})_{A_{\psi}\ln h_{0}}$ are independent of $\hat{\psi}_{0}$ but
$\Delta\phi_{0}$ and its cross correlations,
$(\Gamma^{-1})_{A_{\psi}\phi_{0}}$ and $(\Gamma^{-1})_{\ln h_{0}\phi_{0}}$, do
depend on $\hat{\psi}_{0}$.

Table 2Measurement accuracies of the potential amplitude $A_{\psi}$ (the
fiducial value $=1$) over one year of observation with ${\rm SNR}=100$.

name | $\Delta A_{\psi}$ |  $f\Delta A_{\psi}$ (Hz)
---|---|---
J1745$-$23 |  $0.24$ ($0.026$) | $90$
J1022+1001 |  $0.74$ ($0.077$) | $90$
J1809$-$2332 |  $6.6$ ($0.69$) | $90$
J1730$-$2304 |  $0.37$ ($0.038$) | $91$
J1858-2216 |  $0.11$ ($0.011$) | $92$
J1142+0119 |  $0.23$ ($0.024$) | $92$
J2310-0555 |  $0.12$ ($0.012$) | $92$
J1756-2251 |  $1.3$ ($0.13$) | $92$
J1646-2142 |  $0.27$ ($0.028$) | $93$
J1811-2405 |  $0.12$ ($0.013$) | $93$
J1836-2354B |  $0.16$ ($0.016$) | $93$
J1836-2354A |  $0.15$ ($0.015$) | $93$

Note. — Values in parentheses are calculated without the parameter degeneracy.

Table 2 lists the $\Delta A_{\psi}$ results of all pulsars. Although $\Delta
A_{\psi}$ varies among the pulsars, $f\Delta A_{\psi}$ is almost constant;
i.e., it is independent of the impact parameter to the Sun. $\Delta A_{\psi}$
can then be approximately fitted as

$\Delta A_{\psi}\approx 0.3\left(\frac{f}{300\,{\rm
Hz}}\right)^{-1}\left(\frac{\rm SNR}{100}\right)^{-1},$ (21)

indicating that a higher-frequency pulsar is more promising for detecting
$A_{\psi}$. According to Eq. (21), $A_{\psi}$ can be detected at the $3\sigma$
confidence level when ${\rm SNR}\approx 100\,(f/300{\rm Hz})^{-1}$.

###  4.3. Measurement of the solar density profile

Figure 10.— Division of the solar surface into $N$ annuli. Each annulus has
the same mass in Case (I) and the same radial width in Case (II). Figure 11.—
Measurement accuracies of the solar-density amplitudes in the i-th annulus
$A_{\Sigma{\rm i}}$ for the pulsars crossing behind the Sun. Upper and lower
panels plot the results in Case (I) and (II), respectively. In the upper
panels, $\Delta A_{\Sigma 1}=\Delta A_{\Sigma 2}$. The horizontal axis $T$
denotes the observational period. SNR is set to $10^{4}\,(T/{\rm yr})^{1/2}$
and the results simply scale as $\Delta A_{\Sigma{\rm i}}\propto{\rm
SNR}^{-1}$. Figure 12.— Same as Fig. 11, but for three annuli with ${\rm
SNR}=10^{5}\,(T/{\rm yr})^{1/2}$. Table 3Measurement accuracies of the density
amplitudes in two annuli over a one-year observation period with ${\rm
SNR}=10^{4}$.

| Case (I) | Case (II)
---|---|---
name | $\Delta A_{\Sigma 1,2}$ | $\Delta A_{\Sigma 1},\,\Delta A_{\Sigma 2}$
J1745$-$23 | $0.40$ | $0.11,~{}1.8$
J1022+1001 | $2.0$ | $0.47,~{}8.0$
J1730$-$2304 | $2.9$ | $1.3,~{}22$

Note. — In Case (I), $\Delta A_{\Sigma 1}=\Delta A_{\Sigma 2}$.

Table 4Same as Table 3, but for three annuli and a one-year observation period with ${\rm SNR}=10^{5}$. | Case (I) | Case (II)
---|---|---
name | $\Delta A_{\Sigma 1},\,\Delta A_{\Sigma 2},\,\Delta A_{\Sigma 3}$ | $\Delta A_{\Sigma 1},\,\Delta A_{\Sigma 2},\,\Delta A_{\Sigma 3}$
J1745$-$23 | $0.24,~{}0.33,~{}0.10$ | $0.029,~{}0.24,~{}1.5$
J1022+1001 | $1.3,~{}1.8,~{}0.53$ | $0.17,~{}1.3,~{}6.9$
J1730$-$2304 | $0.43,~{}0.92,~{}0.74$ | $0.22,~{}0.93,~{}6.3$

This subsection presents our main result, namely the measurement accuracy of
the solar density profile. Let the surface density be divided into $N$ annuli
(Fig. 10). We write the surface mass density of the i-th annulus as
$\Sigma_{\rm i}(\theta)$ and the mass as $M_{\rm i}$ (${\rm i}=1,\cdots,N$),
where a smaller i corresponds to an inner annulus. We consider the following
two cases:

Case (I)

Each annulus has the same mass ($=M_{\odot}/N$).

Case (II)

Each annulus has the same radial thickness ($=\theta_{\odot}/N$).

We first consider $N=2$ and $3$. In Case (I), the boundary radius between the
annuli is $0.20\,\theta_{\odot}$ for $N=2$. The boundary radii are
$0.14\,\theta_{\odot}$ and $0.26\,\theta_{\odot}$ between annulus 1 and 2 and
between annulus 2 and 3, respectively, for $N=3$. In Case (II), the mass
ratios are $M_{1}:M_{2}\simeq 17:1$ for $N=2$ and $M_{1}:M_{2}:M_{3}\simeq
57:13:1$ for $N=3$ (i.e., the mass reduces from the innermost to the outermost
annulus).

We change the amplitude of the i-th density as

$\Sigma_{\rm i}(\theta)\rightarrow A_{\Sigma{\rm i}}\Sigma_{\rm i}(\theta),$
(22)

while fixing the total mass, i.e., $\sum_{{\rm i}=1}^{N}A_{\Sigma{\rm
i}}M_{\rm i}=M_{\odot}$. Under these settings, we calculate the measurement
accuracies of $A_{\Sigma{\rm i}}$, $\ln h_{0}$, and $\phi_{0}$. When $N=2$,
the density profile is solely characterized by the fitting parameter
$A_{\Sigma 1}$ because the total mass is fixed, i.e., $M_{1}\Delta A_{\Sigma
1}=M_{2}\Delta A_{\Sigma 2}$. Similarly, when $N=3$, the fitting parameters
are $A_{\Sigma 1}$ and $A_{\Sigma 2}$. The other parameter, $A_{\Sigma 3}$, is
determined via $(M_{3}\Delta A_{\Sigma 3})^{2}=(M_{1}\Delta A_{\Sigma
1})^{2}+(M_{2}\Delta A_{\Sigma 2})^{2}+2M_{1}M_{2}(\Gamma^{-1})_{A_{\Sigma
1}A_{\Sigma 2}}$. To obtain the derivative with respect to $A_{\Sigma{\rm i}}$
in the Fisher matrix, we numerically calculate the derivatives of the
potential in Eq. (4) and the amplification factor in Eq. (6) by changing
$A_{\Sigma{\rm i}}$ by $\pm 1\,\%$.

Figure 11 plots the measurement accuracies of $A_{\Sigma{\rm i}}$ for $N=2$.
As the total mass is fixed, the lensed signal is insensitive to the density
profile when the pulsar is outside the Sun. Therefore, the accuracies are not
improved at $T\gtrsim 10\,{\rm hours}$. J1745-23 gives the best accuracy,
because it has the highest frequency and smallest impact parameter. J1730-2304
has a higher frequency than J1022+1001, but its larger impact parameter
lessens the constraint. In Case (II), $\Delta A_{\Sigma 2}$ is $17$ times
larger than $\Delta A_{\Sigma 1}$, reflecting the mass ratio
$M_{1}/M_{2}\simeq 17$. According to the figure, at least ${\rm SNR}\approx
10^{4}\,(T/{\rm yr})^{1/2}$ is required for probing the solar density profile,
but a lower SNR may be sufficient for probing with J1745-23. We comment that
$\Delta\ln h_{0}$ and $\Delta\phi_{0}$ agree with ${\rm SNR}^{-1}$ outside the
Sun (corresponding to $T>10$ hours; these results are not plotted in the
figure). The numerical values of $\Delta A_{\Sigma{\rm i}}$ are listed in
Table 3.

Figure 12 plots the measurement accuracies of $A_{\Sigma{\rm i}}$ for $N=3$.
In Case (I), $\Delta A_{\Sigma 3}$ shows the highest accuracy for most pulsars
(the exception is J1730-2304), because the third annulus (with radius ranging
from $0.26\,\theta_{\odot}$ to $\theta_{\odot}$) has the largest area.
Therefore, the pulsar remains longer in the third annulus than in the other
annuli. In Case (II), $\Delta A_{\Sigma 1}<\Delta A_{\Sigma 2}<\Delta
A_{\Sigma 3}$ simply because the inner annulus is more massive than the outer
one. The numerical values of $\Delta A_{\Sigma{\rm i}}$ are listed in Table 4.
According to Fig. 12, at least ${\rm SNR}\approx 10^{5}\,(T/{\rm yr})^{1/2}$
is required for measuring the density profile with $N=3$. The SNR must be
further increased for larger $N$ ($>3$).

In geometrical optics, $\varPhi_{F}$ is fully determined by the mass enclosed
within a given radius. Therefore, a pulsar can probe the enclosed mass down to
the impact parameter. In wave optics, $\varPhi_{F}$ also depends on the outer
density profile, owing to diffraction (Subsection 2.3). Therefore, the
accuracy of $A_{\Sigma{\rm i}}$ depends on the frequency and impact parameter
in a complicated manner. Combining the results of different pulsars improves
the measurement accuracy of the solar density and increases the radial range
of the estimation.

Figures 11 and 12 suggest that the accuracies strongly depend on the chosen
boundary. In Case (I), the accuracies are similar in different annuli, but in
Case (II), the accuracy is better in the inner annulus. This implies that the
annulus mass mainly determines the accuracy. The boundary that maximizes the
accuracy will depend on the frequency and impact parameter of the pulsar.

Figure 13.— Measurement accuracies of the projected solar density as functions
of radius. Solid black curve is the fiducial solar model and error bars are
the accuracies for each pulsar. In the left and right panels, the density is
divided into six and ten annuli in Case (I) with ${\rm SNR}=2\times
10^{6}\,(T/{\rm yr})^{1/2}$ and $10^{7}\,(T/{\rm yr})^{1/2}$, respectively.
Vertical dashed lines are the boundaries between the annuli. Down arrows
indicate the ecliptic latitude of each pulsar.

The solar density can (in principle) be plotted as a function of radius in the
limit of large $N$ (but a very high SNR is needed). As a demonstration, Fig.
13 plots the measurement accuracies for $N=6$ and $10$. Here, we follow the
same procedure as done for $N=2$ and $3$. As expected, J1745-23 gives the best
accuracy. Though J1730-2304 has a large impact parameter (denoted by the down
blue arrow), it can constrain the inner density due to the diffraction. The
density of the outermost annulus is determined most precisely, because it has
the largest area and the pulsars remain there for the longest time. The
density of the intermediate annuli ($0.1\lesssim\theta/\theta_{\odot}\lesssim
0.3$) is poorly determined, because the errors of nearby annuli are strongly
correlated. For instance, the cross-correlation coefficient of the accuracies
in between the i- and j-th annuli, $r_{\rm
ij}\equiv(\Gamma^{-1})_{A_{\Sigma{\rm i}}A_{\Sigma{\rm
j}}}/[(\Gamma^{-1})_{A_{\Sigma{\rm i}}A_{\Sigma{\rm
i}}}(\Gamma^{-1})_{A_{\Sigma{\rm j}}A_{\Sigma{\rm j}}}]^{1/2}$, is large for a
close pair of the annuli, $0.6<|r_{\rm ij}|<1$ for $|i-j|=1$.

The strain sensitivities of the planned detectors of the Cosmic Explorer
(CE)999https://cosmicexplorer.org/ and the Einstein Telescope
(ET)101010https://www.et-gw.eu/ will be approximately $100$ times better than
those of the LIGO O3 run. The current upper limits in Table 1 were determined
from observations taken over approximately one year in the O3 run. These
limits correspond to an upper limit of ${\rm SNR}\lesssim 10$ (Appendix C of
Abbott et al., 2019). Therefore, in CE and ET observations, the best SNR for
these known pulsars is estimated as $\approx 10^{3}\,(T/{\rm yr})^{1/2}$ (a
similar discussion is given in Section V of Marchant et al., 2020). Judging
from these estimates, measuring the density profile using these pulsars will
be a difficult task in the near future. To alleviate this problem, we require
more pulsars than those detected by radio telescopes to date (as discussed in
the next section). Moreover, these pulsars must have high frequencies and
small impact parameters.

##  5\. Millisecond pulsars

Table 5Globular clusters crossing behind the Sun

name | $\theta_{\rm{\scriptscriptstyle LAT}}/\theta_{\odot}$ | $\theta_{\rm{\scriptscriptstyle h}}/\theta_{\odot}$ |  $D_{\rm s}\,$(kpc) | mass ($M_{\odot}$)
---|---|---|---|---
NGC 6287 | $0.49$ | $0.046$ | $7.9$ | $1.5\times 10^{5}$
NGC 6717 | $0.50$ | $0.043$ | $7.5$ | $3.6\times 10^{4}$
NGC 6642 | $-0.89$ | $0.046$ | $8.1$ | $3.4\times 10^{4}$

Note. — Listed are the ecliptic latitude of the cluster center
$\theta_{\rm{\scriptscriptstyle LAT}}$, the half-light angular radius
$\theta_{\rm{\scriptscriptstyle h}}$, the distance $D_{\rm s}$, and the mass.
Here, the distance and mass are taken from Baumgardt’s catalog of globular
clusters111111https://people.smp.uq.edu.au/HolgerBaumgardt/globular/ (e.g.,
Baumgardt & Hilker, 2018).

This section roughly estimates the number of MSPs moving behind the Sun that
will be found in the near future.

The present ATNF catalog lists approximately $3300$ pulsars, including around
$500$ highly spinning pulsars with a spin period of $<10\,{\rm ms}$
(corresponding to $f>200\,{\rm Hz}$). These pulsars are the so-called MSPs.
The Milky Way is estimated to host $40000$–$90000$ MSPs (e.g., Bhattacharyya &
Roy, 2022). As the sky fraction of the zodiac belt swept by the Sun is
$2\pi\times 2\theta_{\odot}/(4\pi)\simeq 0.5\%$, the belt probably holds
$200$–$400$ MSPs. Almost none of these MSPs have been found, but tens of them
will be identified in ongoing and upcoming radio surveys. For instance, the
Canadian Hydrogen Intensity Mapping Experiment (CHIME)121212https://chime-
experiment.ca/en, MeerKAT131313https://www.sarao.ac.za/gallery/meerkat/, and
the Five-hundred-meter Aperture Spherical radio Telescope
(FAST)141414https://fast.bao.ac.cn will detect $\approx 1000$ MSPs over the
sky (Lorimer et al., 2019, their Table 1). These performances will be
surpassed by the Square Kilometer Array (SKA) Phase 1, which will detect
$1500$ MSPs (Keane et al., 2015).

At present, approximately $220$ MSPs have been found in $33$ globular
clusters151515https://www3.mpifr-bonn.mpg.de/staff/pfreire/GCpsr.html; $\sim$
seven MSPs per cluster on average. Per unit mass, pulsars are two or three
orders of magnitude more abundant in clusters than in the Galactic disk
(Freire, 2013). Therefore, many MSPs probably exist in clusters. Among the
$157$ clusters listed in the McMaster catalog of Milky Way globular
clusters161616Data may be downloaded from
https://physics.mcmaster.ca/Fac_Harris/mwgc.dat (version 2010, Harris, 1996,
2010), three globular clusters move behind the Sun. The data are summarized in
Table 5. Although no MSPs have been discovered in these
clusters171717https://www3.mpifr-bonn.mpg.de/staff/pfreire/GCpsr.html, tens of
MSPs are expected in each cluster.

##  6\. Lensing by a Galactic star

This section briefly discusses GW lensing by a Milky Way star along the line
of sight to the pulsar. As a star crosses in front of the pulsar, lensing
imprints a time-dependent modulation on the waveform. This event might be
confused with the solar lensing. We comment that these two events can be
distinguished due to the following reasons. First, the stellar lensing is
observed only once, but the solar lensing is observed every year on the same
day. In addition, the probability of the stellar lensing is not high, $\approx
10^{-6}$ ${\rm yr}^{-1}$, on average over the full sky (e.g., Han, 2008). If
the pulsar is in a globular cluster, the probability of lensing by a cluster
star (so-called self-lensing) can reach $\approx 0.1$ ${\rm yr}^{-1}$ (Kıroğlu
et al., 2022). Second, the lensing modulation by the star (modeled as a point
mass) differs from the solar modulation (e.g., Figs. 3 and 5 of Marchant et
al., 2020). For example, the typical event duration is $\approx 10$ days
(e.g., Schneider et al., 2006), which is much longer than the solar lensing
duration $\approx$ half a day. Third, if radio observations are available for
the stellar lensing, double signals (with different amplitudes and arrival
times) from the single pulsar may be detected.

##  7\. Conclusions

We have studied the detectability of solar density profiling using GW lensing
with known pulsars. After selecting suitable pulsars with high frequencies and
small impact parameters to the Sun from the ATNF catalog, we calculated the
measurement accuracy of the overall amplitude of the solar gravitational
potential using a Fisher analysis. The lensing signature can be detected with
$3\sigma$ confidence when ${\rm SNR}\approx 100\,(f/300\,{\rm Hz})^{-1}$
during one year of observation (Table 2 and Eq. (21)). The detection is
therefore improved with high-frequency pulsars. The signature can be detected
even if the pulsar trajectory does not pass behind the Sun (in such cases, the
impact parameter to the Sun does not significantly affect the detectability).
We found that the detectability is degraded by parameter degeneracy with the
constant phase ($\phi_{0}$ in Eq. (16)) of the waveform. Next, we divided the
projected density profile into $N$ annuli ($N=2,3,\cdots$) and calculated the
measurement accuracy of each annulus mass. If three known pulsars move behind
the Sun with $f>100\,{\rm Hz}$, a high SNR ($\gtrsim 10^{4}$) in a year-long
observation is required for measuring the density profile in the two-annulus
case (Fig. 11). The SNR must be raised to $\gtrsim 10^{5}$ in the three-
annulus case (Fig. 12) and $\gtrsim 10^{6}$ in the six-annulus case (Fig. 13).

Currently, only two known MSPs move behind the Sun, but $200$–$400$ can
potentially exist within our Galaxy (Section 5). Ongoing and future radio
surveys such as SKA, CHIME, and FAST will find tens of these MSPs in the
coming decades.

The solar density profile has already been measured to percent-level precision
based on global helioseismology (see the recent review by Basu, 2016). GW
lensing will likely become a complementary and independent method for probing
the Sun in the future.

## Acknowledgments

We thank Hideyuki Hotta for his useful comment on helioseismology. We thank
the CASS (CSIRO (Commonwealth Scientific and Industrial Research Organisation)
Astronomy and Space Science) pulsar group for developing and maintaining the
ATNF pulsar catalog. This work is supported in part by MEXT Japan and JSPS
KAKENHI Grant Numbers of JP22H00130 (RT), JP20H05855 (RT), JP17H06359 (TS),
JP21H05453 (TS), and JP19K03864 (TS).

##  Appendix A Validity of the thin-lens approximation

Suyama et al. (2005) investigated the validity of the thin-lens approximation
in wave optics. They considered plane waves entering a lens and numerically
solved the wave equation to obtain the exact lensed waveform (without the
thin-lens approximation). They considered three spherical lens models (a
uniform density sphere, a singular isothermal sphere truncated at a finite
radius, and the Hernquist profile). They found that the relative error of $F$
is smaller than the lens radius over $D_{\rm L}$ and the maximum error is
achieved when the wavelength is comparable to the Schwarzschild radius of the
lens. This suggests that the error of $F$ is $\lesssim
R_{\odot}/(1\,{\rm{au}})\simeq 0.5\,\%$ in our case.

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2017) —. 2017, ApJ, 839, 12
  * Abbott et al. (2019) —. 2019, ApJ, 879, 10
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, Phys. Rev. D, 103, 064017
  * Abbott et al. (2022a) Abbott, R., Abe, H., Acernese, F., et al. 2022a, ApJ, 935, 1
  * Abbott et al. (2022b) —. 2022b, Phys. Rev. D, 106, 042003
  * Abdo et al. (2009) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Science, 325, 840
  * Backer & Hellings (1986) Backer, D. C., & Hellings, R. W. 1986, ARA&A, 24, 537
  * Bahcall et al. (2005) Bahcall, J. N., Serenelli, A. M., & Basu, S. 2005, ApJ, 621, L85
  * Bartelmann & Schneider (2001) Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291
  * Basu (2016) Basu, S. 2016, Living Reviews in Solar Physics, 13, 2
  * Baumgardt & Hilker (2018) Baumgardt, H., & Hilker, M. 2018, MNRAS, 478, 1520
  * Bhattacharyya & Roy (2022) Bhattacharyya, B., & Roy, J. 2022, Astrophysics and Space Science Library, 465, 1
  * Bliokh & Minakov (1975) Bliokh, P. V., & Minakov, A. A. 1975, Ap&SS, 34, L7
  * Bontz & Haugan (1981) Bontz, R. J., & Haugan, M. P. 1981, Ap&SS, 78, 199
  * Cameron et al. (2020) Cameron, A. D., Champion, D. J., Bailes, M., et al. 2020, MNRAS, 493, 1063
  * Camilo et al. (1996) Camilo, F., Nice, D. J., Shrauner, J. A., & Taylor, J. H. 1996, ApJ, 469, 819
  * Choi et al. (2021) Choi, H. G., Park, C., & Jung, S. 2021, Phys. Rev. D, 104, 063001
  * Cutler (1998) Cutler, C. 1998, Phys. Rev. D, 57, 7089
  * Cutler & Flanagan (1994) Cutler, C., & Flanagan, É. E. 1994, Phys. Rev. D, 49, 2658
  * Cyranski & Lubkin (1974) Cyranski, J. F., & Lubkin, E. 1974, Annals of Physics, 87, 205
  * Dai et al. (2018) Dai, L., Li, S.-S., Zackay, B., Mao, S., & Lu, Y. 2018, Phys. Rev. D, 98, 104029
  * Dalang et al. (2022) Dalang, C., Cusin, G., & Lagos, M. 2022, Phys. Rev. D, 105, 024005
  * Edwards et al. (2006) Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, MNRAS, 372, 1549
  * Ezquiaga et al. (2021) Ezquiaga, J. M., Holz, D. E., Hu, W., Lagos, M., & Wald, R. M. 2021, Phys. Rev. D, 103, 064047
  * Faulkner et al. (2004) Faulkner, A. J., Stairs, I. H., Kramer, M., et al. 2004, MNRAS, 355, 147
  * Ferdman et al. (2014) Ferdman, R. D., Stairs, I. H., Kramer, M., et al. 2014, MNRAS, 443, 2183
  * Freire (2013) Freire, P. C. C. 2013, Neutron Stars and Pulsars: Challenges and Opportunities after 80 years, 291, 243
  * Glampedakis & Gualtieri (2018) Glampedakis, K., & Gualtieri, L. 2018, Astrophysics and Space Science Library, 457, 673
  * Guo & Lu (2020) Guo, X., & Lu, Y. 2020, Phys. Rev. D, 102, 124076
  * Han (2008) Han, C. 2008, ApJ, 681, 806
  * Harris (1996) Harris, W. E. 1996, AJ, 112, 1487
  * Harris (2010) —. 2010, arXiv:1012.3224
  * Hobbs et al. (2006) Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655
  * Hou et al. (2019) Hou, S., Fan, X.-L., & Zhu, Z.-H. 2019, Phys. Rev. D, 100, 064028
  * Jung & Kim (2022) Jung, S., & Kim, S. 2022, arXiv:2210.02649
  * Keane et al. (2015) Keane, E., Bhattacharyya, B., Kramer, M., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 40
  * Keith et al. (2010) Keith, M. J., Jameson, A., van Straten, W., et al. 2010, MNRAS, 409, 619
  * Kıroğlu et al. (2022) Kıroğlu, F., Weatherford, N. C., Kremer, K., et al. 2022, ApJ, 928, 181
  * Liao et al. (2022) Liao, K., Biesiada, M., & Zhu, Z.-H. 2022, Chinese Physics Letters, 39, 119801
  * Lorimer et al. (2019) Lorimer, D., Pol, N., Rajwade, K., et al. 2019, BAAS, 51, 261
  * Lorimer et al. (1995) Lorimer, D. R., Nicastro, L., Lyne, A. G., et al. 1995, ApJ, 439, 933
  * Lynch et al. (2011) Lynch, R. S., Ransom, S. M., Freire, P. C. C., & Stairs, I. H. 2011, ApJ, 734, 89
  * Macquart (2004) Macquart, J. P. 2004, A&A, 422, 761
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993
  * Marchant et al. (2020) Marchant, P., Breivik, K., Berry, C. P. L., Mandel, I., & Larson, S. L. 2020, Phys. Rev. D, 101, 024039
  * Nakamura (1998) Nakamura, T. T. 1998, Phys. Rev. Lett., 80, 1138
  * Nakamura & Deguchi (1999) Nakamura, T. T., & Deguchi, S. 1999, Progress of Theoretical Physics Supplement, 133, 137
  * Ng et al. (2020) Ng, C., Guillemot, L., Freire, P. C. C., et al. 2020, MNRAS, 493, 1261
  * Nice et al. (2015) Nice, D., Demorest, P., Stairs, I., et al. 2015, Tempo: Pulsar timing data analysis, Astrophysics Source Code Library, record ascl:1509.002
  * Oguri (2019) Oguri, M. 2019, Reports on Progress in Physics, 82, 126901
  * Oguri & Takahashi (2020) Oguri, M., & Takahashi, R. 2020, ApJ, 901, 58
  * Ohanian (1974) Ohanian, H. C. 1974, International Journal of Theoretical Physics, 9, 425
  * Park et al. (2017) Park, R. S., Folkner, W. M., Konopliv, A. S., et al. 2017, AJ, 153, 121
  * Patla & Nemiroff (2008) Patla, B., & Nemiroff, R. J. 2008, ApJ, 685, 1297
  * Ray et al. (2011) Ray, P. S., Kerr, M., Parent, D., et al. 2011, ApJS, 194, 17
  * Ray et al. (2012) Ray, P. S., Abdo, A. A., Parent, D., et al. 2012, arXiv:1205.3089
  * Reardon et al. (2021) Reardon, D. J., Shannon, R. M., Cameron, A. D., et al. 2021, MNRAS, 507, 2137
  * Riles (2022) Riles, K. 2022, arXiv:2206.06447
  * Roy & Bhattacharyya (2013) Roy, J., & Bhattacharyya, B. 2013, ApJ, 765, L45
  * Sanpa-Arsa (2016) Sanpa-Arsa, S. 2016, PhD thesis, University of Virginia
  * Schneider et al. (1992) Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses (Springer-Verlag, Berlin)
  * Schneider et al. (2006) Schneider, P., Kochanek, C. S., & Wambsganss, J. 2006, Gravitational Lensing: Strong, Weak and Micro (Springer-Verlag, Berlin)
  * Sonnabend (1979) Sonnabend, D. 1979, PhD thesis, Stanford University, California
  * Stix (2004) Stix, M. 2004, The Sun : an Introduction (Springer-Verlag, Berlin)
  * Suvorov (2022) Suvorov, A. G. 2022, ApJ, 930, 13
  * Suyama et al. (2005) Suyama, T., Takahashi, R., & Michikoshi, S. 2005, Phys. Rev. D, 72, 043001
  * Takahashi (2004) Takahashi, R. 2004, A&A, 423, 787
  * Takahashi (2006) —. 2006, ApJ, 644, 80
  * Takahashi & Nakamura (2003) Takahashi, R., & Nakamura, T. 2003, ApJ, 595, 1039
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021, arXiv:2111.03606
  * The LIGO Scientific Collaboration et al. (2022) —. 2022, arXiv:2201.00697
  * Vinyoles et al. (2017) Vinyoles, N., Serenelli, A. M., Villante, F. L., et al. 2017, ApJ, 835, 202
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29