###### Abstract

The Galactic magnetic field is an integral constituent of the interstellar
medium (ISM), and knowledge of its structure is crucial to understanding
Galactic dynamics. The Rotation Measures (RM) of extragalactic (EG) sources
have been the basis of comprehensive Galactic magnetic field models. Polarised
extended emission (XE) is also seen along lines of sight through the Galactic
disk, and also displays the effects of Faraday rotation. Our aim is to
investigate and understand the relationship between EG and XE RMs near the
Galactic plane, and to determine how the XE RMs, a hitherto unused resource,
can be used as a probe of the large-scale Galactic magnetic field. We used
polarisation data from the Canadian Galactic Plane Survey (CGPS), observed
near 1420 MHz with the Dominion Radio Astrophysical Observatory (DRAO)
Synthesis Telescope. We calculated RMs from a linear fit to the polarisation
angles as a function of wavelength squared in four frequency channels, for
both the EG sources and the XE. Across the CGPS area,
$55^{\circ}<{\ell}<193^{\circ},-3^{\circ}<b<5^{\circ}$, the RMs of the XE
closely track the RMs of the EG sources, with XE RMs about half the value of
EG-source RMs. The exceptions are places where large local HII complexes
heavily depolarise more distant emission. We conclude that there is valuable
information in the XE RM dataset. The factor of 2 between the two types of RM
values is close to that expected from a Burn slab model of the ISM. This
result indicates that, at least in the outer Galaxy, the EG and XE sources are
likely probing similar depths, and that the Faraday rotating medium and the
synchrotron emitting medium have similar variation with galactocentric
distance.

###### keywords:

galaxy: structure; ISM: magnetic fields; polarization; radio continuum: ISM;
techniques: interferometric

xx 1 5 Received: 30 January 2019; Accepted: 22 March 2019 Faraday Rotation of
Extended Emission as a Probe of the Large-Scale Galactic Magnetic Field Anna
Ordog 1,∗, Rebecca A. Booth 1, Cameron L. Van Eck 1,2, Jo-Anne C. Brown 1,∗
and Thomas L. Landecker 3Anna Ordog, Rebecca A. Booth, Cameron L. Van Eck, Jo-
Anne C. Brown and Thomas L. Landecker Correspondence<EMAIL_ADDRESS>(A.O<EMAIL_ADDRESS>(J.C.B.)

##  1 Introduction

The idea that our Galaxy has a magnetic field was first proposed by Alfvén in
1937 Alfvén (1937). Just over ten years later, Fermi corroborated this
proposal, arguing that a Galactic magnetic field (GMF) would be able to
explain the origin and confinement of cosmic rays. This field, he explained,
would have coherence lengths on the order of thousands of parsecs, great
stability, and would be able to prevent cosmic rays from escaping from the
Galaxy Fermi (1949). The first evidence for a large-scale GMF was from
observations of starlight polarisation Hiltner (1949); Hall (1949).

The GMF is now considered to be a fundamental component of the interstellar
medium (ISM). For example, magnetic fields may play a role in star formation
van Loo (2012), and they contribute to the vertical support of gases in the
Galaxy against gravitational collapse Boulares (1990). Magnetic fields may
also influence Galactic evolution Kim (1996).

Most of what we know about the GMF has come from studying the Faraday rotation
of polarised emission from compact, unresolved objects, such as extragalactic
(EG) sources and pulsars Simard-Normandin (1979); Rand (1994); Brown (2003);
Taylor (2009). Modeling the GMF has proven to be challenging, despite various
approaches and use of a variety of datasets and assumptions Jaffe (2010);
VanEck (2011); Jansson (2012); Gressel (2013); Han (2018). Accounting for
features of the GMF, such as the apparent large-scale reversal between the
Local and Sagittarius arms Simard-Normandin (1979); VanEck (2011); Thomson
(1980) (for which no counterpart can be found in external galaxies), the
variation in magnetic spiral pitch angle with galactocentric radius VanEck
(2011), and field components of varying spatial scales and degrees of
coherence, adds complexity to the models being tested. The difficulty in
constraining these models has been compounded by the sparse nature of the EG
source data available.

Faraday rotation of polarised extended emission (XE) is another potential
source of information on the GMF structure. While weaker in polarised
intensity than the EG sources, XE occurs in nearly all directions, and its
Faraday rotation may be useful in conjunction with EG sources. Highly
polarised portions of the sky, for example, the Fan region, have been studied
in terms of Faraday rotation of the diffuse emission Iacobelli (2013), but
very little investigation has been done into utilizing this technique for
studying the large-scale GMF. In this paper we explore that possibility.

For compact sources, it is reasonable to assume that the polarised emission
originates from a single point in space (i.e., that there is no significant
internal Faraday rotation), and that the rotation occurs along a clearly
identified path through the Galaxy. In this case, we calculate the Rotation
Measure (RM) as the slope of a least-squares linear fit to the polarisation
angle ($\tau$) as a function of the square of the wavelength ($\lambda^{2}$),
so that

$\tau=\tau_{\circ}+RM\lambda^{2}.$ (1)

This is subtly different in the case of XE, for which synchrotron emission and
Faraday rotation occur within the same volume, both linked to the same
magnetic field. In this case, we use the term Faraday depth to quantify the
effects of Faraday rotation from a specific volume at a particular distance,
$d$, from the observer located at the origin. Faraday depth is defined by the
line integral:

$\phi=0.812\int_{d}^{0}n_{e}\,B_{\parallel}\,d\ell,$ (2)

where $n_{e}$ is the electron density, $B_{\parallel}$ is the parallel
component of the GMF, and $d\ell$ is the infinitesimal path length along the
line-of-sight (LOS), from the source to the observer. For a Faraday-thin
source without internal Faraday Rotation, as is typically assumed for a
compact EG source, RM and Faraday Depth are identical. In contrast, the
Faraday Depth profile of XE regions can be quite complicated Burn (1966);
Sokoloff (1998), and it is generally argued that wider coverage and higher
resolution in wavelength are required for studying Faraday Rotation of XE than
for compact sources Brentjens (2005).

In general, the RM of a polarised source provides information about the
strength and the direction of the average magnetic field component parallel to
the LOS, scaled by the electron density along the path. A positive [negative]
RM indicates that the average LOS component of the GMF points toward [away
from] the observer.

Recent discussions about the GMF, for example, by Ordog et al. Ordog (2017),
have expanded to use the RMs calculated from four frequency channels as a
proxy for Faraday Depth of XE. From here forward we speak of the “RM of XE”
with this specific meaning. The XE data used in Ordog (2017) were from the
Dominion Radio Astrophysical Observatory (DRAO) Synthesis Telescope without
the addition of single-antenna data, that is, from interferometers of
baselines 13 m to 600 m. It can be argued that XE data observed using radio
interferometry are an incomplete source of information, particularly on
polarisation angle and related quantities, due to a lack of sensitivity to
large-scale structures. However, Ordog et al. Ordog (2017) demonstrated a
remarkable similarity between the RMs seen in the EG and the XE over a large
area using this type of dataset. The reversal in magnetic field direction
between the Local and Sagittarius arms, between longitudes 55∘ and 70∘, can be
seen in the RM patterns of both the EG sources and the XE, occurring in the
same place.

In exploring the relationship between the RMs of EG sources seen through the
Galactic disk and the RMs or Faraday depths of XE, it would be ideal to have
high-resolution data that include information on a wider range of spatial
scales, achieved by combining aperture-synthesis and single-antenna data
(e.g., Landecker (2010)). However, we do not yet have single-antenna data with
multiple frequency channels that can be combined with the synthesis telescope
data. Thus we investigated the utility of the XE data from the DRAO Synthesis
Telescope alone. In particular, we report on our efforts to learn as much as
we can about the GMF by comparing the RMs of EG sources and XE from the
Canadian Galactic Plane Survey (CGPS), with the aim of developing methods for
combining information from the two types of datasets.

Our paper outline is as follows. In Section 2 we describe the datasets used.
In Section 3 we compare the EG and XE sources, commenting on the observed
trends. In Section 4 we discuss possible simple models of the ISM that could
account for the correlations we find between the two types of sources and
comment on anomalous regions. In Section 5 we summarize our results.

##  2 The Data

The CGPS is a radio survey centred at a frequency of 1420 MHz between Galactic
latitudes $-3^{\circ}\leq b\leq 5^{\circ}$ with a Galactic longitude range of
$52^{\circ}\leq\ell\leq 193^{\circ}$. The survey specifications are listed in
Table 1. From the polarisation data collected as part of the CGPS using the
Synthesis Telescope at the DRAO, we have calculated RMs for XE across the
entire dataset.

We derived RM values from the CGPS Stokes $Q$ and $U$ images in four bands,
each 7.5 MHz wide, centred at the frequencies indicated in Table 1 Landecker
(2010). From the Stokes $Q$ and $U$ images we calculated the polarisation
angle, $\tau$, for each frequency using $\tau=\frac{1}{2}\arctan U/Q$. The RM
for each pixel in the image was determined using Equation (1). By using this
algorithm, which is identical to that used in determining the EG source RMs
(Van Eck et al. 2019, in prep., based on the method described in Brown
(2003)), we are treating the two datasets in the same manner. This process
yielded $4.9\times 10^{7}$ RM values spanning the entire CGPS region, stored
in an 1825 $\times$ 29,025 pixel array on a 20 $\times$ 20 arcsecond grid.

Table 1: The CGPS polarisation survey.

Overall coverage | ${53^{\circ}}<{\ell}<{192^{\circ}},{-3^{\circ}}<{b}<{5^{\circ}}$
---|---
| ${101^{\circ}}<{\ell}<{116^{\circ}},{5.0^{\circ}}<{b}<{17^{\circ}}$
Continuum bandwidth | 30 MHz in four bands of
| 7.5 MHz each
Polarisation products | Stokes $I$, $Q$, and $U$
Centre frequencies | 1407.2, 1414.1, 1427.7,
| and 1434.6 MHz
Angular resolution |  $58^{\prime\prime}\times 58^{\prime\prime}$ cosec $\delta$
Sensitivity, $I$ | 200 to 400 $\upmu$Jy/beam rms
Sensitivity, $Q$ and $U$ | 180 to 260 $\upmu$Jy/beam rms
Typical noise in |
mosaicked images |  ${76}\thinspace{\rm{sin{\thinspace}}}\delta$ mK
* Sources of | Effelsberg 100-m Telescope
single-antenna data | and DRAO 26-m Telescope

* These data are part of the full-band CGPS survey, but are not available for the individual 7.5 MHz bands. Therefore we do not use these single-antenna data in our analysis.
---

###  2.1 Data Processing

In order to extract the XE information from the entire CGPS RM array and to
ensure reliable RM data, we examined the corresponding Stokes $I$ image, the
signal to noise (S:N) ratio of the linearly polarised intensity,
$P=\sqrt{Q^{2}+U^{2}}$, and the linear fits of angle as a function of
$\lambda^{2}$ for each pixel. The process by which we filtered the XE RM
dataset involved three steps as described below and illustrated in Figure 1.
111The final filtered XE RM map will be available at the CDS via
http://vizier.u-strasbg.fr/viz-bin/VizieR.

Figure 1: The process by which we filtered the XE RM data, illustrated for a
selected region of the CGPS RM image. (a) The raw RM values from the linear
fit to polarisation angle versus $\lambda^{2}$. (b) The corresponding Stokes
$I$ image used in the first processing step. (c) The RM map with pixels with
Stokes $I$ greater than 3 mJy/beam removed. (d) The signal-to-noise map, used
in the second processing step. (e) The RM map with pixels in polarised
intensity below the signal to noise threshold of 5 removed. (f) The
probability of linear fit used in the third processing step. (g) The RM map
with pixels with a probability of linear fit less than 10% removed.

#### Step 1: Remove pixels with high Stokes I

In addition to the extended synchrotron emission, in the polarisation images
we also see compact sources and discrete objects such as supernova remnants
and HII regions. While the HII regions do not emit polarised radiation, they
can affect the polarised background. Both the compact objects and discrete
objects tend to be brighter in Stokes $I$ than the surrounding XE, making it
relatively straightforward to remove them from the XE RM images. Removing
these objects is necessary, since we are interested in comparing the XE RMs to
the EG RMs. Figure 1a–c demonstrate an example of this stage of the filtering
for a small segment of the RM image.

In order to identify the RMs of XE only, we discarded all pixels from the RM
array with corresponding Stokes $I$ higher than 3 mJy/beam. We experimented
with a variety of thresholds and found that choosing a 3 mJy/beam threshold
successfully removed all EG sources as well as pixels from known bright
Galactic features from the RM data. After this 3 mJy/beam Stokes $I$ threshold
was applied, $4.6\times 10^{7}$ pixels (94% of the original number) remained
in the RM array.

#### Step 2: Removal of pixels with low signal-to-noise ratio

Adopting the practice used in a previous work on the CGPS RM data Brown
(2003), pixels with low S:N, for which the signal detected in polarised
intensity was less than five times the noise, were rejected from the XE RM
analysis, as shown in Figure 1c–e. We determined a S:N map based on the
sensitivity of the Synthesis Telescope (see Table 1) and the weight maps of
the images. The method is presented in detail in Brown (2003) (Appendices 1
and 2) and Brown (2002). We then used this S:N map as a filter for the XE RMs.
Following this step, $4.0\times 10^{6}$ RM pixels (8.2% of the original
number) remained in the array.

#### Step 3: Removal of pixels with fit probability < 10%

In order to ensure reliable RM values, we determined the quality of the linear
fit for each RM using a standard Pearson’s chi-squared test that compares the
calculated chi-squared values to a chi-squared distribution. We only accepted
pixels with reliable linear fit parameters, rejecting those with fit
probability less than 10%, as shown in Figure 1e–g. Following this filtering
step, $2.3\times 10^{6}$ RM pixels (4.7% of the original number) remained.
Although a significant number of the pixels (55% of the original) had a low
probability of fit, more than 75% of these were already discarded due to low
S:N.

The pixels with poor linear fit likely correspond to lines of sight having
increased Faraday complexity. With only four frequency channels, we are not
able to accurately extract useful information for those lines of sight (e.g.,
using RM synthesis Brentjens (2005)). Therefore, we do not include them in
this analysis.

The catalog of EG RMs for the CGPS has a total of 2112 sources in the Galactic
disk, resulting in close to 2 sources per square degree (Van Eck, et al.,
2019, in prep.). However, to ensure that the EG sources and the XE used in our
analysis were probing similar lines of sight, we did not include EG sources
located along lines of sight passing through Stokes $I$ regions brighter than
3 mJy/beam, since those regions were already eliminated from the XE RM map. We
identified such EG sources by their locations within the 3 mJy/beam threshold
mask. Using this criterion, 2084 EG sources remained.

##  3 Analysis

We compared RM trends as a function of longitude in the XE and EG RM datasets
by calculating the average for both within independent bins each spanning 1.5
degrees of longitude and covering the CGPS latitude range. Figure 2 shows an
example of the distribution of XE and EG RMs within one of the bins. There is
a high degree of scatter in the XE RM values, as can also be seen in the
small-scale variability in the map shown in Figure 1g. Nevertheless, the
distributions within the bins do have clearly defined peaks. The use of
sophisticated methods for producing a smoothed EG RM map, such as Oppermann
(2012), is beyond the scope of this paper. For the purpose of observing large-
scale trends as a function of longitude, we have opted to simply average the
EG RMs within each bin.

As described above, the filtered XE RM array contains $2.3\times 10^{6}$
pixels, 4.7% of the original number. Even with this small fraction of the XE
data remaining, the number of lines of sight probed by the XE still greatly
exceeds the number probed by EG sources 222The pixels are not independent,
because the synthesized beam of the telescope consists of between 9 and 24
pixels in the image, depending on the declination. However, even with this
taken into account, we still gain lines of sight in the XE that are missing in
the EG source list.. Figure 3 shows the number of XE RM pixels and the number
of EG sources in each 1.5 degree longitude bin. The 1.5 degree bin width was
selected to ensure that there was a statistically significant number of EG
sources contained in each bin. On average there are 22 EG sources in each bin
with a low of 5 sources/bin at $\ell=79.0^{\circ}$ and $\ell=80.5^{\circ}$.
This is in the Cygnus X region, a local, extended region of thermal emission
which depolarises emission produced behind it. We discuss this region further
in Section 4.

Figure 2: Left: The distribution of XE RM pixels after the filtering process
(dark gray) and EG RM sources (light gray) within one 1.5∘ wide longitude bin,
centred at $\ell=137.5^{\circ}$. Right: The distribution of XE polarisation
angles measured at 1420 MHz within the same longitude bin. While the XE
polarisation angles seem to be randomly distributed within the bin, the XE RM
values have a clear peak. Figure 3: The number of XE RM pixels (dotted line)
and EG sources (solid line) contained in independent bins each spanning 1.5
degrees of longitude and covering the CGPS latitudes, $b=-3^{\circ}$ to
$+5^{\circ}$.

The average number of XE RM pixels in each 1.5 degree bin is $2.3\times
10^{4}$. However, there is a decreasing trend towards the middle of the
longitude range, with only 359 pixels in the bin centred at
$\ell=116.5^{\circ}$. Three effects contribute to this reduced number of
pixels in central longitudes of the CGPS region:

  1. 1.

The intensity of synchrotron emission depends on the component of the magnetic
field perpendicular to the LOS , as well as the relativistic electron density
in the emitting region. The intensity peaks close to the anticentre, in the
Fan region (near $\ell=150^{\circ}$), and falls smoothly toward lower
longitudes.

  2. 2.

The synthesized beam of the telescope varies with declination, with a smaller
beam having less signal for the same polarised brightness temperature. The
beam is smallest at longitude about 120∘, corresponding to the highest
declination, where there are fewer pixels in the beam to survive the various
filters that we have applied.

  3. 3.

The shortest baseline used by the Synthesis Telescope is 12.9 m. At low
declinations foreshortening of the baseline extends the sensitivity of the
telescope to larger angular structures. Conversely, at high declinations the
sensitivity to extended structure is somewhat reduced. The highest
declinations correspond to roughly $\ell=120^{\circ}$. Although it would be
possible to make the sensitivity to extended structures uniform across the
dataset by tapering the visibilities of long baselines, doing so would not
yield significant gains, since it would correspond to an overall loss of high
quality data.

The top panel of Figure 4 is a plot of the averaged RM in 1.5 degree slices in
Galactic longitude for XE and EG, as function of Galactic longitude. The error
bars for the EG RM bins are the standard deviation of the mean within the
bins. The error bars for the XE RM bins were determined by dividing the
standard deviation of the XE RMs within each bin by the square root of the
average number of synthesized beams within the bin, to account for pixels not
being independent. While the average XE RM is smaller in magnitude than the
average EG RM, the most striking result is that the XE and EG RMs seem to be
following the same trend with longitude. The CGPS EG dataset has been shown to
be consistent with other published datasets such as Taylor (2009); VanEck
(2011), as demonstrated in Van Eck et al. 2019, in prep.

In order to evaluate the correlation between the binned XE and EG lists, we
determined the Pearson correlation coefficient and found a result of 0.53 (_p_
= 1$\times 10^{-5}$), indicative of a proportionality between the two types of
RM sources. We experimented with different bin widths, and with different
values of the parameters in the XE data filtering process, and found that the
results shown in Figure 4 were robust with respect to our choice of these
parameters.

Burn (1966) calculated the theoretical RM that would be observed for a LOS
through a slab with uniform electron density, magnetic field, and synchrotron
emission profiles Burn (1966). In this simplified model, the Faraday Depths
average along the LOS resulting in an observed RM for XE that is half that of
the emission from an EG source along the same LOS. Adopting the Burn slab
model as a starting point, we multiplied the XE RM average by two, as shown in
the bottom panel of Figure 4, and were struck by how well the two datasets
then line up. This result is surprising given that the magnetic field,
electron density and synchrotron emissivity are not constant throughout the
Galactic disk. Nonetheless, the agreement in overall trend and detail between
the average XE RM and EG RM curves is undeniable, especially in the portion of
the outer Galaxy spanned by the CGPS ($90^{\circ}\leq\ell\leq 193^{\circ}$).

To further examine the correlation between the XE and EG RMs, we also
calculated average XE RM values within circles centred around each EG source.
Varying the size of the circles, we calculated the resulting correlation
coefficients, and found that the correlation is maximized for a radius of
0.25∘, corresponding to an area of approximately 0.20 square degrees. The
result is shown in Figure 5. Using the size of the XE regions that maximizes
the correletion, we plot the average XE RMs as a function of the corresponding
EG RMs, as shown in Figure 6. We indicate the ratio of 0.5 between the two
types of sources by the dashed line. Although there is a significant degree of
scatter, the points are clustered around this line, with a root-mean-square
difference of 80 rad m-2. The highest degree of scatter occurs for sources at
low longitudes, the reasons for which are discussed in Section 4.3.

Figure 4: Top: The average EG RMs (solid line) and the average XE RMs (dotted
line) within independent bins each spanning 1.5 degrees of longitude and
covering the CGPS latitudes, $b=-3^{\circ}$ to $+5^{\circ}$. Bottom: The
average EG RMs (solid line) and the average XE RMs multiplied by two (dotted
line). Figure 5: The Pearson correlation coefficient between the XE RM values
averaged within circles of varying radii centred at each EG source, and the
corresponding EG RM value. The correlation is strongest for XE areas with a
radius of 0.25∘.

It is important to note that unlike in Figure 4 which utilizes independent
bins, the calculated XE RM values in Figure 6 overlap with some adjacent
points due to the proximity of their locations, which are defined by the
locations of the EG sources. This may, in part, account for the correlation
between XE and EG sources being highest for XE regions significantly smaller
than the independent 1.5${}^{\circ}\times 8^{\circ}$ bins, which result in an
equally robust correlation. Furthermore, due to the variability of the EG
source RMs as well (see Figure 2), averaging both types of sources (Figure 4)
rather than treating the EG sources independently (Figure 6) allows for large-
scale GMF structures to be discerned more easily. Regardless, Figure 6 does
highlight that the lines of sight corresponding to a significant number of EG
sources pass through regions of XE with average RMs of approximately half the
EG RM values.

Figure 6: The XE RM values averaged within circles with radius 0.25∘ centred
at each EG source, plotted as a function of the corresponding EG RM. The
colour scale indicates the locations of the sources in terms of Galactic
longitude, and the dashed line indicates a constant ratio of 0.5 between the
XE and EG RMs. The points are clustered close to this line, highlighting the
correlation between the two types of sources.

##  4 Discussion

Here we comment on the interpretation of the results described above. In
particular, we discuss the relationship between the polarisation angles and
the RM values, and present possible configurations of ISM parameters that may
lead to the observed patterns. Additionally, we comment on regions in this
dataset where the agreement between the XE and EG RMs breaks down.

###  4.1 Polarisation Angles in the XE Data

In addition to analyzing the RM values of the XE pixels, we investigated the
corresponding observed polarisation angles at 1420 MHz (not unwrapped to zero
wavelength). Applying the same filters to the map of polarisation angles as
were applied to the RM map, we found an approximately uniform distribution of
the XE angles within most of the longitude bins. The right panel of Figure 2
shows an example of a bin within which the XE and EG RMs have clear peaks, but
the polarisation angle distribution does not. Less than 25% of the bins
demonstrated a weak peak in the polarisation angle histograms. These bins will
require further investigation beyond the scope of this paper, but the lack of
a dominant angle is likely due to the missing single-antenna data that would
be required in order to detect large-scale structures in angle. Despite the
absence of a coherent pattern in the distribution of angles, the binned RMs
calculated from those angles vary smoothly and systematically with longitude,
indicating that the XE are definitely reliable probes of the large-scale
magnetic field in terms of RM values. The fact that small-scale fluctuations
in the polarisation angle at 1420 MHz do not have a significant effect on bin-
averaged RM values with our chosen bin size is a surprising and useful result.

###  4.2 ISM Configurations

Figure 4 show the strong resemblance between EG and XE RMs across a
significant range of longitudes. This similarity suggests that the two types
of sources probe the same volume of the ISM, and may put into question the
previous assumption that XE RM data probe only the very local ISM (e.g.,
Uyaniker (2003)). If XE and EG probe to a similar depth, this has important
implications for interpreting polarisation data in order to study the GMF. It
may mean that interpreting a Polarisation Horizon Uyaniker (2003) as a
boundary beyond which contributions to the observed polarised intensity are
insignificant, does not apply in the same way to the XE RMs. We investigate
this by considering possible ISM configurations that could give rise to the
observed relationship between the XE and EG RMs.

As described in the previous section, the ratio of EG source to XE source RMs
is close to 2 for a large portion of the CGPS. The simplest configuration for
the XE is the Burn slab model, for which the ratio of EG to XE RMs would
indeed be 2. In the Burn slab model, synchrotron emission and Faraday rotation
occur everywhere along the LOS through a region of uniform magnetic field,
thermal electron density and synchrotron emissivity. Since this is not
accurate for the Galaxy, we explore the extent to which some non-uniform
models of ISM parameters can resemble a slab-like scenario.

For the interstellar medium, the thermal electron density, magnetic field
strength and synchrotron emissivity can be reasonably modelled on large-scales
as exponentials with characteristic scale-lengths: $h_{n_{e}}$ (electron
density), $h_{B}$ (magnetic field) and $h_{s}$ (emissivity). In the simple
model we describe here, the origin of the coordinate system for the
exponential functions is located at the observer (i.e., the Earth) and not at
the Galactic centre. Although this limits the model to lines of sight close to
$\ell=180^{\circ}$, in this way it is possible to determine the Faraday depth
spectrum by means of a straightforward analytical calculation. The resulting
Faraday depth spectrum333We use the tilde to denote a complex phasor encoding
both the polarised intensity and polarisation angle. (derived in the Appendix)
is:

$\tilde{P}(\phi)=\frac{\tilde{P}_{\circ}}{0.812n_{\circ}B_{\circ}}\Big{(}1-\frac{\phi}{\phi_{max}}\Big{)}^{\frac{h_{\phi}}{h_{s}}-1},$
(3)

where $\tilde{P}_{\circ}$, $n_{\circ}$ and $B_{\circ}$ are the polarised
emission, electron density and magnetic field strength at the origin,
$h_{\phi}=(1/h_{n_{e}}+1/h_{B})^{-1}$ is the Faraday rotation scale-length,
and $\phi_{max}=0.812n_{\circ}B_{\circ}h_{\phi}$ is the maximum possible
Faraday depth in this configuration. We determine the Fourier transform
numerically, and calculate observed polarisation angles as a function of
wavelength squared, as shown in Figure 7. This derivation is similar to the
‘asymmetric slab’ scenario presented by Sokoloff et al. 1998 Sokoloff (1998).
Our approach offers new insight in that we formulate the equations in terms of
physical scale-lengths, calculate the Faraday depth spectra, and explicitly
highlight three different regimes: (i) $h_{\phi}=h_{s}$; (ii)
$h_{\phi}>h_{s}$; (iii) $h_{\phi}<h_{s}$.

Figure 7: An example of a polarisation model, comparing the ‘slab-like’ case
(blue) where $h_{\phi}$ and $h_{s}$ are equal ($h_{\phi}/h_{s}=1$), to the
case with $h_{\phi}/h_{s}=1.2$ (orange) and the case with $h_{\phi}/h_{s}=0.8$
(green). (a): the Faraday depth spectrum. (b): polarised intensity as a
function of wavelength squared. (c): polarisation angle as a function of
wavelength squared. (d): the slope of polarisation versus wavelength squared.
Please note that in panels (a) and (b) the units of polarised intensity are
arbitrary and the curves have been normalized to have the same integrated
intensity.

Regime (i): $h_{\phi}=h_{s}$
If Faraday rotation (magnetic field and thermal electrons) and synchrotron
emissivity (magnetic field and cosmic ray electrons) extend the same distance
out along the LOS, the polarised emission has equal contribution from all
Faraday depths out to $\phi_{max}$. In this case, Equation (3) simplifies to
the Burn slab Faraday depth spectrum. We can expect this result to hold for
any distribution of emissivity and Faraday rotation, not necessarily an
exponential, so long as the corresponding parameters are matched. We compare
the other two regimes to this slab-like regime. Henceforth, we use the term
“slab-like” to indicate an ISM configuration that results in the same Faraday
depth spectrum as the Burn-slab, but the parameters of which do not
necessarily have uniform profiles.

Regime (ii): $h_{\phi}>h_{s}$
If Faraday rotation extends further out along the LOS than synchrotron
emissivity, the polarised emission is concentrated at low Faraday depths
(nearby regions). Unlike for the Burn slab, the polarised intensity versus
$\lambda^{2}$ does not have null-points. The corresponding observed
polarisation angles versus $\lambda^{2}$ differ from the slab-like case in
that the slope alternates between positive and negative values. In the limit
as $h_{\phi}/h_{s}$ approaches 1, this approaches a slab-like case: the
Faraday depth spectrum flattens out, and the ranges of $\lambda^{2}$ over
which the RM has a reversed sign decrease in extent.

Regime (iii): $h_{\phi}<h_{s}$
When synchrotron emissivity extends beyond the Faraday rotating region, the
polarised emission is concentrated at higher Faraday depths (more distant
regions). The polarised intensity versus $\lambda^{2}$ deviates from a Burn
slab similarly to Regime (ii), but there is an enhancement of the slope of
polarisation angle versus $\lambda^{2}$ over some ranges of $\lambda^{2}$
rather than reversals in the sign. This leads to the surprising result that
the observed XE RMs can exceed the observed EG RMs for certain wavelengths.

Various models for the GMF have been tested that incorporate exponential
profiles for the cosmic ray electron density and the magnetic field strength.
Examples include the disk magnetic field model by Jaffe et al. Jaffe (2010),
in which the scale-lengths of both the magnetic field and the cosmic ray
electron density are 15 kpc, and the 3-dimensional models by Sun et al. Sun
(2008), in which magnetic field scale-lengths of 6 kpc, 8.5 kpc and 10 kpc are
tested, while the cosmic ray electron density scale-length is 8 kpc . Both of
these models use the NE2001 thermal electron density model Cordes (2002). A 5
kpc scale-length is a reasonable estimate for an exponential fit to the NE2001
model. With this approximation, an estimate for the ratio of $h_{\phi}/h_{s}$
would be 0.8 for the Jaffe et al. Jaffe (2010) model, and 1.1 for the Sun et
al. Sun (2008) model with a magnetic field scale-length of 8.5 kpc. These
values are sufficiently close to 1 to allow for a ratio between XE and EG RMs
of close to 0.5 over significant ranges of wavelengths.

Assuming Galacto-centric exponential models for the parameters discussed, the
spectrum we present is only strictly valid for the LOS directed toward
$\ell=180^{\circ}$. For lines of sight with $\ell<180^{\circ}$, a greater
portion of the LOS is dominated by the local ISM. For lines of sight with
$\ell<90^{\circ}$, the emission and rotation profiles no longer decrease
monotonically with LOS distance, which further complicates the model. In the
scenario of a high Faraday depth dominated regime ($h_{\phi}<h_{s}$), there
will be more emission at smaller Faraday depths in the lower longitude regions
than at $\ell=180^{\circ}$. The effect of this increased emission at lower
Faraday depths would be a shift toward a more slab-like regime for lower
longitudes in the case of $h_{\phi}<h_{s}$ (predicted by the parameters in
Jaffe (2010)). On the other hand, there would be a shift further from a slab-
like regime for lower longitudes in the case of $h_{\phi}>h_{s}$ (predicted by
the parameters in Sun (2008)).

The model presented here is rudimentary and is intended as a toy model only.
Furthermore, it is not possible to constrain the parameters using a dataset
with limited frequency coverage. Nevertheless, the model is illustrative of
the fact that it is possible to construct configurations of ISM parameters
that result in a constant ratio of observed XE to EG RMs. Future work will
include refining these models and incorporating more realistic aspects to
explore exactly how the results would deviate from the basic structure
described here. Another possibility for the ratio of 2 between EG and XE RMs
is an averaging effect arising from the synchrotron emission being
concentrated in a small number of patches with minimal depth, distributed
randomly throughout the Galaxy.

###  4.3 Regions of Disparity between the XE and EG RMs

While there is remarkable agreement between the average EG RMs and twice the
average XE RMs over a significant portion of the CGPS span (as seen in Figure
4), there are a few noteworthy regions where we see disagreement between the
two sources. Since these differences highlight interesting Galactic features,
we discuss them here.

#### Latitudinal variation in HII complex: $\ell=172^{\circ}$

At $\ell=172^{\circ}$ the XE RM $\times 2$ is seen to deviate significantly
from the EG RM curve, and here both the XE and EG RMs suddenly become more
negative than at longitudes $\pm 5^{\circ}$ from this location. For both types
of sources, RM values with such large magnitudes are unexpected as this
longitude is near $\ell=180^{\circ}$, where the LOS becomes perpendicular to
the field lines and as a result we expect RM values that are near zero.

Closer inspection of the region near $\ell=172^{\circ}$, reveals numerous HII
regions at a distance of about 2 kpc Kang (2012), organized into a structure
described by Gao et al. Gao (2010) as “bow-tie”-like, as seen in the Stokes
$I$ image in the left panel of Figure 8. This structure corresponds to a
depolarised bubble, as shown in the polarised intensity image in the middle
panel of Figure 8. As a result, there is very little XE data in this region,
although there is a high density of EG sources. The exception to this is near
the bottom of the structure, where we see unusually large RM magnitudes for
both EG sources and XE. This is a Faraday screen estimated to be nearer than
300 pc Gao (2010). Since the majority of the XE RM pixels in this longitude
slice come from this small patch, the average XE RM magnitude in the bin at
$\ell=172^{\circ}$ is considerably higher than the average EG RM magnitude,
which results from sources at both high and low latitudes within this slice.

Figure 8: An HII region complex with an anomalous ratio of EG to XE RM values.
Left: Stokes $I$ total intensity map. Middle: Polarised intensity map,
revealing the depolarisation in this region. Right: RM values of EG sources
(dots) and the XE (background pixels.) The dashed line box highlights a
grouping of strongly negative EG RM sources.

#### Cygnus X region: $\ell=75^{\circ}$ to $85^{\circ}$

Cygnus X is one of the most complex emission regions in the Galaxy. It
contains hundreds of HII regions, stellar clusters and OB associations. There
are three emission zones along the LOS at distances of 500–800 pc, at 1.7 kpc
(a massive OB association), and other HII regions spread over 1.5 to 2 kpc
Gottschalk (2012). Together, these regions generate a high intensity of radio
emission that is mostly thermal, making Cygnus X extremely bright in Stokes I,
and relatively faint in polarised emission.

In the Cygnus X region there is a series of very high magnitude negative RM EG
sources found at all latitudes. The strong negative dip in RM of EG sources is
likely due to their emission passing though regions of high electron density
in Cygnus X, enhancing their RM values. Polarised XE observed for these
longitudes must have a very local origin, since Cygnus X has very strong
depolarizing effects on more distant emission. Furthermore, nearly all of the
XE pixels are discarded in the filtering process, due to the brightness of
Cygnus X in Stokes $I$, leaving a small sample size of XE RMs. The result is
that in Figure 4 we see a significantly stronger negative average EG RM in the
bins near $\ell=80^{\circ}$ than the average XE RM.

The region near $\ell=108^{\circ}$ also has extremely low number of XE pixels
remaining after filtering. Likely as a consequence of this, the XE and EG RMs
deviate significantly at this longitude. Furthermore, the low number of XE
pixels does not allow any meaningful interpretation of this deviation.

#### Reversal region: $\ell=55^{\circ}$ to $\ell=70^{\circ}$

An interesting feature of the XE and EG RM trends with longitude (Figure 4) is
that the XE RMs appear to change sign at a slightly higher longitude (close to
$\ell=67^{\circ}$) than the EG RMs (close to $\ell=60^{\circ}$). One possible
explanation for this effect may be the following. If the depth to which the XE
RMs probe is limited by a nearby polarisation horizon in this direction, then
it is possible to have the XE RMs dominated by a field region located
relatively nearby, while the EG RMs for the same LOS are dominated by the
field direction on the far side of the reversed region. At lower longitudes,
where a greater portion of the EG source LOS passes through the reversed
region, the EG source RMs will also reverse signs.

An alternative interpretation does not rely on the concept of a polarisation
horizon, but rather is an extension on the slab-like models we presented in
Section 4.2. By including a reversal of the field in those models, we have
found that as a function of the distance to the reversal, the expected ratio
between the EG and XE RMs varies, and can be negative for some distances and
positive for others. Since the distance to the boundary of the reversed region
does vary with longitude, it is possible to construct configurations of the
ISM parameters that lead to the observed result of the XE and EG RMs reversing
sign at slightly different longitudes. Constraining the possible parameters
for this is beyond the scope of this paper, but would be an interesting avenue
for further investigation.

Considering the diagonal orientation of the boundary across which the reversal
occurs, it may be argued that it is artificial to average over vertical bins
in the longitude-latitude maps. Averaging the XE and EG RMs in bins parallel
to the boundary was illustrated in Figure 3 of Ordog (2017). However,
calculating the ratios of those diagonally binned EG and XE RM values still
does not yield a constant value in this region. This may be due to the
possibility suggested in Ordog (2017) that the current sheet defining the
boundary between opposing magnetic field directions is tilted not only in the
plane of the sky-projection, but also along the line of sight, leading to
complicated three-dimensional variations in the ratio of the detected EG to XE
RMs near the reversal.

##  5 Summary

We have analyzed RMs of the polarised extended synchrotron emission (XE) along
lines of sight through the Galactic disk, and compared these to RMs of
extragalactic (EG) sources that are more commonly used as probes of the GMF.
Our aim was to evaluate the reliability of the XE RMs calculated from the
Synthesis Telescope polarisation data of the CGPS as a tool to probe the GMF.
Potential concerns with this particular dataset may be the lack of single-
antenna data in the 4-band observations, which leads to reduced sensitivity to
large-scale structure in the polarisation images, and the limitations in
detecting Faraday complexity that arise from using only four frequency
channels. We have demonstrated that despite these potential limitations, there
does appear to be useful information in this dataset.

We applied stringent filtering criteria to the XE RMs to exclude: (i) pixels
that correspond to compact sources or lines of sight passing through extended
structures that are not part of the diffuse emission, (ii) pixels with low
signal to noise, and (iii) pixels with a poor linear fit to polarisation angle
versus wavelength squared. Although this filtering process excluded a
significant fraction of the available data, the remaining XE data still probe
many more lines of sight than the EG sources.

To compare the CGPS XE and EG RMs, we plotted them as a function of Galactic
longitude and found that over a large portion of the CGPS range they appear to
follow the same pattern. This is promising, because the CGPS EG sources match
the trends observed in other datasets, and have been determined to be a
reliable tool for constraining models of the GMF. The fact that the XE sources
also track the same trend lends credibility to this type of dataset being used
for the same purpose.

Interestingly, the EG RM magnitudes are greater than the XE RM magnitudes by
approximately a factor of 2 over a significant portion of the longitude range
studied, which would be expected in the case of a Burn slab model of the ISM.
A Burn slab model is unlikely to be realistic for most regions in the ISM,
where the Faraday rotation and synchrotron emissivity vary along the line of
sight. However, we have shown that as long as these parameters vary over
similar scales, the resulting Faraday depth profile may be close to uniform,
with the ratio between EG and XE RMs close to 2 over most wavelengths.

One of the regions where the ratio of EG to XE RMs deviates from 2 is near the
reversal in the GMF. We have observed that the change in sign of the XE and EG
RMs does not occur at the same longitude. This could be explained by a
polarisation horizon in this region limiting the depth to which the XE probes
the ISM, or by the asymmetric slab models we explored, combined with the
distance to the reversal varying with longitude.

Based on our results, we conclude that the XE RMs are indeed a useful dataset,
particularly for tracing the GMF in the outer Galaxy, $\ell>90^{\circ}$. The
lines of sight at lower longitudes, probing a portion of the inner Galaxy, are
more complicated due to the apparent large-scale reversal in the GMF in this
region. For a comparison between the EG and XE sources in this region, the
uncertainty in the interpretation of their individual RMs is compounded. At
least for the outer Galaxy, the XE RMs can be incorporated into models of the
GMF in combination with the more typically used compact source RMs.
Improvements to this technique could be made by using combined single-dish and
synthesis telescope data over a wider frequency range, in order to allow for
added sensitivity to large-scale spatial structures and Faraday complexity
along the lines of sight.

Conceptualization, A.O. and J.C.B.; methodology, A.O., R.A.B., C.V.E. and
J.C.B.; formal analysis, A.O., R.A.B. and C.V.E.; investigation, A.O., R.A.B.
and C.V.E.; writing—original draft preparation, A.O. and R.A.B.;
writing—review and editing, C.V.E., J.C.B. and T.L.L.; supervision, J.C.B. and
T.L.L.; funding acquisition, J.C.B. and T.L.L.

The Canadian Galactic Plane Survey was a Canadian project with international
partners. It was supported by the Natural Sciences and Engineering Research
Council. The Dominion Radio Astrophysical Observatory is a national facility
operated by the National Research Council Canada. This research has been
supported by grants from the Natural Sciences and Engineering Research
Council.

###### Acknowledgements.

We thank Bryan Gaensler and Jennifer West for helpful discussions about the
data and ISM models. We also thank the two reviewers whose comments and
suggestions have helped to improve this manuscript. The authors declare no
conflict of interest. The funding sponsors had no role in the design of the
study; in the collection, analyses, or interpretation of data; in the writing
of the manuscript, or in the decision to publish the results. no one

## Appendix A

Here we present the details of the analytical model outlined in Section 4.2
for determining the Faraday Depth spectrum corresponding to exponential models
of ISM parameters. For the Burn slab scenario, the Faraday depth profile is of
the form:

$\phi(x)=0.812n_{e}B_{||}x,$ (4)

where $n_{e}$ is the electron density, $B_{||}$ is the LOS magnetic field
component, and $x$ is the LOS distance between the observer and a specific
location within the emitting slab. The values of $\phi(x)$ range from 0 rad
m-2 to $\phi_{max}$, which is the Faraday depth of the far side of the slab.
Assuming uniform complex synchrotron emission, $\tilde{P_{\circ}}$, the
resulting Faraday depth spectrum Brentjens (2005) is:

$\tilde{P}(\phi)=\begin{cases}\mbox{\Large$\frac{\tilde{P}_{\circ}}{0.812n_{e}B_{||}}$}&0\leq\phi\leq\phi_{max}\\\
0&\text{otherwise}.\end{cases}$ (5)

Applying a Fourier transform to this function from $\phi$ space to
$\lambda^{2}$ space, we can calculate the change in observed polarisation
angles:

$\Delta\tau(\lambda^{2})=\frac{1}{2}\lambda^{2}\phi_{max},$ (6)

which is the well-known result for the expected RM from a Burn slab
configuration Burn (1966); Sokoloff (1998). The Faraday depths from each point
along the LOS average to yield a RM that is half of the value of the RM of a
compact (Faraday thin) source seen through the same volume.

For the non-uniform models presented in Section 4.2, we begin with exponential
functions with different scale-lengths for the thermal electron density,
magnetic field strength and synchrotron emissivity:

$\displaystyle n_{e}(x)$ $\displaystyle=$ $\displaystyle
n_{\circ}e^{-x/h_{n}}$ $\displaystyle B_{||}(x)$ $\displaystyle=$
$\displaystyle B_{\circ}e^{-x/h_{B}}$ $\displaystyle\tilde{P}(x)$
$\displaystyle=$ $\displaystyle\tilde{P}_{\circ}e^{-x/h_{s}}.$ (7)

From the electron density and magnetic field strength, we determine a Faraday
Depth profile along the line of sight:

$\displaystyle\phi(x)$ $\displaystyle=$ $\displaystyle
0.812\int_{0}^{x}n_{\circ}e^{-x/h_{n}}B_{\circ}e^{-x/h_{B}}dx$ (8)
$\displaystyle=$ $\displaystyle
0.812n_{\circ}B_{\circ}\int_{0}^{x}e^{-\frac{x}{h_{\phi}}}dx$ $\displaystyle=$
$\displaystyle 0.812n_{\circ}B_{\circ}h_{\phi}(1-e^{-\frac{x}{h_{\phi}}})$
$\displaystyle=$ $\displaystyle\phi_{max}(1-e^{-\frac{x}{h_{\phi}}}),$

where we have combined the scale-length parameters as
$h_{\phi}=(1/h_{n_{e}}+1/h_{B})^{-1}$. We define
$\phi_{max}=0.812n_{\circ}B_{\circ}h_{\phi}$, which is the maximum possible
Faraday rotation from this configuration. In general, the Faraday Depth
spectrum can be calculated from the emission profile as:

$\tilde{P}(\phi)=\int_{0}^{\infty}\delta(\phi-\phi(x))\tilde{P}(x)dx=\int_{-\infty}^{\infty}\delta(\phi-\phi(x))\tilde{P}(x)dx,$
(9)

which corresponds to summing up the polarised emission from all physical
depths that result in each particular Faraday depth, $\phi$, thus building up
a spectrum. Please note that the integral can be written with either 0 or
$-\infty$ as the lower bound since we only detect emission from one direction
at a time. The delta function can be re-written as:

$\delta(\phi-\phi(x))=\sum\limits_{i=1}^{n}\frac{\delta(x-x_{i})}{|\frac{\partial(\phi-\phi(x))}{\partial
x}|_{x_{i}}}=\sum\limits_{i=1}^{n}\frac{\delta(x-x_{i})}{|\frac{\partial\phi}{\partial
x}|_{x_{i}}},$ (10)

where the $x_{i}$ are the roots of its argument:

$\phi-\phi(x_{i})=0,$ (11)

and $n$ is the number of roots. Combining Equations (9) and Equation (10)
results in:

$\tilde{P}(\phi)=\int_{-\infty}^{\infty}\sum\limits_{i=1}^{n}\frac{\delta(x-x_{i})}{|\frac{\partial\phi}{\partial
x}|_{x_{i}}}\tilde{P}(x)dx=\sum\limits_{i=1}^{n}\frac{\tilde{P}(x_{i})}{|\frac{\partial\phi}{\partial
x}|_{x_{i}}}.$ (12)

The integral can be dropped because the integrand only picks up the locations
where $x=x_{i}$. For the Faraday Depth profile described by Equation (8),
Equation (11) has only one root:

$\displaystyle\phi(x_{i})$ $\displaystyle=$
$\displaystyle\phi=\phi_{max}(1-e^{-\frac{x_{i}}{h_{\phi}}})$ $\displaystyle
x_{i}$ $\displaystyle=$ $\displaystyle-
h_{\phi}\text{ln}\Big{(}1-\frac{\phi}{\phi_{max}}\Big{)}.$ (13)

The derivative evaluated at the root is:

$\displaystyle\Big{|}\frac{\partial\phi}{\partial x}\Big{|}_{x_{i}}$
$\displaystyle=$
$\displaystyle\frac{\phi_{max}}{h_{\phi}}e^{-\frac{x_{i}}{h_{\phi}}}$ (14)
$\displaystyle=$
$\displaystyle\frac{\phi_{max}}{h_{\phi}}\Big{(}1-\frac{\phi}{\phi_{max}}\Big{)}.$

Combining Equation (12) with Equation (13) and Equation (14), results in:

$\displaystyle\tilde{P}(\phi)$ $\displaystyle=$
$\displaystyle\frac{\tilde{P}_{\circ}e^{-x_{i}/h_{s}}}{\frac{\phi_{max}}{h_{\phi}}\Big{(}1-\frac{\phi}{\phi_{max}}\Big{)}}$
(15) $\displaystyle=$
$\displaystyle\frac{\tilde{P}_{\circ}e^{\frac{h_{\phi}}{h_{s}}\text{ln}\Big{(}1-\frac{\phi}{\phi_{max}}\Big{)}}}{\frac{\phi_{max}}{h_{\phi}}\Big{(}1-\frac{\phi}{\phi_{max}}\Big{)}}$
$\displaystyle=$
$\displaystyle\frac{h_{\phi}\tilde{P}_{\circ}}{\phi_{max}}\Big{(}1-\frac{\phi}{\phi_{max}}\Big{)}^{\frac{h_{\phi}}{h_{s}}-1}$
$\displaystyle=$
$\displaystyle\frac{\tilde{P}_{\circ}}{0.812n_{\circ}B_{\circ}}\Big{(}1-\frac{\phi}{\phi_{max}}\Big{)}^{\frac{h_{\phi}}{h_{s}}-1}.$

Applying a Fourier transform to this Faraday depth spectrum allows us to
determine the detected polarisation as a function of wavelength squared that
results from this configuration:

$\tilde{P}(\lambda^{2})=\int_{-\infty}^{\infty}\tilde{P}(\phi)e^{2i\lambda^{2}\phi}d\phi=\int_{0}^{\phi_{max}}\tilde{P}(\phi)e^{2i\lambda^{2}\phi}d\phi.$
(16)

We evaluated this Fourier transform numerically for different values of the
ratio $h_{\phi}/h_{s}$ to yield the results shown in Figure 7.

References

## References

  * Alfvén (1937) Alfvén, H. Cosmic radiation as an intragalactic phenomenon. Ark. Mat. Astron. Fys. 1937, 27A, 1–5.
  * Fermi (1949) Fermi, E. On the Origin of the Cosmic Radiation. Phys. Rev. 1949, 75, 1169–1174, doi:10.1103/PhysRev.75.1169
  * Hiltner (1949) Hiltner, W.A. Polarization of Radiation from Distant Stars by the Interstellar Medium Nature 1949, 163, 283, doi:10.1038/163283a0.
  * Hall (1949) Hall, J.S. Observations of the Polarized Light from Stars. Science 1949, 109, 166–167, doi:10.1126/science.109.2825.166.
  * van Loo (2012) van Loo, S.; Hartquist, T.W.; Falle, S.A.E.G. Magnetic fields and star formation. Astron. Geophys. 2012, 53, 5.31–5.36, doi:10.1111/j.1468-4004.2012.53531.x.
  * Boulares (1990) Boulares, A.; Cox, D.P. Galactic hydrostatic equilibrium with magnetic tension and cosmic-ray diffusion. Astrophys. J. 1990, 365, 544–558, doi:10.1086/169509.
  * Kim (1996) Kim, E.-J.; Olinto, A.V.; Rosner, R. Generation of Density Perturbations by Primordial Magnetic Fields. Astrophys. J. 1996, 468, 28–50, doi:10.1086/177667.
  * Simard-Normandin (1979) Simard-Normandin, M.; Kronberg, P.P. New large-scale magnetic features of the Milky Way. Nature 1979, 279, 115–118, doi:10.1038/279115a0.
  * Rand (1994) Rand, R.J.; Lyne, A.G. New Rotation Measures of Distant Pulsars in the Inner Galaxy and Magnetic Field Reversals. Mon. Not. R. Astron. Soc. 1994, 268, 497–505, doi:10.1093/mnras/268.2.497.
  * Brown (2003) Brown, J.C.; Taylor, A.R.; Jackel, B.J. Rotation Measures of Compact Sources in the Canadian Galactic Plane Survey. Astrophys. J. Suppl. Ser. 2003 145, 213–223, doi:10.1086/346082.
  * Taylor (2009) Taylor, A.R.; Stil, J.M.; Sunstrum, C. A Rotation Measure Image of the Sky. Astrophys. J. 2009, 702, 1230–1236, doi:10.1088/0004-637X/702/2/1230.
  * Jaffe (2010) Jaffe, T.R.; Leahy, J.P.; Banday, A.J.; Leach, S.M.; Lowe, S.R.; Wilkinson A. Modelling the Galactic magnetic field on the plane in two dimensions. Mon. Not. R. Astron. Soc. 2010, 401, 1013–1028, doi:10.1111/j.1365-2966.2009.15745.x.
  * VanEck (2011) Van Eck, C.L.; Brown, J.C.; Stil, J.M.; Rae, K.; Mao, S.A.; Gaensler, B.M.; Shukurov, A.; Taylor, A.R.; Haverkorn, M.; Kronberg, P.P.; McClure-Griffiths, N.M. Modeling the Magnetic Field in the Galactic Disk Using New Rotation Measure Observations from the Very Large Array. Astrophys. J. 2011, 728, 97, doi:10.1088/0004-637X/728/2/97.
  * Jansson (2012) Jansson, R.; Farrar, G.R. A New Model of the Galactic Magnetic Field. Astrophys. J. 2012, 757, 14, doi:10.1088/0004-637X/757/1/14.
  * Gressel (2013) Gressel, O.; Elstner, D.; Ziegler, U. Towards a hybrid dynamo model for the Milky Way. Astron. Astrophys. 2013, 560, A93, doi:10.1051/0004-6361/201322349.
  * Han (2018) Han, J.L.; Manchester, R.N.; van Straten, W.; Demorest, P. Pulsar Rotation Measures and Large-scale Magnetic Field Reversals in the Galactic Disk. Astrophys. J. Suppl. Ser. 2018, 234, 11, doi:10.3847/1538-4365/aa9c45.
  * Thomson (1980) Thomson, R.C.; Nelson, A.H. The interpretation of pulsar rotation measures and the magnetic field of the galaxy. Mon. Not. R. Astron. Soc. 1980, 191, 863–870, doi:10.1093/mnras/191.4.863.
  * Iacobelli (2013) Iacobelli, M.; Haverkorn, M.; Katgert, P. Rotation measure synthesis at the 2 m wavelength of the FAN region: Unveiling screens and bubbles. Astron. Astrophys. 2013, 549, A56, doi:10.1051/0004-6361/201220175.
  * Burn (1966) Burn, B.J. On the depolarization of discrete radio sources by Faraday dispersion. Mon. Not. R. Astron. Soc. 1966, 133, 67–83, doi:10.1093/mnras/133.1.67.
  * Sokoloff (1998) Sokoloff, D.D.; Bykov, A.A.; Shukurov, A.; Berkhuijsen, E.M.; Beck, R.; Poezd, A.D. Depolarization and Faraday effects in galaxies. Mon. Not. R. Astron. Soc. 1998, 299, 189–206, doi:10.1046/j.1365-8711.1998.01782.x.
  * Brentjens (2005) Brentjens, M.A.; de Bruyn, A.G. Faraday rotation measure synthesis. Astron. Astrophys. 2005, 441, 1217–1228, doi:10.1051/0004-6361:20052990.
  * Ordog (2017) Ordog, A.; Brown, J.C.; Kothes, R.; Landecker, T.L. Three-dimensional structure of the magnetic field in the disk of the Milky Way. Astron. Astrophys. 2017, 603, A15, doi:10.1051/0004-6361/201730740.
  * Landecker (2010) Landecker, T.L.; Reich, W.; Reid, R.I.; Reich, P.; Wolleben, M.; Kothes, R.; Uyanıker, B.; Gray, A.D.; Del Rizzo, D.; Fürst, E.; et al. A Survey of the polarized emission from the Galactic plane at 1420 MHz with arcminute angular resolution. Astron. Astrophys. 2010, 520, A80, doi:10.1051/0004-6361/200913921.
  * Brown (2002) Brown, J.C. The Magnetic Field in the Outer Galaxy. Ph.D. Thesis, University of Calgary, Calgary, AB, Canada, 2002.
  * Oppermann (2012) Oppermann, N.; Junklewitz, H.; Robbers, G.; Bell, M.R.; Enßlin, T.A.; Bonafede, A.; Braun, R.; Brown, J.C.; Clarke, T.E.; Feain, I.J.; et al. An improved map of the Galactic Faraday sky. Astron. Astrophys. 2012, 542, A93, doi:10.1051/0004-6361/201118526.
  * Uyaniker (2003) Uyaniker, B.; Landecker, T.L.; Gray, A.D.; Kothes, R. Radio Polarization from the Galactic Plane in Cygnus. Astrophys. J. 2003, 585, 785-800, doi:10.1086/346234.
  * Sun (2008) Sun, X.H.; Reich, W.; Waelkens, A.; Enßlin, T.A. Radio observational constraints on Galactic 3D-emission models. Astron. Astrophys. 2008, 477, 573–592, doi:10.1051/0004-6361:20078671.
  * Cordes (2002) Cordes, J.M.; Lazio, T.J.W. NE2001.I. A New Model for the Galactic Distribution of Free Electrons and its Fluctuations. _arXiv_ 2002, astro-ph/0207156.
  * Kang (2012) Kang, J.-H.; Koo, B.-C.; Salter, C. An Old Supernova Remnant within an H II Complex at l $\approx$173∘: FVW 172.8+1.5. Astron. J. 2012, 143, 75, doi:10.1088/0004-6256/143/3/75.
  * Gao (2010) Gao, X.Y.; Reich, W.; Han, J.L.; Sun, X.H.; Wielebinski, R.; Shi, W.B.; Xiao, L.; Reich, P.; Fürst, E.; Chen, M.Z.; Ma, J. A Sino-German $\lambda$6 cm polarization survey of the Galactic plane. II. The region from 129∘ to 230∘ longitude. Astron. Astrophys. 2010, 515, A64, doi:10.1051/0004-6361/200913793.
  * Gottschalk (2012) Gottschalk, M.; Kothes, R.; Matthews, H.E.; Landecker, T.L.; Dent, W.R.F. The JCMT 12CO(3-2) survey of the Cygnus X region. I. A pathfinder. Astron. Astrophys. 2012, 541, A79, doi:10.1051/0004-6361/201118600.


NuQS Collaboration

# Quantum Simulation of Field Theories Without State Preparation

Siddhartha Harmalkar<EMAIL_ADDRESS>Department of Physics, University of
Maryland, College Park, MD 20742, USA  Henry Lamm<EMAIL_ADDRESS>Fermi
National Accelerator Laboratory, Batavia, Illinois, 60510, USA  Scott Lawrence
<EMAIL_ADDRESS>Department of Physics, University of Maryland, College Park, MD
20742, USA

###### Abstract

We propose an algorithm for computing real-time observables using a quantum
processor while avoiding the need to prepare the full quantum state. This
reduction in quantum resources is achieved by classically sampling
configurations in imaginary-time using standard lattice field theory. These
configurations can then be passed to quantum processor for time-evolution.
This method encounters a signal-to-noise problem which we characterize, and we
demonstrate the application of standard lattice QCD methods to mitigate it.

††preprint: FERMILAB-PUB-20-045-T

Nonperturbative computational methods for quantum field theories are limited
by an inability to access real-time observables, such as conductivities and
viscosities. Stochastic methods (in particular lattice QCD) suffer from a
maximally bad sign problem Alexandru  _et al._ (2016), whereas deterministic
methods are infeasible due to the exponential growth of the Hilbert space with
volume. The creation of quantum processors promises access to these
nonperturbative observables: quantum processors, unlike their classical
counterparts, can efficiently represent and manipulate the exponentially large
Hilbert space.

A natural way to simulate a quantum mechanical system on a quantum processor
is to map the physical Hilbert space to that of the processor, and then
perform unitary operations that mimic time-evolution Lloyd (1996). This
strategy works for lattice-regularized quantum field theories Ortiz  _et al._
(2001); Zohar  _et al._ (2017); Lamm  _et al._ (2019a); Klco  _et al._ . The
primary theoretical difficulty for such simulations is efficiently creating
the initial, strongly-coupled quantum state on the quantum processor. Current
algorithms for this state preparation step are expensive to implement (usually
dominating the circuit depth of simulations), difficult to analyze, and often
unlikely to generalize. In particular, preparing scattering states requires
substantial complexity Jordan  _et al._ (2012, 2011); García-Álvarez  _et al._
(2015); Jordan  _et al._ (2014, 2017); Hamed Moosavian and Jordan (2018);
Moosavian  _et al._ (2019); Gustafson  _et al._ . Consequently, many works
have studied the problem, of which Kokail  _et al._ (2018); Lamm and Lawrence
(2018); Klco and Savage (2019a, b) is a small sample. Here, we show how to
offload this step to a classical computer, reducing the quantum circuit depth.

In Lamm and Lawrence (2018), two of us argued that preparing the full quantum
state on a quantum processor was unnecessary. Instead, the quantum processor
could be used as a black-box oracle for calculating real-time matrix elements
between basis states which are sampled based on a classical Monte Carlo
calculation of the density matrix. However, that method relied on Density
Matrix Quantum Monte Carlo (DMQMC) Blunt  _et al._ (2014) which is unwieldy
for lattice field theory.

In this Letter, we present a simpler algorithm using standard tools from
lattice field theory for the classical portion of the computation. The
algorithm proceeds as a Monte Carlo calculation in Euclidean lattice field
theory, invoking a quantum processor only to compute matrix elements of time-
dependent operators. Our method suffers from a classical signal to noise (StN)
problem of a more general kind than the usual Parisi-Lepage type Parisi
(1984); Lepage (1989). Existing lattice field theory techniques can address
this StN problem. We discuss some of these modifications and their
incorporation into our algorithm if necessary for evaluating real-time
correlators on near-term processors.

Consider a thermal state given by a density matrix $\rho\equiv e^{-\beta
H_{0}}$ with a Hamiltonian $H_{0}$ and inverse temperature $\beta$. Of
interest is the response of this state to time-evolution by $e^{iH_{1}t}$
($H_{1}$ need not be close to $H_{0}$). The expectation value of
$\mathcal{O}(t)\equiv e^{iH_{1}t}\mathcal{O}e^{-iH_{1}t}$ is given by:

$\left<\mathcal{O}(t)\right>=\frac{\operatorname{Tr}e^{-\beta
H_{0}}\mathcal{O}(t)}{\operatorname{Tr}e^{-\beta
H_{0}}}=\frac{\sum_{i,j}\rho_{ji}\mathcal{O}(t)_{ij}}{\sum_{i}\rho_{ii}}=\left(\frac{\sum_{i,j}\rho_{ji}\mathcal{O}(t)_{ij}}{\sum_{i,j}\rho_{ij}}\right)\left(\frac{\sum_{i}\rho_{ii}}{\sum_{i,j}\rho_{ij}}\right)^{-1}\equiv\frac{\langle\mathcal{O}(t)\rangle_{\rho}}{\langle\delta_{ij}\rangle_{\rho}}$
(1)

where $\rho_{ij}$ and $\mathcal{O}(t)_{ij}$ respectively denote matrix
elements $\left<\Psi_{i}\right|\rho\left|\Psi_{j}\right>$ and
$\left<\Psi_{i}\right|\mathcal{O}(t)\left|\Psi_{j}\right>$, with the
$\Psi_{i}$ states in a basis easily prepared on the quantum processor.
Critically, the basis states are cheap to prepare. The notation
$\langle\mathcal{\cdot}\rangle_{\rho}$ denotes expectation values sampled from
the distribution $\rho_{ij}\equiv\langle\Psi_{j}|\rho|\Psi_{i}\rangle$. The
overall normalization $\langle\delta_{ij}\rangle_{\rho}$ measures the weight
of $\sum_{i}\rho_{ii}$; while often unneeded, an efficient classical procedure
for it is in the appendix.

Figure 1:  Two schematic views: (left) the Schwinger-Keldysh contour in the
complex time plane with the $\tau$ path integral ( ) discretized with lattice
spacing $\Delta\tau$ and the $t$ path integral ( ) trotterized with step
$\Delta t$. ( ) are the matching points between the two path integrals, while
( ) correspond to inserting $\mathcal{O}$ (right) cartoon of how the StN can
be improved through smearing, $\Box[\Psi]=\tilde{\Psi}$, and interpolators $P$
permit the preparation of configurations overlapping with non-thermal states

The connection between Eq. (1) and lattice field theory may be seen via the
Schwinger-Keldysh formalism Schwinger (1961); Keldysh (1964), where
$\langle\mathcal{O}(t)\rangle$ is given as a path integral in both real time
$t$ and imaginary time $\tau$:

$\langle\mathcal{O}(t)\rangle=\frac{1}{Z}\int\mathcal{D}\Psi\;e^{iS_{SK}[\Psi]}\;\mathcal{O}(t)$
(2)

where $S_{SK}=\int\mathrm{d}tL_{1}[x]+\int\mathrm{d}\tau L_{0}[x]\equiv
S_{1}[\Psi(t)]+S_{0}[\Psi(\tau)]$ is integrated along a contour $\mathcal{C}$
in both $t$ and $\tau$, shown on the left of Fig. (1). This closed contour can
be decomposed into coupled open contours that are purely Minkowski and
Euclidean, respectively111Since the path integrals have different arguments
$T,H_{1},\Delta t,etc.$ versus $\beta,H_{0},\Delta\tau,etc.$, renormalization
factors will generally differ, and must be computed separately.. These open
contours correspond to the matrix elements $\rho_{ji}$ and
$\mathcal{O}(t)_{ij}$, and the states at the ‘corners’ of the SK contour are
$\Psi_{i}$ and $\Psi_{j}$. To recover Eq. (2), one takes a sum over these
states.

The matrix elements of $\mathcal{O}(t)$ may be efficiently computed on a
quantum processor Low and Chuang (2017); Campbell (2019); Roggero and Carlson
(2019); Zohar and Cirac (2018); Clemente  _et al._ . The elements of $\rho$
correspond to the Euclidean part of the path integral, with state $\Psi_{i}$
at the bottom of the temporal direction and $\Psi_{j}$ at the top. Standard
lattice field theory is an effective technique for calculations involving this
object. For operators $\mathcal{O}$ diagonal in the computational basis, the
the sum over $j$ in Eq. (1) may be eliminated, and the expectation value
$\langle\mathcal{O}\rangle$ is evaluated by sampling with respect to the
diagonal of $\rho$. In contrast, the time-evolved $\mathcal{O}(t)$ is far from
diagonal in our basis, and so we must retain the sums over both $\Psi_{i}$ and
$\Psi_{j}$.

Thus, in Eq. 1, the physical expectation value $\langle\mathcal{O}(t)\rangle$
is written as a ratio of two expectation values taken with respect to the full
density matrix (i.e., the distribution on pairs of states given by
$\rho_{ij}$). When sampling with respect to the diagonal of $\rho$, one
requires that the $\Psi_{i}=\Psi_{j}$ , thus imposing periodic boundary
conditions (PBC) in $\tau$. When sampling from the full $\rho$ instead, both
the $\Psi_{i}$ and $\Psi_{j}$ are summed over, corresponding to open boundary
conditions (OBC). Between these, the operator $\mathcal{O}(t)$ is inserted by
querying the quantum processor for $\mathcal{O}(t)_{ij}$.

OBCs in lattice QCD were first proposed in Luscher and Schaefer (2011) to
study topology and autocorrelation. Other uses were spectroscopy Luscher and
Schaefer (2013), scale setting Höllwieser  _et al._ (2019), and topological
susceptibility Florio  _et al._ . We advocate their use to compute
$\mathcal{O}(t)$.

Our algorithm, then, is to perform a classical lattice Monte Carlo with OBC,
coupled to a quantum processor which evaluates $\mathcal{O}(t)_{ij}$ on each
configuration (with $\Psi_{i}$ and $\Psi_{j}$ the initial and final
$\tau$-slices of the Euclidean lattice). We stress that the imaginary time
evolution is performed on the classical machine, and the quantum processor
only needs enough qubits available to store a single state (i.e. a single
timeslice of the Euclidean lattice).

So far, the discussion has focused on thermal states. Using Euclidean lattice
techniques, we can extend to other states. With configurations sampled from
$\rho$, Euclidean observables are equal to sums over complete set of
eigenstates $|n\rangle$ with weights $e^{-E_{n}\beta}$. If we want the matrix
elements $\langle\psi|\mathcal{O}|\psi\rangle$ of a particular state
$|\psi\rangle$, it is necessary to isolate it from all other states. When
$|\psi\rangle$ is the lowest state with the quantum numbers of some projection
operator $P$, the procedure is simple: insert a source $P$ and a sink
$P^{\dagger}$ between $\rho$ and $\mathcal{O}$, and consider large $\tau$:

$\displaystyle\sum_{n}$ $\displaystyle\langle n|P^{\dagger}\rho
P|n\rangle\langle n|\mathcal{O}|n\rangle$
$\displaystyle=\langle\psi|\mathcal{O}|\psi\rangle\langle\psi|P^{\dagger}P|\psi\rangle
e^{-E_{\psi}\tau}+\mathcal{O}(e^{-\Delta E\tau}).$ (3)

where exponentially-suppressed corrections come from other eigenstates,
$|n\rangle$, that overlap with $P$. To extract
$\langle\psi\mathcal{|}\mathcal{O}|\psi\rangle$, one notes that
$\langle\psi|P^{\dagger}\rho P|\psi\rangle$ can be obtained from a two-point
correlator $\langle P^{\dagger}(\tau)P(0)\rangle$. This method can initialize
the quantum processor with configurations that overlap with the desired state,
avoiding direct preparation.

Figure 2: Average value of the signal-to-noise over early ($t\in[0,0.1]$) and
late ($t\in[5,15]$) times, at two inverse temperatures $\beta$. The StN
decreases exponentially with the size of the system, and is larger for higher
temperatures where $\rho$ is near-diagonal.

While these procedures are theoretically sound, in practice there are StN
problems. This can be seen from Eq. (1) — at nonzero $\beta$, $\rho_{ij}$ is
generally nonzero between any $\Psi_{i}$ and $\Psi_{j}$. Therefore
$\langle\delta_{ij}\rangle_{\rho}=\sum_{i}\rho_{ii}/\sum_{i,j}\rho_{ij}\propto
e^{-V}$ where $V$ is the spatial volume. This is itself not a problem: this
overall normalization can be disregarded or measured efficiently as discussed
in the appendix. However, the physical expectation value
$\langle\mathcal{O}(t)\rangle$ should approach a constant value in the
infinite-volume limit, which implies that
$\langle\mathcal{O}(t)\rangle_{\rho}$ must itself be exponentially small in
the volume. As long as typical $O(t)_{ji}$ are not exponentially small in the
volume, this will constitute a StN problem.

At short times, this StN problem is due to $\mathcal{O}(t)$ having few non-
vanishing matrix elements — an overlap problem between $O(t)_{ij}$ and
$\rho_{ji}$. For an observable chosen such that $\mathcal{O}(0)$ is diagonal,
there are $\sim e^{V}$ non-vanishing elements, out of $\sim e^{2V}$ total. The
nonvanishing $\rho_{ij}$ are thus sampled exponentially infrequently by the
OBC Monte Carlo. At longer times, the generic matrix element is nonzero, and
the StN problem arises instead from fine cancellations between positive and
negative elements.

We demonstrate this by computing the $\mathrm{StN}(t)$ of $\mathcal{O}(t)$,
defined as the ratio of the mean to the standard deviation of the random
variable $\mathcal{O}(t)_{ij}$ sampled from
$\rho_{ij}/\sum_{k,\ell}\rho_{k\ell}$. For this example, the one-dimensional
ferromagnetic Heisenberg model
$H=-J\sum_{i,\alpha}\sigma^{\alpha}_{i}\otimes\sigma^{\alpha}_{i+1}$ with
coupling $J=1$ is simulated with $N\in[4,12]$ sites. The $\mathrm{StN}(t)$ of
the spin-spin correlator $\mathcal{O}(t)=\sigma^{z}_{1}(0)\sigma^{z}_{1}(t)$
acting on one spin, averaged over early and late times, is in Fig. 2.

A variety of methods can reduce the StN in practice. In DMQMC, the StN problem
requires exponentially large samplings to sufficiently populate $\rho$ Blunt
_et al._ (2014). To mitigate these large resource requirements, the authors
reweighted, reducing the probability of sampling far from the diagonal of
$\rho$. This was sufficient for the time-independent $S=1/2$ antiferromagnetic
Heisenberg model. We expect an analogous method is effective in countering the
short-time StN problem, where $\mathcal{O}(t)_{ij}$ vanishes far from the
diagonal.

The severity of the StN problem is basis-dependent, which suggests that
performing these calculations in a different basis can alleviate the issue as
well. In the eigenbasis of $H$, $\rho$ is diagonal and
$\langle\delta_{ij}\rangle_{\rho}=1$. One illustrative example is that of free
fields. For a free theory, the eigenbasis of the Hamiltonian is the momentum
basis, and therefore computations in this basis will have no StN problem. For
strongly-coupled theories, finding the eigenbasis is unlikely, but physical
intuition may motivate bases with better scaling with reduced noise.

Another path to alleviating the StN problem is reducing the dimension of
$\rho$. To do so without affecting the IR, we must wisely remove states from
$\rho$. For gauge theories, the number of group elements $n$ can be infinite
(e.g. $SU(3)$) or merely gigantic (e.g. discrete groups like $S({1080})$,
suggested to approximate $SU(3)$ Alexandru  _et al._ (2019)). In this case, a
pure gauge lattice has $n^{V}$ states in the computational basis and the
dimension of $\rho$ is $n^{2V}$, with some fraction gauge-invariant. The StN
is naively $\langle\delta_{ij}\rangle_{\rho}\approx\frac{n^{V}}{n^{2V}}$.
However, this large density of states is dominated by UV lattice artifacts,
particularly in $(3+1)$d, which do not contribute to IR physics. Thus, these
states can be elided without altering physical measurements.

_Smearing_ methods use projection operators to smooth the field configurations
in lattice QCD. These can be thought of as space or spacetime averaging for
gauge theories which removes UV lattice artifacts, reducing the number of
states to $n_{eff}<n$. Thus $\rho\rightarrow\rho_{eff}$, the size is reduced
to $n_{eff}^{2V}$, and the StN improves to
$\langle\delta_{ij}\rangle_{\rho_{eff}}\approx{n_{eff}^{-V}}$. For gauge
fields, smearing methods include _APE_ Albanese  _et al._ (1987), _HYP_
Hasenfratz and Knechtli (2001), and _Stout_ Morningstar and Peardon (2004).
For fermionic fields, techniques include _Wuppertal_ Gusken  _et al._ (1989)
and _Jacobi_ smearing Allton  _et al._ (1993), which the free fermion kinetic
energy to smooth. Finally, the recent development in momentum smearing for
moving hadrons Bali  _et al._ (2016); Wu  _et al._ (2018) will prove useful
for extracting scattering data from a quantum processor.

Low-lying physical states also contribute to the StN problem. Parisi Parisi
(1984) and Lepage Lepage (1989) recognized that while $|\psi\rangle$ is the
lowest state in $\langle P^{\dagger}(\tau)P(0)\rangle$, it need not be lowest
for the variance, estimated by $\langle
P^{\dagger}(\tau)P(0)P^{\dagger}(0)P(\tau)\rangle$. For example, the proton,
$p$, is the lowest state overlapping with a 3 quark $P_{3q}$ operator. For the
variance, the lowest state is instead 3 pions, $3\pi$. This implies that the
StN for proton observables scales $\propto
e^{-(E_{p}-\frac{3}{2}E_{\pi})\tau}$ at long $\tau$. Together with excited
states, this limits the $\tau$ range to a “golden window” where properties can
be robustly extracted. Practitioners extend this window by another method,
_improved interpolators_ , which removes excited states, variance-overlapping
states, or both.

The idea of improved interpolators is to build projection operators
$\tilde{P}$ from multiple $P$ having the correct quantum numbers and spatial
distributions to improve state overlap. This was done early on variationally
Michael (1985) and recently with more complex methods like distillation
Peardon  _et al._ (2009). Distillation has aided the calculation of a wide
range of observables: multi-hadron states Culver  _et al._ (2019), excited
state spectra Liu  _et al._ (2012), coupled-channel resonances Woss  _et al._
(2019), and nucleon charges Egerer  _et al._ (2019). In a few cases, improved
interpolators directly addressed StN issues Beane  _et al._ (2009); Wagman and
Savage (2017a, b); we highlight the intuitive discussion in Detmold and Endres
(2015). These interpolators are all usable with our method. Quantum processors
may also prove useful in finding improved interpolators Avkhadiev  _et al._
(2019).

Put together, these techniques reduce the StN problem of our method. This is
done by including smearing $\Box$ or improved interpolators $\tilde{P}$
leading to

$\left<\psi|\mathcal{O}(t)|\psi\right>=\frac{\operatorname{Tr}\tilde{P}^{\dagger}\Box^{\dagger}e^{-\beta
H_{0}}\Box\tilde{P}e^{iH_{1}t}\mathcal{O}e^{-iH_{1}t}}{\operatorname{Tr}e^{-\beta
H_{0}}}.$ (4)

Since the StN scales like $e^{-V}$, for large $V$, this may prove necessary,
although for problems amenable to current resources, its use is limited. One
may also imagine smeared $H_{1}$ at the expense of circuit depth to reduce
errors.

As a demonstration, we simulate the $D_{4}$ gauge theory with the Wilson
action of Lamm  _et al._ (2019a) (right of Fig. 3), which needs 12 qubits for
the $8^{4}$ states and 2 ancillary qubits. The Euclidean calculation had
$N_{\tau}=10$ and Wilson coupling $\beta_{W}=0.9$, chosen near the
$\beta_{W}\rightarrow\infty$ ground state with the gauge-invariant projection
of all $U_{i}=\mathbbm{1}$. We sampled $2\times 10^{4}$ configurations
separated by 10 steps. The gauge-invariant states representing the $\Psi_{i}$
and $\Psi_{j}$ were passed to a noise-free quantum simulator. Time evolution
is performed with trotterization step $\delta t=0.5$ and $\beta_{W}=0.8$. The
expectation value of one plaquette is measured versus $t$ (Fig. 3). The real-
time code is written in qiskit Santos (2017) and simulated for 10 shots per
configuration. In additional to the uncertainty from StN, other statistical
uncertainties are the classical Monte Carlo variance and the finite quantum
shots. We estimate one systematic, finite $\delta t$, by the maximum deviation
between the exact and trotterize $\beta_{W}\rightarrow\infty$ results of Lamm
_et al._ (2019a). The error budget is in Table 1.

Table 1: Summary of the fractional uncertainties in determination of $\langle\operatorname{Re}\operatorname{Tr}P_{1}(t)\rangle$ Source | $\alpha=0.0$ | $\alpha=0.7$
---|---|---
Classical Statistics | 2.3% | 2.0%
Quantum Shots | 1.3% | 1.0%
$\langle\delta_{ij}\rangle_{\rho}$ | 1.8% | 1.6%
Trotterization | 2.5% | 2.5%

We also perform $APE$ smearing to improve the StN. Due to the small lattice,
we only smear link $U_{3}$ with $\alpha=0.7$. Agreement is found with the
unsmeared results. The average number of $D_{4}$ elements that $U_{3}$ samples
goes from $n=3.08(2)$ to $n_{eff}=2.98(2)$. $\langle\delta_{ij}\rangle_{\rho}$
increases from 7.8% to 10.0%. From this, we naively estimate an increased StN
of $(3.08/2.98)^{V}$ for larger lattices.

Figure 3: (left) Expectation value of
$\operatorname{Re}\operatorname{Tr}P_{1}=U_{2}^{\dagger}U_{0}^{\dagger}U_{3}U_{0}$
vs. $t$ from a quantum simulator. In blue (red) is the result for
$\beta_{W}=0.9$ with(out) smearing $\alpha=0.7$. (right) Spatial lattice where
dashed lines indicate repeated links due to PBC.

The reduction in circuit depth is a crucial advantage of this procedure.
Adiabatic state preparation requires repeated time evolution, requiring $\sim
1000$ gates _per step_ for the $D_{4}$ model. By contrast, preparing
$\Psi_{i}$ requires only $\sim 1000$ gates, a substantial savings. For larger
lattices, the preparation cost for $\Psi_{i}$ is expected to be $O(V)$,
whereas adiabatic methods are generally $O(V^{4/3})$ for a 3d lattice.

We have proposed a method for computing real-time observables without
performing the expensive step of preparing the highly-entangled state on the
quantum processor. On current quantum processors, this is a worthy trade-off.
To achieve this, we use standard lattice methods to simulate the system at low
temperature, and use the quantum processor as a black box to compute the
matrix elements between basis states. While this method will have issues with
signal to noise, a number of techniques for mitigating can be applied, which
may prove manageable in the NISQ era. Theories with matter fields integrated
out on the Euclidean lattice (i.e. fermions) are the critical next step in
this line of work. Future work could use this method to investigate properties
like PDFs Lamm  _et al._ (2019b) or domain wall dynamics Tan  _et al._ (2019).

###### Acknowledgements.

H.L. is supported by a Department of Energy QuantiSED grant. S.H., H.L, and
S.L., are supported by the U.S. Department of Energy under Contract No. DE-
FG02-93ER-40762. We are grateful to Tom Cohen, Paulo Bedaque, Colin Egerer,
Raju Venugopalan, and Michael Wagman for numerous discussions helpful to
developing this work. Fermilab is operated by Fermi Research Alliance, LLC
under contract number DE-AC02-07CH11359 with the United States Department of
Energy. The authors acknowledge the University of Maryland supercomputing
resources (http://hpcc.umd.edu) made available for conducting the research
reported in this paper.

## References

  * Alexandru  _et al._ (2016) A. Alexandru, G. Basar, P. F. Bedaque, S. Vartak, and N. C. Warrington, Phys. Rev. Lett. 117, 081602 (2016), arXiv:1605.08040 [hep-lat] .
  * Lloyd (1996) S. Lloyd, Science 273, 1073 (1996).
  * Ortiz  _et al._ (2001) G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Phys. Rev. A64, 022319 (2001), arXiv:cond-mat/0012334 [cond-mat] .
  * Zohar  _et al._ (2017) E. Zohar, A. Farace, B. Reznik, and J. I. Cirac, Phys. Rev. A95, 023604 (2017), arXiv:1607.08121 [quant-ph] .
  * Lamm  _et al._ (2019a) H. Lamm, S. Lawrence, and Y. Yamauchi (NuQS), Phys. Rev. D100, 034518 (2019a), arXiv:1903.08807 [hep-lat] .
  * (6) N. Klco, J. R. Stryker, and M. J. Savage, arXiv:1908.06935 [quant-ph] .
  * Jordan  _et al._ (2012) S. P. Jordan, K. S. M. Lee, and J. Preskill, Science 336, 1130 (2012), arXiv:1111.3633 [quant-ph] .
  * Jordan  _et al._ (2011) S. P. Jordan, K. S. M. Lee, and J. Preskill, (2011), [Quant. Inf. Comput.14,1014(2014)], arXiv:1112.4833 [hep-th] .
  * García-Álvarez  _et al._ (2015) L. García-Álvarez, J. Casanova, A. Mezzacapo, I. L. Egusquiza, L. Lamata, G. Romero, and E. Solano, Phys. Rev. Lett. 114, 070502 (2015), arXiv:1404.2868 [quant-ph] .
  * Jordan  _et al._ (2014) S. P. Jordan, K. S. M. Lee, and J. Preskill, (2014), arXiv:1404.7115 [hep-th] .
  * Jordan  _et al._ (2017) S. P. Jordan, H. Krovi, K. S. M. Lee, and J. Preskill, (2017), arXiv:1703.00454 [quant-ph] .
  * Hamed Moosavian and Jordan (2018) A. Hamed Moosavian and S. Jordan, Phys. Rev. A98, 012332 (2018), arXiv:1711.04006 [quant-ph] .
  * Moosavian  _et al._ (2019) A. H. Moosavian, J. R. Garrison, and S. P. Jordan, (2019), arXiv:1911.03505 [quant-ph] .
  * (14) E. Gustafson, P. Dreher, Z. Hang, and Y. Meurice, arXiv:1910.09478 [hep-lat] .
  * Kokail  _et al._ (2018) C. Kokail _et al._ , (2018), arXiv:1810.03421 [quant-ph] .
  * Lamm and Lawrence (2018) H. Lamm and S. Lawrence, Phys. Rev. Lett. 121, 170501 (2018), arXiv:1806.06649 [quant-ph] .
  * Klco and Savage (2019a) N. Klco and M. J. Savage, (2019a), arXiv:1904.10440 [quant-ph] .
  * Klco and Savage (2019b) N. Klco and M. J. Savage, (2019b), arXiv:1912.03577 [quant-ph] .
  * Blunt  _et al._ (2014) N. S. Blunt, T. W. Rogers, J. S. Spencer, and W. M. C. Foulkes, Phys. Rev. B 89, 245124 (2014).
  * Parisi (1984) G. Parisi, Phys. Rept. 103, 203 (1984).
  * Lepage (1989) G. P. Lepage, in  _Boulder ASI 1989:97-120_ (1989) pp. 97–120.
  * Schwinger (1961) J. S. Schwinger, J. Math. Phys. 2, 407 (1961).
  * Keldysh (1964) L. V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964), [Sov. Phys. JETP20,1018(1965)].
  * Low and Chuang (2017) G. H. Low and I. L. Chuang, Phys. Rev. Lett. 118, 010501 (2017).
  * Campbell (2019) E. Campbell, Phys. Rev. Lett. 123, 070503 (2019).
  * Roggero and Carlson (2019) A. Roggero and J. Carlson, Phys. Rev. C100, 034610 (2019), arXiv:1804.01505 [quant-ph] .
  * Zohar and Cirac (2018) E. Zohar and J. I. Cirac, Phys. Rev. B98, 075119 (2018), arXiv:1805.05347 [quant-ph] .
  * (28) G. Clemente _et al._ , arXiv:2001.05328 [hep-lat] .
  * Luscher and Schaefer (2011) M. Luscher and S. Schaefer, JHEP 07, 036 (2011), arXiv:1105.4749 [hep-lat] .
  * Luscher and Schaefer (2013) M. Luscher and S. Schaefer, Comput.Phys.Commun. 184, 519 (2013), arXiv:1206.2809 [hep-lat] .
  * Höllwieser  _et al._ (2019) R. Höllwieser, F. Knechtli, and T. Korzec, (2019), arXiv:1907.04309 [hep-lat] .
  * (32) A. Florio, O. Kaczmarek, and L. Mazur, arXiv:1903.02894 [hep-lat] .
  * Alexandru  _et al._ (2019) A. Alexandru, P. F. Bedaque, S. Harmalkar, H. Lamm, S. Lawrence, and N. C. Warrington (NuQS), Phys.Rev.D 100, 114501 (2019), arXiv:1906.11213 [hep-lat] .
  * Albanese  _et al._ (1987) M. Albanese _et al._ (APE), Phys. Lett. B192, 163 (1987).
  * Hasenfratz and Knechtli (2001) A. Hasenfratz and F. Knechtli, Phys. Rev. D64, 034504 (2001), arXiv:hep-lat/0103029 [hep-lat] .
  * Morningstar and Peardon (2004) C. Morningstar and M. J. Peardon, Phys. Rev. D69, 054501 (2004), arXiv:hep-lat/0311018 [hep-lat] .
  * Gusken  _et al._ (1989) S. Gusken, U. Low, K. H. Mutter, R. Sommer, A. Patel, and K. Schilling, Phys. Lett. B227, 266 (1989).
  * Allton  _et al._ (1993) C. R. Allton _et al._ (UKQCD), Phys. Rev. D47, 5128 (1993), arXiv:hep-lat/9303009 [hep-lat] .
  * Bali  _et al._ (2016) G. S. Bali, B. Lang, B. U. Musch, and A. Schäfer, Phys. Rev. D93, 094515 (2016), arXiv:1602.05525 [hep-lat] .
  * Wu  _et al._ (2018) J. J. Wu, W. Kamleh, D. B. Leinweber, R. D. Young, and J. M. Zanotti, J. Phys. G45, 125102 (2018), arXiv:1807.09429 [hep-lat] .
  * Michael (1985) C. Michael, Nucl. Phys. B259, 58 (1985).
  * Peardon  _et al._ (2009) M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek, R. G. Edwards, B. Joo, H.-W. Lin, D. G. Richards, and K. J. Juge (Hadron Spectrum), Phys. Rev. D80, 054506 (2009), arXiv:0905.2160 [hep-lat] .
  * Culver  _et al._ (2019) C. Culver, M. Mai, R. Brett, A. Alexandru, and M. Döring, (2019), arXiv:1911.09047 [hep-lat] .
  * Liu  _et al._ (2012) L. Liu, G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas, P. Vilaseca, J. J. Dudek, R. G. Edwards, B. Joo, and D. G. Richards (Hadron Spectrum), JHEP 07, 126 (2012), arXiv:1204.5425 [hep-ph] .
  * Woss  _et al._ (2019) A. J. Woss, C. E. Thomas, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Phys. Rev. D100, 054506 (2019), arXiv:1904.04136 [hep-lat] .
  * Egerer  _et al._ (2019) C. Egerer, D. Richards, and F. Winter, Phys. Rev. D99, 034506 (2019), arXiv:1810.09991 [hep-lat] .
  * Beane  _et al._ (2009) S. R. Beane, W. Detmold, T. C. Luu, K. Orginos, A. Parreno, M. J. Savage, A. Torok, and A. Walker-Loud, Phys. Rev. D80, 074501 (2009), arXiv:0905.0466 [hep-lat] .
  * Wagman and Savage (2017a) M. L. Wagman and M. J. Savage, Phys. Rev. D96, 114508 (2017a), arXiv:1611.07643 [hep-lat] .
  * Wagman and Savage (2017b) M. L. Wagman and M. J. Savage, (2017b), arXiv:1704.07356 [hep-lat] .
  * Detmold and Endres (2015) W. Detmold and M. G. Endres,  _Proceedings, 32nd International Symposium on Lattice Field Theory (Lattice 2014): Brookhaven, NY, USA, June 23-28, 2014_ , PoS LATTICE2014, 170 (2015), arXiv:1409.5667 [hep-lat] .
  * Avkhadiev  _et al._ (2019) A. Avkhadiev, P. E. Shanahan, and R. D. Young, (2019), arXiv:1908.04194 [hep-lat] .
  * Santos (2017) A. C. Santos, Revista Brasileira de Ensino de Física 39 (2017), arXiv:1610.06980 [quant-ph] .
  * Lamm  _et al._ (2019b) H. Lamm, S. Lawrence, and Y. Yamauchi (NuQS), (2019b), arXiv:1908.10439 [hep-lat] .
  * Tan  _et al._ (2019) W. L. Tan _et al._ , (2019), arXiv:1912.11117 [quant-ph] .

##  Appendix A OBC Normalization

Naively, the exponentially-small normalization
$\left<\delta_{ij}\right>_{\rho}$ requires $\sim e^{V}$ measurements (in fact,
infinitely many for continuous fields). In this appendix we describe how to
compute this normalization in polynomial time in the volume. The idea is to
interpolate between the OBC and PBC by gradually turning on a term coupling
the first and last time-slices ($\Psi_{i}$ and $\Psi_{j}$). Chose
$f(\alpha)_{ij}$ such that $f(1)_{ij}=\delta_{ij}$ and $f(0)_{ij}=1$. Defining
an interpolating distribution $\rho(\alpha)_{ij}=\rho_{ij}f(\alpha)_{ji}$, we
note that the $\rho(1)$ imposes periodic boundary conditions while $\rho(0)$
has the original open boundary conditions. The relative normalization between
$\rho(\alpha_{i})$ and $\rho(\alpha_{j})$ is given by

$\mathfrak{N}(\alpha_{i},\alpha_{j})=\left<\frac{\rho(\alpha_{j})}{\rho(\alpha_{i})}\right>_{\rho(\alpha_{i})}=\left<\frac{f(\alpha_{j})}{f(\alpha_{i})}\right>_{\rho(\alpha_{i})}\text{.}$
(5)

This normalization is multiplicative in the sense that
$\mathfrak{N}(\alpha,\alpha^{\prime})\mathfrak{N}(\alpha^{\prime},\alpha^{\prime\prime})=\mathfrak{N}(\alpha,\alpha^{\prime\prime})$.
As a result, $\mathfrak{N}(0,1)\equiv\left<\delta_{ij}\right>_{\rho}$ can be
computed by subdividing it into $N$ terms,
$\mathfrak{N}(0,1)=\prod_{i=0}^{N-1}\mathfrak{N}(\alpha_{i},\alpha_{i+1})$.
The performance of this method is sensitive to the choice of interpolating
$\alpha$: if the $\alpha_{i}$ are chosen such that all
$\mathfrak{N}(\alpha_{i},\alpha_{i+1})=e^{-V/N}$ then $\sim Ne^{V/N}$
measurements are required. Choosing $N\sim V$, we see that the measurement of
$\left<\delta_{ij}\right>_{\rho}$ can indeed be performed in polynomial time.