# Local distortion driven magnetic phase switching in pyrochlore
Yb2(Ti1-xSnx)2O7

Yuanpeng Zhang1, Zhiling Dun2, Yunqi Cai2, Chengkun Xing3, Qi Cui2, Naveen
Kumar Chogondahalli Muniraju4,5,6, Qiang Zhang1, Yongqing Li2, Jinguang
Cheng2,†∗, Haidong Zhou3,‡<EMAIL_ADDRESS>
<EMAIL_ADDRESS>1Neutron Scattering Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee 37831, USA
2Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
3Department of Physics and Astronomy, University of Tennessee, Knoxville, TN
37996, USA
4Forschungszentrum Jülich GmbH, Jülich Centre for Neutron Science (JCNS),
Forschungszentrum Jülich, D-52425 Jülich, Germany
5Chemical and Engineering Materials Division, Oak Ridge National Laboratory,
Oak Ridge, Tennessee 37831, USA
6Institute of Physics, Bijenička cesta 46, 10000, Zagreb, Croatia

###### Abstract

While it is commonly accepted that the disorder induced by magnetic ion doping
in quantum magnets usually generates a rugged free-energy landscape resulting
in slow or glassy spin dynamics, the disorder/distortion effects associated
with non-magnetic ion sites doping are still illusive. Here, using AC
susceptibility measurements, we show that the mixture of Sn/Ti on the non-
magnetic ion sites of pyrochlore Yb2(Ti1-xSnx)2O7 induces an antiferromagnetic
ground state despite both parent compounds, Yb2Ti2O7, and Yb2Sn2O7, order
ferromagnetically. Local structure studies through neutron total scattering
reveals the local distortion in the non-magnetic ion sites and its strong
correlation with the magnetic phase switching. Our study, for the first time,
demonstrates the local distortion as induced by the non-magnetic ion site
mixture could be a new path to achieve magnetic phase switching, which has
been traditionally obtained by external stimuli such as temperature, magnetic
field, pressure, strain, light etc.

††preprint: APS/123-QED

The complex interplay between various degree of freedom in frustrated quantum
systems enables rich physics involving complex ground states and exotic
excitations [1, 2]. The balancing between those weak interaction terms in the
Hamiltonian (long-range dipolar interactions, exchange interaction beyond the
nearest neighboring, site disorder, lattice distortions, etc.) and the
frustration could determine the coordination in the phase diagram [3, 4, 5,
6]. Exploring such a complex interplay thus becomes critical in terms of
understanding the formation of those exotic states such as the quantum spin
liquid (QSL) state where the long range magnetic ordering is completely
suppressed by frustration and the quantum fluctuation yields strong dynamics
even at zero temperature. In practice, the QSL state of matter (e.g., the non-
Abelian anyons in 2D QSL) hosts promising application in error-proof quantum
computation [7, 8, 9] and also is believed to be closely related to high
temperature superconductivity [10, 11, 12].

For the study of frustrated quantum systems, of critical concern is the
underlying lattice characteristics, such as the chemical order/disorder,
lattice/sub-lattice distortion, etc. For example, the impact of various types
of disorder has been attracting notable attention – this is especially true
for the disorder effect with respect to the magnetic species, as they are
directly (predominantly through near neighbor interaction) associated with the
magnetic Hamiltonian. Depending on the level of disorder, it can either
perturb spontaneous symmetry breaking [13] or promote magnetic ordering and
break a continuous degeneracy via an order by disorder mechanism [14].
Specifically concerning the pyrochlore systems, on which the current report is
to focus, there have been extensive studies over the magnetic species relevant
disorder effects, including those diluted [15, 16, 17, 18, 19] or stuffed [20,
21, 22, 23] pyrochlores. Following similar pathway, there have been some
recently emerging studies on the disorder effect of the non-magnetic-only
species on the non-magnetic site. Several celebrated examples recently on this
topic are YbMgGaO4 (YMGO) [24, 25, 26, 27, 28], Sr3CuTa2O9 (SCTO) [29], Sr2Cu
(Te1-xWx)O6 (SCTWO) [30, 31, 32], SrLaCuSbO6 (SLCSO) & SrLaCuNbO6 (SLCNO)
[33], and Lu3Sb3Mn2O14 (LSMO) [34]. As compared to the magnetic species
involved chemical disorder, first, the non-magnetic-only disorder has to take
effect through the magnetic interactions beyond the nearest-neighbor. For
example, in SCTWO, substituting W for Te alters the magnetic interactions from
the strong nearest-neighbor type to the strong next-nearest-neighbor type
[31], resulting in strong exchange interaction disorder that is absent in
parent compounds Sr2CuTeO6 [35] and Sr2CuWO6 [36, 37]. Theoretically, the
scenario could be described as the random-singlet (RS) state [32, 30, 38], in
which the randomness in a quantum magnet can induce spin-singlet dimers of
varying strengths with a spatially random manner and therefore account for the
spin liquid like behaviors. Second, for non-magnetic-only order/disorder
systems, it is relatively easier to separate out the lattice distortion effect
from the chemical order/disorder effect, especially for those systems in which
the magnetic interactions beyond the nearest neighbor are non-critical.

Equally important as the order/disorder effect, the local lattice distortion
also plays manifold critical roles in promoting different ground states and
quantum excitations. Both the crystal field environment and the delicate
balance between the anisotropic exchanges could be tuned by the local
distortions [39, 40, 41], and such an interplay could potentially induce
accidental degeneracy in the vicinity of the phase boundary and thus could
lead to the emergence of a QSL [42]. Another typical and intricate scenario
involving local distortion is for those systems presenting strong spin-orbital
coupling – on one side, the local distortion is potentially destroying the
degeneracy of energy levels and thus suppressing the quantum fluctuation that
is necessary for the formation of QSL states. On the other, the local symmetry
breaking may lead to orbital quenching and therefore suppress the spin-orbital
coupling, which is then beneficial for the QSL states formation, since the
spin-orbital coupling is susceptible of lifting the degeneracy of energy
levels and thus is detrimental for the QSL states. As such, several typical
examples were showing the dominant effect of suppressing the quantum
fluctuations and its winning over the orbital quenching effect to give exotic
quantum states were reported for the NiRh2O4 system [43, 44, 45, 46]. Such
impacts of local distortion, along with other potential effects such as the
formation of local dimer or trimer clusters [47, 48, 49], induction of fast
spin fluctuations [50], reduction of the effective dimension of magnetic
coupling [51], and among others, infers the importance of directly probing the
local distortion and constructing a clear picture of the interplay with
magnetic coupling.

With this regard, extensive existing studies focus on the impact of local
distortion on the scheme of exchange interactions, and among other aspects
such as electron localization and charge density wave. Typically for QSL
candidate systems, various studies reported the direct probe of the local
environment and its link to the exotic QSL states, using NMR [52, 53],
Electron Spin Resonance (ESR) [54], X-ray Absorption Spectroscopy (XAS) [55],
and X-ray/neutron total scattering [56, 57]. We realized enormous such
research in the herbertsmithite and barlowite QSL candidate systems. However,
for the pyrochlore systems, direct experimental probing of the local
structural variation and its link to the magnetic coupling is still lacking,
though, excessive theoretical work has provided clear indication for the
impact of local structural distortion, via the spin-lattice coupling effect,
upon the magnetic states and spin ordering [58, 59, 60, 61, 62, 63, 64]. To
the best of our knowledge, there is only limited experimental work focusing on
the spin-lattice coupling effect in pyrochlore systems from the local
perspective, like the report by P. M. Thygesen, et al. demonstrating the local
orbital dimerization of Jahn-Teller active Mo4+ ions instead of random
compositional or site disorder drives the spin-glass state in Y2Mo2O7 [65].
More relevant experimental efforts, though still limited, were mainly from the
average long-range structure perspective [66, 67, 68, 20]. However, while the
local distortion could be extracted from conventional Bragg diffraction, the
local and short-range probe could provide unique pathway and different angle
to such explorations, as is already demonstrated by the work on the
herbertsmithite and barlowite QSL candidate systems [52, 53, 54, 55, 56, 57].
In this report, we were trying to utilize the neutron total scattering
measurement, which incorporates both the Bragg peaks and the diffuse
scattering signal, to study the spin-lattice coupling in the Sn-doped Yb2Ti2O7
(YTSO). For such a system, it is believed only the nearest neighbor coupling
is dominant in the magnetic Hamiltonian [69, 70, 6, 71, 72]. Meanwhile, the
doping in our studied YTSO system is only on the non-magnetic site, which
infers the chemical disorder effect in the YTSO system could be singled out
from the potentially existing distortion effect. Through our modeling for the
neutron total scattering data together with magnetic susceptibility
measurements, a magnetic phase switching could be clearly identified, as
proposed to be induced by the local distortion associated with the Sn-doping
on the non-magnetic Ti sites.

Figure 1: Average structure of Yb2(Ti1-xSnx)2O7 and Yb2(Ti1-xGex)2O7. a.
Neutron powder diffraction pattern (red circles) for Yb2(Ti0.6Sn0.4)2O7
measured at 300 K with the central wavelength of 1.5 Å on POWGEN
diffractometer. The solid black line is the Rietveld refinement using FullProf
[73]. Solid blue line at the bottom of the panel shows the difference curve.
The Bragg peaks are marked with green markers. b. The doping concentration
dependence of the lattice parameter. c. The doping level dependence of the
ratio between the two different Yb-O bond lengths and the Yb-O2-Yb angle.

Neutron diffraction. For comparison, we synthesized both Yb2(Ti1-xSnx)2O7 and
Yb2(Ti1-xGex)2O7 samples and used neutron powder diffraction (NPD) to
characterize their lattice structures. Fig. 1(a) shows the refinement for the
NPD data of Yb2(Ti0.6Sn0.4)2O7 measured at room temperature using the POWGEN
diffractometer. The data could be well fitted by the $Fd\bar{3}m$ pyrochlore
structure. The NPD data for several other Sn and Ge doped samples was also
refined (not shown here), which all exhibits pure pyrochlore structure. As
summarized in Fig. 1(c), the lattice parameter $a$ decreases from Yb2Sn2O7 to
Yb2Ti2O7 and then Yb2Ge2O7 for all doped samples. This is reasonable since the
lattice parameter is strongly related to the ironic radius of the (Sn/Ti/Ge)
site, and therefore the Sn sample has the largest lattice parameter, Ti sample
the second, and Ge sample the smallest. We further used
$\rho=d_{\textrm{Yb-O2}}/d_{\textrm{Yb-O1}}$ and the Yb-O2-Yb angle to
characterize the axial distortion of the YbO8 polyhedra, here
$d_{\textrm{Yb-O1}}$ represents the bond length for the 6 longer Yb-O1 bonds
in the plane perpendicular to the $\langle 111\rangle$ axis and
$d_{\textrm{Yb-O2}}$ represents the bond length for the two shorter Yb-O2
bonds along the $\langle 111\rangle$ axis. As shown in Fig. 1(d), again, both
of them decrease linearly, without abrupt change.

Figure 2: AC susceptibility. a. Real part of AC susceptibility measured under
different DC magnetic fields in the arbitrary unit (a.u.) for Yb2(Ti1-xSnx)2O7
with $x$ = 0.2, 0.4, 0.6, and 0.8. b. Similar data for Yb2(Ti1-xGex)2O7 with
$x$ = 0.2, 0.4, 0.6, and 0.8. The used AC frequency is 3317 Hz for the Sn-
doped system and 331 Hz for the Ge-doped system with a magnitude of 5 Oe.
Dashed arrows indicate the evolution of the peak’s position with increasing DC
fields.

AC susceptibility. Fig. 2 shows the AC susceptibility measured at different DC
magnetic fields for Yb2(Ti1-xSnx)2O7 and Yb2(Ti1-xGex)2O7. For all samples,
the data at zero field exhibits a peak, which represents the long range
magnetic ordering at $T$*, and such a peak shows an obvious shift under
applied DC fields. As demonstrated in Ref. [74], the field dependence of AC
susceptibility can be used as a convenient tool to identify the nature of a
long-range magnetic ordering, i.e., the ferromagnetic (FM) ordering
temperature will shift to higher temperatures with increasing DC field due to
the contribution of domain magnetization while the antiferromagnetic (AFM)
ordering temperature will show an opposite DC field dependence.

The field dependence of $T$* for each doping level was summarized in Fig.
3(a). The data shows (i) for Yb2Ti2O7 and Yb2Sn2O7, the $T$* increases with
increasing field, which is consistent with the fact that both samples have a
splayed ferromagnetic (SF) ground state; (ii) for Yb2Ge2O7, the $T$* decreases
with increasing field, which is consistent with its AFM ground state [75, 76,
74]; (iii) for all Sn- and Ge-doped samples except Yb2(Ti0.8Sn0.2)2O7, the
$T$* decreases first with increasing field and then increases while the field
surpasses a critical value $H_{\text{c}}$. This indicates that as soon as Sn
and Ge are doped, certain volume of AFM phase is introduced. This AFM order
should be in long range nature since it dominates the bulk magnetism at low
fields. With $H$ $>$ $H_{\text{c}}$, the sample comes back to ferromagnetic or
is fully polarized; (iv) for Yb2(Ti0.8Sn0.2)2O7, it has ferromagnetic ground
state since its $T$* monotonically increases with increasing field.

Magnetic Phase Diagram. Accordingly, a magnetic phase diagram of $T$* and
$H_{\text{c}}$ for Yb2(Ti1-xSnx)2O7 and Yb2(Ti1-xGex)2O7 is summarized in Fig.
3(b). For Ge-doped samples, both $T$* and $H_{\text{c}}$ monotonically
increase with increasing Ge-doping level. On the other hand, for Sn-doped
samples, (i) while the $T$* generally decreases with increasing Sn-doping
level, it exhibits a dome around $x$ = 0.5; (ii) the $H_{\text{c}}$ first
increases with increasing Sn-doping level, peaks at $x$ = 0.5, and thereafter
decreases.

Figure 3: Magnetic phase diagram. a. The field dependence of the magnetic
ordering temperature T∗ for Yb2(Ti1-xSnx)2O7 and Yb2(Ti1-xGex)2O7 with
different doping concentrations. b. Magnetic phase diagram as a function of
the doping concentration, $x$, temperature, $T$, and the external DC field,
$\mu_{0}H$. The red and blue regions represent the FM and AFM phases,
respectively. Data points for Yb2Ti2O7, Yb2Sn2O7, and Yb2Ge2O7 are from Dun et
al. [74].

The evolution of $T$* and $H_{\text{c}}$ in the Ge-doped samples is expected.
The magnetic ground states of Yb-pyrochlores are determined by the ratio among
the anisotropic exchange interactions [75, 77, 78]. Unlike Yb2Ti2O7, Yb2Ge2O7
orders antiferromagnetically in the $\Gamma_{5}$ manifold [76]. From the point
view of chemical pressure effect, with increasing the Ge-doping level in
Yb2(Ti1-xGex)2O7, the lattice parameter decreases monotonically and gradually
tunes the balance of anisotropic exchange interactions that drives the system
towards the AFM $\Gamma_{5}$ phase from the SF phase [75, 76]. Alternatively,
if we assume that the phase coexistence in Ge-doped samples is similar to that
of Yb2Ti2O7, then the $H_{\text{c}}$ could be scaled to the volume fraction of
the AFM phase since the larger the $H_{\text{c}}$ is, the more difficult to
polarize the system. Therefore, the evolution of $H_{\text{c}}$ in Fig. 3(b)
means a monotonic increase in volume fraction of the AFM phase in
Yb2(Ti1-xGex)2O7, consistent with the SF and AFM phases in Yb2Ti2O7 and
Yb2Ge2O7, respectively.

However, although we cannot determinate whether it is the $\psi_{2}$ or
$\psi_{3}$ phase of the $\Gamma_{5}$ manifold, the appearance of the long
range AFM order in Yb2(Ti1-xSnx)2O7 is surprising. Since both Yb2Ti2O7 and
Yb2Sn2O7 have the SF ground state, an AFM ground state should not be expected
for Sn-doped samples from the view of chemical pressure effects. Even if there
is still magnetic phase coexistence [79], it is puzzling to observe this non-
monotonic change of the volume fraction of AFM phase (or, $H_{\text{c}}$) in
Yb2(Ti1-xSnx)2O7.

Figure 4: Local distortion as characterized by various local bond angles, as
extracted from the RMC resulted configuration. The inset is an illustration
for the local geometry involving the critical triplet angle – Ti&Sn-O-Ti&Sn.
Yb $\rightarrow$ cyan, Ti&Sn $\rightarrow$ blue, O $\rightarrow$ red. Error
bar was estimated via running 10 repeated RMC modeling followed by the same
bond angle calculation. The actual error bar is smaller than the symbol for
all the data points as presented and here we just take the symbol size as the
representative error bar level.

Total Scattering. The total scattering signal contains both the Bragg peaks
and the diffuse scattering contribution. The diffuse scattering part, which is
usually taken as the background and thus subtracted off in conventional Bragg
peaks analysis, in fact provides unique access to the local structure.
Although both the Bragg peaks and total scattering data would yield
information about local distortions, the former is based on the distance
between average positions whereas the latter is based on the average of
distances ensemble. Therefore, the local probe via total scattering data could
provide access to the local structural variation that is inaccessible through
the conventional Bragg analysis. Here, we collected neutron total scattering
data on POWGEN for the series of Yb2(Ti1-xSnx)2O7 samples and performed the
reverse Monte Carlo (RMC) [80, 81] modeling [82] to extract the local
structure. A 10$\times$10$\times$10 supercell containing more than 100, 000
atoms could be constructed from the RMC modeling and the following statistical
calculation could be performed to extract the key structural aspects.
Typically, here we focused on the various bond angles involving both the Yb,
Ti/Sn, and O atoms – from the RMC resulted configuration, a statistical
distribution of various triplet angles could be obtained [83]. Further, the
width of the triplet angles distribution could be extracted through a Gaussian
peak fitting, which indicates the dispersiveness of the bond angle and thus
infers the significance of the local distortion. The results are presented in
Fig. 4. No abrupt change as the function of the doping level could be observed
except for the Ti&Sn-O-Ti&Sn triplet (see the inset of Fig. 4). This indicates
the non-magnetic-site involved local distortion is enhanced as the result of
the doping. Assuming the nearest neighbor magnetic interaction is dominant in
the YTSO system, the result here infers the potential link between the
magnetic phase switching and the local distortion. As of our knowledge, this
is the first experimental evidence demonstrating the magnetic phase switching
in pyrochlore systems that is potentially coupled with the local structural
distortion beyond the compositional randomness or site disorder. Our finding
here falls in line with various theoretical studies for the spin-lattice
coupling in pyrochlore systems. For example, through many-body quantum-
chemical calculations, N. Bogdanov et al. showed the magnetic interactions and
ordering in Cd2Os2O7 are crucially dependent on the local geometrical features
[62]. H. Shinaoka, et al. showed the coupling of spin-glass transitions to
local lattice distortions on pyrochlore lattices via Monte Carlo simulation
[59]. Also using Monte Carlo simulation, K. Aoyama et al. revealed a lattice
distortion induced spin ordering in the breathing pyrochlore lattice[58].
Based on such theoretical results, we believe the non-magnetic site mixing in
the YTSO system is imposing effect upon the magnetic coupling through the
induced local lattice distortion, which then further onsets the magnetic phase
switching at a certain level of mixing. According to a recent report by A.
Scheie et al., it was revealed that the QSL physics in Yb2Ti2O7 is
fundamentally the FM-AFM phase competition in the low temperature regime [79].
With this regard, the FM-AFM phase switching and its link to the local lattice
distortion is of strong relevance and infers that the non-magnetic site doping
and the accordingly induced local structural variation should be paid serious
attention for the study of QSL physics in pyrochlore systems.

As a final remark, we want to emphasize that the conclusion we arrived at here
about there being a strong correlation between the local distortion and
magnetic phase switching is partly based on the assumption that the magnetic
coupling beyond the nearest neighbor is not taking critical effect in the YTSO
system – in fact, such an assumption was indeed adopted in various previous
reports as mentioned earlier [69, 70, 6, 71, 72]. Without such an assumption,
the induced randomness in the Hamiltonian as the result of chemical disorder
and its effect upon the magnetic ordering cannot be rigorously excluded. As
such, future studies are needed – theoretically, magnetic interaction beyond
the nearest neighbor needs to be inspected and its effect on the magnetic
phase diagram should be studied.

###### Acknowledgements.

We thank Martin Mourigal and Itamar Kimchi for helpful discussion. J.G.C. is
supported by the National Natural Science Foundation of China (12025408,
11874400, 11921004), the Key Research Program of Frontier Sciences of CAS (
QYZDB-SSW-SLH013), the CAS Interdisciplinary Innovation Team (JCTD-2019-01)
and Lujiaxi international group funding of K. C. Wong Education Foundation
(GJTD-2020-01). The work at the University of Tennessee is supported by the U.
S. Department of Energy under grant No. DE-SC0020254. Part of the research
conducted at SNS was sponsored by the Scientific User Facilities Division,
Office of Basic Energy Sciences, US Department of Energy. The following
funding is acknowledged: US Department of Energy, Office of Science (contract
No. DE-AC05-00OR22725). This research used resources of the Computeand Data
Environment for Science (CADES) at the Oak Ridge National Laboratory,which is
supported by the Office of Science of the U.S. Department of Energy under
Contract No. DE-AC05-00OR22725. This research used resources of the Compute
and Data Environment for Science (CADES) at the Oak Ridge National Laboratory,
which is supported by the Office of Science of the U.S. Department of Energy
under Contract No. DE-AC05-00OR22725. This research used resources of the
National Energy Research Scientific Computing Center (NERSC), a U.S.
Department of Energy Office of Science User Facility located at Lawrence
Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231
using NERSC award ERCAP0024340. We thank Dr. Emil S. Bozin for the discussion
about the local structure studies.

## References

  * Gardner  _et al._ [2010] J. S. Gardner, M. J. P. Gingras, and J. E. Greedan, Rev. Mod. Phys. 82, 53 (2010).
  * Shamblin  _et al._ [2017] J. Shamblin, Z. Dun, M. Lee, S. Johnston, E. S. Choi, K. Page, Y. Qiu, and H. Zhou, Phys. Rev. B 96, 174418 (2017).
  * Ross [2012] K. A. Ross, Phd thesis, McMaster University (2012).
  * Ross  _et al._ [2009] K. A. Ross, J. P. C. Ruff, C. P. Adams, J. S. Gardner, H. A. Dabkowska, Y. Qiu, J. R. D. Copley, and B. D. Gaulin, Phys. Rev. Lett. 103, 227202 (2009).
  * Lu  _et al._ [2018] Z. Lu, L. Ge, G. Wang, M. Russina, G. Günther, C. R. dela Cruz, R. Sinclair, H. D. Zhou, and J. Ma, Phys. Rev. B 98, 094412 (2018).
  * Applegate  _et al._ [2012] R. Applegate, N. R. Hayre, R. R. P. Singh, T. Lin, A. G. R. Day, and M. J. P. Gingras, Phys. Rev. Lett. 109, 097205 (2012).
  * Banerjee  _et al._ [2018] M. Banerjee, M. Heiblum, V. Umansky, D. E. Feldman, Y. Oreg, and A. Stern, Nature 559, 205 (2018).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Willett  _et al._ [1987] R. Willett, J. P. Eisenstein, H. L. Störmer, D. C. Tsui, A. C. Gossard, and J. H. English, Phys. Rev. Lett. 59, 1776 (1987).
  * Chamorro  _et al._ [2021] J. R. Chamorro, T. M. McQueen, and T. T. Tran, Chemical Reviews 121, 2898 (2021).
  * Balents [2010] L. Balents, Nature 464, 199 (2010).
  * Anderson [1987] P. W. Anderson, Science 235, 1196 (1987).
  * Villain [1979] J. Villain, Zeitschrift für Physik B Condensed Matter 33, 31 (1979).
  * Henley [1989] C. L. Henley, Phys. Rev. Lett. 62, 2056 (1989).
  * Keren  _et al._ [2004] A. Keren, J. S. Gardner, G. Ehlers, A. Fukaya, E. Segal, and Y. J. Uemura, Phys. Rev. Lett. 92, 107204 (2004).
  * Sheetal  _et al._ [2020] Sheetal, A. Ali, S. Rajput, Y. Singh, T. Maitra, and C. S. Yadav, Journal of Physics: Condensed Matter 32, 365804 (2020).
  * Takahashi  _et al._ [2019] S. K. Takahashi, J. Wang, A. Arsenault, T. Imai, M. Abramchuk, F. Tafti, and P. M. Singer, Phys. Rev. X 9, 031047 (2019).
  * Dantas and Andrade [2022] V. Dantas and E. C. Andrade, Phys. Rev. Lett. 129, 037204 (2022).
  * Kao  _et al._ [2021] W.-H. Kao, J. Knolle, G. B. Halász, R. Moessner, and N. B. Perkins, Phys. Rev. X 11, 011034 (2021).
  * Ross  _et al._ [2012] K. A. Ross, T. Proffen, H. A. Dabkowska, J. A. Quilliam, L. R. Yaraskavitch, J. B. Kycia, and B. D. Gaulin, Phys. Rev. B 86, 174424 (2012).
  * Bowman  _et al._ [2019] D. F. Bowman, E. Cemal, T. Lehner, A. R. Wildes, L. Mangin-Thro, G. J. Nilsen, M. J. Gutmann, D. J. Voneshen, D. Prabhakaran, A. T. Boothroyd, D. G. Porter, C. Castelnovo, K. Refson, and J. P. Goff, Nature Communications 10, 637 (2019).
  * Lau  _et al._ [2006] G. C. Lau, R. S. Freitas, B. G. Ueland, B. D. Muegge, E. L. Duncan, P. Schiffer, and R. J. Cava, Nature Physics 2, 249 (2006).
  * Gardner  _et al._ [2011] J. S. Gardner, G. Ehlers, P. Fouquet, B. Farago, and J. R. Stewart, Journal of Physics: Condensed Matter 23, 164220 (2011).
  * Li  _et al._ [2016] Y. Li, D. Adroja, P. K. Biswas, P. J. Baker, Q. Zhang, J. Liu, A. A. Tsirlin, P. Gegenwart, and Q. Zhang, Phys. Rev. Lett. 117, 097201 (2016).
  * Li  _et al._ [2019] Y. Li, S. Bachus, B. Liu, I. Radelytskyi, A. Bertin, A. Schneidewind, Y. Tokiwa, A. A. Tsirlin, and P. Gegenwart, Phys. Rev. Lett. 122, 137201 (2019).
  * Shen  _et al._ [2016] Y. Shen, Y.-D. Li, H. Wo, Y. Li, S. Shen, B. Pan, Q. Wang, H. C. Walker, P. Steffens, M. Boehm, Y. Hao, D. L. Quintero-Castro, L. W. Harriger, M. D. Frontzek, L. Hao, S. Meng, Q. Zhang, G. Chen, and J. Zhao, Nature 540, 559 (2016).
  * Paddison  _et al._ [2017] J. A. M. Paddison, M. Daum, Z. Dun, G. Ehlers, Y. Liu, M. B. Stone, H. Zhou, and M. Mourigal, Nature Physics 13, 117 (2017).
  * Shen  _et al._ [2018] Y. Shen, Y.-D. Li, H. C. Walker, P. Steffens, M. Boehm, X. Zhang, S. Shen, H. Wo, G. Chen, and J. Zhao, Nature Communications 9, 4138 (2018).
  * Sana  _et al._ [2023] B. Sana, M. Barik, S. Lee, U. Jena, M. Baenitz, J. Sichelschmidt, S. Luther, H. Kuehne, K. Sethupathi, M. S. R. Rao, K. Y. Choi, and P. Khuntia, (2023), arXiv:2304.13116 [cond-mat.str-el] .
  * Hong  _et al._ [2021] W. Hong, L. Liu, C. Liu, X. Ma, A. Koda, X. Li, J. Song, W. Yang, J. Yang, P. Cheng, H. Zhang, W. Bao, X. Ma, D. Chen, K. Sun, W. Guo, H. Luo, A. W. Sandvik, and S. Li, Phys. Rev. Lett. 126, 037201 (2021).
  * Fogh  _et al._ [2022] E. Fogh, O. Mustonen, P. Babkevich, V. M. Katukuri, H. C. Walker, L. Mangin-Thro, M. Karppinen, S. Ward, B. Normand, and H. M. Rønnow, Phys. Rev. B 105, 184410 (2022).
  * Ren  _et al._ [2023] H.-D. Ren, T.-Y. Xiong, H.-Q. Wu, D. N. Sheng, and S.-S. Gong, Phys. Rev. B 107, L020407 (2023).
  * Watanabe  _et al._ [2022] M. Watanabe, N. Kurita, H. Tanaka, W. Ueno, K. Matsui, T. Goto, and M. Hagihala, Phys. Rev. B 105, 054414 (2022).
  * Lee  _et al._ [2023] C. Lee, S. Lee, H.-S. Kim, S. Kittaka, Y. Kohama, T. Sakakibara, K. H. Lee, J. van Tol, D. I. Gorbunov, S.-H. Do, S. Yoon, A. Berlie, and K.-Y. Choi, Phys. Rev. B 107, 214404 (2023).
  * Koga  _et al._ [2016] T. Koga, N. Kurita, M. Avdeev, S. Danilkin, T. J. Sato, and H. Tanaka, Phys. Rev. B 93, 054426 (2016).
  * Vasala  _et al._ [2014a] S. Vasala, M. Avdeev, S. Danilkin, O. Chmaissem, and M. Karppinen, Journal of Physics: Condensed Matter 26, 496001 (2014a).
  * Vasala  _et al._ [2014b] S. Vasala, H. Saadaoui, E. Morenzoni, O. Chmaissem, T.-S. Chan, J.-M. Chen, Y.-Y. Hsu, H. Yamauchi, and M. Karppinen, Phys. Rev. B 89, 134419 (2014b).
  * Kimchi  _et al._ [2018] I. Kimchi, A. Nahum, and T. Senthil, Phys. Rev. X 8, 031028 (2018).
  * Norman [2016] M. R. Norman, Rev. Mod. Phys. 88, 041002 (2016).
  * Yadav  _et al._ [2016] R. Yadav, N. A. Bogdanov, V. M. Katukuri, S. Nishimoto, J. van den Brink, and L. Hozoi, Scientific Reports 6, 37925 (2016).
  * Balz  _et al._ [2016] C. Balz, B. Lake, J. Reuther, H. Luetkens, R. Schönemann, T. Herrmannsdörfer, Y. Singh, A. T. M. Nazmul Islam, E. M. Wheeler, J. A. Rodriguez-Rivera, T. Guidi, G. G. Simeoni, C. Baines, and H. Ryll, Nature Physics 12, 942 (2016).
  * Kermarrec  _et al._ [2017] E. Kermarrec, J. Gaudet, K. Fritsch, R. Khasanov, Z. Guguchia, C. Ritter, K. A. Ross, H. A. Dabkowska, and B. D. Gaulin, Nature Communications 8, 14810 (2017).
  * Chamorro  _et al._ [2018] J. R. Chamorro, L. Ge, J. Flynn, M. A. Subramanian, M. Mourigal, and T. M. McQueen, Phys. Rev. Mater. 2, 034404 (2018).
  * Buessen  _et al._ [2018] F. L. Buessen, M. Hering, J. Reuther, and S. Trebst, Phys. Rev. Lett. 120, 057201 (2018).
  * Das  _et al._ [2019] S. Das, D. Nafday, T. Saha-Dasgupta, and A. Paramekanti, Phys. Rev. B 100, 140408 (2019).
  * Chen [2017] G. Chen, Phys. Rev. B 96, 020412 (2017).
  * Paul  _et al._ [2020] A. Paul, C.-M. Chung, T. Birol, and H. J. Changlani, Phys. Rev. Lett. 124, 167203 (2020).
  * Kelly  _et al._ [2019] Z. A. Kelly, T. T. Tran, and T. M. McQueen, Inorganic Chemistry 58, 11941 (2019).
  * Flint and Lee [2013] R. Flint and P. A. Lee, Phys. Rev. Lett. 111, 217201 (2013).
  * Bonville  _et al._ [2011] P. Bonville, I. Mirebeau, A. Gukasov, S. Petit, and J. Robert, Phys. Rev. B 84, 184409 (2011).
  * Zhang  _et al._ [2018] B. Zhang, P. J. Baker, Y. Zhang, D. Wang, Z. Wang, S. Su, D. Zhu, and F. L. Pratt, Journal of the American Chemical Society 140, 122 (2018).
  * Pasco  _et al._ [2018] C. M. Pasco, B. A. Trump, T. T. Tran, Z. A. Kelly, C. Hoffmann, I. Heinmaa, R. Stern, and T. M. McQueen, Phys. Rev. Mater. 2, 044406 (2018).
  * Malavi  _et al._ [2020] P. Malavi, S. Pal, D. V. S. Muthu, S. Sahoo, S. Karmakar, and A. K. Sood, Phys. Rev. B 101, 214402 (2020).
  * Zorko  _et al._ [2017] A. Zorko, M. Herak, M. Gomilšek, J. van Tol, M. Velázquez, P. Khuntia, F. Bert, and P. Mendels, Phys. Rev. Lett. 118, 017202 (2017).
  * Smaha  _et al._ [2020] R. W. Smaha, I. Boukahil, C. J. Titus, J. M. Jiang, J. P. Sheckelton, W. He, J. Wen, J. Vinson, S. G. Wang, Y.-S. Chen, S. J. Teat, T. P. Devereaux, C. Das Pemmaraju, and Y. S. Lee, Phys. Rev. Mater. 4, 124406 (2020).
  * Shoemaker  _et al._ [2011] D. P. Shoemaker, R. Seshadri, M. Tachibana, and A. L. Hector, Phys. Rev. B 84, 064117 (2011).
  * Bozin  _et al._ [2023] E. S. Bozin, M. Abeykoon, S. Conradson, G. Baldinozzi, P. Sutar, and D. Mihailovic, (2023), arXiv:2301.05670 [cond-mat.str-el] .
  * Aoyama and Kawamura [2019] K. Aoyama and H. Kawamura, Phys. Rev. B 99, 144406 (2019).
  * Shinaoka  _et al._ [2011] H. Shinaoka, Y. Tomita, and Y. Motome, Phys. Rev. Lett. 107, 047204 (2011).
  * Haku  _et al._ [2016] T. Haku, K. Kimura, Y. Matsumoto, M. Soda, M. Sera, D. Yu, R. A. Mole, T. Takeuchi, S. Nakatsuji, Y. Kono, T. Sakakibara, L.-J. Chang, and T. Masuda, Phys. Rev. B 93, 220407 (2016).
  * Huang  _et al._ [2014] Y.-P. Huang, G. Chen, and M. Hermele, Phys. Rev. Lett. 112, 167203 (2014).
  * Savary and Balents [2017] L. Savary and L. Balents, Phys. Rev. Lett. 118, 087203 (2017).
  * Bogdanov  _et al._ [2013] N. A. Bogdanov, R. Maurice, I. Rousochatzakis, J. van den Brink, and L. Hozoi, Phys. Rev. Lett. 110, 127206 (2013).
  * Berg  _et al._ [2003] E. Berg, E. Altman, and A. Auerbach, Phys. Rev. Lett. 90, 147204 (2003).
  * Thygesen  _et al._ [2017] P. M. M. Thygesen, J. A. M. Paddison, R. Zhang, K. A. Beyer, K. W. Chapman, H. Y. Playford, M. G. Tucker, D. A. Keen, M. A. Hayward, and A. L. Goodwin, Phys. Rev. Lett. 118, 067201 (2017).
  * Li  _et al._ [2014] X. Li, W. M. Li, K. Matsubayashi, Y. Sato, C. Q. Jin, Y. Uwatoko, T. Kawae, A. M. Hallas, C. R. Wiebe, A. M. Arevalo-Lopez, J. P. Attfield, J. S. Gardner, R. S. Freitas, H. D. Zhou, and J.-G. Cheng, Phys. Rev. B 89, 064409 (2014).
  * Gao  _et al._ [2021] S. Gao, A. F. May, M.-H. Du, J. A. M. Paddison, H. S. Arachchige, G. Pokharel, C. dela Cruz, Q. Zhang, G. Ehlers, D. S. Parker, D. G. Mandrus, M. B. Stone, and A. D. Christianson, Phys. Rev. B 103, 214418 (2021).
  * Zhao  _et al._ [2016] Z. Y. Zhao, S. Calder, A. A. Aczel, M. A. McGuire, B. C. Sales, D. G. Mandrus, G. Chen, N. Trivedi, H. D. Zhou, and J.-Q. Yan, Phys. Rev. B 93, 134426 (2016).
  * Rau  _et al._ [2019] J. G. Rau, R. Moessner, and P. A. McClarty, Phys. Rev. B 100, 104423 (2019).
  * Ross  _et al._ [2011] K. A. Ross, L. Savary, B. D. Gaulin, and L. Balents, Phys. Rev. X 1, 021002 (2011).
  * Hayre  _et al._ [2013] N. R. Hayre, K. A. Ross, R. Applegate, T. Lin, R. R. P. Singh, B. D. Gaulin, and M. J. P. Gingras, Phys. Rev. B 87, 184423 (2013).
  * Thompson  _et al._ [2011] J. D. Thompson, P. A. McClarty, H. M. Rønnow, L. P. Regnault, A. Sorge, and M. J. P. Gingras, Phys. Rev. Lett. 106, 187202 (2011).
  * Rodríguez-Carvajal [1993] J. Rodríguez-Carvajal, Recent advances in magnetic structure determination by neutron powder diffraction, Phys. B: Condens. Matter 192, 55 (1993).
  * Dun  _et al._ [2014] Z. L. Dun, M. Lee, E. S. Choi, A. M. Hallas, C. R. Wiebe, J. S. Gardner, E. Arrighi, R. S. Freitas, A. M. Arevalo-Lopez, J. P. Attfield, H. D. Zhou, and J. G. Cheng, Phys. Rev. B 89, 064401 (2014).
  * Dun  _et al._ [2015] Z. L. Dun, X. Li, R. S. Freitas, E. Arrighi, C. R. Dela Cruz, M. Lee, E. S. Choi, H. B. Cao, H. J. Silverstein, C. R. Wiebe, J. G. Cheng, and H. D. Zhou, Phys. Rev. B 92, 140407 (2015).
  * Hallas  _et al._ [2016] A. M. Hallas, J. Gaudet, M. N. Wilson, T. J. Munsie, A. A. Aczel, M. B. Stone, R. S. Freitas, A. M. Arevalo-Lopez, J. P. Attfield, M. Tachibana, C. R. Wiebe, G. M. Luke, and B. D. Gaulin, Phys. Rev. B 93, 104405 (2016).
  * Wong  _et al._ [2013] A. W. C. Wong, Z. Hao, and M. J. P. Gingras, Phys. Rev. B 88, 144402 (2013).
  * Taillefumier  _et al._ [2017] M. Taillefumier, O. Benton, H. Yan, L. D. C. Jaubert, and N. Shannon, Phys. Rev. X 7, 041057 (2017).
  * Scheie  _et al._ [2020] A. Scheie, J. Kindervater, S. Zhang, H. J. Changlani, G. Sala, G. Ehlers, A. Heinemann, G. S. Tucker, S. M. Koohpayeh, and C. Broholm, Proc. Natl. Acad. Sci. 117, 27245 (2020).
  * Zhang  _et al._ [2020] Y. Zhang, M. Eremenko, V. Krayzman, M. G. Tucker, and I. Levin, Journal of Applied Crystallography 53, 1509 (2020).
  * Tucker  _et al._ [2007] M. G. Tucker, D. A. Keen, M. T. Dove, A. L. Goodwin, and Q. Hui, Journal of Physics: Condensed Matter 19, 335218 (2007).
  * the [a] See Supplemental Material at [URL to be inserted] for the technical details about the reverse Monte Carlo modeling .
  * the [b] See Supplemental Material at [URL to be inserted] for the statistical distribution plots for the various triplet angles .
  * Rao  _et al._ [2021] X. Rao, G. Hussain, Q. Huang, W. J. Chu, N. Li, X. Zhao, Z. Dun, E. S. Choi, T. Asaba, L. Chen, L. Li, X. Y. Yue, N. N. Wang, J.-G. Cheng, Y. H. Gao, Y. Shen, J. Zhao, G. Chen, H. D. Zhou, and X. F. Sun, Nature Communications 12, 4949 (2021).


# Pulsar-timing measurement of the circular polarization of the stochastic
gravitational-wave background

Gabriela Sato-Polito<EMAIL_ADDRESS>Department of Physics and Astronomy,
Johns Hopkins University,
Baltimore, MD 21218, USA  Marc Kamionkowski  Department of Physics and
Astronomy, Johns Hopkins University,
Baltimore, MD 21218, USA

###### Abstract

Pulsar-timing arrays (PTAs) are in the near future expected to detect a
stochastic gravitational-wave background (SGWB) produced by a population of
inspiralling super-massive black hole binaries. In this work, we consider a
background that can be anisotropic and circularly polarized. We use the
expansion of the intensity and the circular polarization in terms of spherical
harmonics and the overlap reduction functions for each term in this
expansion.We propose an unbiased real-space estimator that can separate the
intensity and circular-polarization contributions of the SGWB to pulsar-
timing-residual correlations and then validate the estimator on simulated
data. We compute the signal-to-noise ratio of a circular-polarization
component that has a dipole pattern under different assumptions about the PTA.
We find that a nearly-maximal circular-polarization dipole may be detectable,
which can aid in determining whether or not the background is dominated by a
handful of bright sources.

##  I Introduction

Pulsar-timing arrays (PTAs) are a network of millisecond pulsars that can be
used as a gravitational-wave detector on galactic scales 1978SvA….22…36S ;
1979ApJ…234.1100D . The measurement is based on the precise determination of
arrival times of radio pulses, which feature spatially correlated fluctuations
in the presence of gravitational waves (GWs). The primary source of GWs in the
PTA band (roughly $1-100$ nHz) is expected to be a stochastic gravitational-
wave background (SGWB) produced by a population of inspiralling super-massive
black hole binaries (SMBHBs)Rajagopal:1994zj ; Jaffe:2002rt . Other sources
include cosmic strings Olmez:2010bi ; Sousa:2013aaa ; Miyamoto:2012ck ;
Kuroyanagi:2012jf , phase transitions Caprini:2010xv , and relic GWs from
inflation Starobinsky:1979ty ; Zhao:2013bba .

Current PTA collaborations include the European PTA (EPTA) Desvignes:2016yex ,
the North American Nanohertz Observatory for Gravitational Waves (NANOGrav)
Brazier:2019mmu , and the Parkes Pulsar Timing Array (PPTA) Kerr:2020qdo ,
which together form the International Pulsar Timing Array (IPTA)
Perera:2019sca . While a detection of a SGWB is yet to be claimed,
increasingly tight upper limits on the SGWB amplitude are being set and
NANOGrav has recently reported evidence for a stochastic red-noise process
across all pulsars. These limits have important consequences for the sources
of the GW background (see, e.g., Refs. Taylor:2016ftv  and Chen:2018znx ).

If the origin of a stochastic background is cosmological or sourced by a large
population of distant objects, the standard assumption of statistical isotropy
is expected to be a good one. In this regime, the correlation between pulsar
timing residuals depends only on their angular separation and is given by the
Hellings and Downs curve 1983ApJ…265L..39H . However, some degree of
statistical anisotropy is expected. For instance, the finite number of SMBHBs
or the presence of a nearby bright source can lead to an anisotropic
background Ravi:2012bz ; Cornish:2013aba ; Sesana:2008mz . Characterizing the
anisotropy of the background can therefore be a powerful probe of this
astrophysical population.

Another relevant feature of the SGWB that can be probed with PTAs is its
polarization. Analogously to electromagnetism, one can define the Stokes
parameters for GWs, and the focus of this work is the circular-polarization
component. While the chirality of a background produced by many distant
sources is expected to vanish, it may be non-negligible if a reduced number of
bright sources dominate, since SMBHBs produce chiral GWs if the orbit is
observed face-on. In this scenario, anisotropies in the circular polarization
would also be produced, as in the case for intensities. Previous work has
addressed the circular polarization of the SWGB in the context of measurements
with ground and space-based interferometers Seto:2006dz ; Seto:2006hf ;
Seto:2007tn ; Seto:2008sr , astrometry Qin:2018yhy , and PTAs Kato:2015bye ;
Belgacem:2020nda .

Ref. Kato:2015bye  studied the effect of intensity and circular-polarization
anisotropies on the correlation between pulsar timing residuals. They propose
a way to infer the spherical-harmonic coefficients of the circular-
polarization anisotropy using linear combinations of two pulsar pairs,
although a priori knowledge of the background is required. By describing the
pulsar positions in the spherical-harmonic basis, Refs. Hotinli:2019tpc  and
Belgacem:2020nda  showed that intensity and circular-polarization anisotropies
correspond to even and odd bipolar spherical harmonic multipoles and can
therefore be measured independently. However, this bipolar-spherical-harmonic
formalism relies on the assumption that one has a complete map of the sky at
hand. Since the expansion is on the angular position of the pulsars, this
would require a dense and uniform distribution of pulsars, which is
incompatible with the current observational capabilities of PTAs.

Here we propose an estimator based on the real-space correlation function,
that can naturally separate the intensity and circular-polarization
contributions. Since circular polarization is defined as the imaginary part of
the cross-correlation between $+$ and $\times$ GW polarizations, it
corresponds to an imaginary contribution to the correlation between pulsars in
the frequency domain. Hence, we construct an estimator for the real and
imaginary parts of the correlation function. We follow a similar approach to
Ref. Anholm:2008wy , but generalize the detection statistic to account for
intensity and circular polarization anisotropies. Furthermore, we extend the
result beyond the weak-signal regime, since any detection of circular-
polarization anisotropy will require a level of sensitivity beyond the
validity of this approximation.

This paper is organized as follows. We begin by outlining the general effect
of GWs on pulsar timing residuals in Sec. II. In Sec. III we discuss the
imprint of the SGWB on the correlation of timing residuals and define the
Stokes parameters for GWs. We present the estimators for intensity and
circular polarization in Sec. IV, apply this estimator for the simple case of
detecting a circular polarization dipole in Sec. V and conclude in Sec. VI.

##  II Timing Residuals

A gravitational wave (GW) propagating between the Earth and a pulsar will
change the observed arrival time of pulses. The fractional frequency shift for
a given pulsar, which we will label by $a$, at a position $\hat{p}_{a}$ on the
sky, induced by a metric perturbation $h_{ij}(t,\hat{\Omega})$ from a GW
propagating in the $\hat{\Omega}$ direction, is given by 1979ApJ…234.1100D ;
Anholm:2008wy

$z_{a}(t,\hat{p}_{a},\hat{\Omega})=\frac{1}{2}\frac{\hat{p}^{i}_{a}\hat{p}^{j}_{a}}{1+\hat{\Omega}\cdot\hat{p}_{a}}\Delta
h_{ij},$ (1)

where $\Delta h_{ij}\equiv
h_{ij}(t_{e},\hat{\Omega})-h_{ij}(t_{p},\hat{\Omega})$ is the difference
between the metric perturbation at the solar system barycenter, with
coordinates $(t_{e},\vec{x}_{e})$, and at the pulsar, with coordinates
$(t_{p},\vec{x}_{p})$. We choose a coordinate system in which the origin is at
the center of the solar system and the pulsar is some distance $L_{a}$ away,
such that $t_{e}=t$, $\vec{x}_{e}=0$, $t_{p}=t-L_{a}$ and
$\vec{x}_{p}=L_{a}\hat{p}_{a}$.

In the transverse-traceless gauge, the metric perturbation at each point can
be written as the following superposition of plane waves Allen:1997ad

$h_{ij}(t,\hat{\Omega})=\sum_{A=+,\times}\int_{-\infty}^{\infty}df\
h_{A}(f,\hat{\Omega})e^{A}_{ij}(\hat{\Omega})e^{2\pi
if(t-\hat{\Omega}\cdot\vec{x})},$ (2)

where the index $A$ labels the $+$ and $\times$ polarizations and $f$ is the
frequency of the GW. The Fourier amplitudes $h_{A}(f,\hat{\Omega})$ are
complex functions that satisfy

$h^{*}_{A}(f,\hat{\Omega})=h_{A}(-f,\hat{\Omega}),$ (3)

and the polarization tensors $e^{A}_{ij}(\hat{\Omega})$ are given by

$\begin{split}e^{+}_{ij}(\hat{\Omega})&=\hat{m}_{i}\hat{m}_{j}-\hat{n}_{i}\hat{n}_{j},\\\
e^{\times}_{ij}(\hat{\Omega})&=\hat{m}_{i}\hat{n}_{j}+\hat{n}_{i}\hat{m}_{j},\end{split}$
(4)

where $\hat{m}$ and $\hat{n}$ are orthogonal unit vectors perpendicular to
$\hat{\Omega}$.

The total frequency shift induced by a stochastic background comprised of GWs
of all frequencies and coming from all directions is then given by

$\begin{split}z_{a}(t)=&\sum_{A=+,\times}\int_{-\infty}^{\infty}df\int
d\hat{\Omega}\ h_{A}(f,\hat{\Omega})F^{A}_{a}(\hat{\Omega})e^{-2\pi ift}\\\
&\times\left(1-e^{2\pi
ifL_{a}(1+\hat{\Omega}\cdot\hat{p}_{a})}\right),\end{split}$ (5)

where $a$ labels each pulsar and the antenna beam pattern is defined as

$F^{A}_{a}(\hat{\Omega})=\frac{\hat{p}^{i}_{a}\hat{p}^{j}_{a}e^{A}_{ij}(\hat{\Omega})}{2(1+\hat{\Omega}\cdot\hat{p}_{a})}.$
(6)

By taking the Fourier transform, we can finally write the timing residual in
the frequency domain as

$z_{a}(f)=\int d^{2}\hat{\Omega}\left(1-e^{2\pi
ifL_{a}(1+\hat{\Omega}\cdot\hat{p}_{a})}\right)\sum_{A}h_{A}(f,\hat{\Omega})F^{A}_{a}(\hat{\Omega}).$
(7)

The explicit expressions for the antenna beam patterns are shown in App. A.

Figure 1: Generalized overlap reduction functions in the computational frame
for intensity and circular polarization as a function of pulsar angular
separation. The monopole is shown on the left panel and the dipole is shown on
the right.

##  III Correlations

The presence of a stochastic gravitational wave background induces
correlations between the timing residuals of pulsars. Here we consider a
background that is Gaussian and stationary, but that can be polarized and
anisotropic. In this case, the correlation function between different GW
polarizations can be written as

$\langle
h^{*}_{A}(f,\hat{\Omega})h_{A^{\prime}}(f^{\prime},\hat{\Omega}^{\prime})\rangle=\delta^{2}(\hat{\Omega},\hat{\Omega}^{\prime})\delta(f-f^{\prime})\mathcal{P}_{AA^{\prime}}(f,\hat{\Omega}),$
(8)

where $\mathcal{P}_{AA^{\prime}}$ is the power spectrum and includes both
frequency and angular dependence.

In analogy with electromagnetism, it is convenient to define the Stokes
parameters for GWs as Kato:2015bye

$\begin{split}I&=\frac{1}{2}\langle|h_{+}|^{2}+|h_{\times}|^{2}\rangle,\\\
Q&=\frac{1}{2}\langle|h_{+}|^{2}-|h_{\times}|^{2}\rangle,\\\
U&=\text{Re}\langle h^{*}_{+}h_{\times}\rangle=\frac{1}{2}\langle
h^{*}_{+}h_{\times}+h^{*}_{\times}h_{+}\rangle,\\\ V&=\text{Im}\langle
h^{*}_{+}h_{\times}\rangle=\frac{1}{2i}\langle
h^{*}_{+}h_{\times}-h^{*}_{\times}h_{+}\rangle.\end{split}$ (9)

With these definitions, the correlation function given in Eq. 8 can be
decomposed as

$\mathcal{P}^{AA^{\prime}}(f,\hat{\Omega})=\begin{pmatrix}I(f,\hat{\Omega})+Q(f,\hat{\Omega})&U(f,\hat{\Omega})-iV(f,\hat{\Omega})\\\
U(f,\hat{\Omega})+iV(f,\hat{\Omega})&I(f,\hat{\Omega})-Q(f,\hat{\Omega})\end{pmatrix}$
(10)

In this work, we will consider only intensity and circular polarization and
therefore assume that $Q=U=0$.

With Eq. 5, we can write the timing residual correlation between pulsars $a$
and $b$. Substituting the definitions of the Stokes parameters in Eq. 9 and
keeping only intensity and circular polarization, we get

$\begin{split}\langle z^{*}_{a}(f)z_{b}(f^{\prime})\rangle=\int
d^{2}\hat{\Omega}\
\kappa_{ab}(f,\hat{\Omega})\delta(f-f^{\prime})\Big{\\{}I(f,\hat{\Omega})(F^{+*}_{a}F^{+}_{b}+F^{\times*}_{a}F^{\times}_{b})+iV(f,\hat{\Omega})(F^{+*}_{a}F^{\times}_{b}-F^{\times*}_{a}F^{+}_{b})\Big{\\}},\end{split}$
(11)

where we have defined

$\kappa_{ab}(f,\hat{\Omega})\equiv\left(1-e^{2\pi
ifL_{a}(1+\hat{\Omega}\cdot\hat{p}_{a})}\right)\left(1-e^{2\pi
ifL_{b}(1+\hat{\Omega}\cdot\hat{p}_{b})}\right).$ (12)

We assume that the frequency and angular dependence of the intensity and
circular polarization power spectra are separable and expand the angular
dependence in spherical harmonics as follows

$X(f,\hat{\Omega})=X(f)\sum_{\ell m}c^{X}_{\ell m}Y_{\ell m}(\hat{\Omega}),$
(13)

where $X=I,V$. Similarly to Ref. Kato:2015bye , we can define the intensity
and circular polarization overlap reduction functions (ORFs) as

${}^{(ab)}\Gamma^{X}=X(f)\sum_{\ell m}c^{X}_{\ell m}{}^{(ab)}\Gamma^{X}_{\ell
m}$ (14)

where we have defined

$\begin{split}{}^{(ab)}\Gamma^{I}_{\ell m}&=\int d^{2}\hat{\Omega}\ Y_{\ell
m}(\hat{\Omega})\kappa_{ab}(f,\hat{\Omega})\left(F^{+*}_{a}F^{+}_{b}+F^{\times*}_{a}F^{\times}_{b}\right),\\\
{}^{(ab)}\Gamma^{V}_{\ell m}&=\int d^{2}\hat{\Omega}\ Y_{\ell
m}(\hat{\Omega})\kappa_{ab}(f,\hat{\Omega})\left(F^{+*}_{a}F^{\times}_{b}-F^{\times*}_{a}F^{+}_{b}\right).\end{split}$
(15)

Note that we have adopted a slightly different convention for the definitions
of the Stokes parameters, but the equations presented here are consistent with
Ref. Kato:2015bye . We take the standard assumption that $L_{a}=L_{b}$ and
that $fL_{a}\gg 1$ 2018JPhCo…2j5002M ; Mingarelli:2014xfa , such that

$\kappa_{ab}(f,\hat{\Omega})\rightarrow(1+\delta_{ab}),$ (16)

and therefore the functions ${}^{(ab)}\Gamma^{X}_{\ell m}$ are independent of
frequency. This approximation is equivalent to neglecting the pulsar term in
the cross-correlations, and only including both Earth and pulsar terms in the
auto-correlations. Notice that including the pulsar term adds a small change
in phase, which corresponds to an imaginary term in the pulsar correlations,
and would therefore contaminate the circular polarization measurement.

Substituting the definitions above into the pulsar timing residual
correlation, we get

$\langle
z^{*}_{a}(f)z_{b}(f^{\prime})\rangle=\delta(f-f^{\prime})\left[I(f)^{(ab)}\Gamma^{I}+iV(f)^{(ab)}\Gamma^{V}\right].$
(17)

We emphasize here that the quantity $V(f)^{(ab)}\Gamma^{V}$ is defined to be
real and that circular polarization therefore induces an imaginary
contribution to the correlation function.

Figure 2: Detectability of the circular polarization dipole as a function of
the pulsar white noise. The panel on the left shows the signal-to-noise ratio
as a function of the noise RMS for PTAs composed of 10, 25, and 50 pulsars
sampled from a uniform spatial distribution. The panel on the right shows the
positions of the pulsars on the sky. Note that the inclusion of pulsars is
additive: the blue curve on the left panel includes only the blue dots in the
map on the right panel, the yellow curve includes the 10 pulsars in blue plus
the 15 pulsars in yellow, and the red curve includes all points in the map
(blue+yellow+red).

##  IV Estimators for intensity and circular polarization

In this section, we introduce the estimators for intensity and circular
polarization of the SGWB. The key result is the signal-to-noise ratio of the
circular polarization statistic when the optimal filter and inverse-variance
weighting of pulsar pairs are employed Allen:1997ad ; Anholm:2008wy . We
consider the covariance between pulsars in the presence of a non-zero SGWB,
thereby accounting for the variance of the background itself Romano:2020sxq .

We assume that the timing residual signal from a given pulsar can be written
as a sum of the contributions from the stochastic background and measurement
noise as follows

$\begin{split}s_{a}(f)=z_{a}(f)+n_{a}(f).\end{split}$ (18)

We consider a simplified noise model in which the intrinsic pulsar red noise
is omitted, and $n_{a}$ is assumed to be an uncorrelated white Gaussian noise,
characterized by the power spectra

$\langle
n^{*}_{a}(f)n_{b}(f^{\prime})\rangle=\frac{1}{2}\delta_{ab}\delta(f-f^{\prime})N_{a}(|f|),$
(19)

where $N_{a}(f)=2\sigma^{2}_{a}\Delta t$ is the one-sided noise power spectrum
of the $a$-th pulsar, $\sigma^{2}_{a}$ is its noise variance, and $\Delta t$
is the sampling period.

The intensity and circular polarization contributions to the background can be
identified through the measurement of the real and imaginary parts of the
cross-correlation between timing residuals. The estimators for the two
components can therefore be written as

$\begin{split}\hat{I}_{ab}=&\frac{1}{2}\int^{\infty}_{-\infty}df\int^{\infty}_{-\infty}df^{\prime}\delta_{T}(f-f^{\prime})\big{[}s^{*}_{a}(f)s_{b}(f^{\prime})\\\
&+s_{a}(f)s^{*}_{b}(f^{\prime})\big{]}Q^{I}_{ab}(f^{\prime}),\\\
\hat{V}_{ab}=&\frac{1}{2i}\int^{\infty}_{-\infty}df\int^{\infty}_{-\infty}df^{\prime}\delta_{T}(f-f^{\prime})\big{[}s^{*}_{a}(f)s_{b}(f^{\prime})\\\
&-s_{a}(f)s^{*}_{b}(f^{\prime})\big{]}Q^{V}_{ab}(f^{\prime}),\end{split}$ (20)

where $\delta_{T}(f)=\sin(\pi fT)/\pi f$, and $Q_{I}$ and $Q_{V}$ are filters
to be determined. We will focus our discussion on the circular polarization
estimator, but note that the optimal estimator for the intensity can be
derived in an identical manner.

Our goal is to derive an expression for the filters that maximize the signal-
to-noise ratio, and combine the measurement from each pulsar pair in an
optimal way. Due to the frequency integration, we require the filter for the
intensity estimator to satisfy $Q_{I}(f)=Q_{I}(-f)$ and for the circular
polarization filter to satisfy $-Q_{V}(f)=Q_{V}(-f)$.

The mean of the circular polarization estimator is given by

$\begin{split}\langle\hat{V}_{ab}\rangle=&T\int_{-\infty}^{\infty}df\
V(f){}^{(ab)}\Gamma^{V}Q^{V}_{ab}(f),\end{split}$ (21)

and its covariance is shown in App. B. In order to shorten the notation, we
write the variance for a pulsar pair given in Eq. 37 as

$\sigma^{2}_{ab}=\frac{T}{2}\int_{-\infty}^{\infty}df\left[Q^{V}_{ab}(f)\right]^{2}\mathcal{C}_{ab}(f),$
(22)

where $\mathcal{C}_{ab}\equiv\mathcal{C}_{ab,ab}$.

The optimal filter can be found by defining the inner product in the space of
complex valued functions 1997rggr.conf..373A

$(A,B)=\int_{-\infty}^{\infty}dfA^{*}(f)B(f)\mathcal{C}_{ab}(f).$ (23)

By writing $\langle\hat{V}_{ab}\rangle$ and $\sigma^{2}_{ab}$ as an inner
product, it can be shown using the Schwarz inequality that the optimal filter
is given by

$Q_{ab}(f)=\chi\frac{V(f){}^{(ab)}\Gamma^{V}}{\mathcal{C}_{ab}(f)},$ (24)

where the constant $\chi$ is chosen so at to set the mean of the circular
polarization estimator. The measurement from each pulsar pair, accounting for
correlations, can then be optimally combined as follows (see, e.g. Ref
Romano:2020sxq )

$\hat{V}_{\text{opt}}=\frac{\sum\limits_{a}\sum\limits_{b<a}\lambda_{ab}\hat{V}_{ab}}{\sum\limits_{a}\sum\limits_{b<a}\lambda_{ab}},$
(25)

where

$\lambda_{ab}=\sum\limits_{c}\sum\limits_{d<c}(C^{-1})_{ab,cd}$ (26)

and the variance of the optimal estimator is given by

$\sigma^{2}_{\text{opt}}=\frac{1}{\sum\limits_{a}\sum\limits_{b<a}\lambda_{ab}}.$
(27)

##  V Case Study: Detecting a circular polarization dipole

To showcase the proposed circular polarization estimator, we consider a simple
scenario in which the gravitational wave background is dominated by a monopole
and a dipole that points in the $\hat{z}$ direction. The background is
therefore assumed to be

$X(f,\hat{\Omega})=X(f)\left(c^{X}_{00}Y_{00}+c^{X}_{10}Y_{10}\right)$ (28)

in the cosmic rest-frame. We note that the cross-correlation of pulsars is
insensitive to an isotropic circularly polarized background, or, in other
words, ${}^{(ab)}\Gamma^{V}_{00}=0$. The inclusion of the $\ell,m=0$ term in
the $V$ component is therefore superfluous. We assume the fiducial value of
the dipole coefficient to be 15% of the monopole for both $I$ and $V$.

We further assume that the frequency spectrum for both intensity and circular
polarization can be written as

$X(f)=\frac{1}{16\pi}A_{\text{X}}^{2}\left(\frac{f}{f_{1\text{yr}}}\right)^{2\alpha_{X}-1}.$
(29)

We choose the fiducial value of $A_{I}=A_{V}=10^{-15}$ and
$\alpha_{I}=\alpha_{V}=-2/3$. We focus on estimating the amplitude of the
circular polarization of the SGWB. In the dipole toy-model we consider here,
the coefficient $c_{10}$ is degenerate with the amplitude of the GW spectrum.
Hence, we choose to normalize the circular polarization estimator such that
$\langle\hat{V}_{ab}\rangle=c^{V}_{10}A^{2}_{V}$. The constant $\chi$ is
therefore given by

$\chi=\frac{c^{V}_{10}A^{2}_{V}}{T}\left(\int_{-\infty}^{\infty}df\frac{\left(V(f)^{(ab)}\Gamma^{V}_{10}\right)^{2}}{\mathcal{C}_{ab}}\right)^{-1}.$
(30)

We generate random pulsar positions uniformly distributed on the sky and
assume that all pulsars have identical white noise. The integral over
frequencies shown throughout this work are, in practice, taken between minimum
and maximum frequencies $f_{\text{min}}$ and $f_{\text{max}}$. The frequency
range is determined by the total observing time $T$ and the cadence time
$\Delta t$ (the time between consecutive pulsar observations). We assume an
observing time of 10 yr, which corresponds to $f_{\text{min}}\sim 3\times
10^{-9}$Hz, and a cadence of 2 weeks, resulting in $f_{\text{max}}\sim 8\times
10^{-7}$Hz.

We obtain the signal-to-noise ratio of the circular polarization amplitude in
two different ways. First, analytically, by computing the covariance matrix
given in Eq. 37 and the variance of the optimal estimator using Eq. 27. We
then validate our results on simulated timing residuals that we generate using
the method described in App. C. The estimator defined in Eq. 20 is then
applied to the simulated data, and the signal and noise are given by the mean
and standard deviation of the optimal estimator across various realizations of
the simulated PTAs. We confirmed that our estimator is unbiased and that the
variance matches our predicted value. We therefore choose to present only the
predicted signal-to-noise ratio in the results shown in this work.

The signal-to-noise ratio of the circular polarization estimator is shown as a
function of noise RMS in Fig. 2, for a network of 10, 25, and 50 pulsars. We
can see three distinct regimes: weak-signal (noise-dominated), intermediate-
signal, and high-signal (SGWB-dominated). The high-signal regime seen here
corresponds to a “cosmic variance” limit of the circular polarization
measurements, which results in a $S/N\sim 3$ for a PTA with 50 pulsars. The
second data release by the IPTA Perera:2019sca  includes 65 pulsars with
values of RMS between $\sim 0.2-14\mu$s, roughly corresponding to the
intermediate-noise portion of Fig. 2 if such a dipolar background were
present.

The results in Fig. 2 correspond to a single realization of pulsars on the
sky. However, the expected signal-to-noise ratio depends on the positions of
the pulsars relative to the SGWB dipole. A different realization of the pulsar
positions would therefore lead to a different signal-to-noise estimate. We
show in Fig. 3 the distributions of the signal-to-noise ratios for different
numbers of pulsars in a PTA. We consider a measurement in the high-signal
regime and assume all pulsars have an RMS of 1ns.

As expected, the distribution is the widest for the smallest PTA considered,
with a median value around $S/N\sim 1$. As we increase the number of pulsars,
the distribution narrows and reaches a signal-to-noise ratio of nearly
$S/N\sim 3$ for the largest array considered. In practice, the distribution of
pulsars on the sky is not uniform and is clustered around the galactic plane.
While the assumption of a uniform spatial distribution of pulsars
overestimates the signal-to-noise ratio of the SGWB, Fig. 3 captures a
sensible range of values for the dipole toy-model of GW anisotropy.

Figure 3: Signal-to-noise ratio of the amplitude of a circular polSarization
dipole as a function of the number of pulsars. The distributions correspond to
the intrinsic variation in $S/N$ due to different realizations of the pulsar
spatial positions. The pulsars were generated assuming a uniform distribution
of angles on the sky.

##  VI Conclusions

Over the last decade, PTA consortia have been gathering increasingly precise
timing data and are beginning to dip into an astrophysically interesting
region of the SGWB parameter space. The prospect of detection in the near
future and the expansion of both telescope-centered programs Bailes:2018azh ;
Ng:2017djg  and new PTAs Susobhanan:2020zmm  motivates us to consider what are
the next targets after an isotropic background is measured.

In this work, we address the circular polarization of the SGWB. By isolating
the real and imaginary parts of the timing residual correlation, we construct
an estimator that can naturally distinguish between the intensity and circular
polarization contributions of the SGWB. The estimator is based on extracting
the real and imaginary parts of the cross-correlation between pulsar timing
residuals. We compute the optimal filter and the combination of pulsar pairs
that maximizes the signal-to-noise ratio of the circular polarization
component.

As a case study, we consider a simple toy-model in which the SGWB is dominated
by a dipole that has 10% of the power of the monopole and that points in the
$\hat{z}$ direction. We then compute the expected signal-to-noise ratio of the
circular polarization amplitude under different assumptions about the number
of pulsars and their white noise. In order to validate our results, we also
apply the proposed estimator on simulated timing residuals and recover the
amplitude of the supposed circular polarization dipole.

By varying the pulsar noise RMS, we show where the three (low-signal,
intermediate-signal, and high-signal) regimes fall in the dipole toy-model.
Under idealized observational assumptions, we show that the dipole can be
recovered with a significance $\sim 3\sigma$ in an array that contains 50
pulsars in the SGWB-dominated high-signal regime. Finally, we show how the
sensitivity of the PTA to the circular polarization amplitude varies with the
locations of the pulsars on the sky. We compute the distribution of signal-to-
noise ratios across multiple realizations for increasing numbers of pulsars.

The circular polarization of the SGWB can offer additional information about
the astrophysical sources of the background and of more exotic parity-
violating physics. In particular, it may be an important tool to distinguish
between scenarios in which the background is dominated by multiple distant
sources and a smaller number of nearby sources. This contribution is naturally
present in the data set of PTAs and therefore extracting the circular
polarization component may be an important step in maximizing the information
recovered from timing data.

###### Acknowledgements.

We thank Tristan L. Smith, Andrea Lommen, Chiara Mingarelli, and José Luis
Bernal for useful discussions. GSP was supported by the National Science
Foundation Graduate Research Fellowship under Grant No. DGE1746891. This work
was supported at Johns Hopkins by NSF Grant No. 1818899 and the Simons
Foundation.

##  Appendix A Overlap Functions

In order to compute the ORFs explicitly, we must specify a coordinate system
for the pulsar positions and the gravitational wave propagation vector.
Following the formalism of Ref. Allen:1996gp , we work in the computational
frame, in which one pulsar is defined to be in the $\hat{z}$ direction and the
second pulsar in the $x-z$ plane. That is, we work in the following coordinate
system:

$\begin{split}\hat{p}&=(0,0,1),\\\ \hat{q}&=(\sin\zeta,0,\cos\zeta),\\\
\hat{\Omega}&=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta),\\\
\hat{m}&=(\sin\phi,-\cos\phi,0),\\\
\hat{n}&=(\cos\theta\cos\phi,\cos\theta\sin\phi,-\sin\theta),\end{split}$ (31)

where $\zeta$ is the angular separation between pulsars $a$ and $b$, and the
angles $\theta$ and $\phi$ define the GW propagation direction. In this
coordinate system, it can then be shown that the antenna beam functions are
given by Mingarelli:2013dsa

$\begin{split}F^{+}_{\hat{p}}(\hat{\Omega})&=-\frac{1}{2}\left(1-\cos\theta\right),\\\
F^{\times}_{\hat{p}}(\hat{\Omega})&=0,\\\
F^{+}_{\hat{q}}(\hat{\Omega})&=\frac{1}{2}\frac{(\sin\phi\sin\zeta)^{2}-\left(\sin\zeta\cos\theta\cos\phi-\sin\theta\cos\zeta\right)^{2}}{1+\cos\cos\zeta+\sin\theta\sin\zeta\cos\phi},\\\
F^{\times}_{\hat{q}}(\hat{\Omega})&=\frac{\left(\sin\phi\sin\zeta\right)\left(\cos\theta\sin\zeta\cos\phi-\sin\theta\cos\zeta\right)}{1+\cos\zeta+\sin\theta\sin\zeta\cos\phi}.\end{split}$
(32)

From the beam functions above, Ref.Mingarelli:2013dsa  computed the ORFs for
intensity, which we quote below up to the dipole terms

$\begin{split}{}^{(ab)}\Gamma^{I}_{00}&=\frac{\sqrt{\pi}}{2}\left[1+\frac{\cos\zeta}{3}+4\left(1-\cos\zeta\right)\log\left(\sin\frac{\zeta}{2}\right)\right],\\\
{}^{(ab)}\Gamma^{I}_{1-1}&=-\frac{1}{2}\sqrt{\frac{\pi}{6}}\left\\{(1+\cos\zeta)+3(1-\cos\zeta)\left[1+\frac{4}{(1+\cos\zeta)}\log\left(\sin\frac{\zeta}{2}\right)\right]\right\\},\\\
{}^{(ab)}\Gamma^{I}_{10}&=-\frac{1}{2}\sqrt{\frac{\pi}{3}}\left\\{1+\cos\zeta+3\left(1-\cos\zeta\right)\left[1+\cos\zeta+4\log\left(\sin\frac{\zeta}{2}\right)\right]\right\\},\\\
{}^{(ab)}\Gamma^{I}_{11}&=-{}^{(ab)}\Gamma^{I}_{1-1}.\end{split}$ (33)

Similarly, the ORFs for circular polarization were computed in Ref.
Kato:2015bye , which we quote up to dipole terms as well:

$\begin{split}{}^{(ab)}\Gamma^{V}_{00}&=0,\\\
{}^{(ab)}\Gamma^{V}_{1-1}&=-\sqrt{\frac{2\pi}{3}}\sin\zeta\left\\{1+3\left(\frac{1-\cos\zeta}{1+\cos\zeta}\right)\log\left(\sin\frac{\zeta}{2}\right)\right\\},\\\
{}^{(ab)}\Gamma^{V}_{10}&=0,\\\
{}^{(ab)}\Gamma^{V}_{11}&={}^{(ab)}\Gamma^{V}_{1-1}.\end{split}$ (34)

The ORFs computed in the computational frame must then be rotated back into
the cosmic-rest frame.

##  Appendix B Covariance

Here we compute the covariance between different pulsar pairs of the intensity
and circular polarization estimators. The estimators in Eq. 20 are defined for
each pulsar pair, which we label $(ab)$ and $(cd)$. The covariance for the
circular polarization estimator is proportional to

$\begin{split}C_{ab,cd}\propto&\langle\left(s^{*}_{a}s_{b}-s_{a}s^{*}_{b}\right)\left(s^{*}_{c}s_{d}-s_{c}s^{*}_{d}\right)\rangle\\\
&-\langle s^{*}_{a}s_{b}-s_{a}s^{*}_{b}\rangle\langle
s^{*}_{c}s_{d}-s_{c}s^{*}_{d}\rangle.\end{split}$ (35)

For the sake of conciseness, we have denoted $s(f,\hat{p}_{a})$ as $s_{a}$.
While the frequency dependence has been omitted, we note that each timing
residual signal corresponds to a certain frequency, which must then be
integrated over.

Under the assumption that the signal has a Gaussian distribution, we can use
Wick’s theorem to expand the expression above. Under the assumption that the
only non-zero Stokes parameters are I and V, we are left with

$\begin{split}C_{ab,cd}=&-\frac{1}{4}\int df_{a}\int df_{b}\int df_{c}\int
df_{d}\ Q^{V}_{ab}(f_{b})Q^{V}_{cd}(f_{d})\times\\\ &\big{[}\langle
s^{*}_{a}s_{d}\rangle\langle s_{b}s^{*}_{c}\rangle+\langle
s_{a}s^{*}_{d}\rangle\langle s^{*}_{b}s_{c}\rangle-\langle
s^{*}_{a}s_{c}\rangle\langle s_{b}s^{*}_{d}\rangle\\\ &-\langle
s_{a}s^{*}_{c}\rangle\langle s^{*}_{b}s_{d}\rangle\big{]}.\end{split}$ (36)

We can now substitute Eq. 18 for the timing residual and take the
correlations, with Eq. 17 describing the GW signal and Eq. 19 describing the
noise. We find that

$\begin{split}C_{ab,cd}=&\frac{T}{2}\int df\
Q^{V}_{ab}(f)Q^{V}_{cd}(f)\big{\\{}V^{2}(f)\left(\Gamma^{V}_{ac}\Gamma^{V}_{bd}-\Gamma^{V}_{ad}\Gamma^{V}_{bc}\right)\\\
&+\left(I(f)\Gamma^{I}_{ac}+\delta_{ac}N_{a}\right)\left(I(f)\Gamma^{I}_{bd}+\delta_{bd}N_{b}\right)\\\
&-\left(I(f)\Gamma^{I}_{ad}+\delta_{ad}N_{a}\right)\left(I(f)\Gamma^{I}_{bc}+\delta_{bc}N_{b}\right)\big{\\}},\end{split}$
(37)

where the pulsar pairs $(ab)$ and $(cd)$ may be 4 distinct pulsars, or may
have one or both pulsars in common. For later convenience, we will call the
term inside the brackets $\mathcal{C}_{ab,cd}$, hence

$C_{ab,cd}\equiv\frac{T}{2}\int df\
Q^{V}_{ab}(f)Q^{V}_{cd}(f)\mathcal{C}_{ab,cd}.$ (38)

##  Appendix C Simulated PTA

As a proof-of-concept, we apply the estimator proposed above to simulated PTA
data 2015PhRvD..91d4048C . Given a set of $N_{p}$ pulsars, we begin by
defining the correlation matrix ${\bf M}$ using Eq. 17, as

$\text{M}_{ab}(f)=\frac{T}{2}\left[I(f)^{(ab)}\Gamma^{I}+iV(f)^{(ab)}\Gamma^{V}\right],$
(39)

where ${\bf M}$ is therefore a $N_{p}\times N_{p}$ Hermitian matrix that is a
function of frequency. Applying the Cholesky decomposition, we can find the
matrix ${\bf H}(f)$ such that

${\bf H}{\bf H}^{\text{T}}={\bf M}$ (40)

for each frequency. The timing residual of a pulsar $a$ can then be written as

$z_{a}(f)=\sum_{b}H_{ab}(f)w_{b}(f),$ (41)

where $w_{a}(f)=x_{a}(f)+iy_{a}(f)$ is random complex number and both $x_{a}$
and $y_{a}$ have zero mean and unit variance. The random variables $w_{a}(f)$
therefore satisfy

$\langle
w^{*}_{a}(f)w_{b}(f^{\prime})\rangle=\frac{2}{T}\delta_{ab}\delta(f-f^{\prime}).$
(42)

We additionally include a source of uncorrelated white Gaussian noise, such
that

$\langle
n^{*}_{a}(f)n_{b}(f^{\prime})\rangle=\frac{1}{2}\delta_{ab}\delta(f-f^{\prime})N_{a},$
(43)

where $N_{a}$ is the noise spectral density. The pulsars simulated in this
manner then have the desired cross-correlation

$\langle
z^{*}_{a}(f)z_{b}(f^{\prime})\rangle=\left[I(f)^{(ab)}\Gamma^{I}+iV(f)^{(ab)}\Gamma^{V}\right].$
(44)

## References

  * (1) M. V. Sazhin, “Opportunities for detecting ultralong gravitational waves,” Sov. Astron.  22 (Feb., 1978) 36–38.
  * (2) S. Detweiler, “Pulsar timing measurements and the search for gravitational waves,” Astrophys. J.  234 (Dec., 1979) 1100–1104.
  * (3) M. Rajagopal and R. W. Romani, “Ultralow frequency gravitational radiation from massive black hole binaries,” Astrophys. J. 446 (1995) 543–549, arXiv:astro-ph/9412038.
  * (4) A. H. Jaffe and D. C. Backer, “Gravitational waves probe the coalescence rate of massive black hole binaries,” Astrophys. J. 583 (2003) 616–631, arXiv:astro-ph/0210148.
  * (5) S. Olmez, V. Mandic, and X. Siemens, “Gravitational-Wave Stochastic Background from Kinks and Cusps on Cosmic Strings,” Phys. Rev. D 81 (2010) 104028, arXiv:1004.0890 [astro-ph.CO].
  * (6) L. Sousa and P. P. Avelino, “Stochastic Gravitational Wave Background generated by Cosmic String Networks: Velocity-Dependent One-Scale model versus Scale-Invariant Evolution,” Phys. Rev. D 88 no. 2, (2013) 023516, arXiv:1304.2445 [astro-ph.CO].
  * (7) K. Miyamoto and K. Nakayama, “Cosmological and astrophysical constraints on superconducting cosmic strings,” JCAP 07 (2013) 012, arXiv:1212.6687 [astro-ph.CO].
  * (8) S. Kuroyanagi, K. Miyamoto, T. Sekiguchi, K. Takahashi, and J. Silk, “Forecast constraints on cosmic strings from future CMB, pulsar timing and gravitational wave direct detection experiments,” Phys. Rev. D 87 no. 2, (2013) 023522, arXiv:1210.2829 [astro-ph.CO]. [Erratum: Phys.Rev.D 87, 069903 (2013)].
  * (9) C. Caprini, R. Durrer, and X. Siemens, “Detection of gravitational waves from the QCD phase transition with pulsar timing arrays,” Phys. Rev. D 82 (2010) 063511, arXiv:1007.1218 [astro-ph.CO].
  * (10) A. A. Starobinsky, “Spectrum of relict gravitational radiation and the early state of the universe,” JETP Lett. 30 (1979) 682–685.
  * (11) W. Zhao, Y. Zhang, X.-P. You, and Z.-H. Zhu, “Constraints of relic gravitational waves by pulsar timing arrays: Forecasts for the FAST and SKA projects,” Phys. Rev. D 87 no. 12, (2013) 124012, arXiv:1303.6718 [astro-ph.CO].
  * (12) G. Desvignes et al., “High-precision timing of 42 millisecond pulsars with the European Pulsar Timing Array,” Mon. Not. Roy. Astron. Soc. 458 no. 3, (2016) 3341–3380, arXiv:1602.08511 [astro-ph.HE].
  * (13) A. Brazier et al., “The NANOGrav Program for Gravitational Waves and Fundamental Physics,” arXiv:1908.05356 [astro-ph.IM].
  * (14) M. Kerr et al., “The Parkes Pulsar Timing Array project: second data release,” Publ. Astron. Soc. Austral. 37 (2020) e020, arXiv:2003.09780 [astro-ph.IM].
  * (15) B. B. P. Perera et al., “The International Pulsar Timing Array: Second data release,” Mon. Not. Roy. Astron. Soc. 490 no. 4, (2019) 4666–4687, arXiv:1909.04534 [astro-ph.HE].
  * (16) S. R. Taylor, J. Simon, and L. Sampson, “Constraints On The Dynamical Environments Of Supermassive Black-hole Binaries Using Pulsar-timing Arrays,” Phys. Rev. Lett. 118 no. 18, (2017) 181102, arXiv:1612.02817 [astro-ph.GA].
  * (17) S. Chen, A. Sesana, and C. J. Conselice, “Constraining astrophysical observables of Galaxy and Supermassive Black Hole Binary Mergers using Pulsar Timing Arrays,” Mon. Not. Roy. Astron. Soc. 488 no. 1, (2019) 401–418, arXiv:1810.04184 [astro-ph.GA].
  * (18) R. W. Hellings and G. S. Downs, “Upper limits on the isotropic gravitational radiation background from pulsar timing analysis.,” Astrophys. J. Lett. 265 (Feb., 1983) L39–L42.
  * (19) V. Ravi, J. S. B. Wyithe, G. Hobbs, R. M. Shannon, R. N. Manchester, D. R. B. Yardley, and M. J. Keith, “Does a ’stochastic’ background of gravitational waves exist in the pulsar timing band?,” Astrophys. J. 761 (2012) 84, arXiv:1210.3854 [astro-ph.CO].
  * (20) N. J. Cornish and A. Sesana, “Pulsar Timing Array Analysis for Black Hole Backgrounds,” Class. Quant. Grav. 30 (2013) 224005, arXiv:1305.0326 [gr-qc].
  * (21) A. Sesana, A. Vecchio, and C. N. Colacino, “The stochastic gravitational-wave background from massive black hole binary systems: implications for observations with Pulsar Timing Arrays,” Mon. Not. Roy. Astron. Soc. 390 (2008) 192, arXiv:0804.4476 [astro-ph].
  * (22) N. Seto, “Quest for circular polarization of gravitational wave background and orbits of laser interferometers in space,” Phys. Rev. D 75 (2007) 061302, arXiv:astro-ph/0609633.
  * (23) N. Seto, “Prospects for direct detection of circular polarization of gravitational-wave background,” Phys. Rev. Lett. 97 (2006) 151101, arXiv:astro-ph/0609504.
  * (24) N. Seto and A. Taruya, “Measuring a Parity Violation Signature in the Early Universe via Ground-based Laser Interferometers,” Phys. Rev. Lett. 99 (2007) 121101, arXiv:0707.0535 [astro-ph].
  * (25) N. Seto and A. Taruya, “Polarization analysis of gravitational-wave backgrounds from the correlation signals of ground-based interferometers: Measuring a circular-polarization mode,” Phys. Rev. D 77 (2008) 103001, arXiv:0801.4185 [astro-ph].
  * (26) W. Qin, K. K. Boddy, M. Kamionkowski, and L. Dai, “Pulsar-timing arrays, astrometry, and gravitational waves,” Phys. Rev. D 99 no. 6, (2019) 063002, arXiv:1810.02369 [astro-ph.CO].
  * (27) R. Kato and J. Soda, “Probing circular polarization in stochastic gravitational wave background with pulsar timing arrays,” Phys. Rev. D 93 no. 6, (2016) 062003, arXiv:1512.09139 [gr-qc].
  * (28) E. Belgacem and M. Kamionkowski, “Chirality of the gravitational-wave background and pulsar-timing arrays,” Phys. Rev. D 102 no. 2, (2020) 023004, arXiv:2004.05480 [astro-ph.CO].
  * (29) S. C. Hotinli, M. Kamionkowski, and A. H. Jaffe, “The search for anisotropy in the gravitational-wave background with pulsar-timing arrays,” Open J. Astrophys. 2 no. 1, (2019) 8, arXiv:1904.05348 [astro-ph.CO].
  * (30) M. Anholm, S. Ballmer, J. D. E. Creighton, L. R. Price, and X. Siemens, “Optimal strategies for gravitational wave stochastic background searches in pulsar timing data,” Phys. Rev. D 79 (2009) 084030, arXiv:0809.0701 [gr-qc].
  * (31) B. Allen and J. D. Romano, “Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities,” Phys. Rev. D 59 (1999) 102001, arXiv:gr-qc/9710117.
  * (32) C. M. F. Mingarelli and A. B. Mingarelli, “Proving the short-wavelength approximation in Pulsar Timing Array gravitational-wave background searches,” Journal of Physics Communications 2 no. 10, (Oct., 2018) 105002, arXiv:1806.06979 [astro-ph.IM].
  * (33) C. M. F. Mingarelli and T. Sidery, “Effect of small interpulsar distances in stochastic gravitational wave background searches with pulsar timing arrays,” Phys. Rev. D 90 no. 6, (2014) 062011, arXiv:1408.6840 [astro-ph.HE].
  * (34) J. D. Romano, J. S. Hazboun, X. Siemens, and A. M. Archibald, “Common-spectrum process versus cross-correlation for gravitational-wave searches using pulsar timing arrays,” Phys. Rev. D 103 no. 6, (2021) 063027, arXiv:2012.03804 [gr-qc].
  * (35) B. Allen, “The Stochastic Gravity-Wave Background: Sources and Detection,” in Relativistic Gravitation and Gravitational Radiation, J.-A. Marck and J.-P. Lasota, eds., p. 373.  Jan., 1997.  arXiv:gr-qc/9604033 [gr-qc].
  * (36) M. Bailes et al., “MeerTime - the MeerKAT Key Science Program on Pulsar Timing,” PoS MeerKAT2016 (2018) 011, arXiv:1803.07424 [astro-ph.IM].
  * (37) CHIME Pulsar Collaboration, C. Ng, “Pulsar science with the CHIME telescope,” IAU Symp. 337 (2017) 179–182, arXiv:1711.02104 [astro-ph.IM].
  * (38) A. Susobhanan et al., “pinta: The uGMRT Data Processing Pipeline for the Indian Pulsar Timing Array,” Publ. Astron. Soc. Austral. 38 (2021) E017, arXiv:2007.02930 [astro-ph.IM].
  * (39) B. Allen and A. C. Ottewill, “Detection of anisotropies in the gravitational wave stochastic background,” Phys. Rev. D 56 (1997) 545–563, arXiv:gr-qc/9607068.
  * (40) C. M. F. Mingarelli, T. Sidery, I. Mandel, and A. Vecchio, “Characterizing gravitational wave stochastic background anisotropy with pulsar timing arrays,” Phys. Rev. D 88 no. 6, (2013) 062005, arXiv:1306.5394 [astro-ph.HE].
  * (41) S. J. Chamberlin, J. D. E. Creighton, X. Siemens, P. Demorest, J. Ellis, L. R. Price, and J. D. Romano, “Time-domain implementation of the optimal cross-correlation statistic for stochastic gravitational-wave background searches in pulsar timing data,” Phys. Rev. D 91 no. 4, (Feb., 2015) 044048, arXiv:1410.8256 [astro-ph.IM].