# Single laser pulse driven thermal limit of the quasi-two dimensional
magnetic ordering in Sr2IrO4

Ruitang Wang1,2,3 J. Sun1 D. Meyers4,5 J. Q. Lin1,2,3 J. Yang6 G. Li1 H.
Ding2,3 Anthony D. DiChiara7 Y. Cao8 J. Liu6 M. P. M. Dean4 Haidan Wen7,8 X.
Liu1<EMAIL_ADDRESS>1School of Physical Science and Technology,
ShanghaiTech University, Shanghai 201210, China.  2 Beijing National
Laboratory for Condensed Matter Physics and Institute of Physics, Chinese
Academy of Sciences, Beijing 100190, China  3 University of Chinese Academy of
Sciences, Beijing 100049, China  4 Condensed Matter Physics and Materials
Science Department, Brookhaven National Laboratory, Upton, New York 11973,
USA.  5 Department of Physics, Oklahoma State University, Stillwater, Oklahoma
74078, USA.  6 Department of Physics and Astronomy, University of Tennessee,
Knoxville, Tennessee 37996, USA.  7 Advanced Photon Source, Argonne National
Laboratory, Argonne, IL, 60439, USA.  8 Materials Science Division, Argonne
National Laboratory, Argonne, Illinois, 60439, USA.

###### Abstract

Upon femtosecond-laser stimulation, generally materials are expected to
recover back to their thermal-equilibrium conditions, with only a few
exceptions reported. Here we demonstrate that deviation from the thermal-
equilibrium pathway can be induced in canonical 3D antiferromagnetically (AFM)
ordered Sr2IrO4 by a single 100-fs-laser pulse, appearing as losing long-range
magnetic correlation along one direction into a glassy condition. We further
discover a ‘critical-threshold ordering’ behavior for fluence above
approximately 12 mJ/cm2 which we show corresponds to the smallest
thermodynamically stable $c$-axis correlation length needed to maintain long-
range quasi-two-dimensional AFM order. We suggest that this behavior arises
from the crystalline anisotropy of the magnetic-exchange parameters in
Sr2IrO4, whose strengths are associated with distinctly different timescales.
As a result, they play out very differently in the ultrafast recovery
processes, compared with the thermal equilibrium evolution. Thus, our
observations are expected to be relevant to a wide range of problems in the
nonequilibrium behavior of low-dimensional magnets and other related ordering
phenomena.

###### pacs:

Valid PACS appear here

##  I Introduction

Understanding the mechanisms behind laser manipulation of magnetism in
materials is indispensable to many problems in both fundamental and applied
science Beaurepaire1996 ; RMP2010 ; KimelNRM2019 ; Bigot2009 ; Nemec2018 .
Phenomenologically, the classic three-temperature modelBeaurepaire1996  works
quite well in explaining the experimental observations on the evolution of the
spins upon external femtosecond-laser stimuli in many systemsRMP2010 ;
Beaurepaire1996 ; QZhang2006 ; Kimling2014 ; Koopmans2010 . In this model, the
spin, electron and lattice baths are coupled through a set of mutual
interaction constants which govern the energy-flow rates when the system is
stimulated into nonthermal equilibrium conditions. Eventually, the spin sector
thermalizes back to equilibrium condition after the externally deposited
energy equilibrates among the charge, spin and lattice
reservoirsBeaurepaire1996 . Here we report an exceptional observation of
single-laser-pulse-induced anisotropic partial recovery of the magnetic
ordering in Sr2IrO4, leading to permanent reduction of the interplane magnetic
correlation. We propose that this unusual recovery behavior could be explained
by the noncooperation of different terms of the exchange interactions
dictating the spin ordering. When these exchanges differ significantly in
strength, they are unfolded along the time axisPeter2019 ; LCLSreview  into
different time scales in the ultrafast recovery process. Because of this
‘time-window mismatch’, the weaker exchange terms could be quenched, leading
to a deviation from the quasiequilibrium relaxation pathway.

Sr2IrO4 is a layered Mott insulator with each Ir site hosting a 1/2
pseudospinBKim2008 . The highly anisotropic crystalline structure leads to
strong magnetic anisotropy(Fig.1(a)). The system starts to develop three-
dimensional anteferromagnetic(AFM) ordering below $T_{N}\sim
240~{}K$Crawford1994 , from the cooperation of a strong intraplane exchange
$J$ of approximately 60 meV and a weak inter-plane exchange $J_{c}$ of
approximately $16~{}\mu$eV JKim2012 ; Porras2019 . The magnetic ordering in
Sr2IrO4 has been well observed with x-ray resonant magnetic scattering (XRMS)
measurements, appearing as Bragg-peaks in the magnetic scattering channel when
the incident x-ray energy is tuned to $\sim$11.216 keV at the Ir-L3 resonance
edgeBKim2009 . This high x-ray resonant energy allows the access to a large
reciprocal space, thus the flexibility in selecting low x-ray incident angle
scattering geometry to best meet the x-ray and laser penetration match
requirementGuzelturk2021 ; Mark2016 . Our experiment was set up such that our
$\sim$100 nm epitaxial thin film sample was stimulated with 1 eV laser
($\sim$100 fs) pulse by pulse, and the AFM ordering peak was continuously
monitored with XRMS at an x-ray incident angle of less than $5^{\circ}$ to
track the recovery of the magnetic ordering at a readout frequency of 1 Hz
(see Ref. suppl ).

Figure 1: Single-laser-pulse-induced suppression of the (1 0 16) magnetic
reflection. (a) The crystal and magnetic structure of Sr2IrO4. J and Jc denote
the in-plane and interplane exchange couplings, respectively. (b) CCD image of
the magnetic peak. (c) Cartoon demonstration of the projection of the CCD
plane in the reciprocal space. (d) Response of the magnetic Bragg-peak height
to a single laser pulse (arrives at t=0) as a function of laser fluence. The
black dots are the average of the last 100 s of each curve. (e) Multiple-shot
evolution of the magnetic-peak height with moderate laser fluence. Multiple
shots arrive in the shaded region and each creates a small intensity drop
until the decrease is fully compensated by an initial recovery.

##  II Experimental observations

Figure 1(b) shows the (1 0 16) 3D AFM ordering peak in Sr2IrO4 observed on the
pixelated area detector at 80 KLupascu2014 , and Fig.1(c) depicts the detector
surface projection in the reciprocal space. In order to catch the true first
laser-shot response from pristine thermal-equilibrium condition, each dataset
is collected after the sample goes through a fresh thermal cycle by warming
above $T_{N}$ and then cooling back again to 80 K with the laser turned off.
Then x-ray is turned on to continuously monitor the magnetic Bragg-peak
intensity. In between, laser stimuli are controlled to come in pulse by pulse.
The arrival time of the first laser pulse is defined as time zero. The
evolution of the experimentally observed magnetic-peak height as a function of
the laser fluence is shown in Fig.1(d).

After a single-laser-pulse stimulus to the magnetic ordering prepared from
thermal-equilibrium evolution, the system does not recover to the initial
state, evidenced by the permanent suppression of the magnetic Bragg-peak
height as shown in Fig.1(d). The degree of the suppression depends on the
laser fluence, i.e., the degree of the damage to the initial ordered spin
network. More details of the response to the pulsed laser stimulation are
shown in Fig.1(e) with a laser fluence of 8.4 mJ/cm2. At this fluence, the
first laser pulse leads to an approximately $33\%$ drop in the magnetic Bragg-
peak height, and the system settles at this condition without further
evolution within an hour of continuous x-ray measurement. A second pulse
arriving after one hour causes further suppression to the measured peak
height, but with a much smaller drop. After 5$\sim$6 pulses, the system
recovery becomes repeatable and more stimulation does not cause further
permanent peak-height reduction.

Figure 2: Evolution of the magnetic ordering upon laser stimulation. (a)
Comparison of the $H$ (in-plane) and $L$ (interplane) scans before and after
single-pulse-laser pumping. The horizontal bars are our instrumental
resolution. Scans are normalized to the peak height determined from $L$ scans
to emphasize the variation in the peak widths. (b) Fluence-dependence
evolution of the $H$\- and $L$\- scans.

These experimental observations clearly demonstrate that the recovery of the
magnetic ordering in Sr2IrO4 depends on the history of laser stimulation, and
the response can be classified into two distinct stages. In the first stage of
the initial few pulses, the degree of magnetic ordering keeps to be
permanently suppressed and the system does not recover to the condition before
pulse arrival. In the second stage, the system does recover but only to a
state prepared by multiple initial laser pulses, which, of course is different
from the thermal-equilibrium evolution. The laser-pulse-induced suppression of
the magnetic ordering can be erased with a full thermal cycle and our
measurements are fully reproducible (see Ref. suppl ), indicating the observed
response is not due to sample damage but intrinsic to the spin system in
Sr2IrO4. Phenomenologically, the observed shot-by-shot dependence is similar
to the photoinduced metastable insulator-to-metal phase transition in
La2/3Ca1/3MnO3JZhang2016 . The laser stimulation drives the system into a
nonthermal-equilibrium condition, and the deviation depends on the history of
laser stimulation, a non-Markovian-type behavior. Such deviation is obviously
beyond the classic three-temperature modelBeaurepaire1996 ; QZhang2006 ;
Koopmans2010 ; Kimling2014 .

To explore the microscopic origin of such evolution, detailed reciprocal space
scans along the in-plane (H scans) and interplane (L scans) directions are
performed across the magnetic peak before and after the first laser-pulse
arrival under various laser fluences, and the results are shown in Fig.2. Note
that the peaks in the $H$ scan, a measure of the in-plane spin ordering, are
very sharp and identical before and after the laser pulses for all laser
fluences applied. From the width of the peaks in $H$ scans, the in-plane spin-
ordering correlation length is estimated to be of a macroscopic scale about
0.9 $\mu m$. Thus, within a layer, the magnetic ordering is always fully
restored to the thermal-equilibrium condition. The observed permanent
suppression in the first stage is due to the incomplete recovery of the spin
correlation along the interplane direction, which is evidenced by the
broadening of the peaks in $L$ scans. With the spin-correlation length
reduced, the scattering intensity is transferred from the center to the tails,
leading to a permanent suppression of the peak height as observed. These
results suggest that the laser pulse leads to highly anisotropic response
between the in-plane and interplane AFM ordering correlations in the recovery
process.

Figure 3: Evolution of the interplane correlation as function of fluence.
(a,b) The fitting results for two representative conditions with weak and
strong laser fluence. (c-e) The evolution of the interplane correlation length
$\xi_{z}$, peak height and the integrated intensity as function of fluence,
which are all normalized to the values before laser stimulation. With fluence
higher than $\sim$12 mJ/cm2, $\xi_{z}$ saturates to $\sim$5 layer thickness.
The exponential lines are guides to the eye. The error bar for fluence is
estimated to be $10\%$. (see Supplementary Materialsuppl )

To quantify such a response, the $L$ scans of the magnetic Bragg peak measured
at thermal-equilibrium condition and after the very first laser pulse stimulus
with varying fluence were analyzed to extract the evolution of the interlayer
spin-ordering texture. It turns out that, for all experimental conditions, the
$L$ scans can be well described by a single Lorentzian function as
$I(q_{z},\xi_{z})\sim\frac{2}{d}\frac{\xi_{z}}{1+q_{z}^{2}\xi_{z}^{2}}$
(Fig.3(a,b), and Ref. suppl ), with $d$ to be the interlayer distance and
$q_{z}$ the interlayer direction relative-momentum transfer. The evolution of
the extracted correlation length $\xi_{z}$ is shown in Fig.3(c). As a function
of the laser fluence, $\xi_{z}$, as well as the peak height (Fig.3(d)),
generally follows an exponential decay. On the other hand, the total area
under the magnetic peak is conserved (Fig.3(e)). The significant data
fluctuation mostly comes from slight misalignment in the H direction where the
peak width corresponds to 0.01∘ in instrument rotational angle. Such
conservation indicates that the size and direction of the local ordered
magnetic moments are the same as those of the pristine conditionCollins1989 ;
Shannon1992 . Thus the incomplete recovery is solely due to the reduction in
the interplane correlation. It is interesting to notice that the exponential
decay is offset from zero. Both the correlation length and the peak height
show a tendency to saturation at higher laser fluence, and the interplane spin
correlation cannot be completely destroyed.

##  III Discussion

Our observations reveal a few interesting aspects of the recovery of the
magnetic-ordering texture in this highly magnetically anisotropic layered
Sr2IrO4. First, starting from the thermal-equilibrium condition, each of the
first few femtosecond-laser pulses induces permanent suppression to the
magnetic ordering, and the degree of the suppression depends on the states
prepared by the preceding pulses. Second, the system does eventually stabilize
into conditions that the recovery of the magnetism becomes repeatable upon
consecutive laser stimulation. This justifies the validity of experiments of
stroboscopic mode in probing the spin dynamicsRMP2010 ; Beaurepaire1996 ;
Mark2016 . Third, although the degree of the partial recovery keeps degrading
with increased laser fluence, it saturates to a robust nonzero value,
indicating a minimum interplane partial recovery is intrinsically protected.
The observed history and laser fluence- dependence hint at the nature of the
magnetic recovery in Sr2IrO4 upon ultrafast laser stimulation.

The history dependence evidences that the demagnetization has certain local
character. In the ultrafast process, the demagnetization does not completely
wipe out all traces of magnetic order, and in the recovery the spins do not
have long enough time to reach a global bath temperature. As a result, there
are unperturbed (or less-perturbed) local spin clusters which preserve the
memory of the prior state. Thus the prevailing global energy-dissipation
picture in the classic three-temperature modelBeaurepaire1996 ; QZhang2006 ;
Koopmans2010 ; Kimling2014  is oversimplified.

Notably, in both the initial and the second stages, the in-plane spin
correlation always recovers to the thermal-equilibrium evolution condition of
micron size. This differentiates our observation from the conventional fast
quenching of a high-temperature state where the correlations along all
directions are expected to changeBray1994 ; Mazenko1990 ; Shah1990 . The full
in-plane recovery and the permanent loss of the interplane spin correlation
are intimately related to the individual terms in the exchange interaction
governing the magnetic dynamics along different directions. Although the 3D
magnetic ordering in Sr2IrO4 is jointly determined by both the intraplane and
interplane magnetic couplings Kastner1998 ; Porras2019 ; Fujiyama2012 , in the
ultrafast recovery process they act differently. The strong in-plane exchange
of tens of meV drives a quick re-establishment of the in-plane correlation
within a few picosecondsMark2016 . We suggest that in a such short duration,
the additional energy introduced by laser pumping into the spin system cannot
be efficiently absorbed by the lattice reservoirJLi2019 ; Yuelin2016 .
Instead, the spin sector is still highly excited in picosecond time scale.
Associated with the weak interlayer exchange $J_{c}$ of approximately
$16~{}\mu$eV, a time window much longer than picoseconds is needed to allow
them to fully dissipate. Such a process is cut off when the macroscopic
intralayer spin correlation is established, since there is an enormous energy
barrier to flip a whole layer.

Figure 4: Theoretical modeling. Calculated magnetic-order parameter at the
middle of the slab, $\langle S_{z}\rangle_{middle}$, as a function of the
layer number from the model Hamiltonian
$H=\frac{J}{2}\sum_{l,ij}\vec{S}_{li}\cdot\vec{S}_{lj}+\Delta\sum_{l,ij}S^{z}_{li}S^{z}_{lj}+\frac{J_{c}}{2}\sum_{ll^{{}^{\prime}},ij}\vec{S}_{li}\cdot\vec{S}_{l^{{}^{\prime}}j}$
at T = 80 K with exchange terms $J=60~{}meV$, $Jc=16.4~{}\mu eV$. Different
curves correspond to different anisotropic exchange $\Delta$ strengths, and
the consequent gap sizes are used as an index.

Following the above arguments, we can understand the saturation plateau at
high laser fluence as shown in Fig.3(c,d). In our earlier reportMark2016 , we
have demonstrated that the laser pulse completely destroys the magnetic order
in the initial hundreds of femtoseconds, when the laser fluence is stronger
than approximately 12 mJ/cm2. Thus, beyond such laser fluence, the system
completely loses its memory of the history and the recovery follows the same
path, regardless of further laser-fluence increasing. This is also consistent
with the observed multishot evolution where, at high fluence, the laser pulses
after the very first shot drive marginal further suppression to the magnetic
Bragg-peak height (see Fig. S9 in Supplemental Material suppl ). Without the
assistance from the remnant order, reestablishment of the global interplane
correlation completely lags behind the intraplane recovery. Thus the system
recovery enters a quasi-two-dimensional regime. For such condition, Mermin and
WagnerMermin1966  prove that spontaneous two-dimensional ferromagnetic or
antiferromagnetic long-range order at finite temperature is highly susceptible
to spin thermal fluctuations, and the nonzero third-dimension correlation is
critical to suppress the thermal fluctuation to realize long-range two-
dimensional ordering.

We confirm such an explanation for our observed plateau by calculating the
magnetic correlation function in a few layers of square lattice from an
effective anisotropic Heisenberg spin-$1/2$ model, whose dynamics is solved
with the equation-of-motion technique and mean-field approximationDiep1991
(see Ref. suppl ). With the realistic parameters from experiment and published
literatureJKim2012 ; Porras2019 , the self-consistent ordered magnetic moment
$\langle S_{z}\rangle$ and the magnetic correlation function $\langle
S^{-}S^{+}\rangle$ are obtained as functions of the model slab thickness. As
the reported exchange anisotropy for Sr2IrO4 is quite controversialPorras2019
; YGim2016 ; Calder2018 ; Pincini2017 , the calculated $\langle
S_{z}\rangle_{middle}$, the ordered magnetic moment at the middle of the slab,
is shown in Fig.4 as a function of the calculated magnon gap size. As
expected, both the anisotropic exchange $\Delta$ strengths and the layer
thickness are critical for true long-range magnetic ordering. When the gap
size approaches approximately 1 meV, the magnetic-order parameter $\langle
S_{z}\rangle_{middle}$ stabilizes around 4-5 layers, which agrees with the
interplane correlation length of the saturation plateau we observe in
experiment (Fig.3(c,d)). Thus, we realize a spin thermal fluctuation limit in
real material with laser-pulse stimulation. Furthermore, this result indicates
that indeed the interplane recovery-time window is set by the in-plane spin
correlation. The longer interplane correlation established at lower fluence is
assisted by memory of the spin network, which is only partially destroyed
below the high-fluence threshold.

##  IV Conclusion

In conclusion, a history and laser-fluence dependence of the partial-recovery
process of the 3D AFM ordering in Sr2IrO4 was observed, which is related to
the distinctly different timescales for the interplane and intraplane
recoveries in the nonthermal-equilibrium ultrafast process. The noncooperation
of the different exchange interactions in the ultrafast process drives the
system to deviate from the quasiequilibrium relaxation pathway. Light-induced
deviation from thermal-equilibrium evolution has been rarely observed, and
previous reports are associated with the latticeStoica2019  or
chargeKiryukhin1997 ; JZhang2016 ; Stojchevska2014 ; SuYang2020  degrees of
freedom. Our results extend the direct observation of deviation from recovery
to thermal-equilibrium condition into the spin sector. Furthermore, we suggest
that such time-window mismatch could generally happen in complex systems
during ultrafast nonthermal evolution, and our observations could be relevant
to a wide range of problems in the nonequilibrium-behavior of low-dimensional
magnets and related ordering phenomena. For example, the laser-induced ‘hidden
quantum state’ in layered 1T-TaS2 could be one special caseStojchevska2014 ;
SLee2019 .

##  V acknowledgments

We thank Yi Zhu for the assistance of experimental setup. The experimental
work by X. L. and R. W. was primarily supported by National Natural Science
Foundation of China under grant No. 11934017. Part of the execution of the
experiment and data interpretation by H. W. and Y. C. were supported by the
U.S. DOE, Office of Science, Office of Basic Energy Sciences, Materials
Sciences and Engineering Division. 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. Work at Brookhaven
National Laboratory was supported by the U.S. DOE, Office of Science, Office
of Basic Energy Sciences, Materials Sciences and Engineering Division under
Contract No. DE-SC0012704. J. L. acknowledges support from the National
Science Foundation under Grant No. DMR-1848269. H. D. acknowledges support
from the National Natural Science Foundation of China (No. 11888101). H. D.
and X. L. acknowledge support from the Ministry of Science and Technology of
China (2016YFA0401000). J. Y. acknowledges funding from the State of Tennessee
and Tennessee Higher Education Commission (THEC) through their support of the
Center for Materials Processing.

## References

  * (1) A. Kirilyuk, A. V. Kimel, and T. Rasing, Ultrafast Optical Manipulation of Magnetic Order, Rev. Mod. Phys. 82, 2731 (2010).
  * (2) A.V. Kimel and M. Li, Writing Magnetic Memory with Ultrashort Light Pulses, Nat. Rev. Mater. 4, 189-200 (2019).
  * (3) J. Bigot, M. Vomir and E. Beaurepaire, Coherent ultrafast magnetism induced by femtosecond laser pulses, Nat. Phys. 5, 515 (2009).
  * (4) P. Némec, M. Fiebig, T. Kampfrath and A. V. Kimel, Antiferromagnetic opto-spintronics, Nat. Phys. 14, 229 (2018).
  * (5) E. Beaurepaire, J.-C. Merle, A. Daunois, and J.-Y. Bigot, Ultrafast Spin Dynamics in Ferromagnetic Nickel, Phys. Rev. Lett. 76, 4250 (1996).
  * (6) Q. Zhang, A. V. Nurmikko, G. X. Miao, G. Xiao and A. Gupta, Ultrafast Spin-Dynamics in Half-Metallic CrO2 Thin Flms. Phys. Rev. B 74, 064414 (2006).
  * (7) B. Koopmans, G. Malinowski, F. Dalla Longa, D. Steiauf, M. Fähnle, T. Roth, M. Cinchetti, and M. Aeschlimann, Explaining the Paradoxical Diversity of Ultrafast Laser-induced Demagnetization, Nat. Mater. 9, 259 (2010).
  * (8) J. Kimling, J. Kimling, R. B. Wilson, B. Hebler, M. Albrecht and D. G. Cahill, Ultrafast Demagnetization of FePt:Cu Thin Films and the Role of Magnetic Heat Capacity, Phys. Rev. B 90, 224408 (2014).
  * (9) M. Mitrano, et al., Ultrafast Time-resolved X-ray Scattering Reveals Diffusive Charge Order Dynamics in La2-xBaxCuO4, Sci. Adv. 5, eaax3346 (2019).
  * (10) C. Bostedt, S. Boutet, D. M. Fritz, Z. Huang, H. J. Lee, H. T. Lemke, A. Robert, W. F. Schlotter, J. J. Turner, and G. J. Williams, Linac Coherent Light Source: The First Five Years, Rev. Mod. Phys. 88, 015007 (2016).
  * (11) B. J. Kim, H. Jin, S. J. Moon, J.-Y. Kim, B.-G. Park, C. S. Leem, J. Yu, T. W. Noh, C. Kim, S.-J. Oh J.-H. Park, V. Durairaj, G. Cao, and E. Rotenberg, Novel Jeff=1/2 Mott State Induced by Relativistic Spin-Orbit Coupling in Sr2IrO4, Phys. Rev. Lett. 101, 076402 (2008).
  * (12) M. K. Crawford, M. A. Subramanian, R. L. Harlow, J. A. Fernandez-Baca, Z. R. Wang, and D. C. Johnston, Structural and Magnetic Studies of Sr2IrO4, Phys. Rev. B 49, 9198 (1994).
  * (13) J. Kim, D. Casa, M. H. Upton, T. Gog, Y.-J. Kim, J. F. Mitchell, M. van Veenendaal, M. Daghofer, J. van den Brink, G. Khaliullin, and B. J. Kim, Magnetic Excitation Spectra of Sr2IrO4 Probed by Resonant Inelastic X-Ray Scattering: Establishing Links to Cuprate Superconductors, Phys. Rev. Lett. 108, 177003 (2012).
  * (14) J. Porras, J. Bertinshaw, H. Liu, G. Khaliullin, N. H. Sung, J.-W. Kim, S. Francoual, P. Steffens, G. Deng, M. Moretti Sala, A. Efimenko, A. Said, D. Casa, X. Huang, T. Gog, J. Kim, B. Keimer, and B. J. Kim, Pseudospin-lattice Coupling in the Spin-Orbit Mott Insulator Sr2IrO4, Phys. Rev. B 99, 085125 (2019).
  * (15) B. J. Kim, H. Ohsumi, T. Komesu, S. Sakai, T. Morita, H. Takagi, and T.-H. Arima, Phase-Sensitive Observation of a Spin-Orbital Mott State in Sr2IrO4, Science 323, 1329 (2009).
  * (16) B. Guzelturk, T. Winkler, T. W. J. Van de Goor, M. D. Smith, S. A. Bourelle, S. Feldmann, M. Trigo, S. W. Teitelbaum, H. G. Steinrück, G. A. de la Pena, R. Alonso-Mori, D. Zhu, T. Sato, H. I. Karunadasa, M. F. Toney, F. Deschler, and A. M. Lindenberg, Visualization of Dynamic Polaronic Strain Fields in Hybrid Lead Halide Perovskites, Nat. Mater. 20, 618 (2021).
  * (17) M. Dean, Y. Cao, X. Liu, S. Wall, D. Zhu, R. Mankowsky, V. Thampy, X. M. Chen, J. G. Vale, D. Casa et al., Ultrafast Energy- and Momentum-Resolved Dynamics of Magnetic Correlations in the Photo-Doped Mott Insulator Sr2IrO4, Nat. Mater. 15, 601 (2016).
  * (18) (see Supplemental Material.)
  * (19) A. Lupascu, J.P. Clancy, H. Gretarsson, Z. Nie, J. Nichols, J. Terzic, G. Cao, S. S. A. Seo, Z. Islam, M. H. Upton, Jungho Kim, D. Casa, T. Gog, A. H. Said, V.M. Katukuri, H. Stoll, L. Hozoi, J. van den Brink, and Y.-J. Kim, Tuning Magnetic Coupling in Sr2IrO4 Thin Films with Epitaxial Strain, Phys. Rev. Lett. 112, 147201 (2014).
  * (20) J. Zhang, X. Tan, M. Liu, S. W. Teitelbaum, K. W. Post, F. Jin, K. A. Nelson, D. Basov, W. Wu, and R. D. Averitt, Cooperative Photoinduced Metastable Phase Control in Strained Manganite Films, Nat. Mater. 15, 956 (2016).
  * (21) M. F. Collins, Magnetic Critical Scattering (Oxford University Press, 1989).
  * (22) R. F. Shannon, S. E. Nagler, C. R. Harkless, and R. M. Nicklow, X-Ray Intensity Fluctuation Spectroscopy Studies on Phase-Ordering Systems, Phys. Rev. B 46, 40 (1992).
  * (23) A. J. Bray, Theory of Phase-Ordering Kinetics, Adv. Phys. 43, 357 (1994).
  * (24) G. F. Mazenko, Theory of Unstable Growth, Phys. Rev. B 42, 4487 (1990).
  * (25) P. J. Shah and O. G. Mouritsen, Dynamics of Ordering Processes in Annealed Dilute Systems: Island Formation, Vacancies at Domain Boundaries, and Compactification, Phys. Rev. B 41, 7003 (1990).
  * (26) M. A. Kastner, R. J. Birgeneau, G. Shirane, and Y. Endoh, Magnetic, Transport, and Optical Properties of Monolayer Copper Oxides, Rev. Mod. Phys. 70, 897 (1998).
  * (27) S. Fujiyama, H. Ohsumi, T. Komesu, J. Matsuno, B. J. Kim, M. Takata, T. Arima, and H. Takagi, Two-Dimensional Heisenberg Behavior of Jeff=1/2 Isospins in the Paramagnetic State of the Spin-Orbital Mott Insulator Sr2IrO4, Phys. Rev. Lett. 108, 247212 (2012).
  * (28) J. Li, R. Wang, H. Guo, Y. Zhu, Y. Cao, J. Liu, H. Ding, H. Wen and X. Liu, Recovery of Photoexcited Magnetic Ordering in Sr2IrO4. J. Phys.: Condens. Matter 31 25580 (2019).
  * (29) Y. Li, R. D. Schaller, M. Zhu, D. A. Walko, J. Kim, X. Ke, L. Miao and Z. Q. Mao, Strong Lattice Correlation of Non-Equilibrium Quasiparticles in a Pseudospin-1/2 Mott Insulator Sr2IrO4, Sci. Rep. 6, 19302 (2016).
  * (30) N. Mermin, and H. Wagner, Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models, Phys. Rev. Lett. 17, 1133 (1966).
  * (31) H. Diep, Quantum Effects in Heisenberg Antiferromagnetic thin Films, Phys. Rev. B 43, 8509 (1991)
  * (32) Y. Gim, A. Sethi, Q. Zhao, J. F. Mitchell, G. Cao, and S. L. Cooper, Isotropic and Anisotropic Regimes of the Field-Dependent Spin Dynamics in Sr2IrO4: Raman Scattering Studies, Phys. Rev. B 93, 024405 (2016).
  * (33) S. Calder, D. M. Pajerowski, M. B. Stone, and A. F. May, Spin-gap and Two-Dimensional Magnetic Excitations in Sr2IrO4, Phys. Rev. B 98, 220402(R) (2018).
  * (34) D. Pincini, J. G. Vale, C. Donnerer, A. de la Torre, E. C. Hunter, R. Perry, M. Moretti Sala, F. Baumberger, and D. F. McMorrow, Anisotropic Exchange and Spin-Wave Damping in Pure and Electron-Doped Sr2IrO4, Phys. Rev. B 96, 075162 (2017).
  * (35) V. A. Stoica, N. Laanait, C. Dai et al., Optical Creation of a Supercrystal with Three-Dimensional Nanoscale Periodicity, Nat. Mater. 18, 377 (2019).
  * (36) V. Kiryukhin, D. Casa, J. Hill, B. Keimer, A. Vigliante, Y. Tomioka and Y. Tokura, An X-Ray-Induced Insulator-Metal Transition in a Magnetoresistive Manganite, Nature 386, 813 (1997).
  * (37) S.-Y. Xu, Q. Ma, Y. Gao, A. Kogar, A. Zong, A. M. M. Valdivia, T. H. Dinh, S.-M. Huang, B. Singh, C.-H. Hsu et al., Spontaneous Gyrotropic Electronic Order in a Transition-Metal Dichalcogenide, Nature (London) 578, 545 (2020).
  * (38) L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S. Brazovskii and D. Mihailovic, Ultrafast Switching to a Stable Hidden Quantum State in an Electronic Crystal, Science 344, 177 (2014).
  * (39) S.-H. Lee, J. S. Goh and D. Cho, Origin of the Insulating Phase and First-order Metal-Insulator Transition in 1T-TaS2, Phys. Rev. Lett. 122, 106404 (2019).


# Grid-VLP: Revisiting Grid Features for Vision-Language Pre-training

Ming Yan, Haiyang Xu, Chenliang Li, Bin Bi, Junfeng Tian, Min Gui and Wei Wang
Alibaba Group
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Existing approaches to vision-language pre-training (VLP) heavily rely on an
object detector based on bounding boxes (regions), where salient objects are
first detected from images and then a Transformer-based model is used for
cross-modal fusion. Despite their superior performance, these approaches are
bounded by the capability of the object detector in terms of both
effectiveness and efficiency. Besides, the presence of object detection
imposes unnecessary constraints on model designs and makes it difficult to
support end-to-end training.

In this paper, we revisit grid-based convolutional features for vision-
language pre-training, skipping the expensive region-related steps. We propose
a simple yet effective grid-based VLP method that works surprisingly well with
the grid features. By pre-training only with in-domain datasets, the proposed
Grid-VLP method can outperform most competitive region-based VLP methods on
three examined vision-language understanding tasks. We hope that our findings
help to further advance the state of the art of vision-language pre-training,
and provide a new direction towards effective and efficient VLP.

##  1 Introduction

Recent studies (Lu et al., 2019; Su et al., 2019; Tan and Bansal, 2019; Chen
et al., 2019; Li et al., 2020c; Yu et al., 2020) on vision-language pre-
training have substantially advanced the state of the art across a variety of
vision-and-language (V+L) tasks. These approaches typically follow a two-stage
pipeline: 1) a pre-trained object detector is first used to encode an image by
identifying a set of salient visual objects from the image; and 2) a cross-
modal fusion model is pre-trained to learn the cross-modal representations.
While most recent VLP methods utilize region-based visual features (Anderson
et al., 2018) extracted by object detection, this paper explores a new way for
vision-language pre-training, by using vanilla grid-based convolutional
features from ConvNets (He et al., 2016; Xie et al., 2017).

Despite their superior performance, region-based VLP solutions suffer from the
three problems: 1) The task-specific object detector has a strong impact on
the performance of the VLP models. Important visual information may be lost
during object detection; 2) In object detection, the extraction and post-
processing of region features is very time-consuming; 3) The feature
extraction process is non-differentiable, which imposes unnecessary
constraints on model designs and makes it difficult to support end-to-end
training.

To address the limitations, we revisit the grid-based convolutional features
for vision-language pre-training. Specifically, we encode an image into a
feature map with convolutional networks such as ResNet (He et al., 2016), and
then conduct cross-modal fusion directly between text and the derived image
feature map. In this way, we step out of the bounding box design of local
region features and make the full use of all visual information for vision and
language learning. On the other hand, computing grid features enables more
flexible and simpler model design and thus leads to inference speedup, by
skipping expensive region-related steps.

In this paper, we propose a simple yet effective grid-based VLP method, namely
Grid-VLP, by directly pre-training with the grid-based convolutional features
instead of the dominant region-based features from bottom-up attention
(Anderson et al., 2018). We follow LXMERT (Tan and Bansal, 2019) and adopt the
typical Masked Language Modeling, Image-Text Matching and Image Question
Answering as our pre-training tasks. To handle the challenge of modeling long
sequences of grid features and accelerate the training process, we utilize a
random grid sampling mechanism as in PixelBERT (Huang et al., 2020), which
helps improve the robustness of visual feature learning. We pre-train Grid-VLP
only with in-domain datasets (Visual Genome (Krishna et al., 2017) and MS-COCO
(Lin et al., 2014)), and then fine-tune it on three popular V+L understanding
tasks: VQA (Antol et al., 2015), NLVR2 (Suhr et al., 2018) and GQA (Hudson and
Manning, 2019). Our results show that with the grid features, Grid-VLP
outperforms the competitive region-based VLP methods across all tasks.
Besides, we also provide in-depth analysis on the influence of different
design choices regarding feature types, image encoder architectures and
resolution of input images.

Figure 1: The overall framework of Grid-VLP.

##  2 Related Work

Existing approaches to VLP (Li et al., 2020a; Su et al., 2019; Chen et al.,
2019; Li et al., 2020c; Tan and Bansal, 2019; Yu et al., 2020; Huang et al.,
2020) mainly take a two-step training pipeline, which consists of extracting
semantic visual features by specific object detector and training the cross-
modal pre-training model to align text and visual features. There are
typically two lines of research about this topic. The first line uses a
single-stream transformer architecture to model both image and text
representations in a unified semantic space such as VLBERT (Su et al., 2019),
UNITER (Chen et al., 2019) and OSCAR (Li et al., 2020c). In contrast, the
other line adopts a two-stream Transformer architecture that first encodes the
image and text modalities separately, and then fuses the cross-modal
representations with another Transformer network, such as LXMERT (Tan and
Bansal, 2019) and ERNIE-ViL (Yu et al., 2020). Recently, VinVL (Zhang et al.,
2021) pre-trained a large-scale object-attribute detection model with much
larger amounts of data on four public object detection datasets for further
improving the performance. In addition to image-text pairs, UNIMO (Li et al.,
2020b) also employed large scale of free text corpus and image collections for
enhancing the cross-modal learning. These methods rely heavily on a task-
specific bounding box (or region) based object detector, which impose
unnecessary constraints on model designs and limit potential applications of
existing vision and language systems. In this paper, Grid-VLP explores a new
way for vision-language pre-training by utilizing the grid-based convolutional
features, which skips the expensive region-related steps and is not restricted
to the necessity of object detector.

##  3 Grid-VLP Approach

###  3.1 Model Architecture

The architecture of Grid-VLP is shown in Figure 1. Given a pair of aligned
image and caption text, we first utilize a pre-trained CNN encoder to extract
the grid features for the image, then a Transformer-based model is used to
conduct cross-modal fusion directly between the token embeddings and image
feature map. Different V+L pre-training tasks are designed to further enhance
the cross-modal learning. To make up for the difficulty of modeling long grid
feature sequence and accelerate training, we adopt a simple random grid
sampling mechanism to select part of image feature map during cross-modal
fusion. The architecture is flexible and no expensive region-related steps are
involved.

###  3.2 Input Representations

The input to Grid-VLP is an image and its related text (e.g. caption text). We
first introduce the way to represent the text sequence and image.

#### Sentence Embeddings

Similar to BERT (Devlin et al., 2018), each sentence is first split into a
sequence of sub-words $\\{w_{1},...,w_{m}\\}$ by WordPiece tokenizer. Then
each token $w_{i}$ is assigned three kinds of learnable embeddings: token,
modal type and position embeddings. The three embeddings are summed and layer-
normalized to represent input sentence as a sequence of embedding vectors
$E_{emb}=\\{e_{CLS},e_{1},...,e_{m},e_{SEP}\\}$, where $[CLS]$ and $[SEP]$ are
special tokens in BERT.

#### Image Grid Features

We use the traditional grid-based convolutional features instead of region-
based features for representing an image. Starting from the initial image
$v_{img}\in R^{3\times H_{0}\times W_{0}}$ with 3 color channels, a CNN-based
image encoder generates a lower-resolution activation map $f_{img}\in
R^{C\times H\times W}$ using the typical values as in Faster R-CNN (Ren et
al., 2015): $C=2048$ and $H=\frac{H_{0}}{32},W=\frac{W_{0}}{32}$. To make the
CNN encoder aware of fine-grained semantics, we add a $1\times 1$ RoIPool
layer on top of the CNN encoder followed by two fully-connected layers and
pre-train the encoder on Visual Genome dataset as in (Jiang et al., 2020).
After pre-training, the CNN encoder is fixed as a grid feature extractor. As
the cross-modal fusion network expects a sequence as input, we collapse the
spatial dimensions of $f_{img}$ into one dimension, resulting in a $HW\times
C$ feature map $s_{img}$. Then we take a linear projection layer to reduce the
channel dimension of the high-level feature map $s_{img}$ from $C$ to a
smaller dimension $d$ so as to match the dimension of token embeddings. To
distinguish between different modalities, we supplement the grid feature map
with a learnable modal type embedding that are added to the output of linear
projection layer. Finally, the sequential image representation
$I_{img}=\\{g_{1},...,g_{HW}\\}$ can be seen as a $HW$ length of
$d$-dimensional vector.

###  3.3 Cross-modal Fusion

Given the embeddings of the tokens for the sentence $\\{e_{i}\\}_{i=1}^{m}$
and the sequential image representations $\\{g_{j}\\}_{j=1}^{n}$, we adopt the
Transformer encoder to learn cross-modal fusion between image grids and
language tokens. To allow a fine-grained feature-level fusion, we concatenate
the derived image features and text embeddings to construct the input
sequence, which is formulated as:
$\\{e_{CLS},e_{1},...,e_{m},e_{SEP},g_{1},...,g_{HW}\\}$. To facilitate cross-
modal understanding, we follow LXMERT (Tan and Bansal, 2019) and conduct three
popular pre-training tasks111Different from LXMERT, we do not have any vision
pre-training task without the concept of region., including Masked Language
Modeling (MLM), Image-Text Matching (ITM) and Image Question Answering (QA).

#### Random Grid Sampling

During pre-training, to further increase the difficulty of cross-modal pre-
training tasks, we adopt a random sampling strategy as in (Huang et al., 2020)
by randomly sampling a fixed number of grids for each image. At each iteration
step, given the sequence of extracted grid features, we will dynamically
sample a part from them and feed it into Transformer. In this way, we
encourage the model to learn cross-modal relation with incomplete visual
input, so as to enhance the robustness. Besides, it can largely accelerate the
training process by reducing the total input sequence length. During fine-
tuning, we still use the completed grid features to keep all the extracted
visual information.

##  4 Experiments

Models | Params | VQA | NLVR2 | GQA
---|---|---|---|---
Test-dev | Test-std | Dev | Test-P | Test-dev | Test-std
Single-stream | VLBERT | 110M | 71.16 | - | - | - | - | -
UNITER | 110M | 72.70 | 72.91 | 77.14 | 77.87 | - | -
OSCAR | 110M | 73.16 | 73.61 | 78.07 | 78.36 | 61.19 | 61.23
UNIMO | 110M | 73.79 | 74.02 | - | - | - | -
VinVL | 110M | 75.95 | 76.12 | 82.05 | 83.08 | 65.05 | 64.65
Two-stream | ViLBERT | 221M | 70.55 | 70.92 | 67.40 | 67.00 | - | -
12-in-1 | 221M | 73.15 | - | - | - | - | 60.65
LXMERT | 183M | 72.42 | 72.54 | 74.90 | 74.50 | 60.00 | 60.33
ERNIE-ViL |  210M | 72.62 | 72.85 | - | - | - | -
End2End | PixelBERT | 170M | 74.45 | 74.55 | 76.50 | 77.20 | - | -
Our Model | Grid-VLP | 110M | 76.05 | 76.18 | 80.03 | 81.11 | 63.50 | 63.81
Table 1: Evaluation Results on VQA, NLVR2 and GQA.

###  4.1 Pre-training Dataset

The same in-domain data as in LXMERT (Tan and Bansal, 2019) is used for pre-
training. It consists of the image caption data from MS-COCO (Lin et al.,
2014), Visual Genome (Krishna et al., 2017), and image question answering data
from VQA v2.0 (Antol et al., 2015), GQA balanced version (Hudson and Manning,
2019) and VG-QA (Zhu et al., 2016). The total amount of the dataset is 9.18M
image-and-sentence pairs on 180K distinct images.

###  4.2 Implementation Details

The maximum sequence length for the sentence is set as 20. The visual encoder
is selected as ResNeXt with different sizes (Xie et al., 2017) as in (Jiang et
al., 2020; Huang et al., 2020). We resize the shorter side of the input image
to 600, and limit the longer side to at most 1000. We select a fixed number of
100 random grids each time during pre-training 222We also test with different
number of selected grids such as 64 and 128, it does not make much
difference.. We pre-train Grid-VLP with 12 layers of Transformer encoder. The
basic settings of the Transformer is the same as BERT (Devlin et al., 2018).
We pre-train Grid-VLP with a total batch size of 512 for 30 epochs on 8 V100
GPUs. We use the AdamW optimizor and set the initial learning rate as
$10^{-4}$.

###  4.3 Main Results on Downstream Tasks

We compare Grid-VLP model against all the state-of-the-art VLP models of the
comparable model size on three popular V+L understanding tasks: VQA (Antol et
al., 2015), NLVR2 (Suhr et al., 2018) and GQA (Hudson and Manning, 2019). The
detailed task descriptions can be referred in LXMERT (Tan and Bansal, 2019).
The results are shown in Table 1. We can see that with only in-domain pre-
training data on MS-COCO and Visual Genome, Grid-VLP can outperform almost all
the other existing VLP methods 333Only one exception is VinVL, which uses
large amounts of objection detection data, out-of-domain pre-training data and
fine-tuning tricks such as using unbalanced data on GQA. by a large margin,
without using any complicated strategies such as data augmentation,
adversarial training and knowledge enrichment. Grid-VLP even outperforms an
end-to-end method PixelBERT, which directly optimizes on pixel level. It shows
the effectiveness of the proposed approach for conducting VLP on a grid
feature level, which provides a brand new way for vision-language pre-
training.

###  4.4 Influence of Image Encoder

Backbone | Avg Time (ms) | VQA | NLVR2
---|---|---|---
x50 | 106 | 74.34 | 78.01
x101 | 153 | 75.18 | 78.95
x152 | 221 | 76.05 | 80.03
x152 (region) | 572 | 75.12 | 79.14
Table 2: Results of different image encoder architectures on VQA and NLVR2
development set.

As the image encoder is critical for extracting grid features, we further
study the importance of the image encoder by changing different ResNet visual
backbone layers. To allow a fair comparison between grid features and region
features, we also add one experiment that using Faster R-CNN (Ren et al.,
2015) with the the same ResNeXt backbone to extract region features of 100
objects. The result is shown in Table 2. We can see that using more
complicated visual backbones can contribute to the final performance gain,
which proves the importance of the image encoder for cross-modal
understanding. Besides, given the same settings with only the difference of
extracted feature type, the performance of grid features can even outperform
region-based features with a much shorter feature extraction cost.

###  4.5 Impact of Input Image Size

Input Size | Speedup | VQA | NLVR2
---|---|---|---
shorter side | longer side
448 | 448 | 3x | 74.71 | 79.22
448 | 746 | 2x | 75.64 | 79.72
600 | 1000 | - | 76.05 | 80.03
Table 3: Impact of input image size on VQA and NLVR2 development set.

As the sequence length of the grid features is determined by the image size
$H$ and $W$. Therefore, the final sequence length of the input to the
transformer also largely depends on the image size, which is important to
control effectiveness and efficiency tradeoff. We further analyze the impact
of input image size to Grid-VLP as shown in Table 3. From the results, we can
see that Grid-VLP benefits from larger image sizes as input, where more fine-
grained information can be extracted. Moreover, resizing the image to a
smaller size can significantly improve the inference speed without decreasing
performance much, e.g., no more than 1% for 2 times speedup.

##  5 Conclusion

In this paper, we revisit grid features as an alternative to the dominant
bottom-up region features for vision-language pre-training. We propose a
rather simple yet effective grid-based VLP method, and find that with larger
and well pre-trained visual backbones, our Grid-VLP method can also outperform
most competitive region-based VLP methods by equipped with the grid-based
convolutional features. We hope our findings can help further advance the
progress of vision-language pre-training and potentially provide new
perspectives to vision-language pre-training.

## References

  * Anderson et al. (2018) Peter Anderson, Xiaodong He, Chris Buehler, Damien Teney, Mark Johnson, Stephen Gould, and Lei Zhang. 2018.  Bottom-up and top-down attention for image captioning and visual question answering.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 6077–6086.
  * Antol et al. (2015) Stanislaw Antol, Aishwarya Agrawal, Jiasen Lu, Margaret Mitchell, Dhruv Batra, C Lawrence Zitnick, and Devi Parikh. 2015.  Vqa: Visual question answering.  In _Proceedings of the IEEE international conference on computer vision_ , pages 2425–2433.
  * Chen et al. (2019) Yen-Chun Chen, Linjie Li, Licheng Yu, Ahmed El Kholy, Faisal Ahmed, Zhe Gan, Yu Cheng, and Jingjing Liu. 2019.  Uniter: Universal image-text representation learning.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. 2016.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778.
  * Huang et al. (2020) Zhicheng Huang, Zhaoyang Zeng, Bei Liu, Dongmei Fu, and Jianlong Fu. 2020.  Pixel-bert: Aligning image pixels with text by deep multi-modal transformers.  _arXiv preprint arXiv:2004.00849_.
  * Hudson and Manning (2019) Drew A Hudson and Christopher D Manning. 2019.  Gqa: A new dataset for real-world visual reasoning and compositional question answering.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 6700–6709.
  * Jiang et al. (2020) Huaizu Jiang, Ishan Misra, Marcus Rohrbach, Erik Learned-Miller, and Xinlei Chen. 2020.  In defense of grid features for visual question answering.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 10267–10276.
  * Krishna et al. (2017) Ranjay Krishna, Yuke Zhu, Oliver Groth, Justin Johnson, Kenji Hata, Joshua Kravitz, Stephanie Chen, Yannis Kalantidis, Li-Jia Li, David A Shamma, et al. 2017\.  Visual genome: Connecting language and vision using crowdsourced dense image annotations.  _International journal of computer vision_ , 123(1):32–73.
  * Li et al. (2020a) Gen Li, Nan Duan, Yuejian Fang, Ming Gong, and Daxin Jiang. 2020a.  Unicoder-vl: A universal encoder for vision and language by cross-modal pre-training.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 11336–11344.
  * Li et al. (2020b) Wei Li, Can Gao, Guocheng Niu, Xinyan Xiao, Hao Liu, Jiachen Liu, Hua Wu, and Haifeng Wang. 2020b.  Unimo: Towards unified-modal understanding and generation via cross-modal contrastive learning.  _arXiv preprint arXiv:2012.15409_.
  * Li et al. (2020c) Xiujun Li, Xi Yin, Chunyuan Li, Xiaowei Hu, Pengchuan Zhang, Lei Zhang, Lijuan Wang, Houdong Hu, Li Dong, Furu Wei, et al. 2020c.  Oscar: Object-semantics aligned pre-training for vision-language tasks.  _arXiv preprint arXiv:2004.06165_.
  * Lin et al. (2014) Tsung-Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick. 2014.  Microsoft coco: Common objects in context.  In _European conference on computer vision_ , pages 740–755. Springer.
  * Lu et al. (2019) Jiasen Lu, Dhruv Batra, Devi Parikh, and Stefan Lee. 2019.  Vilbert: Pretraining task-agnostic visiolinguistic representations for vision-and-language tasks.  In _Advances in Neural Information Processing Systems_ , pages 13–23.
  * Ren et al. (2015) Shaoqing Ren, Kaiming He, Ross Girshick, and Jian Sun. 2015.  Faster r-cnn: Towards real-time object detection with region proposal networks.  _arXiv preprint arXiv:1506.01497_.
  * Su et al. (2019) Weijie Su, Xizhou Zhu, Yue Cao, Bin Li, Lewei Lu, Furu Wei, and Jifeng Dai. 2019\.  Vl-bert: Pre-training of generic visual-linguistic representations.  _arXiv preprint arXiv:1908.08530_.
  * Suhr et al. (2018) Alane Suhr, Stephanie Zhou, Ally Zhang, Iris Zhang, Huajun Bai, and Yoav Artzi. 2018\.  A corpus for reasoning about natural language grounded in photographs.  _arXiv preprint arXiv:1811.00491_.
  * Tan and Bansal (2019) Hao Tan and Mohit Bansal. 2019.  Lxmert: Learning cross-modality encoder representations from transformers.  _arXiv preprint arXiv:1908.07490_.
  * Xie et al. (2017) Saining Xie, Ross Girshick, Piotr Dollár, Zhuowen Tu, and Kaiming He. 2017.  Aggregated residual transformations for deep neural networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 1492–1500.
  * Yu et al. (2020) Fei Yu, Jiji Tang, Weichong Yin, Yu Sun, Hao Tian, Hua Wu, and Haifeng Wang. 2020\.  Ernie-vil: Knowledge enhanced vision-language representations through scene graph.  _arXiv preprint arXiv:2006.16934_.
  * Zhang et al. (2021) Pengchuan Zhang, Xiujun Li, Xiaowei Hu, Jianwei Yang, Lei Zhang, Lijuan Wang, Yejin Choi, and Jianfeng Gao. 2021.  Vinvl: Making visual representations matter in vision-language models.  _arXiv preprint arXiv:2101.00529_.
  * Zhu et al. (2016) Yuke Zhu, Oliver Groth, Michael Bernstein, and Li Fei-Fei. 2016.  Visual7w: Grounded question answering in images.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 4995–5004.


# Red noise in continuous-time stochastic modelling

Andreas Morr∗,1,2, Dörte Kreher3 and Niklas Boers1,2,4

(November 30th, 2022)

###### Abstract

The concept of correlated noise is well-established in discrete-time
stochastic modelling but there is no generally agreed-upon definition of the
notion of red noise in continuous-time stochastic modelling. Here we discuss
the generalization of discrete-time correlated noise to the continuous case.
We give an overview of existing continuous-time approaches to model red noise,
which relate to their discrete-time analogue via characteristics like the
autocovariance structure or the power spectral density. The implications of
carrying certain attributes from the discrete-time to the continuous-time
setting are explored while assessing the inherent ambiguities in such a
generalization. We find that the attribute of a power spectral density
decaying as $S(\omega)\sim\omega^{-2}$ commonly ascribed to the notion of red
noise has far reaching consequences when posited in the continuous-time
stochastic differential setting. In particular, any Itô-differential
$\mathrm{d}Y_{t}=\alpha_{t}\mathrm{d}t+\beta_{t}\mathrm{d}W_{t}$ with
continuous, square-integrable integrands must have a vanishing martingale
part, i.e. $\mathrm{d}Y_{t}=\alpha_{t}\mathrm{d}t$ for almost all $t\geq 0$.
We further argue that $\alpha$ should be an Ornstein-Uhlenbeck process.

††Keywords: red noise, correlated noise, discrete- vs. continuous-time
modelling††∗ Corresponding author. Please contact andreas.morr@pik-
potsdam.de.††1 Potsdam Institute for Climate Impact Research, Potsdam,
Germany.††2 Earth System modelling, School of Engineering and Design,
Technical University Munich, Munich, Germany.††3 Institute for Mathematics,
Humboldt University of Berlin, Berlin, Germany.††4 Department of Mathematics
and Global Systems Institute, University of Exeter, Exeter, UK.††This work has
received funding from the Volkswagen Stiftung, the European Union’s Horizon
2020 research and innovation programme under grant agreement No. 820970 and
under the Marie Sklodowska-Curie grant agreement No. 956170, as well as from
the Federal Ministry of Education and Research under grant No. 01LS2001A.

##  1 Introduction

In many fields of dynamical modelling, it is common practice to introduce an
additive stochastic term to incorporate unresolved dynamics into the otherwise
deterministically defined equation of discrete-time evolution. A simple
example of this may be a stochastic process defined by

$X_{k+1}-X_{k}=f(X_{k},k)+\sigma\varepsilon_{k},\,X_{0}=x_{0}\in\mathbb{R}.$
(1)

In the classical setting, the noise term $\varepsilon_{k}$ will for each time-
step be an independent standard Gaussian random variable [1]. However, in some
cases it may be imperative to discard this assumption of independence, for
instance if the unresolved dynamics are suspected to exhibit persistence in
time [2, 3]. Perhaps the most common example of such correlated noise is so-
called red noise, generated via an AR(1)-process:

$\varepsilon_{k+1}=\varphi\varepsilon_{k}+z_{k},\,\varepsilon_{0}\sim\mathcal{N}\left(0,(1-\varphi^{2})^{-1}\right),$
(2)

where $0<\varphi<1$ and now the $(z_{k})_{k\in\mathbb{N}}$ are i.i.d. standard
Gaussian. The initial distribution of $\varepsilon_{0}$ is chosen specifically
such that the process $\varepsilon$ is stationary in the weak sense. Each
instance $\varepsilon_{k}$ of the discrete-time noise process $\varepsilon$
has exponentially decaying correlation to its neighboring instances

$\displaystyle R_{\varepsilon}(\tau)$
$\displaystyle=\frac{\varphi^{\left|\tau\right|}}{1-\varphi^{2}},\,\tau\in\mathbb{N}$
(3)

and this particular noise process is the basis for many applied discrete-time
stochastic models [4, 5, 6, 7].

When modelling natural systems, however, it is often more rigorous to model
their dynamics in continuous time. Incorporating the concept of noise into
such models requires delicate mathematical constructions. The field of
stochastic analysis offers a wide range of possibilities for introducing
stochasticity into continuous-time dynamics. Discrete-time dynamics with
additive noise, in similar spirit of (1), may be modelled by the following
stochastic differential equation (SDE):

$\mathrm{d}X_{t}=f(X_{t},t)\mathrm{d}t+\sigma\mathrm{d}Y_{t},\,X_{0}=x_{0}\in\mathbb{R}.$
(4)

In the classical setting of uncorrelated noise,
$Y=(Y_{t})_{t\in\mathbb{R}_{+}}$ takes the form of a Wiener process
$W=(W_{t})_{t\in\mathbb{R}_{+}}$ and the equation should be understood as an
equation of Itô-integrals. A wide class of processes $Y$ which may replace $W$
and still result in a well-defined equation of Itô-integrals is the class of
Itô-processes

$Y_{t}=Y_{0}+\int_{0}^{t}\alpha_{s}\mathrm{d}s+\int_{0}^{t}\beta_{s}\mathrm{d}W_{s},$

for suitable processes $\alpha=(\alpha_{t})_{t\in\mathbb{R}_{+}}$ and
$\beta=(\beta_{t})_{t\in\mathbb{R}_{+}}$.

The aim of this work is to bridge the dichotomy between concepts of correlated
noise stemming from discrete-time models of the form (1) and the possible
continuous-time adaptations in the fashion of SDEs like (4). While the Euler-
Mayurama method offers a consistent way of translating the latter to the
former, such a translation is decidedly not unique for the diametric task. We
claim to be able to reduce the number of possible continuous-time
constructions under certain assumptions. In a first step, we will invoke the
Markov-property of the AR(1)-process in (2) to search for a corresponding
continuous-time process (Section 2). Thereafter, we will focus on spectral
characteristics (Section 3). Alternatives will be considered in Sections 4 and
6 and an example will be presented in Section 5.

##  2 Red noise as a stationary Gaussian Markov-process

Starting from the discrete-time formulation of the red noise process in
equation (2), we can make a strong case for a specific continuous-time
analogue: Suppose $(\varepsilon_{k})_{k\in\mathbb{N}}$ is a sub-sample of a
stationary, measurable process $\alpha=(\alpha_{t})_{t\in\mathbb{R}_{+}}$ at
integer time-steps. Then, if one requires that $\alpha$ is a Gaussian Markov
process in continuous time in the same manner as $\varepsilon$ is a Gaussian
Markov process in discrete time, $\alpha$ has to be an Ornstein-Uhlenbeck
process. This is because all stationary, measurable processes which are
simultaneously Gaussian and Markov are of the Ornstein-Uhlenbeck type (Theorem
1.1 in [8]). Therefore, a distinctly motivated continuous-time analogue of (1)
would be

$\mathrm{d}X_{t}=f(X_{t},t)\mathrm{d}t+\alpha_{t}\mathrm{d}t.$

Discretizing this via the Euler method with integration step $\Delta t=1$
would result in the original equation (1). While we will later argue for the
same continuous-time red noise model, this derivation is not entirely
satisfactory. The a priori restriction $\mathrm{d}Y_{t}=\alpha_{t}\mathrm{d}t$
we have made before arguing for any specific process $\alpha$ is rather narrow
considering the wide range of possibilities offered in stochastic calculus and
it is not directly motivated from any characteristic of discrete-time red
noise. In fact, there is an inherent ambiguity which cannot be resolved from
the standpoint of a single discrete-time dynamic equation, i.e. without
knowledge about how the equation scales when the integration time-step $\Delta
t=1$ changes. The assumption we have made here, namely that the discrete-time
noise term for any choice of $\Delta t>0$ should be $\varepsilon_{k}^{(\Delta
t)}=\alpha_{k}\Delta t$, essentially implies that the autocovariance structure
of the discrete-time noise should scale uniformly with a factor of $(\Delta
t)^{2}$. This is decidedly not the case for uncorrelated noise $z_{k}^{(\Delta
t)}$, where the variance classically scales with a factor of $\Delta t$ to
ensure a convergence to the white-noise differential $\mathrm{d}W_{t}$.
However, there exists a characteristic which can be deduced from observations,
enjoys precedent in application, and elegantly implies a certain scaling-
behaviour of the discrete-time differential equation. This characteristic will
be a vanishing power spectral density in the limit of infinitely high
frequencies.

##  3 Characterization via the power spectral density

Throughout this section we will work on a filtered probability space
$(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\in\mathbb{R}_{+}},\mathbb{P})$
supporting a Brownian motion $W=(W_{t})_{t\in\mathbb{R}_{+}}$. For any
$Y_{0}\in\mathcal{F}_{0}$ and predictable processes
$\alpha=(\alpha_{t})_{t\in\mathbb{R}_{+}}$ and
$\beta=(\beta_{t})_{t\in\mathbb{R}_{+}}$ the Itô-process

$Y_{t}=Y_{0}+\int_{0}^{t}\alpha_{s}ds+\int_{0}^{t}\beta_{s}dW_{s},\ t\geq 0,$

is well-defined if e.g. the integrability condition

$\int_{0}^{T}\mathbb{E}\left[\alpha_{t}^{2}+\beta_{t}^{2}\right]\mathrm{d}t<\infty\text{
for all }T>0$ (5)

is satisfied. The power spectral density (PSD) of $\mathrm{d}Y$ is defined as

$S_{\mathrm{d}Y}(\omega):=\lim_{T\rightarrow\infty}\mathbb{E}\left[\frac{1}{T}\left|\int_{0}^{T}\exp(-i\omega
t)\mathrm{d}Y_{t}\right|^{2}\right],$

if the limit exists. One immediate result via the Itô-isometry is
$S_{\mathrm{d}W}\equiv 1$.

If $\mathrm{d}Y_{t}=\alpha_{t}\mathrm{d}t$ for some stationary, centered and
square-integrable process $\alpha$ with absolutely integrable autocovariance
structure $R_{\alpha}(\tau)$, then the Wiener-Khinchin theorem applies:

$S_{\alpha_{t}\mathrm{d}t}(\omega)=\mathcal{F}[R_{\alpha}(\tau)](\omega):=\int_{-\infty}^{\infty}\exp(-i\omega\tau)R_{\alpha}(\tau)\mathrm{d}\tau.$

In most applications of correlated noise, the properties of the PSD are more
prominently featured than the autocovariance counterpart. In fact, the name
red noise stems from the observation that low frequencies exhibit the largest
amplitudes in the PSD. A sufficient but not equivalent condition for this is
to demand an asymptotically vanishing spectral density, i.e.
$S(\omega)\rightarrow 0$ as $\omega\rightarrow 0$. In much of the applied
literature on this topic, a rate of decay of $\mathcal{O}(\omega^{-2})$ is
taken to be the defining characteristic of red noise [9, 10, 11] and the noise
instances are sometimes constructed directly via its PSD [12, 13, 14]. This
particular dependence on $\omega$ can either be derived from observations or
motivated by again assuming that the process $\alpha$ has the autocovariance
structure given in (3) and calculating the implied PSD of this process:

$S_{\alpha_{t}\mathrm{d}t}(\omega)=\mathcal{F}[R_{\alpha}(\tau)](\omega)=\frac{-2\log(\varphi)}{(1-\varphi^{2})(\log(\varphi)^{2}+\omega^{2})}=\mathcal{O}(\omega^{-2}).$

Note that we have not quite resolved the ambiguity of how the autocovariance
structure should scale with the integration time-step but have instead invoked
a stronger definition of red noise on the side of the PSD, namely that it
vanishes in the high-frequencies limit. With this stronger definition, we are
now able to formulate our main result, heavily constraining the possible
choices of continuous-time red noise equivalents that are modelled through an
Itô-differential.

######  Theorem.

Let $\alpha$ and $\beta$ be adapted processes satisfying the integrability
condition (5). Suppose that $\alpha$ has continuous paths and $\beta$ is
predictable. Define the Itô-process

$Y_{t}=\int_{0}^{t}\alpha_{s}\mathrm{d}s+\int_{0}^{t}\beta_{s}\mathrm{d}W_{s},\quad
t\geq 0.$

1\. Finite time horizon. Assume that the finite-time PSD of $\mathrm{d}Y$
vanishes in the limit of infinitely high frequencies, i.e.

$S^{(T)}_{\mathrm{d}Y}(\omega):=\mathbb{E}\left[\frac{1}{T}\left|\int_{0}^{T}\exp(-i\omega
t)\mathrm{d}Y_{t}\right|^{2}\right]\xrightarrow[\omega\to\infty]{}0.$

Then $\beta_{t}=0$ $\mathbb{P}$-a.s. for almost all $t\in[0,T]$.

2\. Infinite time horizon. Assume that $\alpha$ and $\beta$ are stationary
processes and let $\alpha$ be centered around $0$ with an absolutely
integrable autocovariance structure $R_{\alpha}(\tau)$. Assume that the PSD of
$\mathrm{d}Y$ on an infinite time horizon exists. If the PSD of $\mathrm{d}Y$
vanishes in the limit of infinitely high frequencies, i.e.

$S_{\mathrm{d}Y}(\omega)\xrightarrow[\omega\to\infty]{}0,$

then $\beta_{t}=0$ $\mathbb{P}$-a.s. for all $t\geq 0$.

We refer to Appendix A.1 for a proof of this theorem. If one objective of
finding a suitable red noise Itô-differential $\mathrm{d}Y$ is to have it
exhibit a vanishing PSD in the limit of infinitely high frequencies, then all
reasonable choices necessarily have the form
$\mathrm{d}Y_{t}=\alpha_{t}\mathrm{d}t$. The range of reasonable choices
defined through the conditions on $\alpha$ and $\beta$ in the theorem is
narrower if one examines an infinite time horizon, since asymptotic behaviour
needs to be taken into account. The strict stationarity condition for $\alpha$
and $\beta$ may be replaced by other suitable constraints on their asymptotic
behaviour. In the finite time horizon case, the restrictions on $\alpha$ and
$\beta$ reduce to being square-integrable together with the path-continuity of
$\alpha$, which may be motivated from physical principles.

Taking into account the arguments on Gaussian Markov processes from the last
section, we may conclude that the Ornstein-Uhlenbeck process $U$ defined by

$\mathrm{d}U_{t}=-\theta
U_{t}\mathrm{d}t+\mathrm{d}W_{t},\,U_{0}\sim\mathcal{N}\left(0,\frac{1}{2\theta}\right)$

constitutes a unique way of modelling red noise in continuous time via

$\mathrm{d}Y_{t}=\alpha_{t}\mathrm{d}t=U_{t}\mathrm{d}t$

The Ornstein-Uhlenbeck process $U$ exhibits both the property of exponentially
decaying autocovariance and the $\mathcal{O}(\omega^{-2})$ decay of the PSD of
the differential $U_{t}\mathrm{d}t$ that we are accustomed to from the
discrete-time case (see Figure 1b for an illustration and Appendix A.2 for a
derivation):

$\displaystyle R_{U}(\tau)$
$\displaystyle=\frac{1}{2\theta}\exp(-\theta\left|\tau\right|)$ $\displaystyle
S_{U_{t}\mathrm{d}t}(\omega)$ $\displaystyle=\frac{1}{\theta^{2}+\omega^{2}}$

With the condition of a vanishing PSD in the high-frequency limit, there will
be no structurally different reasonable choices of continuous-time red noise
models. This model is featured in some of the applied literature [15, 16, 17,
18], but to the best of our knowledge no comprehensive justification has
hitherto been brought forth in its favour.

##  4 Alternatives under weaker assumptions

In case that red noise should only have to exhibit strong amplitudes in low
frequencies, but not necessarily a vanishing PSD at $\omega\rightarrow\infty$,
we are again left with the aforementioned ambiguity of scaling in the
autocovariance structure. For correlated noise with this relaxed definition,
we propose one notable construction: Given a Brownian motion $W$ we consider
the SDE system

$\displaystyle\mathrm{d}U_{t}$ $\displaystyle=-\theta
U_{t}\mathrm{d}t+\mathrm{d}W_{t},$ $\displaystyle\mathrm{d}Y_{t}$
$\displaystyle=\gamma U_{t}\mathrm{d}t+\mathrm{d}W_{t}.$

In this case, the PSD is (see Appendix A.2 and Fig. 1d)

$S_{\mathrm{d}Y}(\omega)=\frac{(\gamma+\theta)^{2}+\omega^{2}}{\theta^{2}+\omega^{2}}.$

If $\gamma>0$, the condition of decaying PSD toward higher frequencies is
satisfied, but the limit for large frequencies is not 0:

$S_{\mathrm{d}Y}(\omega)\xrightarrow[\omega\to\infty]{}1.$

For $-2\theta<\gamma<0$ however, the results are different. One special case
of this bears mentioning, since it is at times proposed in the literature as a
continuous time red noise model [19, 7, 20]. For $\gamma=-\theta$, we have
$\mathrm{d}Y=\mathrm{d}U$ and it may seem conceptually appealing to label this
red noise. However, the PSD of $\mathrm{d}U$ exhibits the opposite of what we
have so far understood as red noise characteristics:

$S_{\mathrm{d}U}=\frac{\omega^{2}}{\theta^{2}+\omega^{2}}.$

For low frequencies, $S_{\mathrm{d}U}$ tends to 0 and it is monotonically
increasing in $\omega$ (see Figure 1c). The discrete-time noise terms
resulting from discretizing such a differential via the Euler-Mayurama method
at integration time-step $\Delta t$ would also be negatively correlated (here
$\tau\geq\Delta t$):

$\displaystyle\mathrm{Cov}\left(U_{\Delta t}-U_{0},U_{\tau+\Delta
t}-U_{\tau}\right)=\frac{1}{\theta}\exp(-\theta\tau)(1-\cosh(\theta\Delta
t))<0.$

The usage of $\mathrm{d}U$ in this context constitutes a common misconception
about the formulation of continuous-time stochastic models from desired
discrete-time characteristics. If one posits a certain distribution or
correlation in the noise component of a system and encounters the desired
property in a stochastic process $V$, then the differential $\mathrm{d}V$ is
in general not a suitable noise term since it may exhibit entirely different
properties. In [19] and [7], this would instead imply using
$V_{t}\mathrm{d}t$. Even when introducing $\mathrm{d}U$ more generally as
coloured noise [21, 22], one should be aware of the conceptual implications of
negative correlation. If the aim of introducing correlation is to model
temporal persistence in the noise forcing, then this would heuristically
always call for positive correlation.

Figure 1: The theoretically computed (orange) and observed (blue) PSD for each
of the discussed noise differentials, shown on a log-log-scale. The observed
PSD was obtained by taking the squared absolute value of the Fourier-
transformed noise signal and subsequently averaging over $10^{3}$ neighbouring
frequencies respectively. The length of the sample was $2\cdot 10^{7}$, which
implies a total of $10^{7}$ analyzed frequencies and $10^{4}$ after averaging.
(a) PSD of unit white noise $\mathrm{d}W$. (b) PSD of the red noise
differential $U_{t}\mathrm{d}t$. The characteristic $\omega^{-2}$ asymptotic
can be empirically observed until close to the Nyquist frequency, at which
point the discrete nature of the data skews the PSD. (c) PSD of the
differential $\mathrm{d}U$ sometimes erroneously referred to as red noise. (d)
PSD of the red noise alternative $\mathrm{d}Y_{t}=\gamma
U_{t}\mathrm{d}t+\mathrm{d}W_{t}$, $\gamma>0$, when allowing for a non-
vanishing PSD at $\omega\rightarrow\infty$. Parameter values where chosen as
$\theta=0.1$ and $\gamma=0.5$.

##  5 Example: Linearly restoring process

Apart from potentially being a more accurate modelling approach, the
continuous-time realm also allows for the use of more advanced methods from
stochastic analysis in the search of analytical solutions. We can observe this
for the case of an AR(1) process whose noise term is itself an AR(1) process:

$\displaystyle X^{(d)}_{k+1}$ $\displaystyle=\psi
X^{(d)}_{k}+\sigma\varepsilon_{k},\,X^{(d)}_{0}=x_{0}\in\mathbb{R}$
$\displaystyle\varepsilon_{k+1}$
$\displaystyle=\varphi\varepsilon_{k}+z_{k},\,\varepsilon_{0}=0$

Following our previous reasoning, we propose that

$\displaystyle\mathrm{d}X^{(c)}_{t}$ $\displaystyle=-\lambda
X^{(c)}_{t}\mathrm{d}t+\sigma
U_{t}\mathrm{d}t,\,X^{(c)}_{0}=x_{0}\in\mathbb{R}$
$\displaystyle\mathrm{d}U_{t}$ $\displaystyle=-\theta
U_{t}\mathrm{d}t+\mathrm{d}W_{t},\,U_{0}=0$

with $\lambda:=-\log(\psi)$ and $\theta:=-\log(\varphi)$ is a consistent
continuous-time analogue that captures the essence of the discrete-time model
and its noise characteristics. This allows us to solve for the Gaussian
process

$X^{(c)}_{t}=x_{0}\exp(-\lambda
t)+\sigma\int_{0}^{t}\exp(-\lambda(t-s))U_{s}\mathrm{d}s$

and its asymptotically stationary autocorrelation structure

$r_{X^{(c)}}(\tau)\asymp\sigma^{2}\frac{\lambda\exp(-\theta\left|\tau\right|)-\theta\exp(-\lambda\left|\tau\right|)}{\lambda-\theta}.$
(6)

It is in principle possible to obtain similar results in the discrete-time
realm by computing the first and second moment of $X^{(d)}_{k}$ in the
stationary limit via infinite series analysis. Performing the calculations in
the continuous-time realm, however, makes them both more feasible and more
robust. An interesting observation is that the asymptotic distribution of the
resulting process $X^{(c)}$ remains identical after interchanging the values
of $\lambda$ and $\theta$ in the set of the generating differential equations.
We illustrate the efficacy of the continuous-time modelling approach by
sampling both $X^{(d)}$ and $X^{(c)}$ and comparing their observed
autocorrelation structure with the theoretically derived autocorrelation
structure in equation (6) (see Figure 2). When integrating the continuous-time
stochastic differential equations for $X^{(c)}$ it is imperative to minimize
the errors introduced by numerical discretization. This is possible via the
following discrete-time representation of the Ornstein-Uhlenbeck process $U$.
Define for a time-step $\delta t$ the AR(1) process

$q_{k+1}=\exp(-\theta\delta
t)q_{k}+\frac{1}{\sqrt{2\theta}}\left(1-\exp(-2\theta\delta
t)\right)^{1/2}z_{k},\quad
q_{0}\sim\mathcal{N}\left(0,\frac{1}{2\theta}\right),$

where the $z_{k}$ are i.i.d. unit Gaussian. Then we have

$(q_{k})_{k=0,\dots,K-1}\overset{d}{=}(U_{k\delta t})_{k=0,\dots,K-1},$

meaning they have the same $K$-dimensional marginals. This may be used to
integrate continuous-time equations involving the Ornstein-Uhlenbeck process
with minimal loss of accuracy in distribution through discretization. Now
sampling $X^{(c)}$ at time-steps $\Delta t>0$ can be achieved by first
generating $U$ at time-steps $\delta t\ll\Delta t$, integrating
$\mathrm{d}X^{(c)}$ via the Euler-method at steps $\delta t$ and subsequently
sampling the result at steps $\Delta t$.

Figure 2: (a) A sample path of $X^{(d)}$ and (b) a sample path of $X^{(c)}$.
(c) Their respective estimated autocorrelation structure plotted in blue and
orange. The theoretical autocorrelation structure (see eq. (6)) of the
continuous-time model is plotted in red for comparison. The realization of the
discrete-time equations was computed directly from the sampled unit normals
$(z_{k})$. For the continuous-time equations, we sampled $U$ at intervals of
$\delta t=10^{-1}$, integrated the equation for $X^{(c)}$ via the Euler method
and subsequently sampled every $10$th value of the resulting time series as
described in the main text. The length of the considered time series is
$2\cdot 10^{7}$ each. The estimated and theoretical autocorrelation structures
do not exhibit significant differences (within $1$%). Parameter values were
chosen to be $\psi=0.8$, $\varphi=0.9$ and $\sigma=1$.

##  6 Beyond the Itô-framework of stochastic integration

In many applied fields of stochastic modelling, the possibility of using self-
similar processes like fractional Brownian motion $B^{H}_{t}$ with Hurst
parameter $H\neq 1/2$ in lieu of semimartingales for modelling noise has
gained great popularity in applications, cf. e.g. [23, 24]. The non-stationary
autocovariance ($\tau\geq 0$)

$R_{B^{H}}(t,\tau)=\frac{1}{2}(t^{2H}+(t+\tau)^{2H}-\tau^{2H})$

implies stationary increments $B^{H}_{t+\Delta
t}-B^{H}_{t}\sim\mathcal{N}(0,\Delta t^{2H})$ which are positively correlated
if and only if $H>1/2$. This so-called long memory case is interesting to us
because it allows for the modelling of persistence in the noise increments.
Since the paths of $B^{H}$ have vanishing $p$-variation for any $p>H^{-1}$,
the class of admissible integrand processes $\nu$ in the Young-framework of
integration is quite large: If the paths of $\nu$ possess finite
$2$-variation, then the pathwise Riemann-Stieltjes integral is well-defined.
It is therefore sensible to consider stochastic dynamical systems in the form
of integral equations where some terms are integrated with respect to
fractional Brownian motion. Extensive and rigorous theory on this topic may be
found in [25]. A consistent derivation of the PSD associated with the
differential $\mathrm{d}B^{H}_{t}$, often referred to as fractional Gaussian
noise, is given in [26]:

$S_{\mathrm{d}B^{H}_{t}}(\omega)=C_{H}\omega^{1-2H},$

where $C_{H}$ is a constant depending only on $H$. Hence the differential
exhibits a vanishing PSD in the limit of high frequencies and may generally be
considered as an alternative continuous-time modelling approach for the
concept of red noise. However, the persistence of fractional Gaussian noise in
terms of its autocovariance only decays with $\mathcal{O}(\tau^{2H-2})$, in
contrast to the usual exponential decay seen in discrete-time red noise. This
should give pause for concern when we consider to use fractional Gaussian
noise as a continuous-time red noise term. Fractional Gaussian noise is
instead often referred to as a long-memory noise term [27, 28, 29] in analogy
to its continuous-time origin.

## References

  * [1] R. Zwanzig.  Nonequilibrium Statistical Mechanics.  Oxford University Press, 2001.
  * [2] R. Zwanzig.  Memory effects in irreversible thermodynamics.  Phys. Rev., 124:983–992, Nov 1961.
  * [3] A. J. Chorin, O. H. Hald, and R. Kupferman.  Optimal prediction and the Mori-Zwanzig representation of irreversible processes.  Proceedings of the National Academy of Sciences, 97(7):2968–2973, 2000.
  * [4] M. Schwager, K. Johst, and F. Jeltsch.  Does red noise increase or decrease extinction risk? Single extreme events versus series of unfavorable conditions.  The American Naturalist, 167(6):879–888, 2006.
  * [5] M. Rodal, S. Krumscheid, G. Madan, J. Henry LaCasce, and N. Vercauteren.  Dynamical stability indicator based on autoregressive moving-average models: Critical transitions and the Atlantic meridional overturning circulation.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 32(11):113139, 2022.
  * [6] M. E. Mann and J. M. Lees.  Robust estimation of background noise and signal detection in climatic time series.  Climatic Change, 33(3):409–445, 1996.
  * [7] C. Boettner and N. Boers.  Critical slowing down in dynamical systems driven by nonstationary correlated noise.  Phys. Rev. Research, 4:013230, Mar 2022.
  * [8] J. L. Doob.  The Brownian movement and stochastic equations.  Annals of Mathematics, 43(2):351–369, 1942.
  * [9] K. Hasselmann.  Stochastic climate models Part I. Theory.  Tellus, 28(6):473–485, 1976.
  * [10] Z. Liao, K. Ma, S. Sarker, M. S.and Tang, H. Yamahara, M. Seki, and H. Tabata.  Quantum analog annealing of gain-dissipative Ising machine driven by colored Gaussian noise.  Advanced Theory and Simulations, 5(3):2100497, 2022.
  * [11] N. Shibazaki, R. F. Elsner, and M. C. Weisskopf.  The Effect of Decay of the Amplitude of Oscillation on Random Process Models for QPO X-Ray Stars.  The Astrophysical Journal, 322:831, November 1987.
  * [12] J. Timmer and M. Koenig.  On generating power law noise.  Astronomy and Astrophysics, 300:707, August 1995.
  * [13] H. Zhivomirov.  A method for colored noise generation.  Romanian Journal of Acoustics and Vibration, 15:14–19, 2018.
  * [14] L. Kiss, Z. Gingl, Z. Márton, J. Kertész, F. Moss, G. Schmera, and A. Bulsara.  1/f noise in systems showing stochastic resonance.  Journal of Statistical Physics, 70(1):451–462, 1993.
  * [15] P. Hänggi and P. Jung.  Colored Noise in Dynamical Systems.  John Wiley & Sons, Ltd, 1994.
  * [16] P. Hänggi, P. Jung, C. Zerbe, and F. Moss.  Can colored noise improve stochastic resonance?  Journal of Statistical Physics, 70(1):25–47, 1993.
  * [17] Z Liu, P Gu, and T. L. Delworth.  Strong red noise ocean forcing on atlantic multidecadal variability assessed from surface heat flux: Theory and application.  Journal of Climate, pages 1 – 56, 2022.
  * [18] M. Newman, P. D. Sardeshmukh, and C. Penland.  Stochastic forcing of the wintertime extratropical flow.  Journal of the Atmospheric Sciences, 54(3):435 – 455, 1997.
  * [19] B. Bercu, F. Proïa, and N. Savy.  On Ornstein-Uhlenbeck driven by Ornstein-Uhlenbeck processes.  Statistics & Probability Letters, 85, 12 2012.
  * [20] P. Ditlevsen.  Observation of $\alpha$-stable noise induced millenial climate changes from an ice record.  Geophysical Research Letters - GEOPHYS RES LETT, 26:1441–1444, 05 1999.
  * [21] C. Kuehn, K. Lux, and A. Neamţu.  Warning signs for non-Markovian bifurcations: colour blindness and scaling laws.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 478(2259):20210740, 2022.
  * [22] E. Kolářová and L. Brancik.  Stochastic differential equations describing systems with coloured noise.  Tatra Mountains Mathematical Publications, 71:99–107, 12 2018.
  * [23] F. Comte and E. Renault.  Long memory in continuous-time stochastic volatility models.  Mathematical Finance, 8(4):291–323, 1998.
  * [24] M. Rypdal and K. Rypdal.  Long-memory effects in linear response models of earth’s temperature and implications for future global warming.  Journal of Climate, 27(14):5240 – 5258, 2014.
  * [25] J.S. Mišura.  Stochastic Calculus for Fractional Brownian Motion and Related Processes.  Number Nr. 1929 in Lecture Notes in Mathematics. Springer, 2008.
  * [26] M. Li and S. Lim.  A rigorous derivation of power spectrum of fractional Gaussian noise.  Fluctuation and Noise Letters, 6, 12 2006.
  * [27] Y. Chen, Y. Li, and L. Tian.  Moment estimator for an AR(1) model driven by a long memory Gaussian noise.  Journal of Statistical Planning and Inference, 222, 06 2022.
  * [28] D. I. Vyushin and P. J. Kushner.  Power-law and long-memory characteristics of the atmospheric general circulation.  Journal of Climate, 22(11):2890 – 2904, 2009.
  * [29] S. H. Sørbye, E. Myrvoll-Nilsen, and H. Rue.  An approximate fractional Gaussian noise model with $\mathcal{O}(n)$ computational cost.  Statistics and Computing, 29(4):821–833, 2019.

##  Appendix A Appendix

###  A.1 Proof of the Theorem

1\. Finite time horizon. First, we show that the constant term in the finite-
time PSD stemming from the $\beta_{t}\mathrm{d}W_{t}$ term must be matched by
the PSD of $\alpha_{t}\mathrm{d}t$ in order to result in an overall PSD which
is vanishing at $\omega\rightarrow\infty$.

$\displaystyle S^{(T)}_{\mathrm{d}Y}(\omega)$
$\displaystyle=\mathbb{E}\left[\frac{1}{T}\left|\int_{0}^{T}\exp(-i\omega
t)\mathrm{d}Y_{t}\right|^{2}\right]$
$\displaystyle=\mathbb{E}\left[\frac{1}{T}\left|\int_{0}^{T}\exp(-i\omega
t)\alpha_{t}\mathrm{d}t+\int_{0}^{T}\exp(-i\omega
t)\beta_{t}\mathrm{d}W_{t}\right|^{2}\right]$
$\displaystyle=\frac{1}{T}\left\lVert\int_{0}^{T}\exp(-i\omega
t)\alpha_{t}\mathrm{d}t+\int_{0}^{T}\exp(-i\omega
t)\beta_{t}\mathrm{d}W_{t}\right\rVert_{L^{2}(\Omega)}^{2}$
$\displaystyle\geq\left[\frac{1}{\sqrt{T}}\left\lVert\int_{0}^{T}\exp(-i\omega
t)\alpha_{t}\mathrm{d}t\right\rVert_{L^{2}(\Omega)}-\frac{1}{\sqrt{T}}\left\lVert\int_{0}^{T}\exp(-i\omega
t)\beta_{t}\mathrm{d}W_{t}\right\rVert_{L^{2}(\Omega)}\right]^{2}$
$\displaystyle=\left(\sqrt{S^{(T)}_{\alpha_{t}\mathrm{d}t}(\omega)}-\sqrt{\frac{1}{T}\int_{0}^{T}\mathbb{E}\left[\beta_{t}^{2}\right]\mathrm{d}t}\right)^{2},$

where we used the inverse triangle inequality. Since we assumed
$S^{(T)}_{\mathrm{d}Y}(\omega)\xrightarrow[\omega\to\infty]{}0$, we conclude

$S^{(T)}_{\alpha_{t}\mathrm{d}t}(\omega)\xrightarrow[\omega\to\infty]{}\frac{1}{T}\int_{0}^{T}\mathbb{E}\left[\beta_{t}^{2}\right]\mathrm{d}t.$

It remains to prove that this limit is $0$. Since we have uniform
integrability of the random variables $\int_{0}^{T}\exp(-i\omega
t)\alpha_{t}\mathrm{d}t$ for all choices of $\omega\in\mathbb{R}$ via

$\left\lVert\int_{0}^{T}\exp(-i\omega
t)\alpha_{t}\mathrm{d}t\right\rVert_{L^{2}(\Omega)}\leq\int_{0}^{T}\left\lVert\alpha_{t}\right\rVert_{L^{2}(\Omega)}\mathrm{d}t<\infty,$

it suffices to prove that

$\int_{0}^{T}\exp(-i\omega
t)\alpha_{t}\mathrm{d}t\xrightarrow[\omega\rightarrow\infty]{a.s.}0.$ (7)

But each path $a:[0,T]\rightarrow\mathbb{R}$ of $\alpha$ is uniformly
continuous on the compact domain, so for any $\epsilon>0$, we find a
$\delta_{\epsilon}>0$ such that for
$\omega>\omega_{\epsilon}:=2\pi/\delta_{\epsilon}$ we may estimate

$\left|\int_{t_{0}}^{t_{0}+\delta_{\epsilon}}\cos(\omega
t)a(t)\mathrm{d}t\right|\leq\frac{\epsilon}{2}\int_{t_{0}}^{t_{0}+\delta_{\epsilon}}\left|\cos(\omega
t)\right|\mathrm{d}t.$

for all $t_{0}\in[0,T-\delta_{\epsilon}]$. So the real part of the integral in
question becomes arbitrarily small for large $\omega>\omega_{\epsilon}$:

$\left|\int_{0}^{T}\cos(\omega
t)a(t)\mathrm{d}t\right|\leq\frac{\epsilon}{2}\int_{0}^{T}\left|\cos(\omega
t)\right|\mathrm{d}t\leq\frac{\epsilon}{2}T.$

The real part and similarly the imaginary part of (7) converge to $0$ for
every path of $\alpha$. Hence the limit in question is also $0$:

$S^{(T)}_{\alpha_{t}\mathrm{d}t}(\omega)\xrightarrow[\omega\to\infty]{}0=\frac{1}{T}\int_{0}^{T}\mathbb{E}\left[\beta_{t}^{2}\right]\mathrm{d}t.$

We may conclude that for almost all $t\in[0,T]$ we have $\beta_{t}=0$
$\mathbb{P}$-almost surely.

2\. Infinite time horizon. By identical arguments as in the first part of the
proof, we arrive at

$S_{\alpha_{t}\mathrm{d}t}(\omega):=\lim_{T\rightarrow\infty}\mathbb{E}\left[\frac{1}{T}\left|\int_{0}^{T}\exp(-i\omega
t)\alpha_{t}\mathrm{d}t\right|^{2}\right]\xrightarrow[\omega\to\infty]{}\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\mathbb{E}\left[\beta_{t}^{2}\right]\mathrm{d}t=\mathbb{E}\left[\beta_{0}^{2}\right],$

on the condition that the limit of $T\rightarrow\infty$ in the above
definition of $S_{\alpha_{t}\mathrm{d}t}$ exists. Because $\alpha$ is square-
integrable, its autocovariance $R_{\alpha}(\tau)$ is bounded by
$\mathrm{Var}\left(\alpha_{0}\right)$ for all $\tau\in\mathbb{R}$. This,
together with the absolute integrability of $R_{\alpha}$ implies that
$R_{\alpha}\in L^{2}(\mathbb{R})$. By the Wiener-Khinchin theorem and
Plancherel’s theorem, we deduce

$S_{\alpha_{t}\mathrm{d}t}=\mathcal{F}[R_{\alpha}]\in L^{2}(\mathbb{R})$

and thereby confirm that the limit in the definition of
$S_{\alpha_{t}\mathrm{d}t}$ exists. We have established that
$S_{\alpha_{t}\mathrm{d}t}$ converges to a constant as
$\omega\rightarrow\infty$. In order for $S_{\alpha_{t}\mathrm{d}t}$ to be
square-integrable, the constant in question must be $0$. Hence, we must have

$\mathbb{E}\left[\beta_{t}^{2}\right]=\mathbb{E}\left[\beta_{0}^{2}\right]=0$

which implies $\beta_{t}=0$ $\mathbb{P}$-a.s. for all $t\geq 0$.

∎

###  A.2 Derivation of the power spectral densities

For unit white noise $\mathrm{d}W$, the calculation is a straightforward
application of the Itô-isometry:

$S_{\mathrm{d}W}(\omega)=\lim_{T\rightarrow\infty}\mathbb{E}\left[\frac{1}{T}\left|\int_{0}^{T}\exp(-i\omega
t)\mathrm{d}W_{t}\right|^{2}\right]=\lim_{T\rightarrow\infty}\mathbb{E}\left[\frac{1}{T}\int_{0}^{T}1\mathrm{d}t\right]=1$

For the other claims, first define the Ornstein-Uhlenbeck process via the SDE

$\mathrm{d}U_{t}=-\theta U_{t}\mathrm{d}t+\mathrm{d}W_{t},\quad
U_{0}\sim\mathcal{N}\left(0,\frac{1}{2\theta}\right),$

so that $U$ is stationary with autocovariance
$\mathrm{Cov}\left(U_{t},U_{t+\tau}\right)=\frac{1}{2\theta}\exp(-\theta\left|\tau\right|)$
and compute via the Fubini-Tonelli Theorem

$\displaystyle S_{U_{t}\mathrm{d}t}(\omega)$
$\displaystyle=\lim_{T\rightarrow\infty}\mathbb{E}\left[\frac{1}{T}\left|\int_{0}^{T}\exp(-i\omega
t)U_{t}\mathrm{d}t\right|^{2}\right]$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\int_{0}^{T}\exp(-i\omega(t-s))\mathrm{Cov}\left(U_{t},U_{s}\right)\mathrm{d}t\mathrm{d}s$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\int_{0}^{T}\exp(-i\omega(t-s))\frac{1}{2\theta}\exp(-\theta\left|t-s\right|)\mathrm{d}t\mathrm{d}s$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}f_{1}(T,\omega,\theta)=\frac{1}{\theta^{2}+\omega^{2}}.$

The function $f_{1}$ can be computed to be

$f_{1}(T,\omega,\theta)=\frac{T}{\theta^{2}+\omega^{2}}+\frac{(\omega^{2}-\theta^{2})(1-\exp(-\theta
T)\cos(\omega T))-4\theta\omega\exp(-\theta
T))\sin(T\omega)}{\theta(\theta^{2}+\omega^{2})}$

and we will make use of it again in the derivation of $S_{\mathrm{d}Y}$ where
$\mathrm{d}Y_{t}=\gamma U_{t}\mathrm{d}t+\mathrm{d}W_{t}$:

$\displaystyle S_{\mathrm{d}Y}(\omega)$
$\displaystyle=\lim_{T\rightarrow\infty}\mathbb{E}\left[\frac{1}{T}\left|\int_{0}^{T}\exp(-i\omega
t)\mathrm{d}Y_{t}\right|^{2}\right]$
$\displaystyle=\lim_{T\rightarrow\infty}\mathbb{E}\left[\frac{1}{T}\left|\int_{0}^{T}\exp(-i\omega
t)\gamma U_{t}\mathrm{d}t+\int_{0}^{T}\exp(-i\omega
t)\mathrm{d}W_{t}\right|^{2}\right]$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\Bigg{[}\gamma^{2}f_{1}(T,\omega,\theta)+\gamma\mathbb{E}\left[\int_{0}^{T}\exp(-i\omega
t)\mathrm{d}W_{t}\int_{0}^{T}\exp(i\omega t)U_{t}\mathrm{d}t\right]$
$\displaystyle\qquad\qquad\quad+\gamma\mathbb{E}\left[\int_{0}^{T}\exp(-i\omega
t)U_{t}\mathrm{d}t\int_{0}^{T}\exp(i\omega t)\mathrm{d}W_{t}\right]+T\Bigg{]}$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\left[\gamma^{2}f_{1}(T,\omega,\theta)+\gamma
f_{2}(T,\omega,\theta)+\gamma\overline{f_{2}(T,\omega,\theta)}+T\right]$

where we introduced

$\displaystyle f_{2}(T,\omega,\theta)$
$\displaystyle:=\mathbb{E}\left[\int_{0}^{T}\exp(-i\omega
t)\mathrm{d}W_{t}\int_{0}^{T}\exp(i\omega t)U_{t}\mathrm{d}t\right]$
$\displaystyle\>=\mathbb{E}\left[\int_{0}^{T}\int_{0}^{T}\exp(-i\omega(s-t))U_{t}\mathrm{d}W_{s}\mathrm{d}t\right]$
$\displaystyle\>=\mathbb{E}\left[\int_{0}^{T}\int_{0}^{T}\int_{0}^{t}\exp(-i\omega(s-t))\exp(-\theta(t-r))\mathrm{d}W_{r}\mathrm{d}W_{s}\mathrm{d}t\right].$

In the last step we inserted the closed form solution to the Ornstein-
Uhlenbeck SDE, where $U_{0}$ is independent of $W_{t}$ for all $t\geq 0$:

$U_{t}=U_{0}\exp(-\theta t)+\int_{0}^{t}\exp(-\theta(t-r))\mathrm{d}W_{r}$

The integration bounds in the integral with respect to $\mathrm{d}W_{s}$ may
as well be $(0,t)$, since the integral in $(t,T)$ is independent of the
innermost integral and their expectations are $0$ each. Using the Itô isometry
again we get

$\displaystyle f_{2}(T,\omega,\theta)$
$\displaystyle=\int_{0}^{T}\int_{0}^{t}\exp(-i\omega(s-t))\exp(-\theta(t-s))\mathrm{d}s\mathrm{d}t$
$\displaystyle=\frac{T(\theta-i\omega)+e^{-T(\theta-i\omega)}-1}{(\theta-i\omega)^{2}}.$

We see that

$\lim_{T\rightarrow\infty}\frac{1}{T}f_{1}(T,\omega,\theta)=\frac{1}{\theta^{2}+\omega^{2}}=\lim_{T\rightarrow\infty}\frac{1}{T}\frac{1}{2\theta}\left(f_{2}(T,\omega,\theta)+\overline{f_{2}(T,\omega,\theta)}\right)$

and so the spectral density is simply

$\displaystyle S_{\mathrm{d}Y}(\omega)$
$\displaystyle=\lim_{T\rightarrow\infty}\frac{1}{T}\left[\gamma^{2}f_{1}(T,\omega,\theta)+\gamma
f_{2}(T,\omega,\theta)+\gamma\overline{f_{2}(T,\omega,\theta)}+T\right]$
$\displaystyle=\frac{\gamma^{2}+2\theta\gamma}{\theta^{2}+\omega^{2}}+1=\frac{(\gamma+\theta)^{2}+\omega^{2}}{\theta^{2}+\omega^{2}}.$