# No Small Scale Radio Jets Here:
Multi-Epoch Observations of Radio Continuum Structures in NGC 1068 with the
VLBA

Travis C. Fischer AURA for ESA, Space Telescope Science Institute, Baltimore,
MD, USA, 3700 San Martin Drive, Baltimore, MD 21218, USA  Megan C. Johnson
U.S. Naval Observatory, 3450 Massachusetts Ave NW, Washington, DC 20392-5420,
USA  Nathan J. Secrest U.S. Naval Observatory, 3450 Massachusetts Ave NW,
Washington, DC 20392-5420, USA  D. Michael Crenshaw Department of Physics and
Astronomy, Georgia State University, 25 Park Place, Suite 605, Atlanta, GA
30303, USA  Steven B. Kraemer Institute for Astrophysics and Computational
Sciences, Department of Physics, The Catholic University of America,
Washington, DC 20064, USA

###### Abstract

We present recent Very Long Baseline Array (VLBA) 5 GHz radio observations of
the nearby, luminous Seyfert 2 galaxy NGC 1068 for comparison to similar VLBA
observations made on 1997 April 26. By cross-correlating the positions of
emitting regions across both epochs, we find that spatially-resolved extra-
nuclear radio knots in this system have sub-relativistic transverse speeds (v
$\lesssim$ 0.1c). We discuss sources of the observed knots and how the radio
emission relates to additional phases of gas in the central $\sim$150 pcs of
this system. We suggest that the most likely explanation for the observed
emission is synchrotron radiation formed by shocked host media via
interactions between AGN winds and the host environment.

Active galaxies (17), Radio active galactic nuclei (2134), Radio astrometry
(1337)

††facilities: VLBA, VLA, EVLA, ALMA, Gemini, HST; The HST data presented in
this paper were obtained from the Mikulski Archive for Space Telescopes (MAST)
at the Space Telescope Science Institute. The specific observations analyzed
can be accessed via https://doi.org/10.17909/nkex-gm67 (catalog DOI:
10.17909/nkex-gm67). ††software: aips (van Moorsel et al., 1996), Astropy
(Astropy Collaboration et al., 2013), casa (CASA Team et al., 2022)

##  1 Introduction

Active galactic nuclei (AGN) are the central regions of galaxies containing
supermassive black holes that are actively accreting matter. This accretion
generates intense radiation as matter falls toward the black hole and is
rapidly heated up in the surrounding accretion disk. This process is complex
and produces various forms of “feedback” as radiation drives material from the
nuclear region and likely has some effect on how the host galaxy continues to
form stars and evolve. The impact of these feedback processes continues to be
studied, as direct connections between AGN feedback and galaxy evolution have
not yet been fully established.

Currently, two modes of feedback are utilized to describe interactions between
the AGN and its host galaxy; ‘kinetic mode’ and ‘radiative mode’. These
individual modes are often used to describe the dominant form of feedback in
distinctively different types of galaxies. Kinetic-mode feedback is often
associated with radio galaxies and radio-loud quasars, typically massive
elliptical galaxies exhibiting relativistic, collimated jets of synchrotron
radio emission from pc to Mpc scales (Hardcastle & Croston 2020 and references
therein). These jets deposit their energy into the galactic and intergalactic
environments, controlling the cooling of gas in galaxy halos and preventing
star formation. Notably, low-power relativistic jets may also be capable of
clearing out molecular gas in the central parts of galaxies (Morganti et al.,
2013; Mukherjee et al., 2018; Murthy et al., 2022). Radiative-mode feedback
occurs when energetic photons emitted by the AGN accretion disk couple with
gas in the disk or host environment, creating winds that travel at speeds of
hundreds to tens of thousands of km s-1, and are observed to propagate
throughout the host on sub-pc to kpc scales. AGN winds likely impact the
surrounding medium, either by removing significant amounts of gas and dust
from the central regions of their host galaxies, which can presumably prevent
star formation (e.g., Silk & Rees 1998; Fabian 2012), or instead compressing
gas to promote star formation (e.g., Silk 2013; Zubovas et al. 2013; also see
the review by Morganti 2017).

These modes have been defined over the last few decades and, individually,
have been studied in great detail. However, as they are unlikely to be
completely distinct, several studies have begun to examine the interactions
between each mode and their overall impacts on their hosts. These interactions
are often difficult to analyze in radio-loud AGN, as hosts either reside at
distances where spatially resolved studies of AGN winds are difficult, such as
in elliptical galaxies that lack significant amounts of gas and dust for AGN
winds to interact with, or are inconveniently oriented such that evidence of
AGN winds is largely extinguished by dust in the host galaxy (e.g., Centaurus
A; Hamer et al. 2015). Alternatively, numerous radio-quiet AGN exhibit
signatures that might indicate both modes of feedback: collimated, jet-like
radio emission and high-velocity ionized-gas emission (i.e., Falcke et al.
1998; Whittle & Wilson 2004; Gallimore et al. 2006; Storchi-Bergmann et al.
2010; Paggi et al. 2012; Maksym et al. 2017; Fischer et al. 2019; Jarvis et
al. 2019; Girdhar et al. 2022).

Observations of these targets indicate that the brightest radio structures are
intertwined with ionized gas structures, known as the Narrow Line Region
(NLR). In the Unified Model of AGN (Urry & Padovani, 1995), it is predicted
that a radio jet will propagate perpendicular to the accretion disk. The jet
is encompassed by an ionizing radiation field that exhibits a biconical shape
which is constrained by the AGN torus. Therefore, to first order, the
alignment between radio jets and gas ionized by the AGN is to be expected.
However, the AGN-ionized NLR gas is also located on the surface of gas
reservoirs in the host galaxies of these AGN (Fischer et al., 2017; Shimizu et
al., 2019). Therefore, radio emission appearing to be colocated with the NLR
gas suggests that the observed jet not only aligns with the radiation field
but also resides primarily embedded in gas reservoirs along the plane of the
host galaxy. This consistent alignment has led us to suggest that the observed
extended radio continuum emission in radio-quiet AGN is often not a jet
pointed into the plane of the galaxy, but a by-product of radiatively driven
winds originating at smaller radii, which extend to larger distances and drive
into dense reservoirs to produce shocks (Fischer et al., 2019, 2021). These
wind-driven shocks heat and compress the reservoir gas, producing significant
amounts of cosmic rays that flow along enhanced magnetic field lines resulting
from the compression, which in turn allow for the in situ production of
synchrotron emission. In other words, we hypothesize that the radio emission
observed radio-quiet AGN may be both a byproduct of radiative-mode feedback
and an indicator of kinetic-mode feedback.

This hypothesis was supported in a recent study (Fischer et al., 2021) of a
volume-complete sample of 25 nearby AGN ($L_{\mathrm{14-195~{}keV}}>10^{42}$
erg s-1, $D<40$ Mpc) selected from the Swift 105-month Burst Alert Telescope
(BAT) catalog (Oh et al., 2018), where we compared the peak-flux radio
emission of the nuclear sources using the Jansky Very Large Array (VLA) and
Very Long Baseline Array (VLBA) at the same sensitivity.

Assuming that the radio continuum in radio-quiet AGN can be attributed to a
jet emanating from the nucleus, it is likely that a majority of the emission
can be contained in a compact, or potentially unresolved, morphology when the
jet is oriented along our line of sight and observed with the VLA. As such,
the peak-flux measured over the nucleus with the VLA would be approximate to
the total flux of the jet. When observed with the VLBA, a majority of this
emission is then resolved out, creating a large VLA/VLBA flux ratio.
Alternatively, if a jet is pointed into the plane of the sky, a small
percentage of the total jet flux is unresolved over the nucleus when observed
with the VLA. Therefore, when observed with the VLBA, a smaller portion of the
emission from the entire jet is resolved out when compared to the peak-flux
measured with the VLBA, creating a smaller VLA/VLBA flux ratio. From our study
of the nearby BAT sample, we observed the opposite effect such that AGN with
higher ratios of VLA to VLBA nuclear peak-fluxes tended to have extranuclear
radio structures extended into the plane of the sky visible in the VLA
observations.

If we instead assume that the observed radio continuum is not directly due to
a jet but can instead be attributed to shocks through kinetic feedback
processes, the ratio between VLA and VLBA fluxes would instead depend on the
orientation of the AGN relative to its host galaxy. Radiative and kinetic
feedback within the NLR and associated shocks would be generated where the
host galaxy gas density is highest along host dust lanes often located in a
disk. Therefore, targets that exhibit extended radio structures in the VLA
would exhibit more interactions between the AGN and its host medium compared
to targets with compact or unresolved VLA morphologies. We suggest that these
interactions likely translate to parsec-scale levels, where the extranuclear
host interactions may be unresolved into VLA peak-flux measurements and
resolved out of VLBA peak-flux measurements.

Figure 1: Left: Composite image showing the central region of NGC 1068 where
red is _HST_ broadband WFPC2 filter F547M, green is the continuum subtracted
_HST_ [O III], and blue is the VLA A-array 5 GHz continuum map. Right:
Composite image showing the VLBA 6 GHz radio continuum data from Fischer et
al. (2021) for NGC 1068 (yellow) with VLBA 5 GHz radio continuum data from
Gallimore et al. (2004) (cyan).

This extranuclear contamination of VLA fluxes was best highlighted in NGC
1068. Figure 1 includes the employed VLA and VLBA observations from Fischer et
al. (2021), where the measured VLA peak flux is largely due to discrete
extranuclear VLBA emission knot ‘C’ north of the nucleus, without emission
from the AGN itself (S1; Gallimore et al. 2004). We return to the VLBA
observations of this target in order to further study the characteristics of
these extranuclear knots. Namely, as they are discrete and traceable outside
the nucleus, we can study their change in projected distance over time to
understand their relative transverse velocities. This will allow us to test
whether these knots can be attributed to a relativistic jet. In this
manuscript, we discuss the relative positions of several knots across a time
span of 22 years, comparing observations from Fischer et al. (2021) to those
described by Gallimore et al. (2004). In our analysis of NGC 1068, we use the
Tully-Fisher distance to NGC 1068 of D = 10.1 Mpc (Tully et al., 2009) and a
scale of 48.9 pc arcsec-1.

Knot | Position Offset | Projected Velocity |  $i$ at $V=0.1c$aaInclination out of the plane of the sky required for a deprojected velocity of 0.1 $c$ | Velocity PA | Sν(1997)bbIntegrated Flux measurements from Gallimore et al. (2004) | Sν(2019)
---|---|---|---|---|---|---
| (mas) | (pc) | (v/c) | (km s-1 ) | (degrees) | (degrees) | (mJy) | (mJy)
S1 |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | 9.1$\pm$0.8 | 2.6$\pm$0.2
NE | 3.29 | 0.16 | 0.024 | 7124 | 76.1 | 116.6 | 36.7$\pm$4.1 | 7.5$\pm$0.6
C | 6.09 | 0.30 | 0.044 | 13358 | 63.9 | 336.4 | 19.4$\pm$1.6 | 5.7$\pm$0.5
Table 1: VLBA Extranuclear Projected Distance, Velocity, and Integrated Flux
Measurements.

##  2 Data Analysis

The VLBA data shown in the right panel of Figure 1 were observed 22 years
apart. We define Epoch 1 as the data from Gallimore et al. (2004), which we
obtained from J.F. Gallimore through private correspondence and Epoch 2 as the
data from Fischer et al. (2021). For Epoch 1, we rebinned the image pixel
scale to match the 0.8 mas pixels in Epoch 2.

The nuclear emitting region in NGC 1068, labeled S1 in Figure 1, as measured
by Gallimore et al. (2004), exhibits water maser emission and likely marks the
location of the central SMBH (Gallimore et al. 2001 and references therein).
Previous analysis by Gallimore et al. (2004) identified three additional
extranuclear emitting regions, NE, C, and S2. Knots S1, NE, and C are clearly
visible in both epochs. However, knot S2 in Epoch 1 does not have a
corresponding knot at its approximate location in Epoch 2. According to
Gallimore et al. (2004), the integrated flux of source S2 was 4.1 $\pm$ 0.4
mJy at 5 GHz, thus, at the sensitivity of $\sim$22 $\mu$Jy/bm of the
observations from Fischer et al. (2021) for NGC 1068, we would expect to have
detected this source at a signal-to-noise greater than 180. Thus, this source
has either dropped in intensity such that it is below the sensitivity
threshold presented in Fischer et al. (2021), or, it has moved position. A
potential corresponding emitter in Epoch 2 is north of the S2 position in
Epoch 1, as shown in Figure 1, as it is the brightest adjacent feature,
roughly in the same radial path, and no emitting knots are observed south of
the S2 Epoch 1 position in Epoch 2. However, as the relationship between these
structures is unclear, we do not perform analyses for the distance between
these knot positions.

###  2.1 VLBA Correlations

When the data from 1997 are directly compared to the data from 2019, none of
the sources align astrometrically with one another and are offset by a few
milli-arcseconds in the east-west direction. We investigated the possibility
that the observed offset was caused by the phase calibrator used, but we have
ruled this out. This is because we used the same calibrator, IERS
0237$-$027/J0239$-$0234, for both datasets in 1997 and 2019. This target is a
defining source in the most recent International Celestial Reference Frame
(ICRF3), which means it is one of the most astrometrically stable sources in
the catalog. Additionally, the astrometric positional error of this
calibrator, $\sigma_{\alpha}=\pm 2.09\mu$s, $\sigma_{\delta}=\pm 31.9\mu$as,
is lower than the observed offset between the two epochs. Through private
communications with NRAO, we learned that the delay models for the VLBA were
not well-constrained in 1997. These delay models have been improved greatly
over time and thus, the astrometry for the observations taken in 2019 is
likely more robust. Therefore, to begin our analysis of the features probed by
the VLBA, we align the two datasets to the S1 structure as we believe this is
the location of the AGN core and, therefore, likely static in position over
time.

In order to correct for relative offsets in nuclear position between epochs,
we perform auto- and cross-correlation analyses on the position of knot S1.
This analysis was performed using a large, $200\times 200$ pixel ($0\farcs
16\times 0\farcs 16$) base image and a smaller, $25\times 25$ pixel ($0\farcs
02\times 0\farcs 02$) subimage surrounding the S1 knot in both Epoch 1 and
Epoch 2 using identical Epoch 2 ICRS coordinates. The subimage of S1 in Epoch
1 was auto-correlated against its base image to determine its Epoch 1 position
in ICRS coordinates, defined as the auto-correlation distribution peak. We
then cross-correlate the subimage of S1 knot from Epoch 1 against the base
image from Epoch 2 and define the Epoch 2 position of the knot as the location
of the cross-correlation distribution peak. The difference in relative
positions between the auto- and cross-correlated distribution peaks then
provides the projected distance offset between epochs. The S1 knot was aligned
between epochs by altering the ICRS coordinates of the Epoch 1 data such that
the distance between the Epoch 2 auto-correlation peak and cross-correlation
peak is zero.

With the S1 knot from both epochs aligned, we repeated the cross-correlation
measurements for the C and NE knots to measure their differences in position
over time. We provide this difference for each knot in Table 1. From our
measurements of the aligned images, we find angular distances between
correlation peaks for the C and NE knots of 6.09 and 3.29 mas, respectively.
At the given distance for NGC 1068, these differences correspond to a few
tenths of a parsec. Knowing the difference in time between observation epochs,
$t=6.93\times 10^{8}$ s, we can calculate the projected velocities of knots C
and NE to be $\approx$ 13300 and 7100 km s-1, respectively, a few percent of
the speed of light. Using the position offset between the two epochs, we can
also measure the position angle of the velocity vector in the plane of the
sky, with Knot C traveling northwest at a position angle of $336\fdg 4$ east
of north and Knot NE traveling southeast at a position angle of $116\fdg 6$
east of north. We note that these vectors are roughly perpendicular to the
position angle of the overall radio structure in this system, PA $\approx
30\arcdeg$, as illustrated in Figure 1.

##  3 Results and Discussion

At projected values between approximately $7,000-14,000$ km s-1, the
velocities observed for knots NE and C are $<10\%$ the speed of light. The
magnitude of these velocities is relatively unusual in terms of AGN feedback
processes. As described above, jets in kinetic-mode feedback are relativistic,
traveling at a large fraction of the speed of light. This is highlighted in
observations by the MOJAVE program (Lister et al., 2016), which illustrate the
relativistic motion of radio-emitting knots in radio-loud quasars over several
decades. This is also highlighted in theoretical models (Mukherjee et al.,
2018) where jets interacting with and affecting dense, gaseous disc material
also require highly relativistic ($v/c>0.4$) jets within the inner $\sim 200$
pc to couple with the gas and launch outflows. On the other hand, optical
winds in radiative-mode feedback are much slower (V${}_{max}<$ 2000 km s-1;
Fischer et al. 2013, 2018), largely restricted by the increasing enclosed mass
of the host galaxy as a function of outflow radius. As our measurements are
projected values, we also calculate the inclination of the velocity vector
that is required for knots NE and C to reach 0.1c, provided in Table 1, for a
general understanding on the likelihood that we are observing projected,
relativistic motion. Knot C requires an inclination of $\sim 64\arcdeg$ from
the plane of the sky to reach 0.1c and Knot NE would require a larger
inclination of $\sim 76\arcdeg$. Inclinations required to reach velocities of
$\sim$0.2c for either knot are greater still, between $78\arcdeg-84\arcdeg$
away from the plane of the sky, which suggests that the knots would need
trajectories near parallel with our line of sight and perpendicular to the
radiative outflows to be traveling at a significant fraction of the speed of
light.

Further, knots NE and C do not exhibit velocity vectors traveling radially
away from the AGN, illustrated in Figure 2, as one would expect through AGN
feedback processes. Paths to the most extended radio emission to the north and
south of the AGN, visible in Figure 1, are also shown in Figure 2 which are
approximately in the same direction as the normal of the radiation bicone (PA
$\sim 30\arcdeg$; Das et al. 2006). We measure the position angles of the knot
velocity vectors, provided as angles east of north, to be perpendicular to
these axes.

Figure 2: Epoch 2 VLBA 6 GHz radio continuum data with archival VLA A-array
8.49 GHz continuum imaging contours (green) overlaid. Knots NE and C in Epoch
2 are circled in white with arrows pointing in the direction of their measured
velocity vectors (not to scale). The central AGN is circled in yellow. Cyan,
dashed lines point to the apex positions of the extended radio sources shown
in Figure 1. The kinematic major axis of the rotating molecular gas, per
García-Burillo et al. (2019), is shown as blue and red lines, with rotation
traveling counterclockwise on the plane of the sky.

If we define a radio jet to exhibit relativistic velocities, we set 0.1$c$ as
a lower limit for relativistic motion, as this velocity is where one is
typically required to begin accounting for relativistic effects. Comparing the
observed velocities and trajectories to this definition, we do not find the
observed radio emission to exhibit the relativistic velocities present in
radio jets. As described earlier, multi-epoch observations of jetted
structures in radio-loud AGN have shown the velocities to move at relativistic
velocities in a radial motion away from the central AGN. Jet kinematics can be
more complicated at smaller scales. M87, a nearby radio-loud AGN, exhibits an
intricate, polarized double-helix morphology in its radio jet (Pasetto et al.,
2021). One could posit that the velocity vectors observed in NGC 1068 are due
to a similar scenario; however several clues in our observations suggest this
helical, spiral pattern is not occurring. Immediately, the observed motions
are not in a tight and orderly pattern as seen in jets. Also, we know that the
observed knots are not individual point sources or narrow streams that could
be traveling along magnetic field lines. As traced by VLBA, MERLIN, and VLA
radio observatories at decreasing resolution (Figures 1, 2; see also Gallimore
et al. 2004), the observed extranuclear knots are instead peaks of much larger
structures.

###  3.1 Scenarios to Produce Observed Knot Velocities

Assuming the production of radio emission is due to shocks, there are several
possible candidates for the responsible radiative-mode winds. First, the [O
III] outflows themselves. Kinematic analyses (Das et al., 2006; Venturi et
al., 2021), show that the FWHMs near the radio knots are large, greater than
1000 km s-1, but there is no direct evidence for [O III] enhancement
specifically at the knot locations, per Figure 3, described below.
Alternatively, the shocks could be created via X-ray winds. Work by Kraemer et
al. (2020) and Trindade Falcão et al. (2021) revealed X-ray emission-line
regions with masses and mass-outflow rates comparable to, or greatly
exceeding, those of the optical emission-line gas, providing an order of
magnitude larger kinetic power. Finally, ultra-fast outflows (UFOs; Tombesi et
al. 2010) reach velocities $>$ 0.2 $c$, which should have sufficient energy to
cause strong shocks (Maksym et al., 2022).

Although the progenitors of these shocks potentially have relativistic
velocities, we suggest that the observed difference in position between epochs
does not require the individual knots to travel $\sim$10000 km s-1 between the
observed epochs. Instead, as the relationship between the observed radio-
continuum emission and host medium is in a turbulent and dynamic environment,
the impact point at which AGN feedback creates these shocks is likely what
changes over time. In this scenario, the interaction between AGN winds and the
surrounding environment changes the morphology of the dense molecular gas
reservoir observed in each epoch. As gas in this reservoir is ionized and
driven away, the remaining gas is sculpted to possess a new morphology. Winds
then impact and shock different surfaces of the dense host medium due to these
ablative processes to create the apparent motion of the radio emission.

Figure 3: Multi-band imaging of the nuclear region of NGC 1068, centered on
the approximate location of the AGN. First row: Continuum-subtracted ALMA Band
7 observations of CO(3-2) 345.796 GHz emission with contoured comparisons of
[Fe II] 1.644$\mu$m (second column; red), [O III] $\lambda$5007 (third column;
green), and 8.49 GHz radio continuum (fourth column; blue) features. Second
row: Continuum-subtracted Gemini NIFS IFU imaging of [Fe II] 1.644$\mu$m
emission with contoured comparisons of [O III] $\lambda$5007 (second column;
green), and 8.49 GHz radio continuum (third column; blue) features. Third row:
Continuum-subtracted HST WFPC2 imaging of [O III] $\lambda$5007 emission with
a comparison to 8.49 GHz radio continuum (second column; blue). Bottom row:
VLA A-Array imaging of 8.49 GHz radio continuum. Outer contours of [Fe II], [O
III], and radio continuum imaging represent S/N of 12$\sigma$, 100$\sigma$,
and 10$\sigma$, respectively. North is up and east is to the left.

###  3.2 Multiwavelength Nuclear Imaging

To illustrate how radio emission forms within the host environment, we compare
its structure to adjacent structures in other wavebands. Figure 3 compares the
morphologies and locations of different forms of emission surrounding the
nucleus of NGC 1068 across several wavebands and observatories; ALMA Band 7
observations of CO($3-2$) (García-Burillo et al., 2019), Gemini/NIFS H-band
observations of [Fe II] 1.6$\mu$m (Riffel et al., 2014; Barbosa et al., 2014),
HST/WFPC2 narrow-band imaging of [O III] $\lambda$5007 (Schmitt et al., 2003;
Fischer et al., 2018), and VLA A-array observations of 8.49 GHz radio
continuum (observed 8 September 1999, NRAO/VLA Archive Survey), which provides
a large-scale comparison to our VLBA observations, as shown in Figure 2. ALMA
and Gemini observations were aligned using the location of the AGN provided in
previous analyses (Riffel et al., 2014; García-Burillo et al., 2019). HST [O
III] imaging alignment is described in Appendix A.

We have previously analyzed connections between ionized gas and molecular gas
(Fischer et al., 2017) and shown that molecular gas is the reservoir from
which spatially-resolved, ionized-gas structures are formed. We observe this
relationship again in NGC 1068, as [O III] is largely observed north of the
AGN, and peaks of [O III] emission-line knots can be mapped to CO reservoir
walls facing the NLR axis east and west of the nucleus. Kinematic measurements
in these regions suggest the gas is highly turbulent, with complex emission-
line profiles and large mass outflows (Das et al., 2006; Revalski et al.,
2021; Venturi et al., 2021), suggestive of an interaction where radiative
driving by the AGN is carving out significant portions of the molecular gas
reservoir.

Comparing the position of the radio emission knots to the [O III] NLR, we see
that the radio is intertwined with the ionized gas, as observed in other
nearby radio-quiet AGN. However, as seen in Figure 2, if the peak of the most
extended radio structure is used as an axis of feedback, the nuclear radio
emission is not evenly distributed. Specifically, north of the nucleus, we see
knots NE and C are located along the right side of the NLR against the western
CO reservoir and not the left, eastern reservoir. We do see that extended
radio emission to the south of the nucleus is located along the left side of
the NLR, against the eastern CO reservoir, south of the nucleus. Separate from
the VLBA detected knots, there is also a radio knot in Figure 3 south of the
nucleus that also overlaps with [O III] emission and is adjacent to CO
emitting gas on the right, western side of the reservoir.

Comparing these structures to the near-IR [Fe II] emission, we see a conical
morphology similar to that of the [O III] emission to the north of the AGN,
which originates from between the two CO reservoirs east and west of the
nucleus. The flux distribution of the [Fe II] is closer to that of the radio
continuum emission, with the brightest [Fe II] knot aligned with the right,
western [O III] emission knots. Further, south of the AGN, the [Fe II]
emission also wraps around the left, eastern CO reservoir, similar to the
observed [O III] and radio-continuum morphologies. Finally, we note that the
flux distributions across wavebands are fairly consistent; brighter emission
in the north and fainter emission to the south, which is consistent with all
of the observed features being related and adjacent to one another.

In concert, this intertwined emission across the electromagnetic spectrum
suggests a scenario where the observed radio knots in NGC 1068 form within the
optical NLR, as traced by the observed [O III] emission, which represents the
intersection between AGN ionization and host galaxy disk material. This
analysis supports the scenario where the observed radio emission is therefore
a byproduct that occurs only at locations where winds are producing
significant shocks in the radio-quiet AGN host galaxy, similar to what is seen
in supernova remnants (Zakamska & Greene, 2014; Fischer et al., 2019; Smith et
al., 2020; Fischer et al., 2021; Glikman et al., 2022). Analysis of high-
frequency continuum emission mapped by ALMA of these knots finds their
spectral indexes to be negative ($\alpha\leq-1$) throughout, as expected for
synchrotron-like emission (García-Burillo et al., 2019). From our multi-band
comparisons in Figure 3, we see that the distribution of the brightest
synchrotron radio continuum is not filling in the entirety of the NLR and is
concentrated on the western side to the north of the AGN and eastern side to
the south. This is also where we see the most prominent [Fe II] emission,
which is an excellent diagnostic for shock-induced ionization (Sturm et al.,
2002; Inami et al., 2013), and is frequently observed to be co-located with
radio continuum emission in radio-quiet AGN (Storchi-Bergmann et al., 2010;
Riffel & Storchi-Bergmann, 2011; Riffel et al., 2013, 2015). While the direct
relationship is not clearly understood, shocks may create significant amounts
of localized, low-energy ionizing radiation that permeates deeper into gas
reservoirs than photons with ionization potentials greater than 13.6 eV, which
are absorbed by hydrogen at shallower depths.

In NGC 1068, the radio-continuum emission and [Fe II] flux distribution, and
thus shock distribution, are likely uneven across the NLR because these
structures are formed in the densest gas reservoirs that are often observed
entering the radiation field of the AGN. Analysis of the rotating molecular
gas by García-Burillo et al. 2019 shows that the molecular gas reservoir is
rotating counterclockwise, such that the right, western reservoir north of the
AGN and the left, eastern reservoir south of the AGN are rotating into the AGN
radiation field. Therefore, the regions where we observe radio-continuum and
[Fe II] emission in NGC 1068 are likely often the face of the relatively
pristine molecular gas rotating into the AGN radiation field and being
impacted by AGN winds.

The nearby AGN in NGC 4151 exhibits a similar case of gas reservoirs rotating
into AGN radiation fields (Storchi-Bergmann et al., 2010), with H2 2.12 $\mu$m
lanes rotating clockwise into the observed AGN radiation fields and producing
[Fe II] and radio-continuum emission along the impacted face of the gas lane.
Relationships between [Fe II] and radio-continuum observed in NGC 1068, NGC
4151, and others suggest that the radio structures are often aligned with host
disk material, or significant gas lanes, out to distances of $\sim$1000 pc. As
another comparison with radio jet structures, jets do not exhibit any such
preference in orientation and are not required to be adjacent to or
intertwined with host disk material. Therefore, if the radio emission in NGC
1068 and other radio-quiet AGN were due to a radio jet, it is highly
coincidental that they are continuously adjacent to line-emission structures
which aligns their position to the plane of the host galaxy or significant
dust lanes.

###  3.3 Extended Radio Emission in NGC 1068

As both epochs of VLBA observation were focused on the nuclear regions of NGC
1068, analysis of interactions for the most extended radio emission,
particularly the pointed, lobe-like structure to the northeast of the nucleus
visible in Figure 1, are outside the scope of this work. However, we suggest
that the same underlying physics responsible for the nuclear radio structure
produces these features as well.

Kinematic analysis of [O III] across the inner arcminute of NGC 1068 with
VLT/MUSE (Venturi et al., 2021) shows that the line width of AGN-ionized [O
III] gas is large (W${}_{70}\sim$ 1000 km s-1) at small radii and drops
significantly (W${}_{70}\sim$ 150 km s-1) at the base of the lobe-like
structure to larger radii. However, the [O III] emission morphology is
continuous in intensity between the nucleus and the narrow-line emission to
the northeast. This discrepancy between kinematics and morphology does not
suggest a jet interaction, as we would expect to see a peak in FWHM at the
apex of the jet lobe. Instead, we suggest that the more extended radio
emission resides in a shielded or offset region separate from the inner,
turbulent emission region. We hypothesize that the radius where we observe the
difference in line width represents the inner radius of a spiral arm that is
being penetrated by high-velocity outflows. These winds are launched from
smaller radii, which we observe in the regions sampled by our VLBA
observations, and travel radially into the host lane material, ionizing and
excavating it in the process (i.e. Meena et al. (2023)). Material is not
launched at angles that reach above the spiral arm lane. Therefore, the
ionized gas we observe at the larger distances is not disturbed and remains in
rotation, producing the observed narrow emission lines. As the excavation site
is a cavity, the optical ionized gas counterpart of this process is obscured
by the overlying dust lane along our line of sight, preventing us from
observing the correspondingly large optical line-width emission. These winds
also shock the interior surface of the cavity, forming the pointed, lobe-like
structure we observe.

Further analysis is required to examine these hypotheses. Specifically, near-
and mid-infrared observations of the extended radio structures using integral
field spectroscopy instruments such as JWST/NIRSpec, JWST/MIRI, and
Gemini/NIFS would test for the presence of shocks and warm molecular gas
adjacent to the edges of the radio emission as we see near the nucleus.

###  3.4 Radio Loudness at High Redshifts

Fischer et al. (2021) found that radio sources observed with both VLA and VLBA
exhibited discrepancies between peak flux values, suggesting that point
sources observed with larger beams are hosting a significant accumulation of
extranuclear emission at small radii, potentially due to shock processes
similar to what we suggest is occurring in NGC 1068. Unresolved radio sources
at high redshift may be exhibiting similar processes, affecting our
understanding of distant AGN with compact radio sources in two ways. First,
AGN classified as radio-loud may actually be radio-quiet. The hardness ratio
used in this classification, $R_{X}$, compares the X-ray and radio emission
originating in the central engine. AGN emitting as a radio point source, with
current resolution capabilities, may be interacting with their host material,
producing shocks which in turn produce extranuclear radio emission unrelated
to the central engine. Second, variability due to environmental interactions
appears to occur on observable timescales. Knot S2 observed in Epoch 1 of our
VLBA observations either disappears or is replaced by a radially interior
radio knot, suggesting that structures can move in to, or out of, the ionizing
radiation field and affect the integrated radio flux of the circumnuclear
region. At large redshifts, this change in structure morphology would be lost
within the instrument beam, and appear as an increase or drop in flux.
Alternatively, these variations in morphology may also affect the astrometric
position of quasar radio sources with new peak flux positions rising over
time.

We suggest that further study is required on the incidence of significant [Fe
II] emission in nearby AGN with extended radio sources and higher redshift AGN
with unresolved radio sources. Determining if a relationship or specific ratio
exists in nearby AGN where interactions with their host disk occur would
produce a useful diagnostic in determining the source of radio emission in
unresolved sources at higher redshift.

###  3.5 Further Impacts on Radio Survey Studies

Often in recent studies of extended radio structures in radio-quiet AGN,
measurements are performed to quantify the amount of mechanical feedback that
is injected into the host environment (i.e., Jarvis et al. 2019; Villar Martín
et al. 2021; Girdhar et al. 2022). Findings from these studies suggest that
mechanical feedback is not efficient in comparison to radiative-mode feedback.
However, these radio signatures do trace the impact of the most significant
winds on the host galaxy. Further, these locations represent the terminal
radii where AGN winds are no longer able to drive material away from the
central engine due to obstruction via dense host material and may be
compressing gas in these regions, promoting an increase in star formation. As
such, surveying the radio morphologies of radio-quiet AGN in search of
extended structures can probe the incidence of compacting, positive feedback.
Developing a census of how often AGN are significantly interacting with their
hosts and the size of said interactions would then be available for inclusion
in AGN feedback models.

##  4 Conclusions

We compared VLBA 5 GHz radio observations taken 22 years apart for radio-quiet
Seyfert 2 AGN NGC 1068 and measured the projected distance between two
extranuclear flux peaks, knots NE and C, surrounding the central engine. Our
main conclusions are as follows:

  1. 1.

Knots NE and C are well defined across both epochs, exhibit projected distance
offsets $<$ 0.5 pc, and travel in directions perpendicular to the general NLR
structure. Across the timescale between the two observation epochs, the
measured velocities translate to non-relativistic velocities that are a few
percent the speed of light and would require paths along our line of sight to
reach significant fractions of the speed of light.

  2. 2.

The location of the nuclear radio continuum is shown to have an intricate
connection with surrounding molecular and ionized gas, including co-location
with molecular gas reservoirs that are rotating into the NLR ionization field
and low-ionization [Fe II] emission.

  3. 3.

The measured velocities and trajectories of the observed radio flux peaks do
not perform as expected if the emission originated from a relativistic radio
jet. We instead suggest that the flux peaks are not individual, jetted knots
of radio plasma, but peak emission of larger structures whose locations change
with host gas reservoir morphologies over observable time periods.

  4. 4.

These findings continue to be consistent with the explanation that extended
radio structures observed in radio-quiet AGN are not jets, but byproducts
formed via localized shocks, which in turn produce the observed synchrotron
radiation at these locations.

The authors thank the anonymous referee for their helpful comments that
improved the clarity of this work. TCF is thankful for the support of the
European Space Agency (ESA). The authors thank Zo Chapman, Julia Falcone,
Enrique Lopez Rodriguez, Beena Meena, Garrett Polack, Mitchell Revalski, Maura
Shea, and Krista Lynne Smith for helpful discussions while developing this
manuscript. TCF thanks Jack Gallimore for providing the Epoch 1 VLBA
observations and Rogemar Riffel for providing the Gemini/NIFS observations
used in this work.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Barbosa et al. (2014) Barbosa, F. K. B., Storchi-Bergmann, T., McGregor, P., Vale, T. B., & Rogemar Riffel, A. 2014, MNRAS, 445, 2353, doi: 10.1093/mnras/stu1637
  * CASA Team et al. (2022) CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501, doi: 10.1088/1538-3873/ac9642
  * Das et al. (2006) Das, V., Crenshaw, D. M., Kraemer, S. B., & Deo, R. P. 2006, AJ, 132, 620, doi: 10.1086/504899
  * Fabian (2012) Fabian, A. C. 2012, ARA&A, 50, 455, doi: 10.1146/annurev-astro-081811-125521
  * Falcke et al. (1998) Falcke, H., Wilson, A. S., & Simpson, C. 1998, ApJ, 502, 199, doi: 10.1086/305886
  * Fischer et al. (2019) Fischer, T., Smith, K. L., Kraemer, S., et al. 2019, ApJ, 887, 200, doi: 10.3847/1538-4357/ab55e3
  * Fischer et al. (2013) Fischer, T. C., Crenshaw, D. M., Kraemer, S. B., & Schmitt, H. R. 2013, ApJS, 209, 1, doi: 10.1088/0067-0049/209/1/1
  * Fischer et al. (2017) Fischer, T. C., Machuca, C., Diniz, M. R., et al. 2017, ApJ, 834, 30, doi: 10.3847/1538-4357/834/1/30
  * Fischer et al. (2018) Fischer, T. C., Kraemer, S. B., Schmitt, H. R., et al. 2018, ApJ, 856, 102, doi: 10.3847/1538-4357/aab03e
  * Fischer et al. (2021) Fischer, T. C., Secrest, N. J., Johnson, M. C., et al. 2021, ApJ, 906, 88, doi: 10.3847/1538-4357/abca3c
  * Gallimore et al. (2006) Gallimore, J. F., Axon, D. J., O’Dea, C. P., Baum, S. A., & Pedlar, A. 2006, AJ, 132, 546, doi: 10.1086/504593
  * Gallimore et al. (2004) Gallimore, J. F., Baum, S. A., & O’Dea, C. P. 2004, ApJ, 613, 794, doi: 10.1086/423167
  * Gallimore et al. (2001) Gallimore, J. F., Henkel, C., Baum, S. A., et al. 2001, ApJ, 556, 694, doi: 10.1086/321616
  * García-Burillo et al. (2019) García-Burillo, S., Combes, F., Ramos Almeida, C., et al. 2019, A&A, 632, A61, doi: 10.1051/0004-6361/201936606
  * Girdhar et al. (2022) Girdhar, A., Harrison, C. M., Mainieri, V., et al. 2022, MNRAS, 512, 1608, doi: 10.1093/mnras/stac073
  * Glikman et al. (2022) Glikman, E., Lacy, M., LaMassa, S., et al. 2022, ApJ, 934, 119, doi: 10.3847/1538-4357/ac6bee
  * Hamer et al. (2015) Hamer, S., Salomé, P., Combes, F., & Salomé, Q. 2015, A&A, 575, L3, doi: 10.1051/0004-6361/201424808
  * Hardcastle & Croston (2020) Hardcastle, M. J., & Croston, J. H. 2020, New A Rev., 88, 101539, doi: 10.1016/j.newar.2020.101539
  * Inami et al. (2013) Inami, H., Armus, L., Charmandaris, V., et al. 2013, ApJ, 777, 156, doi: 10.1088/0004-637X/777/2/156
  * Jarvis et al. (2019) Jarvis, M. E., Harrison, C. M., Thomson, A. P., et al. 2019, MNRAS, 485, 2710, doi: 10.1093/mnras/stz556
  * Kraemer et al. (2020) Kraemer, S. B., Turner, T. J., Couto, J. D., et al. 2020, MNRAS, 493, 3893, doi: 10.1093/mnras/staa428
  * Lister et al. (2016) Lister, M. L., Aller, M. F., Aller, H. D., et al. 2016, AJ, 152, 12, doi: 10.3847/0004-6256/152/1/12
  * Maksym et al. (2022) Maksym, W. P., Elvis, M., Fabbiano, G., et al. 2022, arXiv e-prints, arXiv:2209.14326, doi: 10.48550/arXiv.2209.14326
  * Maksym et al. (2017) Maksym, W. P., Fabbiano, G., Elvis, M., et al. 2017, ApJ, 844, 69, doi: 10.3847/1538-4357/aa78a4
  * Meena et al. (2023) Meena, B., Crenshaw, D. M., Schmitt, H. R., et al. 2023, ApJ, 943, 98, doi: 10.3847/1538-4357/aca75f
  * Morganti (2017) Morganti, R. 2017, Frontiers in Astronomy and Space Sciences, 4, 42, doi: 10.3389/fspas.2017.00042
  * Morganti et al. (2013) Morganti, R., Fogasy, J., Paragi, Z., Oosterloo, T., & Orienti, M. 2013, Science, 341, 1082, doi: 10.1126/science.1240436
  * Mukherjee et al. (2018) Mukherjee, D., Bicknell, G. V., Wagner, A. Y., Sutherland, R. S., & Silk, J. 2018, MNRAS, 479, 5544, doi: 10.1093/mnras/sty1776
  * Murthy et al. (2022) Murthy, S., Morganti, R., Wagner, A. Y., et al. 2022, Nature Astronomy, 6, 488, doi: 10.1038/s41550-021-01596-6
  * Oh et al. (2018) Oh, K., Koss, M., Markwardt, C. B., et al. 2018, ApJS, 235, 4, doi: 10.3847/1538-4365/aaa7fd
  * Paggi et al. (2012) Paggi, A., Wang, J., Fabbiano, G., Elvis, M., & Karovska, M. 2012, ApJ, 756, 39, doi: 10.1088/0004-637X/756/1/39
  * Pasetto et al. (2021) Pasetto, A., Carrasco-González, C., Gómez, J. L., et al. 2021, ApJ, 923, L5, doi: 10.3847/2041-8213/ac3a88
  * Revalski et al. (2021) Revalski, M., Meena, B., Martinez, F., et al. 2021, ApJ, 910, 139, doi: 10.3847/1538-4357/abdcad
  * Riffel & Storchi-Bergmann (2011) Riffel, R. A., & Storchi-Bergmann, T. 2011, MNRAS, 411, 469, doi: 10.1111/j.1365-2966.2010.17721.x
  * Riffel et al. (2015) Riffel, R. A., Storchi-Bergmann, T., & Riffel, R. 2015, MNRAS, 451, 3587, doi: 10.1093/mnras/stv1129
  * Riffel et al. (2013) Riffel, R. A., Storchi-Bergmann, T., & Winge, C. 2013, MNRAS, 430, 2249, doi: 10.1093/mnras/stt045
  * Riffel et al. (2014) Riffel, R. A., Vale, T. B., Storchi-Bergmann, T., & McGregor, P. J. 2014, MNRAS, 442, 656, doi: 10.1093/mnras/stu843
  * Schmitt et al. (2003) Schmitt, H. R., Donley, J. L., Antonucci, R. R. J., Hutchings, J. B., & Kinney, A. L. 2003, ApJS, 148, 327, doi: 10.1086/377440
  * Shimizu et al. (2019) Shimizu, T. T., Davies, R. I., Lutz, D., et al. 2019, MNRAS, 490, 5860, doi: 10.1093/mnras/stz2802
  * Silk (2013) Silk, J. 2013, ApJ, 772, 112, doi: 10.1088/0004-637X/772/2/112
  * Silk & Rees (1998) Silk, J., & Rees, M. J. 1998, A&A, 331, L1, doi: 10.48550/arXiv.astro-ph/9801013
  * Smith et al. (2020) Smith, K. L., Koss, M., Mushotzky, R., et al. 2020, ApJ, 904, 83, doi: 10.3847/1538-4357/abc3c4
  * Storchi-Bergmann et al. (2010) Storchi-Bergmann, T., Lopes, R. D. S., McGregor, P. J., et al. 2010, MNRAS, 402, 819, doi: 10.1111/j.1365-2966.2009.15962.x
  * Sturm et al. (2002) Sturm, E., Lutz, D., Verma, A., et al. 2002, A&A, 393, 821, doi: 10.1051/0004-6361:20021043
  * Tombesi et al. (2010) Tombesi, F., Cappi, M., Reeves, J. N., et al. 2010, A&A, 521, A57, doi: 10.1051/0004-6361/200913440
  * Trindade Falcão et al. (2021) Trindade Falcão, A., Kraemer, S. B., Fischer, T. C., et al. 2021, MNRAS, 505, 3054, doi: 10.1093/mnras/stab1549
  * Tully et al. (2009) Tully, R. B., Rizzi, L., Shaya, E. J., et al. 2009, AJ, 138, 323, doi: 10.1088/0004-6256/138/2/323
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
  * van Moorsel et al. (1996) van Moorsel, G., Kemball, A., & Greisen, E. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 37
  * Venturi et al. (2021) Venturi, G., Cresci, G., Marconi, A., et al. 2021, A&A, 648, A17, doi: 10.1051/0004-6361/202039869
  * Villar Martín et al. (2021) Villar Martín, M., Emonts, B. H. C., Cabrera Lavers, A., et al. 2021, A&A, 650, A84, doi: 10.1051/0004-6361/202039642
  * Whittle & Wilson (2004) Whittle, M., & Wilson, A. S. 2004, AJ, 127, 606, doi: 10.1086/380940
  * Zakamska & Greene (2014) Zakamska, N. L., & Greene, J. E. 2014, MNRAS, 442, 784, doi: 10.1093/mnras/stu842
  * Zubovas et al. (2013) Zubovas, K., Nayakshin, S., King, A., & Wilkinson, M. 2013, MNRAS, 433, 3079, doi: 10.1093/mnras/stt952

##  Appendix A HST [O III] Alignment

To investigate the relationship between the AGN in NGC 1068 and its
surrounding environment in exquisite high spatial resolution, we utilized a
myriad of archival data spanning a wide range of wavelengths across the
electromagnetic spectrum, each of which probe various aspects of the AGN and
its feedback on its surrounding medium. In this work, we used Hubble Space
Telescope optical [O III] narrow line emission as part of a suite of emission
structures to suss out the source of the observed offsets of the peak emission
features seen on the highest resolution VLBA spatial scales. In this section,
we describe our treatment of the data to account for any astrometry offsets
that might exist between the various datasets.

###  A.1 HST [O III] Data

We obtained an [O III] emission line image from the HST data archive as well
as a broad spectrum continuum emission image that we used to subtract
continuum emission sources from the [O III] image, thus leaving only the
narrow band [O III] line emission in the central region of NGC 1068. However,
before applying the continuum subtraction, we used the unsubtracted [O III]
image to astrometrically align the continuum point sources in the central
region of the galaxy to Gaia EDR3. Once the astrometry was corrected, we then
subtracted off the continuum emission from the [O III] image and applied the
astrometric solution to the resulting narrow band [O III] emission line image.

####  A.1.1 HST Astrometric Alignment

The native World Coordinate System (WCS) solution of the HST data is not
sufficient for comparison with the high angular resolution and astrometric
precision of the VLA/VLBA data, with astrometric offsets of 2–3 pixels
(100–150 mas) visible by eye. We corrected the WCS solution of the HST data
using the following procedure:

  1. 1.

Query the Gaia EDR3 catalog for all sources overlapping with the HST image and
plot them as region files.

  2. 2.

For Gaia sources associated with a compact knot of emission, manually center
the region file on the knot. Delete remaining Gaia sources associated with
extended emission. This produced 46 Gaia sources with compact knots in the HST
image.

  3. 3.

Fit two-dimensional Gaussians to the Gaia associations to determine the
precise centroid coordinates $(x,y)$ and convert to $(\alpha,\delta)$ using
the original WCS solution.

  4. 4.

Perform a least-squares 6-parameter fit (translation, rotation, scale) to the
Gaia coordinates using the set of $(\alpha,\delta)$ from (3). To mitigate the
effect of astrometric outliers, the 6-parameter fit was carried out in an
iterative fashion, at each iteration using a random subset (without
replacement) of the 46 Gaia sources from (2). We used subsets of size 16,
which was chosen as a balance between maximizing the degrees of freedom of the
fit and minimizing outliers. New fit coordinates $(\alpha,\delta)$ are
calculated at each iteration, as well as the astrometric residuals
$\Delta\alpha=(\alpha-\alpha_{\mathrm{Gaia}})\cos{\delta}$,
$\Delta\delta=\delta-\delta_{\mathrm{Gaia}}$. If the absolute mean values of
$\Delta\alpha$, $\Delta\delta$ are less than their standard errors and the
maximum value of the total offset
$\sqrt{(\Delta\alpha)^{2}+(\Delta\delta)^{2}}$ is less than that of the
previous best WCS solution, the WCS solution of the iteration is stored as the
new best.

By using the iterative fitting strategy in (4), We were able to achieve high
astrometric precision relative to Gaia EDR3, with a mean offset in $\alpha$
and $\delta$ of $-0.6\pm 1.5$ mas and $-0.9\pm 1.7$ mas, where the
uncertainties are the standard errors of the means. These 16 sources
circumscribe NGC 1068 (Figure 4) so the nuclear region is not affected by
errors of extrapolation, and the 6-parameter fit does not include any tilt or
higher-order distortion terms, which can introduce errors of interpolation.

Figure 4: Locations of Gaia EDR3 astrometric reference points used to correct
the HST astrometry, overlaid on the F547M image used for the counterparts.


$\displaystyle\left|\frac{\|\theta_{*}\|^{2}}{\|\theta_{\sharp}\|^{2}}\frac{2}{n}\sum_{i=1}^{n/2}R_{i}-\rho\right|$
$\displaystyle=\left|\left(1+\frac{\|\theta_{*}\|^{2}-\|\theta_{\sharp}\|^{2}}{\|\theta_{\sharp}\|^{2}}\right)\frac{2}{n}\sum_{i=1}^{n/2}R_{i}-\rho\right|$
(236)
$\displaystyle\leq\left|\frac{2}{n}\sum_{i=1}^{n/2}R_{i}-\rho\right|+\frac{\left|\|\theta_{*}\|^{2}-\|\theta_{\sharp}\|^{2}\right|}{\|\theta_{\sharp}\|^{2}},$
(237)

where the last inequality follows since
$\left|\frac{2}{n}\sum_{i=1}^{n/2}R_{i}\right|\leq 1$. In this last equation,
the last term can be considered a bias of the estimator due to the mismatch
between $\theta_{*}$ and $\theta_{\sharp}$. Now, by Bernstein’s inequality for
a sum of independent, zero-mean, bounded random variables
$R_{i}-\mathbb{E}[R_{i}]\in[-2,2]$ (cf. (323) from [Wainwright, 2019,
Propostion 2.14])

$\mathbb{P}\left[\left|\frac{2}{n}\sum_{i=1}^{n/2}R_{i}-\mathbb{E}[R_{i}]\right|\geq
t\right]\leq
2\exp\left(-\frac{\frac{n}{2}t^{2}}{2\left(4\delta+\frac{2t}{3}\right)}\right),$
(238)

and so with probability larger than $1-\epsilon$

$\displaystyle\left|\frac{2}{n}\sum_{i=1}^{n/2}R_{i}-\mathbb{E}[R_{i}]\right|$
$\displaystyle\leq\frac{4}{3n}\log\left(\frac{2}{\epsilon}\right)+\frac{4}{3n}\sqrt{9n\delta\log\left(\frac{2}{\epsilon}\right)+\log^{2}\left(\frac{2}{\epsilon}\right)}$
(239)
$\displaystyle\leq\sqrt{\frac{16\delta\log\left(\frac{2}{\epsilon}\right)}{n}}+\frac{8\log\left(\frac{2}{\epsilon}\right)}{3n}$
(240) $\displaystyle\leq
7\log\left(\frac{2}{\epsilon}\right)\sqrt{\frac{\delta}{n}},$ (241)

since $\delta\geq\frac{1}{n}$.

The second term in (233) (and similarly, the third term in (233)) is

$\frac{2}{\|\theta_{\sharp}\|^{2}n}\sum_{i=1}^{n/2}S_{2i}\theta_{\sharp}^{\top}Z_{2i-1}\sim
N\left(0,\frac{2}{\|\theta_{\sharp}\|^{2}n}\right).$ (242)

Thus, by the standard Chernoff bound for Gaussian random variables

$\mathbb{P}\left[\left|\frac{2}{\|\theta_{\sharp}\|^{2}n}\sum_{i=1}^{n/2}S_{2i}\theta_{\sharp}^{\top}Z_{2i-1}\right|>t\right]\leq
2e^{-\frac{\|\theta_{\sharp}\|^{2}nt^{2}}{4}},$ (243)

and so with probability larger than $1-\epsilon$

$\left|\frac{2}{\|\theta_{\sharp}\|^{2}n}\sum_{i=1}^{n/2}S_{2i}\theta_{\sharp}^{\top}Z_{2i-1}\right|\leq\sqrt{\frac{4}{n\|\theta_{\sharp}\|^{2}}\log\left(\frac{2}{\epsilon}\right)}.$
(244)

The fourth term in (233) satisfies

$\frac{2}{\|\theta_{\sharp}\|^{2}n}\sum_{i=1}^{n/2}Z_{2i}^{\top}Z_{2i-1}\stackrel{{\scriptstyle
d}}{{=}}\frac{2}{\|\theta_{\sharp}\|^{2}n}\sum_{i=1}^{nd/2}W_{i}\tilde{W}_{i}$
(245)

where $\\{W_{i}\\}_{i\in[nd/2]}$ and $\\{\tilde{W}_{i}\\}_{i\in[nd/2]}$ are
i.i.d. and $W_{i},\tilde{W}_{i}\sim N(0,1)$ are independent. It is well known
[Vershynin, 2018, Lemma 2.7.7] that the product of two subGaussian random
variables (even if they are not independent) is sub-exponential. Here, we have
a simpler and exact characterization. Letting

$W_{i}\tilde{W}_{i}=\left(\frac{W_{i}+\tilde{W}_{i}}{2}\right)^{2}-\left(\frac{W_{i}-\tilde{W}_{i}}{2}\right)^{2}\stackrel{{\scriptstyle
d}}{{=}}V_{i}^{2}-\tilde{V}_{i}^{2}$ (246)

and since $V_{i}=\frac{1}{2}(W_{i}+\tilde{W}_{i})\sim N(0,1)$ and
$\tilde{V}_{i}=\frac{1}{2}(W_{i}-\tilde{W}_{i})\sim N(0,1)$ are uncorrelated,
they are independent. Letting $\chi_{nd/2}^{2}$ and $\tilde{\chi}_{nd/2}^{2}$
be a pair of independent chi-square random variables with $nd/2$ degrees of
freedom, it then holds that

$\displaystyle\frac{2}{\|\theta_{\sharp}\|^{2}n}\sum_{i=1}^{n/2}Z_{2i}^{\top}Z_{2i-1}$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\frac{2}{\|\theta_{\sharp}\|^{2}n}\left(\chi_{nd/2}^{2}-\tilde{\chi}_{nd/2}^{2}\right)$
(247)
$\displaystyle=\frac{2}{\|\theta_{\sharp}\|^{2}n}\left(\chi_{nd/2}^{2}-\mathbb{E}[\chi_{nd/2}^{2}]+\mathbb{E}[\tilde{\chi}_{nd/2}^{2}]-\tilde{\chi}_{nd/2}^{2}\right).$
(248)

From the chi-square tail bound in (330) and (331) it holds that

$\mathbb{P}\left[\left|\chi_{nd/2}^{2}-\mathbb{E}[\chi_{nd/2}^{2}]\right|\geq
2\sqrt{nd/2t}+2t\right]\leq 2e^{-t}$ (249)

and so it holds with probability $1-\epsilon$ that

$\left|\chi_{nd/2}^{2}-\mathbb{E}[\chi_{nd/2}^{2}]\right|\leq\sqrt{2nd\log\left(\frac{2}{\epsilon}\right)}+2\log\left(\frac{2}{\epsilon}\right)\leq
4\sqrt{nd}\log\left(\frac{2}{\epsilon}\right).$ (250)

Hence, by the union bond

$\frac{2}{\|\theta_{\sharp}\|^{2}n}\sum_{i=1}^{n/2}Z_{i}^{\top}Z_{i-1}\leq\frac{16\log\left(\frac{2}{\epsilon}\right)}{\|\theta_{\sharp}\|^{2}}\sqrt{\frac{d}{n}}.$
(251)

with probability $1-\epsilon$. The claim (20) follows from the analysis of the
terms above, the choice of $\epsilon=\frac{2}{n}$ and a union bound. ∎

###  C.2 Proof of Proposition 5: Impossibility lower bound

#### The effect of knowledge of $\theta_{*}$

To begin, it is apparent that if $\|\theta_{*}\|^{2}=\|\theta_{\sharp}\|^{2}$,
the first term in (20) vanishes, and the loss is bounded by the remaining
terms. Furthermore, if $\|\theta_{*}\|\leq 1$ then the dominant term in the
brackets is $\frac{1}{\|\theta_{*}\|^{2}}\sqrt{\frac{d}{n}}$, and evidently,
dominant error term suffers from a penalty of $\sqrt{d}$, even though,
essentially, $\rho$ is a one-dimensional parameter. If $\theta_{*}$ is known
exactly up to a sign, that is $\theta_{\sharp}=\pm\theta_{*}$ (and it is not
just their norms which are equal), then the estimation error can be reduce to
the case of $d=1$. This is a simple consequence of the rotational invariance
of the distribution of the Gaussian noise, which implies that the projections
$(\pm\theta_{*}^{\top}X_{i})_{i=1}^{n}$ are sufficient statistics for the
estimation of $\rho$. The estimator constructs the projections

$U_{i}:=\frac{\theta_{\sharp}^{\top}X_{i}}{\|\theta_{\sharp}\|}=\pm\|\theta_{*}\|\cdot
S_{i}+\frac{\pm\theta_{*}^{\top}Z_{i}}{\|\theta_{*}\|}\stackrel{{\scriptstyle
d}}{{=}}\|\theta_{*}\|\cdot S_{i}+W_{i},$ (252)

where $W_{i}\sim N(0,1)$. This is effectively a one-dimensional model with
parameter given by $\|\theta_{*}\|$, and so Corollary 4 immediately follows
from Theorem 3.

#### Tightness of the impossibility lower bound

Evidently, Corollary 4 and Proposition 5 match in their dependence on the
number of samples $\Theta(\frac{1}{\sqrt{n}})$, but it is not clear what the
optimal dependence on $\|\theta_{*}\|$ is. The estimator we propose is based
on the moment $\rho=\mathbb{E}[U_{i}U_{i-1}]$ and can be contrasted with
likelihood based methods as follows. Letting $p_{\delta}(s_{0}^{n})$ denote
the probability of the sign sequence $s_{0}^{n}$ with flip probability
$\delta$, and letting $\varphi(x)$ denote the standard Gaussian density, the
likelihood of $u_{1}^{n}$ is given by

$P_{\delta,\theta_{*}}(u_{1}^{n})=\sum_{s_{0}^{n}\in\\{\pm
1\\}^{n+1}}p_{\delta}(s_{0}^{n})\cdot\prod_{i\in[n]}\varphi(u_{i}-s_{i}\|\theta_{*}\|).$
(253)

This function is a large degree polynomial in $\delta$ on the order of $n$.
Even if one sums only over $s_{0}^{n}$ with the typical number of flips, then
this degree is $\Theta(\delta n)$, which means this polynomial has a degree
which blows up with $n$. Thus, the MLE may indeed be sensitive to empirical
errors. The update of the Baum-Welch algorithm (or EM) can also be easily
computed and contrasted with our proposed estimator. Letting
$P_{\delta,\theta_{*}}$ denote the probability distribution of the
corresponding model with flip probability $\delta$ and mean parameter
$\theta_{*}$, the Baum-Welch estimate $\hat{\delta}^{(j)}$ at iteration $j$ is
given by

$\hat{\delta}^{(j)}=\frac{1}{n}\sum_{i=1}^{n}P_{\hat{\delta}^{(j-1)},\theta_{*}}(S_{i}\neq
S_{i-1}\mid u_{1}^{n}).$ (254)

Evidently, the inner estimate is the probability that $S_{i}\neq S_{i-1}$
conditioned on the entire sample $u_{1}^{n}$, rather than just $u_{i},u_{i-1}$
as in our proposed estimator in (19). However, deriving sharp error rate
bounds for this estimator seems to be a challenging task.

In terms of minimax lower bounds, the difficulty arises since a reduction to
an (almost) memoryless GMM as in lower bound for the estimation of
$\theta_{*}$ in Theorem 2 does not seem fruitful. A standard application of
Le-Cam’s method requires bounding the total variation between models
$P_{\delta_{1},\theta_{*}}$and $P_{\delta_{2},\theta_{*}}$ for some “close”
$\delta_{1},\delta_{2}$ (e.g., $\delta_{1}=0$ and $\delta_{2}=\epsilon$). The
total variation is then typically bounded by the KL divergence. However, the
KL divergence does not tensorize (due to the memory), and using the chain rule
requires evaluating the KL divergence for the process $U_{1}^{n}$ which is not
a Markov process (but rather an HMM). Alternatively, further bounding the
$n$-dimensional KL divergence with a chi-square divergence, which is a
convenient choice for mixture models, leads to an excessively large bound.

#### The proof’s main ideas

The proof is based on Le-Cam’s method that is applied to a genie-aided model
in which the estimator knows every other sign
$S_{0},S_{2},S_{4},\ldots,S_{n}$. The main technical challenge is then to
bound the total variation for this genie-aided model. A complete proof is
presented below.

###### Proof of Proposition 5..

For a given $\delta\in[0,1]$ let $P_{\delta}^{(n)}$ denote the probability
distribution of $U_{1}^{n}$ under the model (252), that is $U_{i}=t\cdot
S_{i}+W_{i}$, $i\in[n]$, where $t\equiv\|\theta_{*}\|$ and $S_{i}$ is a binary
symmetric Markov chain with flip probability $\delta$, and
$\mathbb{P}[S_{0}=1]=\frac{1}{2}$. It should be noted that $P_{\delta}^{(n)}$
is a distribution on $\mathbb{R}^{n}$ which is _not_ a product distribution.
To lower bound the estimation error of an estimator $\hat{\delta}(U_{1}^{n})$,
we consider a genie-aided estimator which is informed with the values of
$S_{0},S_{2},S_{4},\ldots S_{n}$ (assuming for simplicity that $n$ is even).
We then set $\epsilon\in(0,\frac{1}{2})$ and use Le-Cam’s two point method
[Wainwright, 2019, Section 15.2] with $\delta_{0}=\frac{1}{2}$ and
$\delta_{1}=\frac{1}{2}-\epsilon$. Let us denote by $Q_{U_{1}^{n}S_{0}^{n}}$
(resp. $P_{U_{1}^{n}S_{0}^{n}}$) the joint probability distribution of
$U_{1}^{n}$ and $S_{0}^{n}$ under $\delta=\frac{1}{2}$ (resp.
$\delta=\frac{1}{2}-\epsilon$), and marginals and conditional versions by
standard notation, e.g., $Q_{U_{2}U_{4}\mid S_{0}S_{2}}$. The proof of the
proposition follows from Le-Cam’s two point method, which states that for any
estimator $\hat{\delta}(U_{1}^{n},S_{0},S_{2},S_{4},\ldots,S_{n})$, and thus
also for any less informed estimator $\hat{\delta}(U_{1}^{n})$,

$\mathbb{E}\left[\left|\hat{\delta}-\delta\right|\right]\geq\frac{\epsilon}{2}\cdot\left(1-\mathrm{d_{{TV}}}\left(P_{U_{1}U_{2}\cdots
U_{n}S_{0}S_{2}S_{4}\cdots S_{n}},Q_{U_{1}U_{2}\cdots
U_{n}S_{0}S_{2}S_{4}\cdots S_{n}}\right)\right).$ (255)

We next obtain a bound on the total variation distance in Lemma 16 of
$4\sqrt{n}\epsilon$ for $t\leq\frac{1}{\sqrt{2}}$ and choosing
$\epsilon=\frac{1}{8\sqrt{n}}$ in (255) completes the proof of the lower
bound. ∎

######  Lemma 16.

If $t\leq 1/\sqrt{2}$ then

$\mathrm{d_{{TV}}}\left(P_{U_{1}U_{2}\cdots U_{n},S_{0}S_{2}S_{4}\cdots
S_{n}},Q_{U_{1}U_{2}\cdots U_{n},S_{0}S_{2}S_{4}\cdots
S_{n}}\right)\leq\sqrt{\frac{5}{2}n}\epsilon+\sqrt{8n}t\epsilon.$ (256)

###### Proof.

By Pinsker’s inequality (e.g., [Tsybakov, 2008, Lemma 2.5]), it holds for any
pair of probability measures $P,Q$ that
$\mathrm{d_{{TV}}}(P,Q)\leq\sqrt{\frac{1}{2}\mathrm{D_{KL}}(Q,P)}.$ We next
upper bound the KL divergence. To this end, recall the chain rule, that for
any joint distributions $P_{XY}$ and $Q_{XY}$ (with conditional distributions
$P_{Y\mid X},Q_{Y\mid X})$ the chain rule for the KL divergence states that
$\mathrm{D_{KL}}(P_{XY}\mid\mid Q_{XY})=\mathrm{D_{KL}}(P_{X}\mid\mid
Q_{X})+\mathrm{D_{KL}}(P_{Y\mid X}\mid\mid Q_{Y\mid X}\mid P_{X})$ where
$\mathrm{D_{KL}}(P_{Y\mid X}\mid\mid Q_{Y\mid X}\mid
P_{X}):=\int\mathrm{D_{KL}}(P_{Y\mid X}\mid\mid Q_{Y\mid X})\mathrm{d}P_{X}$
is the conditional KL divergence [Cover and Thomas, 2006, Theorem 2.5.3].
Thus,

$\displaystyle\mathrm{D_{KL}}\left(P_{U_{1}U_{2}\cdots
U_{n},S_{0}S_{2}S_{4}\cdots S_{n}}\mid\mid Q_{U_{1}U_{2}\cdots
U_{n},S_{0}S_{2}S_{4}\cdots S_{n}}\right)$
$\displaystyle=\mathrm{D_{KL}}\left(P_{S_{0}S_{2}S_{4}\cdots S_{n}}\mid\mid
Q_{S_{0}S_{2}S_{4}\cdots S_{n}}\right)$ (257)
$\displaystyle\hphantom{==}+\mathrm{D_{KL}}\left(P_{U_{1}U_{2}\cdots U_{n}\mid
S_{0}S_{2}S_{4}\cdots S_{n}}\mid\mid Q_{U_{1}U_{2}\cdots U_{n}\mid
S_{0}S_{2}S_{4}\cdots S_{n}}\mid P_{S_{0}S_{2}S_{4}\cdots S_{n}}\right).$
(258)

We next bound each of the two KL divergences appearing in (257) and (258).
First,

$\displaystyle\mathrm{D_{KL}}\left(P_{S_{0}S_{2}S_{4}\cdots S_{n}}\mid\mid
Q_{S_{0}S_{2}S_{4}\cdots S_{n}}\right)$
$\displaystyle\overset{{\scriptstyle(a)}}{\leq}\mathrm{D_{KL}}\left(P_{S_{0}^{n}}\mid\mid
Q_{S_{0}^{n}}\right)$ (259)
$\displaystyle=\sum_{s_{0}^{n}\in\\{-1,1\\}^{n+1}}P_{S_{0}^{n}}(s_{0}^{n})\log\frac{P_{S_{0}^{n}}(s_{0}^{n})}{Q_{S_{0}^{n}}(s_{0}^{n})}$
(260) $\displaystyle\overset{{\scriptstyle(b)}}{=}(n+1)\log
2-\sum_{s_{0}^{n}\in\\{-1,1\\}^{n+1}}P_{S_{0}^{n}}(s_{0}^{n})\log\frac{1}{P_{S_{0}^{n}}(s_{0}^{n})}$
(261) $\displaystyle\overset{{\scriptstyle(c)}}{=}(n+1)\log
2-H(S_{0},S_{1},\ldots,S_{n})$ (262)
$\displaystyle\overset{{\scriptstyle(d)}}{=}(n+1)\log
2-H(S_{0})-\sum_{i=1}^{n}H(S_{i}\mid S_{0}^{i-1})$ (263)
$\displaystyle\overset{{\scriptstyle(e)}}{=}n\cdot\log
2-\sum_{i=1}^{n}H(S_{i}\mid S_{i-1})$ (264)
$\displaystyle\overset{{\scriptstyle(f)}}{=}n\cdot\left[\log
2-h_{b}\left(\frac{1}{2}-\epsilon\right)\right],$ (265)

where $(a)$ follows from the convexity of the KL divergence (recall
$S_{0}^{n}=S_{0},S_{1},S_{2},\ldots,S_{n}$), $(b)$ since
$Q_{S_{0}^{n}}(s_{0}^{n})=2^{-(n+1)}$ for any $s_{0}^{n}\in\\{\pm 1\\}^{n+1}$,
$(c)$ follows by defining the entropy $H$ of $S_{0}^{n}$ (under the
probability measure $P$), $(d)$ follows from the chain rule of entropy and the
definition of conditional entropy [Cover and Thomas, 2006, Theorem 2.2.1],
$(e)$ follows from Markovity and $H(S_{0})=\log 2$, $(f)$ follows from
$H(S_{i}\mid S_{i-1})=h_{b}(\delta_{1})=h_{b}(\frac{1}{2}-\epsilon)$ where
$h_{b}(\delta):=-\delta\log\delta-(1-\delta)\log(1-\delta)$ is the binary
entropy function. Now, for $\epsilon\in(0,\frac{1}{2})$ the power series
expansion of the binary entropy function results the bound

$\displaystyle h_{b}\left(\frac{1}{2}-\epsilon\right)$ $\displaystyle=\log
2-\sum_{k=1}^{\infty}\frac{(2\epsilon)^{2k}}{2k(2k-1)}$ (266)
$\displaystyle\geq\log 2-(2\epsilon)^{2}\sum_{k=1}^{\infty}\frac{1}{2k(2k-1)}$
(267) $\displaystyle\geq\log
2-4\epsilon^{2}\sum_{k=1}^{\infty}\frac{1}{(2k-1)^{2}}$ (268)
$\displaystyle=\log 2-\frac{\pi^{2}}{2}\epsilon^{2}$ (269)
$\displaystyle\geq\log 2-5\epsilon^{2}.$ (270)

Inserting this bound into (265) results that the first term in (257) is upper
bounded as

$\mathrm{D_{KL}}\left(P_{S_{0}S_{2}S_{4}\cdots S_{n}}\mid\mid
Q_{S_{0}S_{2}S_{4}\cdots S_{n}}\right)\leq 5n\epsilon^{2}.$ (271)

We now move on to bound the second term in (258). To this end, we note that
the distribution of $U_{1}^{n}$ conditioned on $S_{0},S_{2},S_{4},\ldots
S_{n}$ can be decomposed in a simple way. First, under $Q$, the signs
$S_{0}^{n}$ are i.i.d., and so $U_{1}^{n}$ is a vector of independent samples
from a Gaussian mixture model. Furthermore, $U_{i}$ depends on the sign
$S_{i}$ but otherwise is independent of all other
$S_{1}^{n}\backslash\\{S_{i}\\}$. Hence,

$Q_{U_{1}U_{2}\cdots U_{n}\mid S_{0}S_{2}S_{4}\cdots S_{n}}=Q_{U_{1}}\cdot
Q_{U_{2}\mid S_{2}}\cdot Q_{U_{3}}\cdot Q_{U_{4}\mid S_{4}}\cdots
Q_{U_{n-1}}\cdot Q_{U_{n}\mid S_{n}},$ (272)

that is, a model of independent samples, where the odd samples are drawn from
a Gaussian mixture model, and the even samples from a Gaussian location model
with known sign $S_{2i}$. Second, under $P$, an application of Bayes rule and
the Markovity assumption results the decomposition to pairs of samples given
by

$\displaystyle P_{U_{1}U_{2}\cdots U_{n}\mid S_{0}S_{2}S_{4}\cdots S_{n}}$
$\displaystyle=P_{U_{1}U_{2}\mid S_{0}S_{2}S_{4}\cdots S_{n}}\cdot
P_{U_{3}U_{4}\mid S_{0}S_{2}S_{4}\cdots S_{n}U_{1}U_{2}}\cdots
P_{U_{n-1}U_{n}\mid S_{0}S_{2}S_{4}\cdots S_{n}U_{1}^{n-2}}$ (273)
$\displaystyle=P_{U_{1}U_{2}\mid S_{0}S_{2}}\cdot P_{U_{3}U_{4}\mid
S_{2}S_{4}}\cdots P_{U_{n-1}U_{n}\mid S_{n-2}S_{n}}.$ (274)

The chain rule for the KL divergence therefore implies that

$\displaystyle\mathrm{D_{KL}}\left(P_{U_{1}U_{2}\cdots U_{n}\mid
S_{0}S_{2}S_{4}\cdots S_{n}}\mid\mid Q_{U_{1}U_{2}\cdots U_{n}\mid
S_{0}S_{2}S_{4}\cdots S_{n}}\mid P_{S_{0}S_{2}S_{4}\cdots S_{n}}\right)$
$\displaystyle\overset{{\scriptstyle(a)}}{=}\mathrm{D_{KL}}\left(P_{U_{1}U_{2}\mid
S_{0}S_{2}}\mid\mid Q_{U_{1}}\cdot Q_{U_{2}\mid S_{2}}\mid
P_{S_{0}S_{2}}\right)+\mathrm{D_{KL}}\left(P_{U_{3}U_{4}\mid
S_{2}S_{4}}\mid\mid Q_{U_{3}}\cdot Q_{U_{4}\mid S_{4}}\mid
P_{S_{2}S_{4}}\right)$
$\displaystyle\hphantom{==}+\cdots+\mathrm{D_{KL}}\left(P_{U_{n-1}U_{n}\mid
S_{n-2}S_{n}}\mid\mid Q_{U_{n-1}}\cdot Q_{U_{n}\mid S_{n}}\mid
P_{S_{n-2}S_{n}}\right)$ (275)
$\displaystyle\overset{{\scriptstyle(b)}}{=}\frac{n}{2}\cdot\mathrm{D_{KL}}\left(P_{U_{1}U_{2}\mid
S_{0}S_{2}}\mid\mid Q_{U_{1}}\cdot Q_{U_{2}\mid S_{2}}\mid
P_{S_{0}S_{2}}\right)$ (276)
$\displaystyle\overset{{\scriptstyle(c)}}{=}\frac{n}{2}\cdot\mathrm{D_{KL}}\left(P_{U_{1}\mid
S_{0}S_{2}}\mid\mid Q_{U_{1}}\mid
P_{S_{0}S_{2}}\right)+\frac{n}{2}\cdot\mathrm{D_{KL}}\left(P_{U_{2}\mid
S_{0}S_{2}}\mid\mid Q_{U_{2}\mid S_{2}}\mid P_{S_{0}S_{2}}\right)$ (277)
$\displaystyle\overset{{\scriptstyle(d)}}{=}\frac{n}{2}\cdot\mathrm{D_{KL}}\left(P_{U_{1}\mid
S_{0}S_{2}}\mid\mid Q_{U_{1}}\mid P_{S_{0}S_{2}}\right)$ (278)
$\displaystyle\overset{{\scriptstyle(e)}}{=}\frac{n}{2}\cdot\mathrm{D_{KL}}\left(P_{U_{1}\mid
S_{0}=1,S_{2}}\mid\mid Q_{U_{1}}\mid P_{S_{2}\mid S_{0}=1}\right)$ (279)

where $(a)$ follows from (274) and the chain rule for KL divergence, and $(b)$
follows from the stationarity of the Markov chain, $(c)$ follows again from
the chain rule, $(d)$ follows since $\mathrm{D_{KL}}(P_{U_{2}\mid
S_{0}S_{2}}\mid\mid Q_{U_{2}\mid S_{2}}\mid
P_{S_{0}S_{2}})=\mathrm{D_{KL}}(P_{U_{2}\mid S_{2}}\mid\mid Q_{U_{2}\mid
S_{2}}\mid P_{S_{2}})=0$, and $(e)$ follows since by symmetry, we may
condition on $S_{0}=1$. Thus, the last KL divergence in (279) should be
averaged over the cases $S_{2}=-1$ and $S_{2}=1$. For the first case, it holds
that $P_{S_{1}\mid S_{0}S_{2}}(\cdot\mid 1,-1)$ is a uniform distribution on
$\\{\pm 1\\}$. Hence, conditioned on $S_{0}=1,S_{2}=-1$, under $P$, $U_{1}$ is
a sample from a balanced Gaussian mixture, just as under $Q$. So,

$\mathrm{D_{KL}}\left(P_{U_{1}\mid S_{0}=1,S_{2}=-1}\mid\mid
Q_{U_{1}}\right)=0.$ (280)

So, the KL divergence is only comprised of the term in the second case
$S_{2}=1$. Continuing to evaluate the KL divergence in (279), we get

$\displaystyle\mathrm{D_{KL}}\left(P_{U_{1}\mid S_{0}=1,S_{2}}\mid\mid
Q_{U_{1}}\mid P_{S_{2}\mid S_{0}=1}\right)$ $\displaystyle=P_{S_{2}\mid
S_{0}=1}(1)\cdot\mathrm{D_{KL}}\left(P_{U_{1}\mid S_{0}=1,S_{2}=1}\mid\mid
Q_{U_{1}}\right)$ (281)
$\displaystyle\overset{{\scriptstyle(a)}}{\leq}\mathrm{D_{KL}}\left(P_{U_{1}\mid
S_{0}=1,S_{2}=1}\mid\mid Q_{U_{1}}\right)$ (282)
$\displaystyle\overset{{\scriptstyle(b)}}{\leq}\chi^{2}\left(P_{U_{1}\mid
S_{0}=1,S_{2}=1}\mid\mid Q_{U_{1}}\right)$ (283)

where $(a)$ follows since $P_{S_{2}\mid S_{0}}(1\mid
1)=(1-\delta_{1})^{2}+\delta_{1}^{2}\leq 1$, and $(b)$ follows from the chi-
square divergence bound on the KL divergence (e.g., [Tsybakov, 2008, Eq.
(2.27)]). Now, recall that under $Q_{U_{1}}$, it holds that
$U_{1}\sim\frac{1}{2}N(t,1)+\frac{1}{2}N(-t,1)$, and under $P_{U_{1}\mid
S_{0}=1,S_{2}=1}$ it holds that $U_{1}\sim(1-\alpha)\cdot N(t,1)+\alpha
N(-t,1)$ where

$\alpha=P_{S_{1}\mid S_{0}S_{2}}(-1\mid
1,1)=\frac{\delta_{1}^{2}}{(1-\delta_{1})^{2}+\delta_{1}^{2}}.$ (284)

Thus, the chi-square divergence in (283) is between a balanced Gaussian
mixture (with probability that $S_{1}=1$ being $\frac{1}{2}$) and an
unbalanced Gaussian mixture (with probability that $S_{1}=1$ being
$1-\alpha>\frac{1}{2}$). Let $\psi_{\beta}(u)$ denote the probability density
function of $U=t\cdot S+W$ with $\mathbb{P}[S=-1]=1-\mathbb{P}[S=1]=\beta$ and
$W\sim N(0,1)$. Then, if $\varphi(u)\equiv\frac{1}{\sqrt{2\pi}}e^{-u^{2}/2}$
is the standard Gaussian density function, it holds that

$\displaystyle\psi_{\beta}(u)$
$\displaystyle=\beta\cdot\varphi(u+t)+(1-\beta)\cdot\varphi(u+t)$ (285)
$\displaystyle=\beta\cdot\frac{1}{\sqrt{2\pi}}e^{-(u+t)^{2}/2}+(1-\beta)\cdot\frac{1}{\sqrt{2\pi}}e^{-(u-t)^{2}/2}$
(286) $\displaystyle=e^{-t^{2}/2}\cdot\varphi(u)\cdot\left[\beta\cdot
e^{-ut}+(1-\beta)\cdot e^{ut}\right].$ (287)

Specifically, for $\beta=\frac{1}{2}$ it holds that
$\psi_{1/2}(u)=e^{-t^{2}/2}\cdot\varphi(u)\cdot\cosh(ut)$. Therefore, the chi-
square divergence from (283) is upper bounded as

$\displaystyle\chi^{2}\left(P_{U_{1}\mid S_{0}=1,S_{2}=1}\mid\mid
Q_{U_{1}}\right)$
$\displaystyle=\int_{-\infty}^{\infty}\frac{\left[P_{U_{1}\mid
S_{0}S_{2}}(u\mid 1,1)-Q_{U_{1}}(u)\right]^{2}}{Q_{U_{1}}(u)}\mathrm{d}u$
(288)
$\displaystyle=e^{-t^{2}/2}\cdot\int_{-\infty}^{\infty}\varphi(u)\frac{\left[(\alpha-\frac{1}{2})\cdot
e^{-ut}+(1-\alpha-\frac{1}{2})\cdot e^{ut}\right]^{2}}{\cosh(ut)}\mathrm{d}u$
(289)
$\displaystyle\overset{{\scriptstyle(a)}}{\leq}e^{-t^{2}/2}\cdot\int_{-\infty}^{\infty}\varphi(u)\left[\left(\alpha-\frac{1}{2}\right)\cdot
e^{-ut}+\left(1-\alpha-\frac{1}{2}\right)\cdot e^{ut}\right]^{2}\mathrm{d}u$
(290)
$\displaystyle=e^{-t^{2}/2}\cdot\int_{-\infty}^{\infty}\varphi(u)\left[\left(\alpha-\frac{1}{2}\right)^{2}\cdot
e^{-2ut}+2\left(\alpha-\frac{1}{2}\right)\left(1-\alpha-\frac{1}{2}\right)+\left(1-\alpha-\frac{1}{2}\right)^{2}\cdot
e^{2ut}\right]\mathrm{d}u$ (291)
$\displaystyle=e^{-t^{2}/2}\cdot\left[\left(\alpha-\frac{1}{2}\right)^{2}e^{2t^{2}}+2\left(\alpha-\frac{1}{2}\right)\left(1-\alpha-\frac{1}{2}\right)+\left(1-\alpha-\frac{1}{2}\right)^{2}e^{2t^{2}}\right]$
(292)
$\displaystyle\overset{{\scriptstyle(b)}}{\leq}\left[\left(\alpha-\frac{1}{2}\right)^{2}(1+4t^{2})+2\left(\alpha-\frac{1}{2}\right)\left(1-\alpha-\frac{1}{2}\right)+\left(1-\alpha-\frac{1}{2}\right)^{2}(1+4t^{2})\right]$
(293)
$\displaystyle=4t^{2}\cdot\left[\left(\alpha-\frac{1}{2}\right)^{2}+\left(1-\alpha-\frac{1}{2}\right)^{2}\right]$
(294) $\displaystyle=2t^{2}\cdot\left[1-4\alpha(1-\alpha)\right]$ (295)
$\displaystyle\overset{{\scriptstyle(c)}}{=}2t^{2}\cdot\left[1-\left(\frac{2-8\epsilon^{2}}{2+8\epsilon^{2}}\right)^{2}\right]$
(296)
$\displaystyle=2t^{2}\cdot\left[1-\frac{2-8\epsilon^{2}}{2+8\epsilon^{2}}\right]\left[1+\frac{2-8\epsilon^{2}}{2+8\epsilon^{2}}\right]$
(297)
$\displaystyle=2t^{2}\cdot\left[\frac{64\epsilon^{2}}{(2+8\epsilon^{2})^{2}}\right]$
(298) $\displaystyle\leq 32t^{2}\epsilon^{2},$ (299)

where $(a)$ follows since $\cosh(x)\geq 1$ for all $x\in\mathbb{R}$, $(b)$
follows using the assumption $t\leq\frac{1}{\sqrt{2}}$, and so $e^{2t^{2}}\leq
1+4t^{2}$, and $(c)$ is by setting

$4\alpha(1-\alpha)=4\frac{\delta_{1}^{2}(1-\delta_{1})^{2}}{\left[(1-\delta_{1})^{2}+\delta_{1}^{2}\right]^{2}}=4\frac{\left(\frac{1}{2}-\epsilon\right)^{2}\left(\frac{1}{2}+\epsilon\right)^{2}}{\left[\left(\frac{1}{2}+\epsilon\right)^{2}+\left(\frac{1}{2}-\epsilon\right)^{2}\right]^{2}}=\left(\frac{2-8\epsilon^{2}}{2+8\epsilon^{2}}\right)^{2}.$
(300)

Using the bound (299) in (283), We thus conclude that

$\mathrm{D_{KL}}\left(P_{U_{1}\mid S_{0}=1,S_{2}}\mid\mid Q_{U_{1}}\mid
P_{S_{2}\mid S_{0}=1}\right)\leq 32t^{2}\epsilon^{2}$ (301)

and then in (279) that the second KL term of (258) is upper bounded as

$\mathrm{D_{KL}}\left(P_{U_{1}U_{2}\cdots U_{n}\mid S_{0}S_{2}S_{4}\cdots
S_{n}}\mid\mid Q_{U_{1}U_{2}\cdots U_{n}\mid S_{0}S_{2}S_{4}\cdots S_{n}}\mid
P_{S_{0}S_{2}S_{4}\cdots S_{n}}\right)\leq 16nt^{2}\epsilon^{2}.$ (302)

Combining this bound with the bound on the first KL term of (258) given in
(271) we obtain

$\mathrm{D_{KL}}\left(P_{U_{1}U_{2}\cdots U_{n},S_{0}S_{2}S_{4}\cdots
S_{n}}\mid\mid Q_{U_{1}U_{2}\cdots U_{n},S_{0}S_{2}S_{4}\cdots
S_{n}}\right)\leq 5n\epsilon^{2}+16nt^{2}\epsilon^{2}.$ (303)

The aforementioned Pinsker bound on the total variation distance, and the
relation $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ then completes the proof. ∎

##  Appendix D Proofs for Section 3: Analysis of Algorithm 1

###### Proof of Theorem 6..

Before presenting the analysis of Algorithm 1, we first specify the choices of
the constants $\lambda_{\theta}\geq 1$ and $\lambda_{\delta}\geq 1$. Recall
that from Theorem 1 for a Markov model $(\theta_{*},\delta)$ in low dimension,
$d\leq\delta n$, there exists a numerical constant $\lambda_{\theta}>0$,333The
constants can be deduced from the proof, though they were not optimized. We
reiterate that both $\lambda_{\theta}$ and $\lambda_{\delta}$ below are
universal numerical constants. In particular, they do _not_ depend on
$\theta_{*},\delta$ or any other problem parameters. The subscripts are merely
to emphasize that they are obtained in the estimation procedure for
$\theta_{*}$ and $\delta$, respectively.  so that the estimator
$\hat{\theta}\equiv\hat{\theta}_{\text{cov}}(X_{1}^{n};k=\frac{1}{8\delta})$,
which assumes a perfect knowledge of $\delta$, achieves with probability
$1-O(\frac{1}{n})$,

$\mathsf{loss}(\hat{\theta},\theta_{*})\leq\lambda_{\theta}\cdot\log(n)\cdot\begin{cases}\|\theta_{*}\|,&\|\theta_{*}\|\leq\left(\frac{\delta
d}{n}\right)^{1/4}\\\ \frac{1}{\|\theta_{*}\|}\sqrt{\frac{\delta
d}{n}},&\left(\frac{\delta
d}{n}\right)^{1/4}\leq\|\theta_{*}\|\leq\sqrt{\delta}\\\
\sqrt{\frac{d}{n}},&\|\theta_{*}\|\geq\sqrt{\delta}\end{cases}$ (304)

if $d\leq\delta n$, and

$\mathsf{loss}(\hat{\theta},\theta_{*})\leq\lambda_{\theta}\cdot\log(n)\cdot\begin{cases}\|\theta_{*}\|,&\|\theta_{*}\|\leq\sqrt{\frac{d}{n}}\\\
\sqrt{\frac{d}{n}},&\|\theta_{*}\|\geq\sqrt{\frac{d}{n}}\end{cases}$ (305)

if $\delta n\leq d\leq n$. We assume that $\lambda_{\theta}\geq 1$, and
otherwise replace $\lambda_{\theta}$ by $1$.

In addition, it holds from Theorem 3 that for the estimator $\hat{\delta}$ for
$\delta$, which is based on a mismatched mean $\theta_{\sharp}$, i.e.,

$\hat{\delta}\equiv\hat{\delta}_{\text{corr}}(X_{1}^{n};\theta_{\sharp}):=\frac{1}{2}\left(1-\hat{\rho}_{\text{corr}}(X_{1}^{n};\theta_{\sharp})\right),$
(306)

there exists a numerical constant $\lambda_{\delta}>0$, for which with
probability larger than $1-O(\frac{1}{n})$,

$\mathsf{loss}(\hat{\delta},\delta)\leq\lambda_{\delta}\left[\frac{\left|\|\theta_{*}\|^{2}-\|\theta_{\sharp}\|^{2}\right|}{\|\theta_{\sharp}\|^{2}}+\frac{\log(n)}{\|\theta_{\sharp}\|^{2}}\sqrt{\frac{d}{n}}\right]$
(307)

assuming that $\|\theta_{\sharp}\|\leq 2$ in order to simplify the bound to
the regime of interest. We assume here too that $\lambda_{\delta}\geq 1$, and
otherwise replace $\lambda_{\delta}$ by $1$.

For the analysis we assume that all three steps of the algorithm are
successful estimation events. By the union bound, this occurs with probability
$1-O(\frac{1}{n})$.

Analysis of Step A:

First, suppose that

$\|\theta_{*}\|\leq\lambda_{\theta}\log(n)\cdot\left(\frac{d}{n}\right)^{1/4}.$
(308)

Then, it holds by (304) and (305) that

$\|\hat{\theta}^{(A)}\|\leq\|\theta_{*}\|+\lambda_{\theta}\log(n)\cdot\left(\frac{\delta
d}{n}\right)^{1/4}\leq 2\lambda_{\theta}\log(n)\left(\frac{d}{n}\right)^{1/4}$
(309)

and thus the algorithm will stop and output $\hat{\theta}=0$, for which it
holds that $\mathsf{loss}(\theta_{*},\hat{\theta}^{(A)})=\|\theta_{*}\|$. This
agrees with (26) and (27).

Second, suppose that $\|\theta_{*}\|\geq 1$. Now, it holds by (304) and (305)

$\|\hat{\theta}^{(A)}\|\geq\|\theta_{*}\|-\lambda_{\theta}\log(n)\sqrt{\frac{d}{n}}\geq\frac{1}{2},$
(310)

by the assumption $d\leq\frac{1}{4\lambda_{\theta}^{2}\log^{2}(n)}\cdot n$.
Thus the algorithm will stop and output $\hat{\theta}=\hat{\theta}^{(A)}$, for
which it also holds that

$\mathsf{loss}(\theta_{*},\hat{\theta}^{(A)})\leq\lambda_{\theta}\cdot\log(n)\cdot\sqrt{\frac{d}{n}}$
(311)

(irrespective of the value of $\delta$). This also agrees with (26) and (27).

Thus, the algorithm will continue to Step B only if

$\lambda_{\theta}\log(n)\cdot\left(\frac{d}{n}\right)^{1/4}\leq\|\theta_{*}\|\leq
1$ (312)

which we henceforth assume. As a preparation for Step B, in which
$\hat{\theta}^{(A)}$ will play the rule of $\theta_{\sharp}$ (the mismatched
mean estimator of Section 3) in an estimation of $\delta$, we bound the
absolute value $\left|\|\hat{\theta}^{(A)}\|-\|\theta_{*}\|\right|$ and also
show that $\|\hat{\theta}^{(A)}\|\geq\frac{1}{2}\|\theta_{*}\|$. To this end,
we may repeat the arguments of the proof of Theorem 1, and specifically,
arguments similar to (126) to show that

$\displaystyle\left|\|\hat{\theta}^{(A)}\|-\|\theta_{*}\|\right|$
$\displaystyle\lesssim\log(n)\cdot\psi\left(n,d,\delta=\frac{1}{2},k=1\right)$
(313)
$\displaystyle=\log(n)\cdot\left[2\sqrt{\frac{\delta}{n}}\cdot\|\theta_{*}\|^{2}+2\sqrt{\frac{d}{n}}\cdot\|\theta_{*}\|+13\sqrt{\frac{d}{n}}+10\frac{d}{n}\right].$
(314)

Note that the probability of this event is $1-O(\frac{1}{n})$, and is included
in the successful estimation event mentioned in the beginning of the proof.
Under the assumption $\|\theta_{*}\|\leq 1$ it then holds that

$\left|\|\hat{\theta}^{(A)}\|-\|\theta_{*}\|\right|\leq\lambda_{\theta}\log(n)\cdot\sqrt{\frac{d}{n}}.$
(315)

Note that we take $\lambda_{\theta}>0$ to be large enough so that this holds.
The assumed case
$\|\theta_{*}\|\geq\lambda_{\theta}\log(n)\cdot\left(\frac{d}{n}\right)^{1/4}$,
and the assumption of the theorem $d\leq$$\frac{n}{16}$ then imply that

$\frac{1}{2}\|\theta_{*}\|\leq\|\hat{\theta}^{(A)}\|\leq 2\|\theta_{*}\|\leq
2.$ (316)

Analysis of Step B:

Setting $\theta_{\sharp}=\hat{\theta}^{(A)}$ and utilizing the two properties
just derived in (315) and (316), (307) implies that

$\displaystyle\mathsf{loss}(\hat{\delta}^{(B)},\delta)$
$\displaystyle\leq\lambda_{\delta}\left[\frac{\left|\|\theta_{*}\|^{2}-\|\theta_{\sharp}\|^{2}\right|}{\|\theta_{\sharp}\|^{2}}+\frac{\log(n)}{\|\theta_{\sharp}\|^{2}}\sqrt{\frac{d}{n}}\right]$
(317) $\displaystyle\leq
4\lambda_{\delta}\lambda_{\theta}\frac{\log(n)}{\|\theta_{*}\|^{2}}\sqrt{\frac{d}{n}}.$
(318)

There are two cases. First suppose that
$\|\theta_{*}\|\geq\sqrt{8\lambda_{\delta}\lambda_{\theta}\log(n)}(\frac{d}{\delta^{2}n})^{1/4}$.
In this case,

$\mathsf{loss}(\hat{\delta}^{(B)},\delta)\leq\frac{\delta}{2}$ (319)

and so

$\frac{1}{2}\delta\leq\hat{\delta}^{(B)}\leq 2\delta.$ (320)

In addition, from (316) it holds that

$64\lambda_{\delta}\lambda_{\theta}\frac{\log(n)}{\|\hat{\theta}^{(A)}\|^{2}}\sqrt{\frac{d}{n}}\geq
16\lambda_{\delta}\lambda_{\theta}\frac{\log(n)}{\|\theta_{*}\|^{2}}\sqrt{\frac{d}{n}}\geq
2\delta\geq\hat{\delta}^{(B)}$ (321)

and so the algorithm will continue to Step C. Now, suppose that
$\|\theta_{*}\|\leq\sqrt{8\lambda_{\delta}\lambda_{\theta}\log(n)}(\frac{d}{\delta^{2}n})^{1/4}$.
If the algorithm stops then $\hat{\theta}=\hat{\theta}^{(A)}$, and the
estimation rates of the Gaussian mixture model are achieved, which agrees with
(26) and (27). If the algorithm does not stop and proceeds to Step C then it
holds that

$\displaystyle\hat{\delta}^{(B)}$ $\displaystyle\geq
64\lambda_{\delta}\lambda_{\theta}\frac{\log(n)}{\|\hat{\theta}^{(A)}\|^{2}}\sqrt{\frac{d}{n}}\geq
16\lambda_{\delta}\lambda_{\theta}\frac{\log(n)}{\|\theta_{*}\|^{2}}\sqrt{\frac{d}{n}}\geq\delta,$
(322)

by (316) and the assumption
$\|\theta_{*}\|\leq\sqrt{8\lambda_{\delta}\lambda_{\theta}\log(n)}(\frac{d}{\delta^{2}n})^{1/4}$.

To conclude, there are cases in which the algorithm proceeds to Step C. If
$\|\theta_{*}\|\geq\sqrt{8\lambda_{\delta}\lambda_{\theta}\log(n)}(\frac{d}{\delta^{2}n})^{1/4}$
then it proceeds Step C and (320) holds. Otherwise, the algorithm might
proceed to Step C, yet now only $\hat{\delta}^{(B)}\geq\frac{1}{2}\delta$ is
assured.

Analysis of Step C:

If the algorithm has proceeded to Step C, then it is guaranteed that
$\hat{\delta}^{(B)}\geq\frac{1}{2}\delta$ (in any event). Recall that,
ideally, had $\delta$ was known, the choice of the block length $k$ for the
estimator $\hat{\theta}_{\text{cov}}(X_{2n+1}^{3n};k)$ is
$k_{*}:=\frac{1}{8\delta}$ (Section 2). It can be readily verified that the
analysis of Section 2 is valid when using any _smaller_ blocklength $k$, and
that the error rate improve as $k$ increases from $k=1$ to $k=k_{*}$. In
accordance to the analysis of Step B, there are two possible cases. If
$\|\theta_{*}\|\geq\sqrt{8\lambda_{\delta}\lambda_{\theta}\log(n)}(\frac{d}{\delta^{2}n})^{1/4}$
then it holds that $\frac{1}{2}\delta\leq\hat{\delta}^{(B)}\leq 2\delta$.
Using $k=\frac{1}{16\hat{\delta}^{(B)}}$ then assures that $k\leq k_{*}$. In
addition, since the error bound scales linearly with $\psi(n,d,\delta,k)$ (see
(92)), the error can increase by a factor at most $2$. Thus,
$\hat{\theta}^{(C)}$ achieves the error rates in (304) and (305), with a
factor of $2$. There are again two cases to consider. If
$d\geq\frac{1}{64\lambda_{\delta}^{2}\lambda_{\theta}^{2}\log^{2}(n)}\delta^{4}n$,
then the interval
$(\sqrt{8\lambda_{\delta}\lambda_{\theta}\log(n)}(\frac{d}{\delta^{2}n})^{1/4},\sqrt{\delta})$
is empty (its right end point is smaller than its left end point), and
plugging $k$ in the error rates of (304) and (305), the resulting error rates
are as in the Gaussian mixture model. If
$d\leq\frac{1}{64\lambda_{\delta}^{2}\lambda_{\theta}^{2}\log^{2}(n)}\delta^{4}n$
then the aforementioned interval is non-empty, and the resulting error rates
agree with (26) and (27). Otherwise, if
$\|\theta_{*}\|\leq\sqrt{8\lambda_{\delta}\lambda_{\theta}\log(n)}(\frac{d}{\delta^{2}n})^{1/4}$,
and the estimator $\hat{\theta}_{\text{cov}}(X_{2n+1}^{3n};k)$ operates with
$k\leq k_{*}$, but $k$ can be as low as $1$. Thus, error rates of the Gaussian
mixture are again achieved, and this agrees with (26) and (27). ∎

##  Appendix E Useful results

#### Bernstein’s inequality

Let $X_{1},\cdots,X_{\ell}$ be independent random variables and $|X_{i}|\leq
b$ almost surely for every $i\in[\ell]$. Then [Wainwright, 2019, Proposition
2.14]444In the original form in [Wainwright, 2019, Proposition 2.14], on the
right hand side of the inequality, $\mathbb{V}[X_{i}]$ is replaced with
$\mathbb{E}[X_{i}^{2}]$ which leads to a seemingly worse concentration bound.
However, applying this form of Bernstein’s inequality to
$Y_{i}=X_{i}-\mathbb{E}[X_{i}]$ allows us to conclude (323).  states

$\displaystyle\mathbb{P}\left[\left|\frac{1}{\ell}\sum_{i=1}^{\ell}(X_{i}-\mathbb{E}[X_{i}])\right|\geq\delta\right]$
$\displaystyle\leq
2\exp\left(-\frac{\ell\delta^{2}/2}{\frac{1}{\ell}\sum_{i=1}^{\ell}\mathbb{V}[X_{i}]+b\delta/3}\right).$
(323)

#### Norm of subGaussian random vectors

A random vector $X\in\mathbb{R}^{d}$ is said to be $\sigma^{2}$-subGaussian if
$\mathbb{E}[X]=0$ and

$\displaystyle\mathbb{E}\left[\exp(tv^{\top}X)\right]$
$\displaystyle\leq\exp\left(\frac{\sigma^{2}t^{2}}{2}\right)$ (324)

for every $t\in\mathbb{R}$ and every $v\in\mathbb{S}^{d-1}$. Let
$X\in\mathbb{R}^{d}$ be a $\sigma^{2}$-subGaussian random vector. Then from
[Rigollet and Hütter, 2019, Theorem 1.19], we have for any $\delta>0$,

$\displaystyle\mathbb{P}\left[\max_{v\in\mathbb{B}^{d}}v^{\top}X\leq
4\sigma\sqrt{d}+2\sigma\sqrt{2\log\left(\frac{1}{\delta}\right)}\right]=\mathbb{P}\left[\|X\|\leq
4\sigma\sqrt{d}+2\sigma\sqrt{2\log\left(\frac{1}{\delta}\right)}\right]\geq
1-\delta.$ (325)

#### Gaussian covariance estimation

Let $W\in\mathbb{R}^{\ell\times d}$ be a matrix with i.i.d. $N(0,1)$ entries.
Then [Wainwright, 2019, Example 6.2] implies that for any $\delta>0$,

$\displaystyle\mathbb{P}\left[\left\|\frac{1}{\ell}W^{\top}W-I_{d}\right\|_{\text{op}}\leq
2\left(\sqrt{\frac{d}{n}}+\delta\right)+\left(\sqrt{\frac{d}{n}}+\delta\right)^{2}\right]$
$\displaystyle\geq 1-2e^{-n\delta^{2}/2}.$ (326)

#### Davis-Kahan’s perturbation bound

Let $\Sigma,\hat{\Sigma}$ be symmetric matrices with the same dimensions. Let
$\lambda_{i}(\Sigma)$ and $v_{i}(\Sigma)$ denote the $i$th largest eigenvalue
and the associated eigenvector (of unit norm) of $\Sigma$. Fix $i$ and assume
that $\lambda_{i}(\Sigma)$ is well-separated from the rest of the spectrum of
$\Sigma$:

$\displaystyle\min_{j\neq i}|\lambda_{i}(\Sigma)-\lambda_{j}(\Sigma)|$
$\displaystyle=\delta>0.$ (327)

Then from [Vershynin, 2018, Theorem 4.5.5], we have

$\displaystyle\mathsf{loss}(v_{i}(\Sigma),v_{i}(\hat{\Sigma}))$
$\displaystyle\leq 4\frac{\|\Sigma-\hat{\Sigma}\|_{\text{op}}}{\delta},$ (328)

where $\mathsf{loss}(\cdot,\cdot)$ was defined in (3).

#### Packing number

Let $M(\delta;\mathbb{B}^{d},\|\cdot\|)$ and
$N(\delta;\mathbb{B}^{d},\|\cdot\|)$ denote the $\delta$-packing number and
$\delta$-covering number of $\mathbb{B}^{d}$ w.r.t. $\|\cdot\|$, respectively.
From [Wainwright, 2019, Lemma 5.7 and Example 5.8], we have

$\displaystyle M(\delta;\mathbb{B}^{d},\|\cdot\|)$ $\displaystyle\geq
N(\delta;\mathbb{B}^{d},\|\cdot\|)\geq\left(\frac{1}{\delta}\right)^{d}.$
(329)

#### Chi-square tail bounds

From the chi-square tail bound in [Boucheron et al., 2013, Remark 2.11]

$\mathbb{P}\left[\chi_{\ell}^{2}-\ell\geq 2\sqrt{\ell t}+2t\right]\leq e^{-t}$
(330)

and

$\mathbb{P}\left[\chi_{\ell}^{2}-\ell\leq-2\sqrt{\ell t}\right]\leq e^{-t}.$
(331)

##  Appendix F Numerical validation

In this section, we provide numerical validation of the performance of the
estimators $\hat{\theta}_{\text{cov}}(X_{1}^{n};k)$ (cf. (14)),
$\hat{\delta}_{\text{corr}}(X_{1}^{n};\theta_{\sharp})=\frac{1}{2}(1-\hat{\rho}_{\text{corr}}(X_{1}^{n};\theta_{\sharp}))$
(cf. (19)) and Algorithm 1 proposed and analyzed in Theorem 1, Theorem 3 and
Corollary 4, and Theorem 6, respectively. Numerical results for these
estimators/algorithm are plotted in Figures 3, 4 and 5, respectively.

[width=0.5]plot-theta-estimation.pdf

Figure 3: Numerical validation of the performance of the
$\theta_{*}$-estimator $\hat{\theta}_{\text{cov}}(X_{1}^{n};k)$ (cf. (14))
used in the proof of Theorem 1 which assumes a known $\delta$. We take the
minimum loss (cf. (3)) achieved by
$\hat{\theta}_{\text{cov}}\left(X_{1}^{n};k=\frac{1}{8\delta}\right)$ and the
trivial estimator $\hat{\theta}_{0}(X_{1}^{n})=0$. We plot the loss as a
function of $t=\|\theta_{*}\|$ for $t\in[0,5]$ with step size $0.05$. We take
$n=5000,d=250,\delta=0.05$. Therefore, $\delta>\frac{1}{n}$ and $d<\delta n$.

[width=0.48]plot-delta-estimation-mismatched.pdf [width=0.48]plot-delta-
estimation-matched.pdf

Figure 4: Numerical validation of the performance of the $\delta$-estimator
$\hat{\delta}_{\text{corr}}(X_{1}^{n};\theta_{\sharp})=\frac{1}{2}(1-\hat{\rho}(X_{1}^{n};\theta_{\sharp}))$
where $\hat{\rho}_{\text{corr}}(X_{1}^{n};\theta_{\sharp})$ is given by (19).
This estimator is used in Theorem 3 (which assumes $\theta_{\sharp}$ is a
mismatched estimate of $\theta_{*}$), Corollary 4 (which assume
$\theta_{\sharp}=\pm\theta_{*}$) and also Theorem 6 (which concerns estimating
$\theta_{*}$ without the knowledge of $\delta$). We take the loss
$|\hat{\delta}-\delta|$ achieved by $\hat{\delta}_{\text{corr}}$ and the
trivial estimators
$\hat{\delta}_{0}(X_{1}^{n})=0,\hat{\delta}_{1}(X_{1}^{n})=1,\hat{\delta}_{1/2}(X_{1}^{n})=\frac{1}{2}$.
We plot the loss as a function of $t=\|\theta_{*}\|$ for $t\in[0,1]$ with step
size $0.05$. We take $n=500,d=250,\delta=0.1$. In the left panel, we assume
the estimator has access to $\theta_{\sharp}$ with
$\|\theta_{\sharp}\|=1.2\cdot\|\theta_{*}\|$. In the right panel, we assume
$\theta_{\sharp}=\theta_{*}$ and therefore the model (2) is equivalent to the
model (252).

[width=0.5]plot-theta-estimation-joint.pdf

Figure 5: Numerical validation of the performance of Algorithm 1 in Theorem 6.
We take the minimum loss (cf. (3)) achieved by Algorithm 1 and the trivial
estimator $\hat{\theta}_{0}(X_{1}^{3n})=0$. We plot the loss as a function of
$t=\|\theta_{*}\|$ for $t\in[0,4]$ with step size $0.05$. We take
$n=100,d=5,\delta=0.1$.