# A $\gamma$-Ray Emitting Blazar at Redshift 3.64: Fermi-LAT and OVRO
Observations of PKS 0201+113

Hai Lei  Department of Physics and Astronomy, College of Physics, Guizhou
University, Guiyang 550025, China  Ying-Kang Zhang  Shanghai Astronomical
Observatory, Key Laboratory of Radio Astronomy, 80 Nandan Road, 200030
Shanghai, China.  Xiong Jiang  Department of Physics and Astronomy, College of
Physics, Guizhou University, Guiyang 550025, China  S. Kiehlmann  Institute of
Astrophysics, Foundation for Research and Technology-Hellas, GR-71110
Heraklion, Greece  Department of Physics, Univ. of Crete, GR-70013 Heraklion,
Greece  A. C. S. Readhead  Owens Valley Radio Observatory, California
Institute of Technology, Pasadena, CA 91125, USA  Liang Chen  Shanghai
Astronomical Observatory, Key Laboratory of Radio Astronomy, 80 Nandan Road,
200030 Shanghai, China.  Neng-Hui Liao Department of Physics and Astronomy,
College of Physics, Guizhou University, Guiyang 550025, China  Tao An Shanghai
Astronomical Observatory, 80 Nandan Road, 200030 Shanghai, P. R. China  Key
Laboratory of Radio Astronomy and Technology, Chinese Academy of Sciences, A20
Datun Road, Beijing, 100101, P. R. China

###### Abstract

High-redshift ($z>3$) $\gamma$-ray blazars are rare, but they are crucial for
our understanding of jet evolution, $\gamma$-ray production and propagation,
and the growth of supermassive black holes in the early universe. A new
analysis of Fermi-LAT data reveals a significant (5$\sigma$), spectrally soft
($\Gamma\simeq$ 3.0) $\gamma$-ray source in a specific 4-month epoch,
cospatial with PKS 0201+113 ($z$ = 3.64). Monitoring of PKS 0201+113 at 15 GHz
by the Owens Valley Radio Observatory 40 m Telescope from 2008 to 2023 shows a
prominent flare that dominates the radio light curve. The maximum of the radio
flare coincides with the $\gamma$-ray flare, strongly suggesting an
association ($\textrm{p-value}=0.023$) between the $\gamma$-ray and the radio
sources. PKS 0201+113 is only the third $\gamma$-ray blazar to be identified
with $z>3.5$, and it is the first such object to be identified by the
detection of quasi-simultaneous $\gamma$-ray and radio flares. The jet
properties of this peculiar blazar have been investigated. A detailed study of
a two-zone leptonic model is presented that fits the broadband spectral energy
distribution. An alternative scenario is also briefly discussed.

galaxies: active – galaxies: high redshift – galaxies: jets – gamma rays:
galaxies – quasars: individual (PKS 0201+113)

##  1 Introduction

High-redshift blazars, with Doppler-boosted emission from relativistic jets
aligned close to the line of sight (Rees, 1967; Readhead et al., 1978;
Blandford & Rees, 1978; Cohen et al., 1979; Readhead, 1980; Blandford et al.,
2019), provide luminous beacons from early cosmic epochs (e.g., Romani et al.,
2004; Romani, 2006; Sbarrato et al., 2012; Yi et al., 2014; Ghisellini et al.,
2014). Studies of high-redshift blazars are central to astrophysics, since
they provide key insights into the evolution of relativistic jets, jet
formation, and the dynamics of AGN activity in the early universe (Ghisellini
et al., 2010; Volonteri, 2010). Among thousands of extragalactic $\gamma$-ray
emitters, only ten blazars beyond redshift 3 have been detected by Fermi-LAT
(e.g., Ackermann et al., 2017; Liao et al., 2018), making them particularly
significant. Their detected $\gamma$-ray emission not only confirms their
blazar nature but also provides critical information on the jet properties.
Moreover, they are valuable probes of extragalactic background light (EBL)
across cosmic time (e.g., $Fermi$-LAT Collaboration et al., 2018).

Blazars are notable for broadband violent variability on timescales ranging
from minutes to years (Wagner & Witzel, 1995; Ulrich et al., 1997).
Contemporaneous $\gamma$-ray variations together with those in other windows
of the electromagnetic spectrum have frequently been detected (e.g., Abdo et
al., 2010; Liao et al., 2014; Fuhrmann et al., 2014). Given the relatively low
angular resolution of Fermi-LAT compared to radio and optical bands,
identifying the low-energy counterpart of a $\gamma$-ray source can be
challenging. Particularly for high-redshift blazars, which exhibit soft
$\gamma$-ray spectra, and given that the spatial resolution of sub-GeV photons
is lower than that of GeV photons, model localization constraints from
$\gamma$-ray observations are rather loose (e.g., Liao et al., 2018).
Therefore, the detection of synchronized variations in multiple bands is
essential in order to claim a discovery of a high-redshift $\gamma$-ray
emitting blazar. However, only a few such cases have been reported (e.g., Liao
et al., 2019, 2020; Li et al., 2020).

PKS 0201+113 is a well studied high-redshift ($z$ = 3.64, Jin et al. 2023)
blazar (Massaro et al., 2009). It exhibits high radio and hard X-ray
luminosities (Healey et al., 2008; Marcotulli et al., 2020). Because of its
bright radio emission, PKS 0201+113 has been included in several radio
monitoring programs (Teräsranta et al., 2005; Richards et al., 2011; Lister et
al., 2018). Previous attempts to detect its $\gamma$-ray emission have been
unsuccessful (Ackermann et al., 2017; Marcotulli et al., 2020; Paliya et al.,
2020a). Considering that $\gamma$-ray emitting blazars at high-redshifts (z
$>$ 3.5) are extremely rare, in this study, we carry out a thorough
investigation of the multiwavelength variability properties of PKS 0201+113
(in Section 2), using Fermi-LAT data and data in lower-energy bands, and
discuss the astrophysical implications (in Section 3). Here we take a
$\Lambda$CDM cosmology with $H_{0}=67~{}{\rm km~{}s^{-1}~{}Mpc^{-1}}$,
$\Omega_{\rm m}=0.32$, and $\Omega_{\Lambda}=0.68$ (Planck Collaboration et
al., 2014).

##  2 Data Analysis and Results

###  2.1 Fermi-LAT data

We have collected the archival Fermi-LAT data (P8R3_SOURCE_V3) spanning 2008
August 4 to 2023 February 25 (i.e., MJD 54683 to 60000), covering an energy
range of 100 MeV to 500 GeV. We employed Fermitools software (version 2.0.8),
along with Fermitools-data (version 0.18) for this study. During the
procedure, filtering cuts, such as zenith angle $<90^{\circ}$ and DATA_QUAL>0
&& LAT_CONFIG==1, were applied to the photon data by the gtselect and gtmktime
tasks, respectively. The $\gamma$-ray flux extraction was then executed using
the gtlike task, applying an Unbinned likelihood approach. The test statistic
(TS, TS = 2$\Delta\log\zeta$, Mattox et al. 1996) method was used to assess
the significance of the $\gamma$-ray source, where $\zeta$ are the maximum
likelihood values of different models with and without the source. Background
4FGL-DR3 sources (Abdollahi et al., 2022) within a 15$\arcdeg$ radius of PKS
0201+113, together with the galactic and extragalactic $\gamma$-ray diffuse
components (i.e., gll_iem_v07.fits and iso_P8R3_SOURCE_V3_v1.txt), were
considered during the flux extraction. We only adjusted the parameters of
inner region background sources (i.e., within a radius of $10^{\circ}$ from
the target), and the normalization of the diffuse templates, while keeping the
other parameters fixed at the default values. Residual TS maps were generated
to check for additional background sources not in 4FGL-DR3. Because of the
variability of most extragalactic $\gamma$-ray sources, and Fermi-LAT’s
sensitivity dependence on exposure length, weak background sources with TS $<$
10 were excluded from the $\gamma$-ray temporal analysis.

An initial analysis of the entire dataset suggested the existence of a
$\gamma$-ray source (TS = 53), not included in 4FGL-DR3, located at
$0.43^{\circ}$ from the radio position of PKS 0201+113. Given the $\gamma$-ray
localization radius of $0.11^{\circ}$ at 95% confidence level (C.L.), PKS
0201+113 is unlikely the low-energy counterpart of this particular
$\gamma$-ray source. Additionaly, several other possible $\gamma$-ray sources
were also identified, but they are relatively faint (TS $\leq$ 100) and
distant (more than $5^{\circ}$ away). After including these sources, we added
a test $\gamma$-ray source, representative of PKS 0201+113 and modeled with a
power-law function assumed as its spectral template, into the analysis source
model for a subsequent likelihood analysis. The resulting TS value of $\simeq$
1 for the test source indicates no evidence for the detectable $\gamma$-ray
emission from PKS 0201+113 under these circumstances. This null result is
consistent with those reported in the literature (Ackermann et al., 2017;
Marcotulli et al., 2020; Paliya et al., 2020a).

We then extracted a 3-month time bin $\gamma$-ray light curve to search for
any activity potentially obscured by background emissions. Upper limits at 95%
C.L. level were calculated for low TS values ($<$ 10) using the pyLikelihood
UpperLimits tool. As shown in Figure 1, one time bin is distinguished by a
significant TS value ($\simeq$ 25). Since the angular resolution of Fermi-LAT
for sub-GeV photons is as large as $\gtrsim~{}5^{\circ}$, possible
contamination from nearby bright background sources has been investigated.
$\gamma$-ray light curves of these sources have been extracted as well. The
$\gamma$-ray light curves of these sources show no concurrent TS value
increases, suggesting the $\gamma$-ray brightening observed in direction of
PKS 0201+113 is not attributable to nearby sources. A subsequent zoomed-in
1-month time bin light curve, shown in Figure 1, marks a 4-month period (i.e.,
MJD 55980 to 56100) of potential $\gamma$-ray activity.

Analysis focusing on this specific 4 months of data reveals a significant
$\gamma$-ray source cospatial with PKS 0201+113, see Figure 2. The optimized
localization of this $\gamma$-ray source is R.A. 30.892∘ and DEC. 11.54∘, with
a 95% C. L. error radius of 0.29∘. The angular separation between the
$\gamma$-ray location and radio position of PKS 0201+113 is 0.07∘. No
additional blazar candidates (Healey et al., 2007; Massaro et al., 2009) are
found within the $\gamma$-ray error radius.

Adopting the updated $\gamma$-ray position, a single power-law function
provides an acceptable description of the source. The photon flux is estimated
to be $\rm(5.6\pm 1.1)\times 10^{-8}$ ph $\rm cm^{-2}$ $\rm s^{-1}$ and the
corresponding TS value is 43 (5.7$\sigma$). Considering that the 4-month epoch
was selected from a span of 14.5 years (i.e., 44 trials), the global
significance after trial factor correction is 5.0$\sigma$. Our analysis shows
that the $\gamma$-ray source is spectrally soft (i.e., $\Gamma$ = 2.98 $\pm$
0.23), which supports the relatively loose localization constraints on our
model. If this $\gamma$-ray source is indeed associated with PKS 0201+113, the
corresponding $K$-corrected apparent isotropic $\gamma$-ray luminosity between
0.1 and 500 GeV is $10^{49}$ erg $\rm s^{-1}$.

The influence on the transient source, of the nearby neighbor from the
analysis of the entire dataset has been investigated as follows. In addition
to the source TS Map shown in Figure 2, a residual TS map is shown in Figure
3, were the contributions of the transient source and the background sources
apart from the nearby neighbor have been subtracted. No significant excess (TS
= 1) due to the nearby neighbor source is seen, suggesting that, at least in
this specific epoch, the contribution of the nearby neighbor is not
significant. In addition, since the source corresponding to PKS 0201+113 is
highly time sensitive, we have analyzed the entire 14-yr Fermi-LAT data set,
leaving out this 4-month interval, to check whether there is a reduction of
the TS of the nearby neighbor. The analysis yields a TS value of 48, which is
comparable with that derived from the intact dataset (i.e. TS = 53). This
confirms that the existence of the neighbor is irrelevant to the detection of
the $\gamma$-ray source spatially coincident with PKS 0201+113.

We divided the entire dataset into three parts to investigate the long-term
$\gamma$-ray variability. The analysis of the data, prior to the 4-month epoch
(totally 3.5-year, see Figure 2), suggests a tentative detection in this
direction (i.e., TS = 19). The limited statistics preclude precise
determination of the spectral index. The flux level during this period is
moderate, $\rm(6.3\pm 1.8)\times 10^{-9}$ ph $\rm cm^{-2}$ $\rm s^{-1}$, about
one order of magnitude lower than that of the high flux state. These results
are consistent with a nearby $\gamma$-ray source, 3FGL J0203.6+1148, as
reported in the 3FGL catalog (Acero et al., 2015). However, no low-energy
counterpart for this 3FGL source is identified. Our updated localization
likely benefits from the latest galactic $\gamma$-ray diffuse template.
Meanwhile, analysis of the data after the flare epoch shows no sign of an
excess above the background emission (TS $\leq$ 1), see Figure 2. The target
remains quiescent for more than 10 years, and a 95% C. L. upper limit of $\rm
2.1\times 10^{-9}$ ph $\rm cm^{-2}$ $\rm s^{-1}$ is found. These findings
strongly suggest that the $\gamma$-ray source exhibits substantial long-term
variability. In addition, a weekly time bin light curve focused on the flare
epoch was also analyzed, but reveals no evidence of variability on timescales
of several days.

###  2.2 Zwicky Transient Facility Light-curve Data

We collected the latest ZTF data (Graham et al., 2019; Bellm et al., 2019;
Masci et al., 2019; PTF Team, 2020)111https://www.ztf.caltech.edu/ztf-public-
releases.html to analyze the target. Due to its faintness in $g$-band and
sparse sampling in $i$ band, we extracted $r$-band light curve data only for
objects within a 5″ radius of the target position based on the co-added
reference images. Only frames with labels of catflags = 0 and chi < 4 were
chosen. The snapshots were then aggregated into a daily time bin light curve,
containing 342 nights of observations from MJD 58301 to 58803, shown in Figure
2. No significant optical activity was found, and the standard deviation of
the magnitudes ($\simeq$ 0.11 mag) is comparable to the averaged measurement
uncertainty ($\simeq$ 0.10 mag).

###  2.3 WISE data

The Wide-field Infrared Survey Explorer (_WISE_ ; Wright et al. 2010) has
performed repetitive all-sky surveys in the infrared since 2010 (WISE Team,
2019; Mainzer et al., 2014; WISE Team, 2020), except for a gap from 2011 to
2013. Time-resolved WISE/NEOWISE co-added magnitudes in the _W1_ and _W2_
bands, centered at 3.4, 4.6 $\mu$m (observed), were
adopted222https://github.com/fkiwy/unTimely_Catalog_explorer (Lang, 2014;
Meisner et al., 2018, 2023). The dataset comprises 17 measurements in each
infrared band between MJD 55215 and 59056, sampling every half a year except
during the gap, see Figure 2. We applied a $\chi^{2}$-test to investigate the
infrared variability, following the method of Liao et al. (2019). The analysis
reveals no significant variability ($<$ 5$\sigma$) in either of the infrared
light curves.

###  2.4 Radio Data

####  2.4.1 Owens Valley Radio Observatory light curve

The OVRO 40m telescope has been monitoring over 1800 objects, on a $\sim$
3-day cadence at a frequency of 15 GHz (Richards et al., 2011) since 2008.
OVRO provides an extensive radio light curve for PKS 0201+113 spanning 15
years, see Figure 2. The remarkable feature in the OVRO light curve is a
prominent radio flare peaking around MJD 56000. At the beginning of the OVRO
observations, the radio flux density level is $\simeq$ 0.55 Jy. This is
followed by a drop to 0.40 Jy at MJD 55000, followed by a major radio flare.
The flare lasts roughly 1000 days and reaches a maximum of 0.82 Jy at MJD
55974. Despite the sparse observations during the highest flux state, it is
clear that the flare flux density decline starts around MJD 56109 (0.83 Jy),
and it takes roughly 185 days to drop to 0.55 Jy at MJD 56294. The flux
density then decays gradually over the next $\sim$ 2300 days, reaching 0.33 Jy
by MJD 58582. Subsequently, moderate activity, with peak fluxes lower than 0.5
Jy, is seen.

The 15 GHz light curve around the time of the peak of the giant radio flare
($\gtrsim$ 0.5 Jy, between MJD 55600 and 56200) was decomposed into a series
of sub-flares. The sub-flares are assumed to have an exponential shape, and
they are superposed on a constant baseline flux density (Savolainen et al.,
2002) of 0.34 Jy. This corresponds to the observed minimum value at OVRO. As
shown in Figure 5, a set of four sub-flares, together with the baseline
component, fits the radio data. The peak times of the sub-flares are MJD
55743, MJD 55833, MJD 55972 and MJD 56117, respectively, with uncertainties of
several days. The rise times of the sub-flares range from a few months to
nearly one year (observed), consistent with the results from studies on radio
variability of Fermi blazars (e.g., León-Tavares et al., 2011).

####  2.4.2 VLBI radio images

We selected the X-band Astrogeo data333http://astrogeo.org (Petrov et al.,
2009), focusing on epochs starting in 2011 to match the burst period around
2012. Hybrid imaging processing of the calibrated (u,v) visibility data and
subsequent model fitting were performed using DIFMAP (Shepherd, 1997). The
corresponding parameters are presented in Table 1. Analyses of both the images
and the visibility data fitting reveal that the source exhibits a typical
blazar-like compact core-jet structure, similar to other high-redshift blazars
(e.g. J0906+6930 from Zhang et al., 2017). The VLBI data indicate the source
comprises a bright core, accounting for most of the flux density, and a
compact jet extending towards the northwest. Due to the insufficient (u,v)
sampling from the geodetic scans and the small number of available epochs
around 2012, tracking the motion of the jet components based on model-fitting
results from these four epochs ws not possible. Thus only total flux densities
and those of the core component are given. The uncertainties are estimated by
combining the error propagation Fomalont (1999) with an additional 5% of the
respective flux densities.

###  2.5 Implications of $\gamma$-ray and radio variability

A close connection between centimeter radio (i.e., observed at 8.4 and 15 GHz)
and $\gamma$-ray emission of Fermi-LAT detected blazars has been recognized
(e.g., Kovalev et al., 2009; Ackermann et al., 2011). The relationship holds
for mm/sub-mm radio emission (León-Tavares et al., 2012). These correlations
imply that the $\gamma$-ray emission likely originates from the parsec scale
core-jet. Benefiting from extensive radio monitoring programs (e.g., Aller et
al., 1999; Teräsranta et al., 2005; Richards et al., 2011; Fuhrmann et al.,
2016), together with the wide field of view of Fermi-LAT, comparisons of
temporal behaviors in the two bands have been investigated. A non-zero time
delay between VLBA 15 GHz radio emission and $\gamma$-ray emission for Fermi
blazars, otherwise known as $\gamma$-ray leading radio, has been detected in
other objects (e.g., Pushkarev et al., 2010). The observed time lag is
typically months to years, taking into account the width of the cross-
correlation functions (e.g., León-Tavares et al., 2011; Max-Moerbeck et al.,
2014). Interestingly, a tendency that the lag (source rest frame) decreases
from cm to mm/sub-mm bands is reported, and no significant lag has been found
in stacked 3mm/$\gamma$-ray emissions (Fuhrmann et al., 2014).

In PKS 0201+113, the time of emergence of the $\gamma$-ray source (i.e., MJD
55980 and 56100) coincides, to within the uncertainties, with the peak of the
giant radio flare (i.e., MJD 55970 and 56109). Given that the timespan of the
$\gamma$-ray and radio light curves is $>$ 5000 days, the chance probability
of flares occurring in both bands in the same 4-month interval is $\sim$ 1/44
= 0.023, which suggests a causal connection, and hence an association between
the $\gamma$-ray source and the radio source. We therefore identify PKS
0201+113 as a $\gamma$-ray emitting blazar. Due to the relatively limited
statistics of the $\gamma$-ray data, cross correlation analysis of the
$\gamma$-ray and OVRO light curves cannot be carried out. Nevertheless, the
time bin with the largest TS value, in the monthly $\gamma$-ray light curve
(Figure 2), is centered at MJD 56055 when the radio emission remains at a high
flux level. We note that, during the period when we have optical data, the
optical emission does not vary and the radio variation is also small.

The most significant observational result of this paper is the apparent near-
coincidence in time between the $\gamma-$ray flare and the radio flare. This
establishes the association between these two flares, i.e. that they both
originate in the same blazar. Furthermore, it shows that the emission at these
two frequencies is physically linked. The data are insufficient to determine
the exact relationship between these flares, i.e. whether the peak of the
radio flare lags or precedes the $\gamma$-ray flare, although it is clear that
the radio flux starts to rise around MJD 55000, well before the peak of the
$\gamma-$ray flare.

The location of the jet dissipation region in blazars, where the bulk energy
of the jet is converted into relativistic particles, remains an unresolved
aspect of blazar physics. Rapid $\gamma$-ray variability, with timescales as
short as minutes, indicates an extremely compact emitting region (e.g.,
Aharonian et al., 2007; Ackermann et al., 2016). Assuming a conical jet
geometry, this region is likely to be in close proximity to the central
supermassive black hole (SMBH). In that case, significant absorption of
$\gamma$-ray emission via the $\gamma\gamma$ process could occur (Blandford &
Levinson, 1995; Begelman et al., 2008; Poutanen & Stern, 2010). On the other
hand, the location of the $\gamma-$ray emission region relative to the radio
emission region can be determined via the correlated radio/$\gamma$-ray
emission. The observed time lag between the observed emission in these two
energy ranges is directly related to the distance between the $\gamma$-ray
production site (i.e., $r_{\gamma}$) and the radio core (i.e., $r_{c}$), where
the opacity of self-synchrotron absorption (SSA) is unity. The latter is
inferred by the core shift measure (Lobanov, 1998). Based on this approach,
the location of the jet dissipation region could be several parsecs away from
the SMBH (e.g., Kramarenko et al., 2022).

For PKS 0201+113, unfortunately, no simultaneous VLBI images corresponding to
the radio/$\gamma$-ray flares are available. Hence the size of the radio core
is poorly determined. In addition, the time lag between the OVRO and
$\gamma$-ray light curves is not well determined. This makes precise
localization of the jet dissipation region impossible. Nevertheless,
contemporaneous brightening between the radio and $\gamma$-ray emission
provides a critical clue. If the dissipation region is deeply embedded within
the broad line region, a compact dissipation region leading to dense
synchrotron photons and severe SSA would weaken the relationship between the
emission in the $\gamma-$ray and radio bands. Conversely, if the region is
significantly distant from the SMBH, beyond the edge of the dust torus, the
efficiency of the external Compton (EC) process would decrease rapidly,
complicating the explanation of the luminous $\gamma$-ray emission.

As mentioned above, and in spite of the papers cited above, it remains the
case that the phenomenology of the correlation of $\gamma$-ray and radio
variations in blazars is still largely unknown. We know, from radio VLBI maps,
that the radio emission originates in extended regions along the blazar jet,
and we also know that any distance between the $\gamma$-ray emission site and
the radio emission site that is deduced from time delays must be multiplied by
$2\gamma^{2}$ (Rees, 1966, 1967) for blazars aligned within an angle
$\theta<1/\gamma$ to the line of site. In the two cases (i.e., AO 0235+164 and
B2 2308+34) where clear correlations were seen by Max-Moerbeck et al. (2014),
we think that an increase of activity in the blazar, caused by for example a
sudden influx of new fueling (Falcke & Biermann, 1995; Gutiérrez et al.,
2023), leads to an increase of activity and of emission which is seen first at
lower energies (i.e., radio), and then, as the activity builds up it
eventually leads to $\gamma$-ray emission. Thus parts of the radio flare can
precede the $\gamma$-ray flare.

##  3 Discussions and Summary

Among the thousands of sources listed in Data Release 3 of the Fourth Catalog
of AGN Detected by the Fermi-LAT (Ajello et al., 2022), only three are located
beyond a redshift of 3.5, including GB 1508+5714 ($z$ = 4.3), PKS 1351-018
($z$ = 3.72) as well as SDSS J163547.23+362930.0 ($z$ = 3.62). Additional
sources identified in individual studies are PMN J2219-2719 ($z$ = 3.63, Li et
al. 2020) and NVSS J121915+365718 ($z\sim$ 3.59, Paliya et al. 2020b). The
identification of PKS 0201+113 ($z$ = 3.64) as a $\gamma$-ray emitter places
it among the top five distant $\gamma$-ray sources known so far. Since the
energy of horizon photons (i.e., $\rm\tau_{\gamma\gamma}^{EBL}$ = 1) at such a
redshift is $\simeq$ 30 GeV (e.g., Finke et al., 2010), while the most
energetic photons from PKS 0201+113 are $\sim$ 1 GeV, current EBL models are
not inconsistent with our findings. The observed rather soft $\gamma$-ray
spectrum of PKS 0201+113 is consistent with those of other known high-redshift
$\gamma$-ray blazars (e.g., Li et al., 2018). Since PKS 0201+113 is not
detected by Fermi-LAT when considering the entire dataset, its flux level is
obviously lower than most known high-redshift $\gamma$-ray blazars. A similar
case is that of PMN J2219-2719, for which a Fermi-LAT detection is only
possible when choosing a specific 5-month flaring epoch (Li et al., 2020).
Large $\gamma$-ray variability amplitudes, concurrent with low-energy
brightening can play a crucial role in proving an association with a
$\gamma$-ray source. Another noteworthy case is B3 1428+422 at z = 4.72, where
a co-spatial, spectrally soft $\gamma$-ray source was reported during a
10-month epoch, yet no detections were found outside this period (Liao et al.,
2018). In this case, the absence of simultaneous observations in other energy
bands prevents the confirmation of a $\gamma$-ray association.

Multiwavelength campaigns for high-redshift blazars, particularly those above
a redshift of 3.5, are challenging due to the faintness of the sources when
they are located at high redshifts. Beyond such a distance, so far
$\gamma$-ray emissions of GB 1508+5714 and PMN J2219-2719 have been identified
by contemporaneous brightenings between emissions in $\gamma$-ray and
optical/near-infrared regimes (Liao et al., 2020; Li et al., 2020). In fact,
PKS 0201+113 is the first high-redshift blazar whose $\gamma$-ray emission has
been identified with the aid of a radio monitoring program, which is
remarkable, given that it is one of the brightest radio sources in the early
universe ($\simeq$ 1 Jy observed at 5 GHz). The identification of $\gamma$-ray
and radio variability in PKS 0201+113 provides a unique perspective on high-
redshift blazars. Radio and $\gamma$-ray emission in blazars originates almost
solely from the plasma jet without significant contributions from other AGN
components. However, the connection between these emission regions is complex.
First: radio emission at long wavelengths tends to suffer significant SSA,
which is also influenced by the distance of radiating region from the SMBH,
and by the strength of the magnetic field. In addition, unlike the high energy
emission, which varies rapidly, the radio emission is believed to emanate from
an extended region, and it varyies over months to years (e.g., Ciaramella et
al., 2004; Richards et al., 2014; Fuhrmann et al., 2016). Nevertheless, as
seen in PKS 0201+113 and a few other blazars, quasi-simultaneous radio and
$\gamma$-ray flares (Max-Moerbeck et al., 2014; Karamanavis et al., 2016; Liao
et al., 2016; Wang et al., 2023; Cheong et al., 2024) suggest a possible close
connection between the radio and the $\gamma$-ray emission. Several other
high-redshift objects are noteworthy. CGRaBS J0733+0456 ($z$ = 3.0, Liao et
al. 2019), which is in the OVRO monitoring program, exhibits unprecedented
quasi-simultaneous IR and $\gamma$-ray flares during the ascent phase of the
radio emission. QSO J0906+6930, once proposed as the farthest blazar ($z$ =
5.48, Romani et al. 2004; Romani 2006; An et al. 2020), shows no significant
$\gamma$-ray emission, but a tentative signal in a specific 11-day epoch has
been reported (An & Romani, 2018), coinciding with its peak OVRO radio flux
density (Zhang et al., 2017). In addition, in the case of the blazar candidate
PSO J047.4478+27.2992 ($z$ = 6.1), the RATAN-600 telescope has observed mild
radio variability (Mufakharov et al., 2021). These cases show the important
role that the radio monitoring programs play in exploring faint non-thermal
transients in the early Universe.

The broadband spectral energy distribution (SED) of PKS 0201+113 is presented
in Figure 5. The SED of a blazar exhibits a characteristic two-bump shape
(Blandford et al., 2019). The first bump is attributed to synchrotron emission
from relativistic electrons spiraling in a magnetic field, while the second,
peaking in the $\gamma$-ray range, is generally attributed to inverse Compton
(IC) emission. In the “single zone” model these processes involve the same
electron population either scattering virtual photons, in the synchrotron case
(Blumenthal & Gould, 1970), or soft photons originating from within or from
outside the jet, in the IC case (e.g. Maraschi et al., 1992; Dermer &
Schlickeiser, 1993; Sikora et al., 1994; Błażejowski et al., 2000). We do not
have multi-frequency radio, infrared, optical, or X-ray observations covering
the flares seen at the OVRO at 15 GHz (0.81 Jy at MJD 55990) and on the Fermi-
LAT (between MJD 55980 and 56100), so it is not possible to construct a
simultaneous SED. We therefore use archival multi-frequency radio measurements
taken from NED444https://ned.ipac.caltech.edu/, and data from NuSTAR and
Swift-XRT observations at MJD 57700 (Marcotulli et al., 2020), ALLWISE and
SDSS magnitudes (WISE Team, 2019; Ahumada et al., 2020) to construct the SED
shown in Figure 5. Based on these data, the synchrotron peak frequency of PKS
0201+113 in its quiescent state likely falls in the sub-mm/infrared domain, a
characteristic typical of Flat Spectrum Radio Quasars (FSRQs) (Madejski &
Sikora, 2016). The optical/UV emission is dominated by the accretion disk,
indicative of a central SMBH with $\rm\sim$ 3$\times 10^{9}$ $\rm M_{\sun}$
(Marcotulli et al., 2020).

Since the duration of the $\gamma-$ray flare in PKS 0201+113 is $\lesssim 30$
days, whereas the ascent phase of 15 GHz radio flare alone lasts $\sim 3$ yr,
it is clear that the full extent of the radio emission regions is much more
extended than that of the $\gamma-$ray emission region. We therefore adopt a
two-zone leptonic model, based on the central sub-flare of the radio light
curve, to describe the broadband SED. In this two-zone model an inner compact
region radiating at $\gamma-$ray energies is embedded in an outer extended
radio emission region, as shown in Figure 6. The radius of the outer region is
assumed to be $\sim 2\times 10^{18}$ cm, while that inner region is assumed to
be $\sim 2\times 10^{16}$ cm. The inner region could, for example, be a mini-
jet, or a narrow beam, which might correspond to a plasmoid in relativistic
reconnection scenarios. This model could explain any ultra-fast $\gamma$-ray
variability (e.g., Giannios et al., 2009; Petropoulou et al., 2016). Multiple
mini-jets have been evoked to describe complex variability behavior (Marscher,
2014). For a review of these scenarios see Marscher (2016). A different
structured jet model, the spine-jet model, successfully explains the luminous
GeV/TeV $\gamma$-ray emission of radio galaxies (e.g., Ghisellini et al.,
2005; Chen, 2017). In these structured jet models the inner region moves
faster than the outer region (Marscher, 2016). In our two-zone leptonic model,
it is similarly assumed that the inner region has a higher bulk Lorentz factor
than the outer region. In each region, the emitting electrons are assumed to
be homogeneously distributed and to follow a broken power-law spectrum,

$N(\gamma)\propto\left\\{\begin{array}[]{ll}\gamma^{-p_{1}}&\mbox{
$\gamma_{\rm min}\leq\gamma\leq\gamma_{br}$}\\\ \gamma_{\rm
br}^{p_{2}-p_{1}}\gamma^{-p_{2}}&\mbox{ $\gamma_{\rm br}<\gamma\leq\gamma_{\rm
max}$, }\end{array}\right.$ (1)

where $\rm\gamma_{br}$, $\rm\gamma_{min}$ and $\rm\gamma_{max}$ are the break,
the minimum and maximum energies of the electrons, and $p_{1,2}$ are indices
of the broken power-law particle distribution. In calculating the total
radiation output of PKS 0201+113, both synchrotron and IC emission has been
taken into account for both the inner and the outer emission regions. Our IC
analysis includes both synchrotron self-Compton (SSC) and EC processes within
each region. Additionally, IC scattering of synchrotron photons originating
from one region by electrons in the other region is also taken into account.
In the analysis, two important physical effects are properly incorporated: SSA
and the Klein–Nishina effect in the IC scattering region. The former is
crucial for understanding the radio emission, especially at lower frequencies,
while the latter is significant at high energies, affecting the IC scattering
cross-section. The transformations of frequency and luminosity between the jet
rest frame and the observers frame are $\nu=\delta\nu^{\prime}/(1+z)$ and $\nu
L_{\nu}=\delta^{4}\nu^{\prime}L^{\prime}_{\nu^{\prime}}$. These
transformations are essential for reproducing the observed spectrum, as they
account for the relativistic effects due to the jet motion and the expansion
of the Universe. The relative Lorentz gamma factor, ($\Gamma_{rel}$), between
the inner and outer emission regions is
$\Gamma_{rel}=\Gamma_{in}\Gamma_{out}(1-\beta_{in}\beta_{out})$ (Ghisellini et
al., 2005). where the individual Lorentz gamma factors of the two regions are
($\Gamma_{in},\Gamma_{out}$) and the corresponding velocities are
($\beta_{in},\beta_{out}$). This aspect of the model is vital for
understanding the dynamics within the jet, and how this affects the observed
emission.

The leptonic scenario gives an explanation of the high flux state SED of PKS
0201+113, assuming the parameters listed in Table 2. Based on the radius and
the Doppler factor of the outer emission region, the variability timescale in
radio band is $t_{obs,var}\gtrsim(1+z)R_{out}/c\delta$, which is consistent
with the OVRO light curve. Assuming a conical jet geometry, the distance
between the emission region and the SMBH is approximated as
$r_{\gamma}\sim\delta R_{out}\simeq$ 2 pc, in agreement with findings in the
literature (e.g., Max-Moerbeck et al., 2014; Kramarenko et al., 2022). The
infrared emission from the dust torus supplies the external photon field, with
an energy density of 2.5 $\times 10^{-4}$ erg $\rm cm^{-3}$ (Ghisellini &
Tavecchio, 2009). On this model the observed OVRO 15 GHz emission is due
primarily to synchrotron radiation in the outer emission region. Based on the
magnetic field strength of 0.09 Gauss, the opacity of SSA at 70 GHz in source
rest frame (17 GHz in the jet rest frame) is found to be $\lesssim 1$, which
allows for a tight connection between the radio and $\gamma-$ray emission.
Note that, at the observed frequency of 15 GHz the spectral index is
$\alpha\sim-0.6$, confirming that this is optically thin emission, and we
therefore expect a tight connection between the radio and $\gamma-$ray
emission. On the other hand, the contribution from the inner emission region
to the observed OVRO emission is negligible due to severe SSA effects. The
Doppler factor of 28.7, and radius of 2 $\times 10^{16}$ cm, implies an
observed variability timescale of $\simeq$ 1 day. IC scattering by electrons
of the inner emission region, of the infrared radiation from dust, is
responsible for the observed $\gamma$-ray emission. Modeling of such an inner
emission region of PKS 0201+113 in a quiescent state, when no $\gamma$-ray
signal is detected, has been carried out by Marcotulli et al. (2020). The
major difference between these results and our model is our significantly
higher value of the Doppler factor and $\rm\gamma_{br}$. This discrepancy may
be attributed to the ejection of a new emission region, leading to the
injection of fresh emitting particles and consequently to $\gamma$-ray
brightening. This scenario has also been proposed in multi-wavelength studies
for other high-redshift blazars (Liao et al., 2019, 2020).

While the above two-zone model provides a satisfactory explanation of all the
known facts about PKS 0201+113, and therefore deserves serious consideration,
it is important to note that the apparent near-simultaneity of the OVRO 15 GHz
peak flux density and the $\gamma-$ray, peak flux needs to be interpreted with
caution. It is entirely possible that the emission at radio and $\gamma-$ray
energies originates in widely-separated regions, in which case, although it is
highly likely that the flares seen at radio and $\gamma-$ray energies are
causally connected, the above nested two-zone model would not apply. It should
be noted that we cannot say that the observed flares at $\gamma-$ray and radio
energies are exactly “simultaneous” because we are considering 30-day averages
of the $\gamma-$ray flux. In the two-zone model we have broken down the large
flare into sub-flares to address this issue. But the sub-flare scenario, while
plausible, may not be correct. All we can say with certainty is that the peaks
in the $\gamma-$ray and radio light curves occur within the same 30-day
period. Thus the time difference between the $\gamma-$ray and the radio peaks
could be as large as 30 days. This means that the distance between the
$\gamma-$ray emission region and the radio emission region could be
$30/(1+3.64)\times 2\gamma^{2}c$ lt dy. For not unreasonable values of
$\gamma$ this would place the $\gamma-$ray emission site outside the radio
emission site.

A more general concern is that the assumption that the peak $\gamma-$ray
emission and the peak radio emission should occur at the same time is an
assumption that is not straightforward to justify. One could have a model, for
example, in which increased activity is brought about by a slow and steady
increase in fueling. On such a model, one could easily imagine that this would
be seen first at low (radio) energies, so that the radio flux begins to
increase first. As more fuel is added, the activity level increases, adding
steadily to the radio flux, and initiating higher energy emission. On this
model, it may be that the level of activity only becomes high enough to cause
$\gamma-$ray emission some time after the increase in fueling begins. There is
no guarantee on such a model that the peaks in the $\gamma-$ray and radio
emission would be simultaneous. The possibility of variations in fueling
leading to enhanced radio activity (Falcke & Biermann, 1995; Gutiérrez et al.,
2023) is frequently discussed in the case of tidal disruption events, and has
recently been invoked to explain the relatively short lifetimes of Type 2
compact symmetric objects (CSO-2s) (Readhead et al., 2024). While this is a
very different case to that of the standard blazar jets discussed in this
paper, it provides a clear instance in which variations in fueling the SMBH
has a dramatic impact on the emission from relativistic jets.

Thus, although we have shown that the observations are consistent with the
two-zone model we have adopted, the 30-day uncertainty in the time of the peak
of the $\gamma-$ray emission allows the possibility that the radio and
$\gamma-$ray emission regions are well separated and are not nested, as is
assumed in our two-zone model.

In summary, our comprehensive multi-wavelength analysis of the high-redshift
blazar PKS 0201+113 reveals several key findings. We identified a transient
$\gamma$-ray source in a specific 4-month epoch by re-examining Fermi-LAT
data. This $\gamma$-ray source, located at the radio position of PKS 0201+113
and with no other blazar candidates in the vicinity, has a TS value of 43,
corresponding to a global $5\sigma$ detection significance. The $\gamma$-ray
source exhibits a spectrally soft spectrum, consistent with the behavior of
other known high-redshift $\gamma$-ray emitting blazars. Subsequent Fermi-LAT
data over a decade following this epoch showed no detectable signals, which
indicates $\gamma$-ray variability exceeding a factor 20. A prominent radio
flare was detected in the OVRO monitoring program in which the radio flux
density reaches the maximum value approximately when the $\gamma$-ray flare is
seen. The detection of quasi-contemporaneous radio/$\gamma$-ray flares
strongly suggests a causal connection between the emission in these two energy
bands, thus providing crucial evidence of an association between the
$\gamma$-ray source and PKS 0201+113. This is the third source identified with
$\gamma$-ray emission above a redshift of 3.5, and the first to be identified
with the aid of a radio monitoring program, as opposed to optical/infrared
band observations. A two-zone leptonic model effectively reproduces the
broadband SED of PKS 0201+113, with reasonable input parameters. Future
multiwavelength campaigns, initiated by radio time-domain surveys and
supplemented by simultaneous observations from infrared to $\gamma$ rays, will
provide further constraints on model jet properties of this intriguing high-
redshift blazar.

We appreciate the instructive suggestions from the anonymous referee that led
to a substantial improvement of this work. Talvikki Hovatta is appreciated for
the helpful comments. Lea Marcotulli is appreciated for sharing the
observational SED data of PKS 0201+113. This research has made use of data
obtained from the High Energy Astrophysics Science Archive Research Center
(HEASARC), provided by NASA’s Goddard Space Flight Center. This research has
made use of the NASA/IPAC Infrared Science Archive, which is funded by the
NASA and operated by the California Institute of Technology. This research has
made use of data from the OVRO 40-m monitoring program (Richards, J. L. et al.
2011, ApJS, 194, 29), supported by private funding from the California
Institute of Technology and the Max Planck Institute for Radio Astronomy, and
by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants
AST-0808050 and AST-1109911. This work was supported in part by the National
SKA Program of China (No. 2022SKA0120102, 2022SKA0130103) and NSFC under
grants U2031120 and 11703093. This work was also supported in part by the
Special Natural Science Fund of Guizhou University (grant No. 201911A) and the
First-class Physics Promotion Programme (2019) of Guizhou University. SK
acknowledges support from the European Research Council (ERC) under the
European Unions Horizon 2020 research and innovation programme under grant
agreement No. 771282.

## References

  * Abdo et al. (2010) Abdo, A. A., et al. 2010, Nature, 463, 919, doi: 10.1038/nature08841
  * Abdollahi et al. (2022) Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53, doi: 10.3847/1538-4365/ac6751
  * Acero et al. (2015) Acero, F., et al. 2015, Astrophys. J. Suppl., 218, 23, doi: 10.1088/0067-0049/218/2/23
  * Ackermann et al. (2011) Ackermann, M., Ajello, M., Allafort, A., et al. 2011, ApJ, 741, 30, doi: 10.1088/0004-637X/741/1/30
  * Ackermann et al. (2016) Ackermann, M., Anantua, R., Asano, K., et al. 2016, ApJ, 824, L20, doi: 10.3847/2041-8205/824/2/L20
  * Ackermann et al. (2017) Ackermann, M., Ajello, M., Baldini, L., et al. 2017, ApJ, 837, L5, doi: 10.3847/2041-8213/aa5fff
  * Aharonian et al. (2007) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, ApJ, 664, L71, doi: 10.1086/520635
  * Ahumada et al. (2020) Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3, doi: 10.3847/1538-4365/ab929e
  * Ajello et al. (2022) Ajello, M., Baldini, L., Ballet, J., et al. 2022, ApJS, 263, 24, doi: 10.3847/1538-4365/ac9523
  * Aller et al. (1999) Aller, M. F., Aller, H. D., Hughes, P. A., & Latimer, G. E. 1999, ApJ, 512, 601, doi: 10.1086/306799
  * An & Romani (2018) An, H., & Romani, R. W. 2018, ApJ, 856, 105, doi: 10.3847/1538-4357/aab435
  * An et al. (2020) An, T., Mohan, P., Zhang, Y., et al. 2020, Nature Communications, 11, 143, doi: 10.1038/s41467-019-14093-2
  * Begelman et al. (2008) Begelman, M. C., Fabian, A. C., & Rees, M. J. 2008, MNRAS, 384, L19, doi: 10.1111/j.1745-3933.2007.00413.x
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002, doi: 10.1088/1538-3873/aaecbe
  * Blandford et al. (2019) Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467, doi: 10.1146/annurev-astro-081817-051948
  * Blandford & Levinson (1995) Blandford, R. D., & Levinson, A. 1995, ApJ, 441, 79, doi: 10.1086/175338
  * Blandford & Rees (1978) Blandford, R. D., & Rees, M. J. 1978, in BL Lac Objects, ed. A. M. Wolfe, 328–341
  * Błażejowski et al. (2000) Błażejowski, M., Sikora, M., Moderski, R., & Madejski, G. M. 2000, ApJ, 545, 107, doi: 10.1086/317791
  * Blumenthal & Gould (1970) Blumenthal, G. R., & Gould, R. J. 1970, Reviews of Modern Physics, 42, 237, doi: 10.1103/RevModPhys.42.237
  * Chen (2017) Chen, L. 2017, ApJ, 842, 129, doi: 10.3847/1538-4357/aa7744
  * Cheong et al. (2024) Cheong, W. Y., Lee, S.-S., Kim, S.-H., et al. 2024, MNRAS, 527, 882, doi: 10.1093/mnras/stad3250
  * Ciaramella et al. (2004) Ciaramella, A., Bongardo, C., Aller, H. D., et al. 2004, A&A, 419, 485, doi: 10.1051/0004-6361:20035771
  * Cohen et al. (1979) Cohen, M. H., Pearson, T. J., Readhead, A. C. S., et al. 1979, ApJ, 231, 293, doi: 10.1086/157192
  * Dermer & Schlickeiser (1993) Dermer, C. D., & Schlickeiser, R. 1993, ApJ, 416, 458, doi: 10.1086/173251
  * Falcke & Biermann (1995) Falcke, H., & Biermann, P. L. 1995, A&A, 293, 665, doi: 10.48550/arXiv.astro-ph/9411096
  * $Fermi$-LAT Collaboration et al. (2018) $Fermi$-LAT Collaboration, Abdollahi, S., Ackermann, M., et al. 2018, Science, 362, 1031, doi: 10.1126/science.aat8123
  * Finke et al. (2010) Finke, J. D., Razzaque, S., & Dermer, C. D. 2010, ApJ, 712, 238, doi: 10.1088/0004-637X/712/1/238
  * Fomalont (1999) Fomalont, E. B. 1999, in Astronomical Society of the Pacific Conference Series, Vol. 180, Synthesis Imaging in Radio Astronomy II, ed. G. B. Taylor, C. L. Carilli, & R. A. Perley, 301
  * Fuhrmann et al. (2014) Fuhrmann, L., Larsson, S., Chiang, J., et al. 2014, MNRAS, 441, 1899, doi: 10.1093/mnras/stu540
  * Fuhrmann et al. (2016) Fuhrmann, L., Angelakis, E., Zensus, J. A., et al. 2016, A&A, 596, A45, doi: 10.1051/0004-6361/201528034
  * Ghisellini et al. (2014) Ghisellini, G., Sbarrato, T., Tagliaferri, G., et al. 2014, MNRAS, 440, L111, doi: 10.1093/mnrasl/slu032
  * Ghisellini & Tavecchio (2009) Ghisellini, G., & Tavecchio, F. 2009, MNRAS, 397, 985, doi: 10.1111/j.1365-2966.2009.15007.x
  * Ghisellini et al. (2005) Ghisellini, G., Tavecchio, F., & Chiaberge, M. 2005, A&A, 432, 401, doi: 10.1051/0004-6361:20041404
  * Ghisellini et al. (2010) Ghisellini, G., Della Ceca, R., Volonteri, M., et al. 2010, MNRAS, 405, 387, doi: 10.1111/j.1365-2966.2010.16449.x
  * Giannios et al. (2009) Giannios, D., Uzdensky, D. A., & Begelman, M. C. 2009, MNRAS, 395, L29, doi: 10.1111/j.1745-3933.2009.00635.x
  * Graham et al. (2019) Graham, M. J., Kulkarni, S. R., Bellm, E. C., et al. 2019, PASP, 131, 078001, doi: 10.1088/1538-3873/ab006c
  * Gutiérrez et al. (2023) Gutiérrez, E. M., Combi, L., Romero, G. E., & Campanelli, M. 2023, arXiv e-prints, arXiv:2301.04280, doi: 10.48550/arXiv.2301.04280
  * Healey et al. (2007) Healey, S. E., Romani, R. W., Taylor, G. B., et al. 2007, ApJS, 171, 61, doi: 10.1086/513742
  * Healey et al. (2008) Healey, S. E., Romani, R. W., Cotter, G., et al. 2008, ApJS, 175, 97, doi: 10.1086/523302
  * Jin et al. (2023) Jin, J.-J., Wu, X.-B., Fu, Y., et al. 2023, ApJS, 265, 25, doi: 10.3847/1538-4365/acaf89
  * Karamanavis et al. (2016) Karamanavis, V., Fuhrmann, L., Krichbaum, T. P., et al. 2016, A&A, 586, A60, doi: 10.1051/0004-6361/201527225
  * Kovalev et al. (2009) Kovalev, Y. Y., Aller, H. D., Aller, M. F., et al. 2009, ApJ, 696, L17, doi: 10.1088/0004-637X/696/1/L17
  * Kramarenko et al. (2022) Kramarenko, I. G., Pushkarev, A. B., Kovalev, Y. Y., et al. 2022, MNRAS, 510, 469, doi: 10.1093/mnras/stab3358
  * Lang (2014) Lang, D. 2014, AJ, 147, 108, doi: 10.1088/0004-6256/147/5/108
  * León-Tavares et al. (2012) León-Tavares, J., Valtaoja, E., Giommi, P., et al. 2012, ApJ, 754, 23, doi: 10.1088/0004-637X/754/1/23
  * León-Tavares et al. (2011) León-Tavares, J., Valtaoja, E., Tornikoski, M., Lähteenmäki, A., & Nieppola, E. 2011, A&A, 532, A146, doi: 10.1051/0004-6361/201116664
  * Li et al. (2020) Li, S., Sun, L.-M., Liao, N.-H., & Fan, Y.-Z. 2020, ApJ, 900, 72, doi: 10.3847/1538-4357/abab03
  * Li et al. (2018) Li, S., Xia, Z.-Q., Liang, Y.-F., Liao, N.-H., & Fan, Y.-Z. 2018, ApJ, 853, 159, doi: 10.3847/1538-4357/aaa3fb
  * Liao et al. (2014) Liao, N. H., Bai, J. M., Liu, H. T., et al. 2014, ApJ, 783, 83, doi: 10.1088/0004-637X/783/2/83
  * Liao et al. (2019) Liao, N.-H., Dou, L.-M., Jiang, N., et al. 2019, ApJ, 879, L9, doi: 10.3847/2041-8213/ab2893
  * Liao et al. (2018) Liao, N.-H., Li, S., & Fan, Y.-Z. 2018, ApJ, 865, L17, doi: 10.3847/2041-8213/aae20d
  * Liao et al. (2020) Liao, N.-H., Li, S., Sheng, Z.-F., & Fan, Y.-Z. 2020, ApJ, 898, L56, doi: 10.3847/2041-8213/aba82e
  * Liao et al. (2016) Liao, N.-H., Xin, Y.-L., Fan, X.-L., et al. 2016, ApJS, 226, 17, doi: 10.3847/0067-0049/226/2/17
  * Lister et al. (2018) Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJS, 234, 12, doi: 10.3847/1538-4365/aa9c44
  * Lobanov (1998) Lobanov, A. P. 1998, A&A, 330, 79, doi: 10.48550/arXiv.astro-ph/9712132
  * Madejski & Sikora (2016) Madejski, G. G., & Sikora, M. 2016, ARA&A, 54, 725, doi: 10.1146/annurev-astro-081913-040044
  * Mainzer et al. (2014) Mainzer, A., Bauer, J., Cutri, R. M., et al. 2014, ApJ, 792, 30, doi: 10.1088/0004-637X/792/1/30
  * Maraschi et al. (1992) Maraschi, L., Ghisellini, G., & Celotti, A. 1992, ApJ, 397, L5, doi: 10.1086/186531
  * Marcotulli et al. (2020) Marcotulli, L., Paliya, V., Ajello, M., et al. 2020, ApJ, 889, 164, doi: 10.3847/1538-4357/ab65f5
  * Marscher (2014) Marscher, A. P. 2014, ApJ, 780, 87, doi: 10.1088/0004-637X/780/1/87
  * Marscher (2016) —. 2016, Galaxies, 4, 37, doi: 10.3390/galaxies4040037
  * Masci et al. (2019) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003, doi: 10.1088/1538-3873/aae8ac
  * Massaro et al. (2009) Massaro, E., Giommi, P., Leto, C., et al. 2009, A&A, 495, 691, doi: 10.1051/0004-6361:200810161
  * Mattox et al. (1996) Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396, doi: 10.1086/177068
  * Max-Moerbeck et al. (2014) Max-Moerbeck, W., Hovatta, T., Richards, J. L., et al. 2014, MNRAS, 445, 428, doi: 10.1093/mnras/stu1749
  * Meisner et al. (2023) Meisner, A. M., Caselden, D., Schlafly, E. F., & Kiwy, F. 2023, AJ, 165, 36, doi: 10.3847/1538-3881/aca2ab
  * Meisner et al. (2018) Meisner, A. M., Lang, D., & Schlegel, D. J. 2018, AJ, 156, 69, doi: 10.3847/1538-3881/aacbcd
  * Mufakharov et al. (2021) Mufakharov, T., Mikhailov, A., Sotnikova, Y., et al. 2021, MNRAS, 503, 4662, doi: 10.1093/mnras/staa3688
  * Paliya et al. (2020a) Paliya, V. S., Ajello, M., Cao, H. M., et al. 2020a, ApJ, 897, 177, doi: 10.3847/1538-4357/ab9c1a
  * Paliya et al. (2020b) Paliya, V. S., Domínguez, A., Cabello, C., et al. 2020b, ApJ, 903, L8, doi: 10.3847/2041-8213/abbc06
  * Petropoulou et al. (2016) Petropoulou, M., Giannios, D., & Sironi, L. 2016, MNRAS, 462, 3325, doi: 10.1093/mnras/stw1832
  * Petrov et al. (2009) Petrov, L., Gordon, D., Gipson, J., et al. 2009, Journal of Geodesy, 83, 859, doi: 10.1007/s00190-009-0304-7
  * Planck Collaboration et al. (2014) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A16, doi: 10.1051/0004-6361/201321591
  * Poutanen & Stern (2010) Poutanen, J., & Stern, B. 2010, ApJ, 717, L118, doi: 10.1088/2041-8205/717/2/L118
  * PTF Team (2020) PTF Team. 2020, Palomar Transient Factory Level 1, IPAC, doi: 10.26131/IRSA155
  * Pushkarev et al. (2010) Pushkarev, A. B., Kovalev, Y. Y., & Lister, M. L. 2010, ApJ, 722, L7, doi: 10.1088/2041-8205/722/1/L7
  * Readhead (1980) Readhead, A. C. S. 1980, in Objects of High Redshift, ed. G. O. Abell & P. J. E. Peebles, Vol. 92, 165–175
  * Readhead et al. (1978) Readhead, A. C. S., Cohen, M. H., Pearson, T. J., & Wilkinson, P. N. 1978, Nature, 276, 768, doi: 10.1038/276768a0
  * Readhead et al. (2024) Readhead, A. C. S., Ravi, V., Blandford, R. D., et al. 2024, ApJ, 961, 242, doi: 10.3847/1538-4357/ad0c55
  * Rees (1966) Rees, M. J. 1966, Nature, 211, 468, doi: 10.1038/211468a0
  * Rees (1967) —. 1967, MNRAS, 135, 345, doi: 10.1093/mnras/135.4.345
  * Richards et al. (2014) Richards, J. L., Hovatta, T., Max-Moerbeck, W., et al. 2014, MNRAS, 438, 3058, doi: 10.1093/mnras/stt2412
  * Richards et al. (2011) Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29, doi: 10.1088/0067-0049/194/2/29
  * Romani (2006) Romani, R. W. 2006, AJ, 132, 1959, doi: 10.1086/508216
  * Romani et al. (2004) Romani, R. W., Sowards-Emmerd, D., Greenhill, L., & Michelson, P. 2004, ApJ, 610, L9, doi: 10.1086/423201
  * Savolainen et al. (2002) Savolainen, T., Wiik, K., Valtaoja, E., Jorstad, S. G., & Marscher, A. P. 2002, A&A, 394, 851, doi: 10.1051/0004-6361:20021236
  * Sbarrato et al. (2012) Sbarrato, T., Ghisellini, G., Nardini, M., et al. 2012, MNRAS, 426, L91, doi: 10.1111/j.1745-3933.2012.01332.x
  * Shepherd (1997) Shepherd, M. C. 1997, in Astronomical Society of the Pacific Conference Series, Vol. 125, Astronomical Data Analysis Software and Systems VI, ed. G. Hunt & H. Payne, 77
  * Sikora et al. (1994) Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153, doi: 10.1086/173633
  * Teräsranta et al. (2005) Teräsranta, H., Wiren, S., Koivisto, P., Saarinen, V., & Hovatta, T. 2005, A&A, 440, 409, doi: 10.1051/0004-6361:20053356
  * Ulrich et al. (1997) Ulrich, M.-H., Maraschi, L., & Urry, C. M. 1997, ARA&A, 35, 445, doi: 10.1146/annurev.astro.35.1.445
  * Volonteri (2010) Volonteri, M. 2010, A&A Rev., 18, 279, doi: 10.1007/s00159-010-0029-x
  * Wagner & Witzel (1995) Wagner, S. J., & Witzel, A. 1995, ARA&A, 33, 163, doi: 10.1146/annurev.aa.33.090195.001115
  * Wang et al. (2023) Wang, A., An, T., Guo, S., et al. 2023, ApJ, 944, 187, doi: 10.3847/1538-4357/acaf02
  * WISE Team (2019) WISE Team. 2019, AllWISE Source Catalog, IPAC, doi: 10.26131/IRSA1
  * WISE Team (2020) —. 2020, NEOWISE 2-Band Post-Cryo Single Exposure (L1b) Source Table, IPAC, doi: 10.26131/IRSA124
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Yi et al. (2014) Yi, W.-M., Wang, F., Wu, X.-B., et al. 2014, ApJ, 795, L29, doi: 10.1088/2041-8205/795/2/L29
  * Zhang et al. (2017) Zhang, Y., An, T., Frey, S., et al. 2017, MNRAS, 468, 69, doi: 10.1093/mnras/stx392

Figure 1: $\gamma$-ray light curves of PKS 0201+113 and its neighbors. Panels
(a): 3-month time bin light curves for the entire dataset; Panels (b): zoomed-
in 1-month time bin light curves. Blue circles and red triangles are flux
estimations and upper limits, respectively. Pink vertical bars are the
corresponding TS values. Two dashed vertical lines are used to mark the epoch
that a possible $\gamma$-ray source is seen.

Figure 2:  Upper panels: smoothed $\gamma$-ray TS maps in units of degrees for
different epochs (4$\arcdeg\times 4\arcdeg$ scale with 0.1$\arcdeg$ per pixel,
PKS 0201+113 not included in the analysis source model). Panel (b) marks the
emergence of the transient $\gamma$-ray source between MJD 55980 and 56100, in
which the green circle is the 95% C.L. error radius and the green X-shaped
symbol is the radio position of PKS 0201+113. Panel (a) corresponds to the
first 3.5-year Fermi-LAT data, while panel (c) corresponds to the last 10
years of data. Note that all the TS maps are centered at the optimized
location of the transient $\gamma$-ray source. Bottom panels: 3-month
$\gamma$-ray time bin light curve, also shown in Figure 1, OVRO light curve
spanning 15 years, ZTF $r$-band light curve and NEOWISE light curves. The
filled yellow region marks the epoch when the transient $\gamma$-ray source
appears in Figure 2. Red vertical lines mark the start and end times of the
epoch. Figure 3: A smoothed $\gamma$-ray TS map with the same scale and source
of data as Figure 2, but here PKS 0201+113 is included in the subtracted model
to check whether there is any excess in direction of the neighbor, found in
the analysis of the entire dataset. The green circle and the green X-shaped
symbol are the same as in Figure 2. The yellow circle is the 95% C.L. error
radius of the nearby source obtained from the analysis of the entire 14-yr
dataset, where no significant excess is found. Figure 4: A zoomed-in view of
the OVRO light curve, focusing on the high radio flux state. The total flux
density curve observed at 15 GHz has been decomposed using a series of four
exponential sub-flares together with a constant baseline flux density,
following the method described in Savolainen et al. (2002). The filled yellow
region is as same as that shown in Figure 2, while the vertical dashed yellow
line marks the central time of the bin with the hightest TS value in the
monthly $\gamma$-ray light curve (see Figure 1) Figure 5: Broadband SED of PKS
0201+113, along with the predictions of the two-zone leptonic model. Dotted,
dashed dotted and dashed lines correspond to the synchrotron, SSC and EC
components for the inner (colored in red) and outer emission regions (colored
in green), respectively. The dashed line with purple stars is responsible for
IC scattering of synchrotron emission of outer emission region by electrons of
inner emission region. Radiation by interactions the other way round is too
low to be seen in this diagram. The gray line represents the accretion disk
emission by a central SMBH with $\rm\sim$ 3$\times 10^{9}$ $\rm M_{\sun}$
(Marcotulli et al., 2020).  Figure 6: A schematic illustration (not to scale)
of geometry of the jet. The grey regions are responsible for generations of
the radio emission, at different evolutionary stages. The left one is at the
edge of the radio core where the SSA opacity is unity, and the radio flux
density starts to rise. Since the emergence of the $\gamma$-ray source is at
the ascent phase of a radio sub-flare peaking around MJD 56117 (see Figure 4),
the $\gamma$-ray dissipation region (a tiny peacock green blob), embedded in
an extended region (the middle grey blob), might locate downstream. The right
grey region, rather distant from the central SMBH, has been plotted for an
early radio activity.

Table 1: VLBI information of 0203+1134 Epoch | Freq | $\sigma$ | Sp | Stot | Score | Bmaj | Bmin | Bpa
---|---|---|---|---|---|---|---|---
| (GHz) | (mJy b-1) | (mJy b-1) | (mJy) | (mJy) | (mas) | (mas) | ($\arcdeg$)
2011-05-16 | 8.4 | 1.1 | 546.7$\pm$27.4 | 741.8$\pm$37.5 | 543.0$\pm$27.2 | 2.3 | 1.0 | -4.7
2013-05-07 | 7.6 | 1.0 | 649.3$\pm$32.5 | 833.1$\pm$43.8 | 537.9$\pm$27.0 | 3.0 | 1.5 | -9.2
2014-04-01 | 8.6 | 1.0 | 489.6$\pm$24.5 | 671.4$\pm$33.9 | 409.7$\pm$20.5 | 1.7 | 0.7 | -6.9
2015-06-12 | 8.6 | 0.6 | 483.6$\pm$24.2 | 671.1$\pm$34.3 | 397.7$\pm$16.9 | 2.5 | 1.0 | -5.1
Notes for columns: Col. 1 – observing date; Col. 2 – Observing frequency; Col.
3 – Peak flux density of the image; Col. 4 – Total flux density of the fitted
VLBI components; Col. 5 – The flux density of the fitted VLBI core component;
Col. 6-9 – The FWHM size of the restoring beam of the image, which is also
demonstrated at the bottom left corner of each image.
Table 2: List of parameters for the theoretical jet SEDs in Figure 5

Component | $p_{1}$ | $p_{2}$ | $\gamma_{br}$ | $\gamma_{min}$ | $\gamma_{max}$ |  $B$[Gauss] | $\delta$ |  $R_{j}^{\prime}$[cm]
---|---|---|---|---|---|---|---|---
Inner emission region | 1.6 | 4.5 | 700 | 1 | $\rm 10^{4}$ | 1 | 28.7 | $2\times 10^{16}$
Outer emission region | 1.4 | 4.0 | 250 | 1 | 500 | 0.09 | 4 | $2\times 10^{18}$

Note. — $p_{1,2}$ are the indexes of the broken power-law radiative electron
distribution; $\gamma_{br}$, $\gamma_{min}$, and $\gamma_{max}$ are the break,
minimum, and maximum energies of the electron distribution, respectively; B is
the magnetic field strength; $\delta$ is the Doppler boosting factor; and
$R^{\prime}_{j}$ is the radius of the emission emission region at the jet rest
frame.


.tifpng.pngconvert #1 .tif

# Effect of Buffer Termination on Intermixing and Conductivity in
LaTiO3/SrTiO3 Heterostructures Integrated on Si(100)

Tongjie Chen  Department of Physics, North Carolina State University, Raleigh,
North Carolina 27695, USA  Kamyar Ahmadi-Majlan  Department of
Physics,University of Texas-Arlington, Arlington, TX, 76019 USA  Zheng Hui Lim
Department of Physics,University of Texas-Arlington, Arlington, TX, 76019 USA
Zhan Zhang  Advanced Photon Source, Argonne, IL, 76019 USA  Joseph H. Ngai
Department of Physics,University of Texas-Arlington, Arlington, TX, 76019 USA
Divine .P. Kumah<EMAIL_ADDRESS>Department of Physics, North Carolina State
University, Raleigh, North Carolina 27695, USA

###### Abstract

The control of chemical exchange across heterointerfaces formed between ultra-
thin functional transition-metal oxide layers provides an effective route to
manipulate the electronic properties of these systems. We show that cationic
exchange across the interface between the Mott insulator, LaTiO3(LTO) grown
epitaxially on SrTiO3(STO)-buffered Silicon by molecular beam epitaxy depends
strongly on the surface termination of the strained STO buffer. Using a
combination of temperature-dependent transport and synchrotron X-ray crystal
truncation rods and reciprocal space mapping, an enhanced conductivity in
STO/LTO/SrO- terminated STO buffers compared to heterostructures with
TiO2-terminated STO buffers is correlated with La/Sr exchange and the
formation of metallic La1-xSrxTiO3. La/Sr exchange effectively reduces the
strain energy of the system due to the large lattice mismatch between the
nominal oxide layers and the Si substrate.

##  I Introduction

Structural, electronic and chemical interactions at the interfaces between
ultra-thin complex perovskite oxides can lead to exciting physical properties
which are not found in the bulk analogue materials including interfacial
magnetism, two-dimensional electron gases, metal-insulator transitions and
superconductivity.Hellman  _et al._ (2017); Chakhalian, Millis, and Rondinelli
(2012); Ngai, Walker, and Ahn (7 01) The ability to control these interactions
at heterointerfaces using atomic layer-by-layer synthesis techniques such as
molecular beam epitaxy allows for the unprecedented tailoring of electronic
and magnetic ground states. Structural coupling including strain, oxygen
octahedral distortions and interfacial polar distortions can be induced by
epitaxial constraints provided by appropriate substrate and buffer
layers.Pentcheva and Pickett (2009); Jia  _et al._ (2009); Zubko  _et al._ (1
03); Rondinelli, May, and Freeland (2012) Additionally, thermodynamic and
kinetic effects can lead to chemical interdiffusion across heterointerfaces
leading to significant changes to composition away from nominal values which
can significantly alter their physical properties.Nakagawa, Hwang, and Muller
(2006); Qiao  _et al._ (2013) For example, intermixing between the nominally
insulating LaTiO3 (LTO) and SrTiO3 (STO) layers can lead to the formation of a
conducting LaxSr1-xTiO3 (LSTO) interface in addition to the charge transfer
mechanism proposed to alleviate the divergent electric field which arises due
to the polar discontinuity at the LTO/STO interface.Tokura  _et al._ (1993);
Veit  _et al._ (2019) Interfacial intermixing at polar/non-polar interfaces
has been proposed as a contributing mechanism to the high-mobility two
dimensional electron gas formed at the interface between LaAlO3 and SrTiO3
Willmott  _et al._ (2007); Qiao  _et al._ (2010); Lee  _et al._ (2016); Ohtomo
and Hwang (2004); Nakagawa, Hwang, and Muller (2006); Chambers (2011);
Pentcheva and Pickett (2010) and interfacial conductivity in LaCrO3/SrTiO3
superlattices.Comes  _et al._ (2017)

Of particular interest is understanding how the interfacial structure and
electronic behavior of epitaxial oxides evolve when integrated on Si. Oxides
exhibit a variety of properties that can potentially be exploited in device
applications, provided integration onto Si is achieved.Kumah, Ngai, and
Kornblum (0 05); Hu  _et al._ (2003); Woicik  _et al._ (2006); Mi  _et al._
(2008); Chambers  _et al._ (2001); Kornblum  _et al._ (5 18) The surface unit-
cell of Si(100) has a lattice constant of 3.84 Å, which in many cases will
impart compressive strain on an epitaxial oxide. Intermixing across a
heterointerface may be energetically favorable if it leads to a reduction of
the strain energy of a system.Sankara Rama Krishnan  _et al._ (2014); Lipinski
_et al._ (2000) Developing ways to either promote or minimize intermixing,
depending on its desirability, is crucial for controlling the functionality of
epitaxial oxides on Si.

In this letter, we explore the effects of deposition sequence and chemical
termination of a buffer layer on intermixing at polar/non-polar interfaces
within STO/LTO/STO trilayers grown on Si(100). We find that deposition
sequence and chemical termination have pronounced effects on intermixing and
electronic behavior in heterostructures that are nominally identical in terms
of composition. To determine the degree of chemical intermixing at the LTO/STO
interface and the effect on the transport properties of the heterostructures,
we investigate the atomic-scale structures of 1.5 uc STO cap/ 3 uc LTO/ N uc
STO buffer/(001) Si samples grown by molecular beam epitaxy where the STO
buffer thickness N is either 6 unit cells or 4 unit cells. (1 unit cell (uc)
=1 monolayer (ML) La(Sr)O + 1 ML TiO2) We find that in addition to strain-
driven chemical intermixing, La-Sr exchange is more favorable if the LTO is
grown by co-depositing LaO and TiO2 on SrO-terminated STO buffers leading to
higher conductivity compared to samples where the STO is TiO2-terminated. The
structural profiles are obtained by synchrotron-diffraction based crystal
truncation rod measurements and high-resolution reciprocal space maps.

Bulk LTO is a Mott insulator with a pseudo-cubic lattice constant of 3.96
Å.Cwik  _et al._ (2003) The Mott insulating state can be broken resulting in
metallicity by epitaxial strain, over-oxidation or doping with divalent
Sr.Zhang  _et al._ (2019); He  _et al._ (2012); Wong  _et al._ (2010); Ishida
and Liebsch (2008); Hays  _et al._ (1999) Bulk STO is a band insulator with a
pseudocubic lattice constant of 3.905 Å. A high mobility two-dimensional
electron gas (2DEG) forms at the LTO/STO interface exhibiting
superconductivity.Seo  _et al._ (2007); Biscaras  _et al._ (2010) Interest in
integrating the unique properties of the oxide 2DEG system with semiconductor-
based technologies has led to the fabrication of LTO/STO heterostructures on
Si and Ge.Jin  _et al._ (2014) STO grows epitaxially on (001)-oriented Si with
the epitaxial relationship given by STO[110]/Si[100] and STO[001]//Si[001].
The c-axis of STO and Si are aligned in the out-of-plane direction with the
perovskite STO lattice rotated in-plane by 45∘ with respect to the Si
surface.McKee, Walker, and Chisholm (1998); Reiner  _et al._ (2010), The in-
plane lattice mismatch between STO and Si is given by
$\epsilon=\frac{a_{\text{STO}}-a_{\text{Si surface}}}{a_{\text{Si}}}$ where
the in-plane lattice spacing of STO, $a_{\text{STO}}$=3.905 Å and the in-plane
lattice spacing of the Si (001) surface $a_{\text{Si surface}}$=3.84 Å. Given
these values, the calculated STO/Si lattice mismatch is 1.66 %. Due to the
large lattice mismatch and the step-structure of the Si surface, strain
relaxation of the STO lattice is known to occur within a critical thickness,
$t_{\text{critical}}$= 5-10 unit cells.Kumah  _et al._ (2010); Warusawithana
_et al._ (2009); Chen  _et al._ (2018); Woicik  _et al._ (2006); Segal  _et
al._ (2010).

By controlling the thickness of the STO buffer, $t_{STO_{b}}$, the strain
state of the LTO layers can be effectively tuned. For LTO films grown on thin
STO buffer layers with $t_{STO}$ less than$t_{critical}$, the lattice mismatch
between the LTO films and the STO buffer coherently strained to Si is 2.5 %
and thus, it becomes energetically favorable to reduce the strain energy by
relaxing the strain through the formation of dislocations, and/or by chemical
intermixing with the STO layers to form LaxSr1-xTiO3 which has a smaller bulk
lattice constant than LTO.Hays  _et al._ (1999)

##  II Experiment

The nominal 1.5 uc STO cap/ 3 uc LTO/ N uc STO buffer/(001) Si
heterostructures were grown by molecular beam epitaxy. 2” diameter Si
(100)-oriented wafers (Virginia Semiconductor) were loaded into the MBE
chamber and cleaned by exposing to activated oxygen generated by a radio
frequency source to remove residual organics from the surface. To desorb the
native oxide layer which formed at the surface of the Si substrate, 2
monolayers of Sr were deposited at a substrate temperature of 550 ∘C, and the
sample was heated to 870 ∘C to remove the native layer of SiOx through the
formation and desorption of SrO.

Following the appearance of a 2 $\times$ 1 reconstruction in the reflection
high energy electron diffraction (RHEED) pattern, indicative of a clean
reconstructed Si surface, half a monolayer of Sr was deposited at 660 ∘C to
form a template for subsequent layers of STO. The substrate was then cooled to
room temperature, and 3 ML of SrO and 2 ML of TiO2 were co-deposited at room
temperature and then heated to 500 ∘C to form 2.5 ucs (two complete STO ucs
terminated by an SrO layer) of crystalline STO, as shown in Figure 1.
Subsequent layers of STO and LTO of various thicknesses were grown at a
substrate temperature of 500 ∘C.

Figure 1: Schematics of 1.5 uc STO cap/ 3 uc LTO/ 4 uc STO/ Si samples grown
by molecular beam epitaxy. Though identical in composition, the sequence of
deposition and the terminating layer of the STO buffer (highlighted in red)
prior to LTO deposition differ.

Two sets of samples were then grown for each STO buffer thickness as shown in
the schematic in Figure 1 for the N=4 heterostructure. While identical in
composition, the two sets differ in the sequence of deposition, and thereby
the terminating layer of the STO buffer prior to LTO deposition. For the so-
called TiO2-terminated buffer, a single uc of STO was deposited on the 2.5 uc
base layer of crystallized STO through co-deposition of Sr and Ti, followed by
1 monolayer of TiO2 to form the buffer, as indicated in Figure 1. Three ucs of
LTO followed by 1 uc of STO were then deposited by co-deposition of La/Sr and
Ti. A single monolayer of SrO was then deposited to cap the heterostructure.
In contrast, for the so-called SrO-terminated buffer, 1 uc of STO, 3 ucs of
LTO and 2 ucs of STO were deposited sequentially all through co-deposition of
Sr/La and Ti fluxes on top of the 2.5 uc crystallized base, as indicated in
Figure 1. For both methods, the STO buffer was briefly annealed in vacuum at
580 ∘C immediately prior to deposition of LTO, to enhance crystallinity. For
the N = 6 heterostructures, samples with TiO2-terminated and SrO-terminated
buffers were deposited in the manner outlined above, except that the buffer
has 2 additional ucs of STO in comparison to the N = 4 samples.

The transport properties of the trilayers were measured in the Van-der-Pauw
configuration with Au contacts sputtered on the sample corners. The
temperature-dependent sheet resistance was measured using a Keithley 2400
sourcemeter and a Keithley 2700 multiplexer in a Quantum Design Physical
Property Measurements System (PPMS).

To determine the relationship between the structural and transport properties
of the films, the atomic-scale structure of the samples were determined by
high-resolution synchrotron X-ray diffraction crystal truncation rod (CTR)
measurements.Robinson (1986); Kumah  _et al._ (2010) X-ray diffraction
measurements were performed at the 33ID beamline at the Advanced Photon
Source. The samples were mounted in a Be-dome chamber evacuated to a base
pressure of 5$\times 10^{-5}$ Torr. The incident photon energy was fixed at 16
keV ($\lambda=0.7749\AA{}$). The diffracted intensities were measured using a
Pilatus 100K 2D detector.Schlepütz  _et al._ (2005)

##  III Results and Discussion

The temperature-dependent transport properties of the samples as a function of
the STO buffer thickness (4 uc and 6 uc) and termination (SrO or TiO2) are
given in Figure 2. The N=4 samples show metallic behavior for both buffer
terminations, however, the N=4 sample with an SrO-terminated buffer where
La/Sr intermixing is enhanced as discussed below, exhibits a lower sheet
resistance compared with the sample with a TiO2-terminated buffer. The sheet
resistance of the N=6 sample with the TiO2-terminated buffer increases as the
temperature decreases.Ahmadi-Majlan  _et al._ (2018) For the corresponding N=6
sample grown on an SrO-terminated buffer, the resistivity at room-temperature
is reduced by 1 order of magnitude at room temperature and decreases as the
temperature decreases indicative of a metallic ground state. The reduced
resistivity for the SrO-terminated buffer is postulated to be due to enhanced
La-Sr intermixing leading to the formation of metallic LaxSr1-xTiO3.Tokura
_et al._ (1993)

Figure 2: Comparison of the sheet resistance as a function of temperature for
nominal 1.5 uc STO cap/ 3 uc LTO/ N=4,6 uc STO buffer/ (001) Si substrate
heterostructures with SrO and TiO2 terminated STO buffers.

Crystal truncation rods along the Si substrate-defined reciprocal lattice
vectors (1 Si reciprocal lattice unit (r.l.u.) = 1/5.43 Å${}^{-1})$were
measured to determine atomic structure of the coherently strained fractions of
the oxide heterostructures. The diffraction data of the coherently strained
fraction of the film was analyzed using the genetic-algorithm based GenX X-ray
fitting program.Björck and Andersson (2007) In addition to the CTRs, relaxed
film peaks were observed at non-integer in-plane Si reciprocal lattice
vectors. The in-plane lattice vectors of the incoherent fraction of the films
do not coincide with the lattice vectors of the Si substrate, thus,
reciprocal-space mapping measurements were performed to determine the lattice
parameters of the strain-relaxed portions of the films.

Figure 3 shows a comparison of the measured crystal truncation rods along the
Si H=1 K=1 direction for the N=4 heterostructures with SrO and TiO2 terminated
buffer layers. Due to the rotation of the perovskite unit cell by 45o relative
to the Si lattice, the perovskite (20Lfilm) peaks are present along the Si
(11L) crystal truncation rod.

A significant difference observed in the measured data for the 2 buffer
terminations is the position of the film Bragg peak. The film Bragg peaks of
the SrO-terminated sample are shifted to higher L values as compared to the
TiO2-terminated sample. This indicates a smaller (larger) average out-of-plane
lattice spacing for the SrO (TiO2) terminated buffer samples. The calculated
average film lattice parameters for the SrO and TiO2 terminated buffer samples
are 3.95 Å and 4.02 Å, respectively.

Figure 3: Comparison of the measured and fitted diffraction intensities along
the Si 11L crystal truncation rods for nominal 1.5 uc STO cap/ 3 uc LTO/ 4 uc
STO buffer/ (001) Si substrate heterostructures with SrO and TiO2 terminated
STO buffers. The dashed vertical lines indicate the locations of the Bragg
peaks for the TiO2-terminated buffer sample.

To determine the layer-resolved structural profile of the fractions of the
sample coherently strained (i.e. in-plane lattice constant is the same as the
substrate) to the Si substrate, the measured data were fit using the GenX
X-ray fitting program.Björck and Andersson (2007) The fit parameters are the
lattice parameters of each layer and the La and Sr occupations of the A-site
of the perovskite unit cell to account for La-Sr intermixing across the
STO/LTO interfaces. The simulated CTR’s for the best fit structures are shown
as solid lines in Figure 3.

Figure 4 shows a comparison of the layer-resolved out-of-plane lattice
spacings and La/Sr chemical profiles for the 2 N=4 samples obtained from the
CTR analysis. For both samples, the lattice spacings of the STO buffer layers
(layers 1-4) adjacent to the Si substrate are measured to be 3.96 Å $\pm 0.01$
Å corresponding to a c/a ratio of 1.03 (a=3.84 Å). The layer spacings in the
nominal LTO layers (layers 5-7) for the TiO2-terminated (SrO-terminated)
buffer sample have an average value of 4.05 $\pm 0.02$ (3.97 $\pm 0.02$) Å
corresponding to a c/a of 1.06 (1.034).

The composition profiles along the growth direction for the 2 samples are
shown in the lower panel of Figure 4. While the La fractional occupation of
the LTO layers in the TiO2-terminated buffer sample are close to the expected
value of 1, a significant reduction in the La content and a corresponding
increase in the Sr content is observed within the LTO layers for the SrO-
terminated buffer sample. The Sr-incorporation into the nominal LTO layer
leads to the formation of metallic LaxSr1-xTiO3 where the lattice volume
decreases with the Sr concentration Hays  _et al._ (1999) effectively reducing
the lattice mismatch with the buffer layer. Thus, the reduced lattice constant
observed for the SrO-terminated buffer sample is consistent with the measured
chemical profiles.

Figure 4: Comparison of the layer-resolved (a)lattice spacings and (b) La
occupation profile for nominal 1.5 uc STO cap/ 3 uc LTO/ 4 uc STO buffer/
(001) Si substrate heterostructures with SrO and TiO2 terminated STO buffers.
The purple and green dashed vertical lines indicate the locations of the
nominal STO buffer/ LTO and LTO/STO cap interfaces, respectively.  Figure 5:
Reciprocal space maps around the (1 1 2) film peak for nominal 1.5 uc STO cap/
3 uc LTO/ 4 uc STO buffer/ (001) Si substrate heterostructures with (a) TiO2
and (b) SrO terminated STO buffers. (c) Line profiles in the in-plane
H-direction through the maps in (a) and (b) for L=2.72 Si r.l.u. (TiO2
termination) and L=2.76 Si r.l.u. (SrO termination). (d) Line profiles for
fixed H in the L-direction through the strained and relaxed film peaks.

The reciprocal space maps (RSM) around the H=2 K=0 L=2.7 Si r.l.u. are
compared for the TiO2 and SrO buffer terminated N=4 samples in Figure 5(a) and
(b). Due to the epitaxial relationship between the Si and the perovskite film
unit cell, the (2 0 2.7) Si peak corresponds to the (1 1 2) perovskite film
peak. For each sample, 2 peaks are observed. The narrow peak along H=2 Si
r.l.u. corresponds to the fraction of the film coherently strained to the Si
substrate with the in-plane lattice parameter values of 3.84 Å. The broader
peak at lower H values corresponds to relaxed fractions of the film with in-
plane lattice parameters larger than 3.84 Å. The observation of relaxed film
peaks indicates that lateral strain relaxation occurs independent of the STO
buffer termination. The average lattice parameters of the relaxed and strained
fractions of the film can be determined from the peak positions in Figure 5(a)
and (b). Table 1 summarizes the average lattice parameters calculated from the
RSMs for the two samples.

Figure 5(c) shows a cut in the H-direction along K=0 and L=2.72 Si r.l.u.
(TiO2 termination) and L=2.76 Si r.l.u. (SrO termination). The peaks for the
relaxed fractions are located at H=1.973 Si r.l.u. and H=1.962 Si r.l.u. for
the TiO2 and SrO terminated buffers, respectively. The line profiles along the
L direction for the strained and relaxed fractions are shown in Figure 5(d).
Here, the peak for the TiO2-terminated sample is shifted to a higher L value
than the SrO sample indicating a smaller out-of-plane spacing for the TiO2
sample.

Parameter |  SrO-term |  TiO2-term
---|---|---
Hstrained(Si r.l.u.) |  2.00 |  2.00
a${}_{strained}(\AA{})$ |  3.84 |  3.84
Hrelaxed (Si r.l.u.) |  1.962 |  1.973
a${}_{relaxed}(\AA{})$ |  3.915 |  3.893
Lstrained (Si r.l.u.) |  2.7308 |  2.7085
c${}_{strained}(\AA{})$ |  3.977 |  4.010
Lrelaxed (Si r.l.u.) |  2.7564 |  2.7181
c${}_{relaxed}(\AA{})$ |  3.9406 |  3.9962
V${}_{strained}(\AA{}^{3})$ |  58.643 |  59.13
V${}_{relaxed}(\AA{}^{3})$ |  60.398 |  60.56
Table 1: Comparison of lattice parameters for the relaxed and strained Bragg
peaks for the TiO2 terminated buffer and the SrO-terminated buffer sample.

When LaO and TiO2 are co-deposited on a SrO-terminated STO buffer, La species
can directly react with SrO layer to form an alloy which reduces the lattice
volume and minimizes the strain energy of the system. The suppressed La/Sr
exchange observed for LTO deposited on TiO2 terminated STO buffer layers
suggest that the TiO2 layer serves as an effective barrier layer for La-Sr
interdiffusion between LTO and STO. For the samples with the TiO2 terminated
buffers, as the STO buffer thickness is increased from N=4 to N=6, the
epitaxial strain in the STO relaxes leading to a reduced lattice mismatch
between the STO and LTO adlayer. Thus, strain-driven intermixing is expected
to be further suppressed in the N=6 sample leading to more chemically abrupt
LTO/STO interfaces and insulating behavior observed for the N=6 sample in
Figure 2.Ahmadi-Majlan  _et al._ (2018)

##  IV Conclusion

In conclusion, we have demonstrated how strain-driven chemical intermixing at
the LTO/STO interface is strongly dependent on the chemical termination of the
STO buffer. Enhanced La/Sr intermixing for LTO co-deposited on SrO-terminated
STO strained to Si leads to enhanced conductivity due to the formation of
metallic LSTO. We find that the La/Sr exchange is suppressed for
TiO2-terminated STO buffers leading to a significant increase in the sheet
resistance for samples with nominally identical composition. These results
highlight the critical importance of chemical interactions driven by epitaxial
strain and the composition of the interface terminal layer on the physical
properties of functional oxide materials. We demonstrate that deposition
sequence and terminating layers can be exploited to promote or minimize cation
intermixing in layered heterostructures integrated on Si(100).

##  V Acknowledgements

This work was supported by the National Science Foundation (NSF) under Awards
No. DMR-1751455. Synthesis and transport measurements were supported by the
National Science Foundation (NSF, Award No. DMR-1508530). The use of the
Advanced Photon Source was supported by the U. S. Department of Energy, Office
of Science, Office of Basic Energy Sciences, under Contract No. DE-
AC02-06CH11357.

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Hellman  _et al._ (2017) F. Hellman, A. Hoffmann, Y. Tserkovnyak, G. S. Beach, E. E. Fullerton, C. Leighton, A. H. MacDonald, D. C. Ralph, D. A. Arena, H. A. Dürr, _et al._ , Reviews of modern physics 89, 025006 (2017).
  * Chakhalian, Millis, and Rondinelli (2012) J. Chakhalian, A. Millis, and J. Rondinelli, Nature Materials 11, 92 (2012).
  * Ngai, Walker, and Ahn (7 01) J. Ngai, F. Walker, and C. Ahn, Annual Review of Materials Research 44, 1 (2014-07-01), publisher: Annual Reviews.
  * Pentcheva and Pickett (2009) R. Pentcheva and W. E. Pickett, Physical Review Letters 102, 107602 (2009).
  * Jia  _et al._ (2009) C. Jia, S. Mi, M. Faley, U. Poppe, J. Schubert, and K. Urban, Physical Review B 79, 081405 (2009).
  * Zubko  _et al._ (1 03) P. Zubko, S. Gariglio, M. Gabay, P. Ghosez, and J.-M. Triscone, Annual Review of Condensed Matter Physics 2, 141 (2011-03).
  * Rondinelli, May, and Freeland (2012) J. M. Rondinelli, S. J. May, and J. W. Freeland, MRS bulletin 37, 261 (2012).
  * Nakagawa, Hwang, and Muller (2006) N. Nakagawa, H. Y. Hwang, and D. A. Muller, Nature Materials 5, 204 (2006).
  * Qiao  _et al._ (2013) L. Qiao, K. Zhang, M. E. Bowden, T. Varga, V. Shutthanandan, R. Colby, Y. Du, B. Kabius, P. V. Sushko, M. D. Biegalski, _et al._ , Advanced Functional Materials 23, 2953 (2013).
  * Tokura  _et al._ (1993) Y. Tokura, Y. Taguchi, Y. Okada, Y. Fujishima, T. Arima, K. Kumagai, and Y. Iye, Physical Review Letters 70, 2126 (1993).
  * Veit  _et al._ (2019) M. Veit, M. Chan, B. Ramshaw, R. Arras, R. Pentcheva, and Y. Suzuki, Physical Review B 99, 115126 (2019).
  * Willmott  _et al._ (2007) P. R. Willmott, S. A. Pauli, R. Herger, C. M. Schlep tz, D. Martoccia, B. D. Patterson, B. Delley, R. Clarke, D. Kumah, C. Cionca, and Y. Yacoby, Physical Review Letters 99, 155502 (2007).
  * Qiao  _et al._ (2010) L. Qiao, T. C. Droubay, V. Shutthanandan, Z. Zhu, P. Sushko, and S. A. Chambers, Journal of Physics: Condensed Matter 22, 312201 (2010).
  * Lee  _et al._ (2016) P. W. Lee, V. N. Singh, G. Y. Guo, H.-J. Liu, J.-C. Lin, Y.-H. Chu, C. H. Chen, and M.-W. Chu, Nature Communications 7, 12773 (2016).
  * Ohtomo and Hwang (2004) A. Ohtomo and H. Y. Hwang, Nature 427, 423 (2004).
  * Chambers (2011) S. A. Chambers, Surface Science 605, 1133 (2011).
  * Pentcheva and Pickett (2010) R. Pentcheva and W. E. Pickett, Journal of Physics: Condensed Matter 22, 043001 (2010).
  * Comes  _et al._ (2017) R. B. Comes, S. R. Spurgeon, D. M. Kepaptsoglou, M. H. Engelhard, D. E. Perea, T. C. Kaspar, Q. M. Ramasse, P. V. Sushko, and S. A. Chambers, Chemistry of Materials 29, 1147 (2017).
  * Kumah, Ngai, and Kornblum (0 05) D. P. Kumah, J. H. Ngai, and L. Kornblum, Advanced Functional Materials 30, 1901597 (2020-05).
  * Hu  _et al._ (2003) X. Hu, H. Li, Y. Liang, Y. Wei, Z. Yu, D. Marshall, J. Edwards Jr, R. Droopad, X. Zhang, A. Demkov, _et al._ , Applied physics letters 82, 203 (2003).
  * Woicik  _et al._ (2006) J. C. Woicik, H. Li, P. Zschack, E. Karapetrova, P. Ryan, C. R. Ashman, and C. S. Hellberg, Physical Review B 73, 024112 (2006).
  * Mi  _et al._ (2008) S.-B. Mi, C.-L. Jia, V. Vaithyanathan, L. Houben, J. Schubert, D. G. Schlom, and K. Urban, Applied physics letters 93, 101913 (2008).
  * Chambers  _et al._ (2001) S. Chambers, Y. Liang, Z. Yu, R. Droopad, and J. Ramdani, Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films 19, 934 (2001).
  * Kornblum  _et al._ (5 18) L. Kornblum, E. N. Jin, D. P. Kumah, A. T. Ernst, C. C. Broadbridge, C. H. Ahn, and F. J. Walker, Applied Physics Letters 106, 201602 (2015-05-18), publisher: American Institute of Physics.
  * Sankara Rama Krishnan  _et al._ (2014) P. Sankara Rama Krishnan, A. N. Morozovska, E. A. Eliseev, Q. M. Ramasse, D. Kepaptsoglou, W.-I. Liang, Y.-H. Chu, P. Munroe, and V. Nagarajan, Journal of Applied Physics 115, 054103 (2014).
  * Lipinski  _et al._ (2000) M. Lipinski, H. Schuler, O. Schmidt, K. Eberl, and N. Jin-Phillipp, Applied Physics Letters 77, 1789 (2000).
  * Cwik  _et al._ (2003) M. Cwik, T. Lorenz, J. Baier, R. M ller, G. Andr , F. Bour e, F. Lichtenberg, A. Freimuth, R. Schmitz, E. M ller-Hartmann, and M. Braden, Physical Review B 68, 060401 (2003).
  * Zhang  _et al._ (2019) T. Zhang, C. Gu, Z. Mao, X. Chen, Z. Gu, P. Wang, Y. Nie, and X. Pan, Applied Physics Letters 115, 261604 (2019).
  * He  _et al._ (2012) C. He, T. Sanders, M. Gray, F. Wong, V. Mehta, and Y. Suzuki, Physical Review B 86, 081401 (2012).
  * Wong  _et al._ (2010) F. J. Wong, S.-H. Baek, R. V. Chopdekar, V. V. Mehta, H.-W. Jang, C.-B. Eom, and Y. Suzuki, Physical Review B 81, 161101 (2010).
  * Ishida and Liebsch (2008) H. Ishida and A. Liebsch, Physical Review B 77, 115350 (2008).
  * Hays  _et al._ (1999) C. C. Hays, J.-S. Zhou, J. T. Markert, and J. B. Goodenough, Physical Review B 60, 10367 (1999).
  * Seo  _et al._ (2007) S. S. A. Seo, W. S. Choi, H. N. Lee, L. Yu, K. W. Kim, C. Bernhard, and T. W. Noh, Physical Review Letters 99, 266801 (2007).
  * Biscaras  _et al._ (2010) J. Biscaras, N. Bergeal, A. Kushwaha, T. Wolf, A. Rastogi, R. Budhani, and J. Lesueur, Nature Communications 1, 89 (2010).
  * Jin  _et al._ (2014) E. N. Jin, L. Kornblum, D. P. Kumah, K. Zou, C. C. Broadbridge, J. H. Ngai, C. H. Ahn, and F. J. Walker, APL Materials 2, 116109 (2014).
  * McKee, Walker, and Chisholm (1998) R. A. McKee, F. J. Walker, and M. F. Chisholm, Physical Review Letters 81, 3014 (1998).
  * Reiner  _et al._ (2010) J. W. Reiner, A. M. Kolpak, Y. Segal, K. F. Garrity, S. Ismail-Beigi, C. H. Ahn, and F. J. Walker, Advanced Materials 22, 2919 (2010).
  * Kumah  _et al._ (2010) D. P. Kumah, J. W. Reiner, Y. Segal, A. M. Kolpak, Z. Zhang, D. Su, Y. Zhu, M. S. Sawicki, C. C. Broadbridge, C. H. Ahn, and F. J. Walker, Applied Physics Letters 97, 251902 (2010).
  * Warusawithana  _et al._ (2009) M. P. Warusawithana, C. Cen, C. R. Sleasman, J. C. Woicik, Y. Li, L. F. Kourkoutis, J. A. Klug, H. Li, P. Ryan, L.-P. Wang, M. Bedzyk, D. A. Muller, L.-Q. Chen, J. Levy, and D. G. Schlom, Science 324, 367 (2009).
  * Chen  _et al._ (2018) T. Chen, K. Ahmadi-Majlan, Z. H. Lim, Z. Zhang, J. H. Ngai, A. F. Kemper, and D. P. Kumah, Applied Physics Letters 113, 201601 (2018).
  * Segal  _et al._ (2010) Y. Segal, J. W. Reiner, Z. Zhang, C. H. Ahn, and F. J. Walker, Journal of Vacuum Science & Technology B, Nanotechnology and Microelectronics: Materials, Processing, Measurement, and Phenomena 28, C5B1 (2010).
  * Robinson (1986) I. K. Robinson, Physical Review B 33, 3830 (1986).
  * Schlepütz  _et al._ (2005) C. Schlepütz, R. Herger, P. Willmott, B. Patterson, O. Bunk, C. Brönnimann, B. Henrich, G. Hülsen, and E. Eikenberry, Acta Crystallographica Section A: Foundations of Crystallography 61, 418 (2005).
  * Ahmadi-Majlan  _et al._ (2018) K. Ahmadi-Majlan, T. Chen, Z. H. Lim, P. Conlin, R. Hensley, M. Chrysler, D. Su, H. Chen, D. P. Kumah, and J. H. Ngai, Applied Physics Letters 112, 193104 (2018).
  * Björck and Andersson (2007) M. Björck and G. Andersson, Journal of Applied Crystallography 40, 1174 (2007).