# Realistic Kilonova Up Close

Alexandra Ruth Stewart  Li-Ta Lo  Oleg Korobkin  Irina Sagert  Julien Loiseau
Hyun Lim  Mark Alexander Kaltenborn  Christopher Michael Mauney  Jonah Maxwell
Miller  CCS-2, Computational Physics and Methods, Los Alamos National
Laboratories, Los Alamos, NM 87545, USA  CCS-3, Information Sciences, Los
Alamos National Laboratories, Los Alamos, NM 87545, USA  CCS-7, Applied
Computer Science, Los Alamos National Laboratories, Los Alamos, NM 87545, USA
HPC Environments, Los Alamos National Laboratories, Los Alamos, NM 87545, USA
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA  Department of Physics, The George Washington University,
Washington, DC 20052, USA

###### Abstract

Neutron star mergers are cosmic catastrophes that produce some of the most
energetic observed phenomena: short gamma-ray bursts, gravitational wave
signals, and kilonovae. The latter are optical transients, powered by
radioactive nuclides which are synthesized when the neutron-rich ejecta of a
disrupted neutron star undergoes decompression. We model this decompression
phase using data from simulations of post-merger accretion disk winds. We use
smoothed particle hydrodynamics with realistic nuclear heating to model the
expansion over multiple scales, from initially several thousand km to billions
of km. We then render a realistic image of a kilonova ejecta as it would
appear for a nearby observer. This is the first time such a visualization is
performed using input from state-of-the-art accretion disk simulations,
nuclear physics and atomic physics. The volume rendering in our model computes
an opacity transfer function on the basis of the physical opacity, varying
significantly with the inhomogeneity of the neutron richness in the ejecta.
Other physical quantities such as temperature or electron fraction can be
visualized using an independent color transfer function. We discuss several
difficulties with the ParaView application that we encountered during the
visualization process, and give descriptions of our solutions and workarounds
which could be used for future improvements.

###### keywords:

scientific visualization , high performance computing , astrophysics , neutron
star merger , kilonova , smoothed particle hydrodynamics

††journal: Parallel Computing

##  1 Introduction

In recent years, coalescences of binary neutron stars (BNS) and black hole
(BH)-NS systems have taken central stage in astrophysics due to their multi-
messenger character and wide-ranging opportunities to advance areas of
physics, including gravitational waves, high-density nuclear matter, extreme
gravity, and the origin of heavy elements. The first confirmed BNS merger was
a multi-messenger event which manifested itself via a gravitational wave
signal [1], a short gamma-ray burst [2, 3], and a range of transients spanning
the entire electromagnetic spectrum, from X-rays [4] to radio emissions [5].
The latter, originates from heavy-element nucleosynthesis which takes place in
the ejecta of NS mergers [6]. The heat from radioactive decay following the
rapid neutron capture process (the $r$-process) powers a unique astrophysical
transient known as a “kilonova”. This transient lasts for about a week, is
characterized by rapid reddening, and becomes a thousand times more luminous
than a conventional nova; hence its name. Kilonovae are essentially the
$r$-process observed in-the-making with the first one, AT2017gfo, being
discovered in coincidence with the first BNS gravitational wave detection,
GW170817 [7, 8].

Multiple kilonova models have been proposed to interpret the signal of
AT2017gfo (see e.g. review by Metzger [9]). However, the complexity of the
problem precludes definitive conclusions [10]. The most accurate observations
of AT2017gfo consist of time-dependent ultraviolet, optical, and infrared
spectra [11, 12], currently unmatched by theoretical models. This prevents
scientists from answering some of the critical questions about kilonovae, such
as: how much $r$-process is produced? What is the detailed composition of the
ejecta? Are BNS mergers the main $r$-process sites in the universe?

Most current models assume that the ejecta contains two well-separated
components with uniform composition: very neutron-rich “red” ejecta, producing
the heaviest elements, and moderately neutron-rich “blue” ejecta. The colors
come from the fact that heavy nuclei include lanthanides and actinides, which
have notoriously high opacity in the opital bands [13], emitting light mostly
in the near infrared. The “blue” ejecta is lanthanide-free and may thus
produce a short-lived “blue” optical transient. Progress can be made by more
comprehensive and advanced models, for example via the inclusion of more
complex compositional and morphological elements of the ejecta. In this
project, we go beyond this assumption and simulate more general, non-uniform
inhomogeneous expansion of the ejected material in 3D.

There are several mechanisms that can unbind neutron-rich matter from a BNS
merger [14]. Here, we focus on the wind ejecta from a post-merger accretion
disk [15]. Such a disk is left after the central hypermassive merger remnant
collapses to a BH. Visualizing the kilonova helps to gain insights in the ways
the morphology and composition are featured in the observed spectra and is a
step towards an accurate interpretation and characterization of kilonova
spectra for AT2017gfo and future detections.

Our pipeline incorporates several state-of-the art tools developed at Los
Alamos National Laboratory (LANL) to tackle different facets of the problem.
The most realistic, state-of-the art morphologies for accretion disk wind and
ejecta to date have been recently provided by Miller et al. [15] who used the
$\nu$bhlight code [16] to simulate the disk dynamics in a fixed general-
relativistic background including Monte-Carlo neutrino transport. This
simulation required approximately 400 cores for 3 weeks of walltime. As a
hybrid simulation, it used over one million finite volume zones, two million
tracer particles, and ten million Monte Carlo packets. An example of the
results can be seen in Figure 1.

Kilonova ejecta expansion over many scales is then simulated with FleCSPH
[17], which is a smoothed particle hydrodynamics (SPH) code. SPH is an
explicit numerical mesh-free Lagrangian method designed to solve the equations
of hydrodynamics by representing the fluid with a set of moving and
interacting particles [18]. Each particle is equipped with a smoothing length
and a kernel, measuring the range and intensity of hydrodynamic interactions
with other particles, respectively. SPH possesses excellent conservational
properties, adaptive resolution across multiple scales, and natural treatment
of vacuum. The FleCSPH code is based on FleCSI [19]—a compile-time
configurable framework, which supports parallel and distributed computing for
multiple data topologies, such as mesh, n-array, or ntree. FleCSPH’s ntree is
a parallel binary, quad, and octree used for neighbor search and gravitational
interactions in one, two, and three spatial dimensions, respectively. Features
currently implemented in FleCSPH include fast multipole method for long-range
gravitational interactions [20], variety of effective potentials, tabulated
equation of state (EOS), tabulated radioactive heating source, material
strength, and both Newtonian and relativistic SPH [21], making it perfectly
suitable for the task.

##  2 Astrophysical model

In this work, we focus on the expanding ejecta and the kilonova aspect of a NS
merger. It is believed that in most cases the heavy remnant of the merger
collapses to a BH, leaving behind an accretion disk [9]. The latter ejects
neutron-rich material which becomes a production site for $r$-process
elements. It has recently been demonstrated that the morphology of the ejecta
plays a substantial role in shaping the spectra and light curves of the
kilonova [10, 22]. We captured the output flux from the accretion disk
simulations with $\nu$bhlight [16] as initial data for the kilonova expansion.

Figure 1:  Initial distribution of the electron fraction (number of electrons
per nucleon) for accretion disk ejecta in velocity space [15]. The velocity is
sampled at the times when particles exit the outflow boundary of the spherical
grid. The ejecta is discretized with 27M particles. The red color corresponds
to high neutron richness with strong $r$-process nucleosynthesis. The blue
color indicates high electron fraction and lower opacity, marking regions that
are potentially responsible for the bright optical emission and the “blue”
kilonova.

The expansion of the kilonova covers multiple orders of magnitude in scale,
density, and temperature. With that, a Lagrangian particle method, such as
SPH, is better suited than a grid-based code to model this phase. To map the
ejecta data from $\nu$bhlight into FleCSPH, we sample the outgoing flux from
the accretion disk wind simulation and convert it into particles using the
Wendland $C_{6}$ kernel for the SPH interpolation (see reviews of the SPH
method for details [18]). We discretize the flux emerging from the spherical
boundary using regular time intervals $\Delta t$. The total mass of the ejecta
$m_{\rm ej}$ is computed by integrating the flux over the entire accretion
disk simulation time (about 1 second). Assuming the ejecta will be represented
by $N_{p}$ particles of equal mass $m_{p}=m_{\rm ej}/N_{p}$, we estimate how
many particles will emerge from the boundary during $\Delta t$, and then
randomly sample that many particles using the normalized flux as the
probability distribution over the sphere. To avoid unphysical inwardly
directed pressure gradients at the inflow boundary, we prepare the
compactified configuration of all the particles inside the spherical inflow
boundary using an inversion map. During the simulation, particles inside the
sphere are advected using the time-dependent inversion, passing through the
boundary at times such that their emergent configuration correctly mimics the
outgoing flux. The equation of advection for a particle $a$ inside the initial
sphere (with radius $R_{\rm in}$) is:

${r_{a}(t)=\frac{R_{\rm in}}{1+v_{r}(t-t_{\rm ex})/R_{\rm in}}}$ (1)

where $r_{a}$ is the radial distance of the particle from the origin, $v_{r}$
is its radial velocity upon exit, and $t_{\rm ex}$ is the expected exit time.
Because the density in our SPH formulation is computed using the local density
of particles, this setup provides continuity of the flux and ensures accurate
pressure gradients at the boundary. As soon as particles exit the injection
sphere ($R_{\rm in}=3500$ km in our case), they undergo regular hydrodynamic
evolution. Figure 2 demonstrates this initial compactified particle
distribution.

Figure 2:  Sample initial compactified particle configuration, color-coded
according to the internal energy that particles will have when they cross the
boundary. Inside the spherical region, the particles move radially according
to the formula (1). When the particles cross the boundary and enter the SPH
domain, their density is in accordance with the simulated outgoing flux from
the grid-based code $\nu$bhlight, and other physical quantities are
interpolated from the grid. Thus, the outer boundary of the grid-based code
serves as the inner flux boundary of the SPH code.

Since $r$-process nucleosynthesis happens on timescales of a second and it is
generally very challenging to include an inline $r$-process network in
hydrodynamic simulations, we apply precomputed radioactive heating
contributions. The latter are given in form of tables and were computed with
WinNet [23] using the Helmholtz EOS [24] and asumming that $40\%$ of the
generated energy escapes in the form of neutrinos. For validity checks, we
also compared with results from SkyNet [25] and PRISM [26] for a few
representative trajectories.

At early times of kilonova expansion, the internal energy of the outflow is
heavily dominated by radiation, followed by a subdominant contribution from
electrons and an almost negligible contribution from ions. For this reason, we
find it sufficient to use the ideal fluid equation of state
$P=(\Gamma-1)\rho\varepsilon$, with $P$, $\rho$, and $\varepsilon$
representing the pressure, density, and specific internal energy,
respectively, and $\Gamma=4/3$. At early times, the optical depth of the
outflow is extremely high, such that radiation is fully trapped. Simple
estimates for density and temperature of the ejecta show that throughout the
morphology-forming times (the first minute) the contribution to energy and
pressure from ions is completely negligible. The electron contribution can be
neglected as well as it is heavily subdominant. Contributions from electron-
positron pairs can be estimated following e.g. [27], but would only make
significant differences at high temperatures which exceed $5\times 10^{9}$ K.
This is not the case in our simulations.

The ejecta is expanded in a static gravitational potential of a central point
mass, mimicking the $2.58$-solar-mass BH from the original accretion disk
simulation, with added softening of the $1/r$ potential at a length scale of
$\varepsilon=10^{3}\;{\rm km}$. This simulation used 128 cores for 48 hours.
After about 30 seconds of simulation time, some of the particles form dynamic
fallback. Since, in this work, we are interested in the expanding part of the
flow with progressively larger timesteps and sizes, we truncate these
particles to avoid complications due to potentially small timesteps near the
origin. A similar method was employed in Rosswog et al. [14]. Nuclear heating
and gravity reshape the ejecta: the former makes it more isotropic [18], while
the latter slows it down and takes up some of the material, forming dynamical
fallback. Without these two factors, the output flux would just expand as
shown in Figure 1. The figure shows the initial electron fraction in the
ejecta, with red indicating more neutron-rich and blue less neutron-rich
material. The “red” ejecta would subsequently synthesize lanthanides while the
“blue” part would be lanthanide-free and produce the “blue” kilonova. As can
be seen, there is a broad equatorial belt of very neutron-rich material, while
the polar and medium-latitude regions are occupied by the “blue” component.

##  3 Visualization

ParaView offers several options to visualize particles from a SPH simulation.
One example is the Surface representation which we can use to obtain a quick
overview of the simulation results. We can also apply either the Point
Gaussian representation or the Glyph filter to plot particles as spheres with
sizes defined by the smoothing length. However, since the objective of SPH is
to represent the continuum with particles, we are most interested in applying
volume rendering to our simulation results.

###  3.1 Volume rendering of SPH particles

Paraview includes a SPH Volume Interpolation filter to turn particles into a
uniform grid which can then be visualized with volume rendering. As shown in
Figure 3a, the Interpolator takes two input datasets, a uniform grid provided
by the vtkBoundedVolumeSource and the SPH particles. Through the graphical
user interface (GUI), a user can set the bounding box of the particle dataset
and configure the ${xyz}$-dimensions, origin and spacing for the uniform grid.
A user can also specify various aspects of the SPH particles, for example,
arrays for particle mass, density, and cutoff radius (proportional to
smoothing length) as well as the SPH smoothing kernel to match the one used in
the simulation. At first sight, this filter seemed to be the perfect tool for
our work. However, we soon encountered several major difficulties which made
it unsuitable for our purposes.

The first issue lies in the design of the current ParaView pipeline
architecture. Although ParaView provides a widget for users to interactively
specify the bounding box of the particle dataset, the Interpolator cannot
automatically calculate the bounding box from its particles input and use the
result to configure the vtkBoundedVolumeSource input. This would require
_information_ to flow in the inverse direction of the pipleline. Over the
course of our simulation, the accretion disk ejecta expands from the size of a
small planet to the size of the solar system. The lack of automatic adjustment
of the bounding box poses a major difficulty. We were therefore required to
use an external interpolation to manually compute and set the bounding box
through the GUI when producing animations.

The second and more severe issue lies in the implementation of the
interpolation algorithm in the SPH Volume Interpolation filter. Here, the
Interpolator uses a user-specified parameter called spatial step to accelerate
the search of particles. As shown in Figure 3b, for each grid point G in the
uniform grid, the Interpolator performs a spatial range search within the
radius of spatial step. It then iterates through all the particles _within_
the radius, calculates, and accumulates contributions from each of them. This
approach works well when particle densities and thus cutoff radii don’t vary
by a large amount. However, in our case they span several orders of magnitude.
Using a single fixed search radius therefore leads to inaccurate interpolation
results. In our simulation, particles near the center of the ejecta have much
higher densities compared to the particles at the periphery. In lower-density
regions, particles are also spaced more sparsely. As depicted in Figure 3b, a
uniform spatial step will be able to find the smaller, higher density
particles H and J near the grid point G but will miss the larger but lower
density particle L whose center is outside of the search radius and completely
ignores the contribution from L to G. As a result, particles in the outer
regions of the ejecta are mostly unaccounted for and are invisible in the
volume rendered image.

We also encountered some further unidentified issues with the SPH Volume
Interpolator, resulting in some un-physical interpolation. For example, the
electron fraction $Y_{e}$ is physically constrained in the range of $[0,0.5]$.
However, the Interplator sometimes gives us either negative values or values
much larger than 1. All the described difficulties prompted us to implement
our own SPH particle-to-mesh interpolator which we will describe in the
following. Fortunately, it is straightforward to generalize the standard SPH
interpolation formula so that it can be applied to grid points [18]. For a
physical quantity $F_{b}$ specified on particles, the interpolated value
$F_{ijk}$ can be computed using the following expression:

$\displaystyle F_{ijk}$ $\displaystyle=\sum_{b\in
O(\vec{r}_{ijk})}F_{b}\frac{m_{b}}{\rho_{b}}W(|\vec{r}_{ijk}-\vec{r}_{b}|,h_{ijk}),$
(2)

where $\rho_{b}$ and $m_{b}$ are the density and the mass of the particle $b$,
$h_{ijk}$ is the “smoothing length” for the grid point $\vec{r}_{ijk}$,
$O(\vec{r}_{ijk})$ is the $h_{ijk}$-neighborhood of $\vec{r}_{ijk}$, and
$W(r,h)$ is the SPH kernel. To apply these formulae, we need to first compute
the smoothing length $h_{ijk}$ for each grid point. Our interpolation
algorithm is therefore performed in two passes over the particles. In the
first pass, we determine the smoothing length on a grid by simply borrowing
the smoothing length of the nearest particle for each grid point. For each
particle, we traverse the grid points within its interaction radius and update
their nearest particle information. Taking Figure 3b as example, the algorithm
first marks L as the nearest particle of grid point G and records the distance
${|G-L|}$. This is later updated to particle H and the distance ${|G-H|}$. The
first pass ensures that valid contributions from distant particles are
accounted for. In the second pass, we use expression (2) with the known
smoothing length on the grid to compute values of the physical quantities.

Even with properly interpolated grid data, we still faced severe numerical
precision issues with the volume renderer in ParaView. Our simulations adopt
the centimetre–gram–second system of units (CGS) which are widely used in the
astrophysics community. Since we are simulating astronomical events, many of
our data values are also astronomical. For example, the $xyz$-coordinates
might range from $-10^{10}$ to $10^{10}$ cm. At any reasonable mesh
resolution, we will have a cell size (called _spacing_ in ParaView) in the
order of $10^{8}$. ParaView seems to have difficulty propagating rays along
such _astronomically large_ cells and shows either nothing or all black pixels
as a result. As suggested by Kitware Inc., we scaled the spacing for the grid
to be in the units of astronomical unit (AU) which seems to resolve the
numerical precision issue.

|
---|---
a | b
Figure 3: Left: ParaView SPH Volume Interpolator pipeline. Right: An example
for uniform grid interpolation of data from SPH particles with varying
smoothing radius, i.e. size, represented by circles.

###  3.2 Physics-based transfer function

The main goal of this project is to answer the question: _what would a
kilonova realistically look like up close?_ To answer this, we need an opacity
transfer function that is physics-based to apply in our volume rendering. In
general, both the color and opacity transfer functions depend on the electron
fraction $Y_{e}$. Regions with low $Y_{e}$ correspond to a neutron-rich
environment where the heaviest elements, including lanthanides and actinides,
are synthesized. Due to half-filled atomic $f$-orbitals, the latter possess an
extremely complex level structure with tens of thousands of energy levels and
up to a few million lines in the optical bands. As a consequence, ejecta that
contains lanthanides and actinides is associated with a dense line blanketing
of the bound-bound opacity in the optical range of the spectrum [13].
Lanthanide-free material, on the other hand, is less opaque in the optical and
appears as “blue” to the observers. With that, it is convenient to adopt a
red-blue colormap on $Y_{e}$ and use red color to indicate more neutron-rich
material where lanthanides are present, and blue color to show less neutron-
rich and lanthanide-free matter. We map the physical opacity to the
absorptance of the grid cells using the following formula:

$\displaystyle A=1-e^{-\tau},$ (3)

where $\tau$ is the optical depth across the cell, estimated as
$\tau=\varkappa\ \rho\ \Delta x.$ Here, $\rho$ is the density in the grid
cell, $\Delta x$ is its size, and $\varkappa$ is the material opacity. In this
work, we use $\varkappa=10\ \textrm{cm}^{2}\ \textrm{g}^{-1}$ for lanthanide-
rich ejecta, and $\varkappa=1\ \textrm{cm}^{2}\ \textrm{g}^{-1}$ for the
lighter $r$-process, lanthanide-free ejecta, which is a simple approximation
that has been successfully applied for constructing kilonova models in the
past [28]. An electron fraction of $Y_{e}=0.25$ is used as threshold value
between lanthanide-rich and lanthanide-free ejecta, similar to previous
studies [23]. Figure 4 shows the effectiveness of such an opacity transfer
function to reveal the complex structure of lanthanide-rich material
responsible for the blue kilonova, in comparison with the naïve use of density
in that role. Because density spans several orders of magnitude, when its
normalized value is used as the opacity transfer function, it looks
significantly more transparent and redder than expected.

Figure 4: Volume rendering of the electron fraction using two different
opacity transfer functions: simple density (left) vs the physics-based
absorptance (right) from Eq.(3). On the left, the opacity is strongly
underestimated, creating an image of the ejecta that is almost completely
transparent. The right half shows correct physical picture for this epoch
(2.23 days) with complex structures in the electron fraction.

##  4 Results

|
---|---
Figure 5:  Left: kilonova rendering at $t=1$ day. Right: same at $t=20$ days.
Opacity transfer function is computed according to Eq.(3), while the color
maps out distribution of the electron fraction.

Our most realistic visualization of the kilonova is given in Figure 5. The
left panel shows the ejecta at the epoch $t=1$ day, while the right panel
illustrates the view at $t=20$ days.

At day 1, we can see various ”blob-like” structures of lanthanide-rich (red)
material. They are mostly located around the central equatorial region of the
kilonova but also reach into the polar lobes. These blobs are seeded by
turbulence in the accretion disk from the magneto-rotational instabilities.
The inhomogeneities are then exaggerated in the disk wind. Even in artificial
colors and the presence of red material blobs, the kilonova looks blue at this
point in time. By day 20 the ejecta expanded by one order of magnitude in size
and we can see a significant color change from blue to red. As the kilonova
expands, a photospheric recession advances the photosphere into the core
uncovering more lanthanide-rich “red” material. Our visualization thus finds a
rich structure of regions where the lighter $r$-process takes place. For once,
there no is simple concentration e.g. only in the polar lobes; instead regions
span a band in the middle latitudes, which would show up in observations as a
“blue” kilonova for some orientations.

##  5 Conclusion

We have produced, for the first time, a physically realistic rendering of a
kilonova, resulting from the coalescence of a neutron star binary. The realism
of our model is based on the fact that we used a state-of-the art model for
the accretion disk wind, realistic nucleosynthesis, plausible models for
physical opacities, and correctly expanded hydrodynamical flow. We use flux
compactification to interface the grid-based hydrodynamics with the particle-
based one (Section 2). For visualization, the expanded flow is interpolated on
the grid, and the inhomogeneous physical opacity is mapped to the absorptance
using Eq.(3). The final result, a synthesized image of the disk wind ejecta,
illustrates the complexity of the kilonova structure. The nature of the rich
features and their influence on the observed kilonova spectra will need
further investigation. Current astrophysical models use simplistic assumptions
such as spherical symmetry and single- or two-component morphology. Our work
clearly demonstrates how little may actually be captured by such simple
models. Future work may include more advanced mapping between the electron
fraction and the opacity ($\kappa(Y_{e})$). Furthermore, the ejecta and
associated kilonova are expected to differ depending on the final outcome of
the NSM. Depending on the neutron stars’ masses, the result of the merger can
be a stable neutron star, a short-lived hypermassive neutron star, or a black
hole. Visualization and analysis of such ejecta data and kilonova for
different merger configurations would be an interesting future undertaking.

During our work, we encountered several issues when visualizing SPH particle
data with highly variable smoothing length in ParaView. Tracking down the
origin of the problem, we give a suggestion for mapping SPH data onto a grid
in order to take full advantage of Paraview’s volume rendering capabilities.

## Acknowledgments

The authors thank Ryan T. Wollaeger and Soumi De for their helpful feedback.
Funding for this project was provided by LANL Ristra Project and Laboratory
Directed Research and Development (LDRD). Alexandra Stewart would like to
thank LANL’s Data Science at Scale Summer School for the opprotunity to work
on this project. This work used resources provided by the LANL Institutional
Computing Program. LANL is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of the U.S. DOE (Contract No.
89233218CNA000001). This work is authorized for unlimited release under LA-
UR-21-28322.

## References

  * [1] B. P. Abbott, et al., GW170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett.119 (2017) 161101.  doi:10.1103/PhysRevLett.119.161101.
  * [2] A. Goldstein, et al., An ordinary short gamma-ray burst with extraordinary implications: Fermi-GBM detection of GRB 170817a, ApJ848 (2) (2017) L14.  doi:10.3847/2041-8213/aa8f41.
  * [3] V. Savchenko, et al., INTEGRAL detection of the first prompt gamma-ray signal coincident with the gravitational-wave event GW170817, ApJ848 (2) (2017) L15.  doi:10.3847/2041-8213/aa8f94.
  * [4] E. Troja, et al., The X-ray counterpart to the gravitational-wave event GW170817, Nature551 (7678) (2017) 71–74.  doi:10.1038/nature24290.
  * [5] K. D. Alexander, et al., The electromagnetic counterpart of the binary neutron star merger LIGO/virgo GW170817. VI. radio constraints on a relativistic jet and predictions for late-time emission from the kilonova ejecta, ApJ848 (2) (2017) L21.  doi:10.3847/2041-8213/aa905d.
  * [6] J. M. Lattimer, D. N. Schramm, Black-Hole-Neutron-Star Collisions, ApJL192 (1974) L145.  doi:10.1086/181612.
  * [7] B. P. Abbott, et al., Multi-messenger observations of a binary neutron star merger, ApJ848 (2) (2017) L12.  doi:10.3847/2041-8213/aa91c9.
  * [8] B. P. Abbott, et al., Gravitational waves and gamma-rays from a binary neutron star merger: GW170817 and GRB 170817a, ApJ848 (2) (2017) L13.  doi:10.3847/2041-8213/aa920c.
  * [9] B. D. Metzger, Kilonovae, Living Reviews in Relativity 23 (1) (2019) 1.  doi:10.1007/s41114-019-0024-0.
  * [10] O. Korobkin, et al., Axisymmetric Radiative Transfer Models of Kilonovae, ApJ910 (2) (2021) 116.  doi:10.3847/1538-4357/abe1b5.
  * [11] R. Chornock, et al., The electromagnetic counterpart of the binary neutron star merger LIGO/virgo GW170817. IV. detection of near-infrared signatures of r-process nucleosynthesis with gemini-south, ApJ848 (2) (2017) L19.  doi:10.3847/2041-8213/aa905c.
  * [12] M. Nicholl, et al., The electromagnetic counterpart of the binary neutron star merger LIGO/virgo GW170817. III. optical and UV spectra of a blue kilonova from fast polar ejecta, ApJ848 (2) (2017) L18.  doi:10.3847/2041-8213/aa9029.
  * [13] C. J. Fontes, C. L. Fryer, A. L. Hungerford, R. T. Wollaeger, O. Korobkin, A line-binned treatment of opacities for the spectra and light curves from neutron star mergers, MNRAS493 (3) (2020) 4143–4171.  doi:10.1093/mnras/staa485.
  * [14] S. Rosswog, O. Korobkin, A. Arcones, F. K. Thielemann, T. Piran, The long-term evolution of neutron star merger remnants - I. The impact of r-process nucleosynthesis, MNRAS439 (1) (2014) 744–756.  doi:10.1093/mnras/stt2502.
  * [15] J. M. Miller, B. R. Ryan, J. C. Dolence, $\nu$bhlight: Radiation GRMHD for neutrino-driven accretion flows, ApJS241 (2) (2019) 30.  doi:10.3847/1538-4365/ab09fc.
  * [16] J. M. Miller, et al., Full transport model of GW170817-like disk produces a blue kilonova, Phys. Rev. D100 (2019) 023008.  doi:10.1103/PhysRevD.100.023008.
  * [17] J. Loiseau, H. Lim, M. A. Kaltenborn, O. Korobkin, C. M. Mauney, I. Sagert, W. P. Even, B. K. Bergen, FleCSPH: The next generation FleCSIble parallel computational infrastructure for smoothed particle hydrodynamics, SoftwareX 12 (2020) 100602.  doi:https://doi.org/10.1016/j.softx.2020.100602.
  * [18] S. Rosswog, SPH Methods in the Modelling of Compact Objects, Living Reviews in Computational Astrophysics 1 (1) (2015) 1.  doi:10.1007/lrca-2015-1.
  * [19] B. K. Bergen, M. R. J. Charest, Flexible computational science infrastructure (FleCSI): Overview & application progress, Tech. rep., Los Alamos National Lab.(LANL), Los Alamos, NM (United States) (2017).
  * [20] O. Korobkin, H. Lim, I. Sagert, J. Loiseau, C. Mauney, M. A. R. Kaltenborn, B.-J. Tsao, W. P. Even, Conservation of Angular Momentum in the Fast Multipole Method, in: Proc. of 2021 international SPHERIC workshop (virtual), New Jersey Institute of Technology, 2021, pp. 112–119.  arXiv:2107.07166.
  * [21] B.-J. Tsao, I. Sagert, O. Korobkin, I. Tews, H. Lim, G. Dilts, J. Loiseau, Modeling Neutron Star Oscillations in a Fixed General Relativistic Background Including Solid Crust Dynamics, in: Proc. of 2021 international SPHERIC workshop (virtual), New Jersey Institute of Technology, 2021, pp. 297–304.  arXiv:2107.07315.
  * [22] J. Heinzel, et al., Comparing inclination-dependent analyses of kilonova transients, MNRAS502 (2) (2021) 3057–3065.  doi:10.1093/mnras/stab221.
  * [23] O. Korobkin, S. Rosswog, A. Arcones, C. Winteler, On the astrophysical robustness of the neutron star merger r-process, MNRAS426 (3) (2012) 1940–1949.  doi:10.1111/j.1365-2966.2012.21859.x.
  * [24] F. X. Timmes, D. Arnett, The Accuracy, Consistency, and Speed of Five Equations of State for Stellar Hydrodynamics, ApJS125 (1) (1999) 277–294.  doi:10.1086/313271.
  * [25] J. Lippuner, L. F. Roberts, SkyNet: A modular nuclear reaction network library, ApJS233 (2) (2017) 18.  doi:10.3847/1538-4365/aa94cb.
  * [26] M. R. Mumpower, et al., $\beta$-delayed fission in r-process nucleosynthesis, ApJ869 (1) (2018) 14.  doi:10.3847/1538-4357/aaeaca.
  * [27] K. Farouqi, et al., Charged-particle and neutron-capture processes in the high-entropy wind of core-collapse supernovae, ApJ712 (2) (2010) 1359–1377.  doi:10.1088/0004-637x/712/2/1359.
  * [28] D. Grossman, O. Korobkin, S. Rosswog, T. Piran, The long-term evolution of neutron star merger remnants - II. Radioactively powered transients, MNRAS439 (1) (2014) 757–770.  arXiv:1307.2943, doi:10.1093/mnras/stt2503.


# On Lamperti transformation and characterisations of discrete random fields

Marko Voutilainen111Turku School of Economics, Department of Accounting and
Finance, FI-20014 University of Turku, Finland, Lauri Viitasaari222Uppsala
University, Department of Mathematics, Box 480, 751 06 Uppsala, Sweden,
Pauliina Ilmonen333Aalto University, Department of Mathematics and Systems
Analysis, P.O. Box 11100, FI-00076 Aalto, Finland

###### Abstract

In this article we characterise discrete time stationary fields by difference
equations involving stationary increment fields and self-similar fields. This
gives connections between stationary fields, stationary increment fields and,
through Lamperti transformation, self-similar fields. Our contribution is a
natural generalisation of recently proved results covering the case of
stationary processes.

AMS 2010 Mathematics Subject Classification: 60G60, 60G10, 60G18

Keywords: random fields, stationary fields, self-similar fields, Lamperti
transformation, fractional Ornstein-Uhlenbeck fields

##  1 Introduction

Stationary processes $X=(X_{t})_{t\in T}$ have numerous applications in many
different fields, and they are a topic of active research. Similarly, self-
similar processes and stationary increment processes have many applications in
various disciplines of science. For details on self-similar processes, we
refer to the monograph [4] and the references therein.

All of these three classes are intimately connected. Indeed, it was already
observed by Lamperti in [8] that there exists a one-to-one correspondence
between stationary processes and self-similar processes. Later on, this
connection was used in [13] to obtain relation between stationary processes
and stationary increment processes in continuous time, i.e. $T=\mathbb{R}$,
through Langevin equation

$dX_{t}=-\theta X_{t}dt+dG_{t},$ (1)

where $\theta>0$ is a parameter, $X$ is stationary, and $G$ has stationary
increments (along with certain other properties). Most notably, this gives
rise to the well-known Ornstein-Uhlenbeck process when one plugs in $G=W$, the
Brownian motion. Connection (1) was later extended for discrete time processes
$X$, i.e. $T=\mathbb{Z}$, in [15], where the authors proved discrete analogue

$\Delta X_{t}=-\theta X_{t-1}+\Delta G_{t}$ (2)

of (1), and studied the estimation of the unknown parameter $\theta$. A
vector-valued version was later provided in [16] and [14], covering both
continuous and discrete time cases together with estimation of the unknown
parameter matrix.

Similarly to stationary processes, stationary fields form an important
subclass of random objects. In this case, $X=(X_{t})_{t\in T}$ with
$T=\mathbb{R}^{N}$ in continuous time or $T=\mathbb{Z}^{N}$ in discrete time
(naturally, $T$ can be a more general parameter space). Also, self-similarity
and stationarity of the increments are wanted features in many applications.
However, while the notion of stationarity is essentially unchanged in the
context of random fields, the notion of self-similarity and stationarity of
the increments become more complicated when $t$ is multidimensional. For
notion of self-similarity for fields, one typically introduces componentwise
self-similarity and considers $H$-self-similarity with $H$ as an
$N$-dimensional vector of componentwise self-similarity indices, see e.g. [11,
5]. In [5] a version of the Lamperti theorem was proved (in continuous time)
for fields, providing a connection between stationary fields and self-similar
fields. For other notions of self-similarity for fields, see for example [3]
and [1].

The notion of stationary increments for fields is even more complicated due to
the fact that the definition of increment is not obvious. One approach is to
consider _rectangular increments_ , where increments are taken over
$N$-dimensional rectangulars. Gaussian self-similar fields and Gaussian
rectangular increment fields were studied, again in continuous time, e.g. in
[10] and [9].

In this article we extend the Characterisation (2) provided in [15] to
discrete time fields, providing a connection between stationary fields, self-
similar fields, and stationary increment fields. More precisely, we provide a
characterisation of the type (see Theorem 2.20)

$X_{t}=\langle\hat{\Theta},\hat{X}^{-}_{t}\rangle+\Delta_{t}G.$

Here $\hat{\Theta}$ is a vector of parameters, and $\hat{X}_{t}^{-}$ is a
vector of ”previous value” consisting of previous values in different
coordinate directions. Our notion of increment $\Delta_{t}G$ corresponds to
the notion of stationary rectangular increments of [10], cf. Remark 2.9. As
such and exactly as in [13, 15] in the case of processes, we obtain
correspondence between stationary fields, stationary rectangular increment
fields, and self-similar fields.

The rest of the article is organised as follows. In Section 2 we introduce and
prove our main results. We introduce our notation and main definitions in
Section 2.1, while our main results and their proofs are presented in Section
2.2. In Section 2.3 we briefly illustrate how our characterisation can be used
to construct discrete time fractional Ornstein-Uhlenbeck fields, extending
notions of (generalized) Ornstein-Uhlenbeck processes of [13]. We end the
paper with conclusions, Section 3, describing future directions, in particular
to cover continuous time parameter space and statistical inference.

##  2 Connections between stationary, self-similar, and stationary increment
fields

###  2.1 Preliminaries and notations

We begin with by introducing some definitions and notations.

######  Definition 2.1 (Stationarity).

A random field $X=(X_{t})_{t\in\mathbb{Z}^{N}}$ is stationary if

$(X_{t+s})_{t\in\mathbb{Z}^{N}}\overset{\text{law}}{=}(X_{t})_{t\in\mathbb{Z}^{N}}$

for every $s\in\mathbb{Z}^{N}$ in the sense of finite dimensional
distributions.

######  Definition 2.2 (Self-similarity).

Let
$Y=(Y_{e^{t}})_{t\in\mathbb{Z}^{N}}=(Y_{e^{t_{1}},\dots,e^{t_{N}}})_{t\in\mathbb{Z}^{N}}$
be a random field. In addition, let
$\Theta=(\theta_{1},\dots,\theta_{N})\in(0,\infty)^{N}$ be a positive multi-
index. If

$(Y_{e^{t+s}})_{t\in\mathbb{Z}^{N}}\overset{\text{law}}{=}(e^{\langle
s,\Theta\rangle}Y_{e^{t}})_{t\in\mathbb{Z}^{N}}$

for every $s\in\mathbb{Z}^{N}$, where $\langle s,\Theta\rangle$ is the
standard inner product of vectors, then $Y$ is a $\Theta$-self-similar random
field.

######  Remark 2.3.

The exponential terms in Definition 2.2 are introduced in order to take into
account the discrete nature of the field. Definition 2.2 is analogous to the
classical definition in continuous time. See e.g. [8] for the definition in
the one parameter continuous case and [13] for the definition in the one
parameter discrete case.

The following definition provides a notion of Lamperti transformation in our
setting.

######  Definition 2.4 (Lamperti).

Let $\Theta=(\theta_{1},\dots,\theta_{N})\in(0,\infty)^{N}$. The Lamperti
transformation $\mathcal{L}_{\Theta}$ and its inverse
$\mathcal{L}^{-1}_{\Theta}$ for discrete random fields are defined by

$\displaystyle(\mathcal{L}_{\Theta}X)_{e^{t}}$ $\displaystyle=e^{\langle
t,\Theta\rangle}X_{t},\quad t\in\mathbb{Z}^{N},$
$\displaystyle(\mathcal{L}^{-1}_{\Theta}Y)_{t}$ $\displaystyle=e^{-\langle
t,\Theta\rangle}Y_{e^{t}},\quad t\in\mathbb{Z}^{N}.$

######  Remark 2.5.

The formulae in Definitions 2.2 and 2.4 differ slightly from the standard
forms in continuous settings. In the case of random fields, the definition of
component-wise self-similarity and the corresponding Lamperti transformation
together with a one-to-one correspondence between self-similar and stationary
fields was presented in [5]. In comparison, our definitions are obtained via
change of variables from the standard ones, ensuring that we stay within our
discrete parameter set when applying the transformation.

######  Definition 2.6 (Increments).

The square increment $\Delta_{t}X$ of a field $(X_{t})_{t\in\mathbb{Z}^{N}}$
at a point $t=(t_{1},\dots,t_{N})$ $\in\mathbb{Z}^{N}$ is given by

$\Delta_{t}X=\sum_{(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}}(-1)^{\sum_{l=1}^{N}i_{l}}X_{t_{1}-i_{1},\dots,t_{N}-i_{N}}.$
(3)

######  Example 2.7.

In the particular case $N=2$, we have

$\Delta_{t}X=X_{t_{1},t_{2}}+X_{t_{1}-1,t_{2}-1}-X_{t_{1}-1,t_{2}}-X_{t_{1},t_{2}-1}.$

######  Definition 2.8 (Stationary increment field).

A field $X=(X_{t})_{t\in\mathbb{Z}^{N}}$ has stationary increments if the
increment field $(\Delta_{t}X)_{t\in\mathbb{Z}^{N}}$ is stationary. That is

$(\Delta_{t+s}X)_{t\in\mathbb{Z}^{N}}\overset{\text{law}}{=}(\Delta_{t}X)_{t\in\mathbb{Z}^{N}}$

for every $s\in\mathbb{Z}^{N}$ in the sense of finite dimensional
distributions.

######  Remark 2.9.

The authors in [10] introduced a continuous time analogous notion of _strictly
stationary rectangular increments_ by assuming stationarity of the increments
over arbitrary rectangular increments. In comparison, in our definition, we
consider stationary increments over unit rectangulars. However, due to the
discrete nature of our index space $t\in\mathbb{Z}^{N}$, one can show that our
definition is equivalent to assuming stationarity over arbitrary (discrete)
rectangular increments.

We also need a notion of _previous value_ that is not so straightforward in a
multi-parameter setting.

######  Definition 2.10 (Previous value).

The previous value of the field $X=(X_{t})_{t\in\mathbb{Z}^{N}}$ at a point
$t=(t_{1},\dots,t_{N})\in\mathbb{Z}^{N}$ is given by

$X_{t}^{-}=X_{t}-\Delta_{t}X=\sum_{\begin{subarray}{c}(i_{1},\ldots,i_{N})\in\\{0,1\\}^{N}\\\
(i_{1},\ldots,i_{N})\neq\bf{0}\end{subarray}}(-1)^{1+\sum_{l=1}^{N}i_{l}}X_{t_{1}-i_{1},\dots,t_{N}-i_{N}}.$

######  Remark 2.11.

Note that in our definition of the previous value, we take into account the
terms in (3) that have a smaller index in at least one of the coordinate
directions. Indeed, in the two-dimensional case, we have

$\Delta_{t}X=X_{t_{1},t_{2}}+X_{t_{1}-1,t_{2}-1}-X_{t_{1}-1,t_{2}}-X_{t_{1},t_{2}-1}$

while the previous value is given by

$X^{-}_{t}=X_{t_{1}-1,t_{2}}+X_{t_{1},t_{2}-1}-X_{t_{1}-1,t_{2}-1}.$

######  Definition 2.12 (Inner product
$\langle\hat{\Theta},\hat{X}^{-}_{t}\rangle$).

Let $\Theta=(\theta_{1},\dots,\theta_{N})\in(0,\infty)^{N}$ and
$X=(X_{t})_{t\in\mathbb{Z}^{N}}$ be a field. We define vectors $\hat{\Theta}$
and $\hat{X}^{-}_{t}$ of length $2^{N}-1$ having elements of the forms

$(-1)^{1+\sum_{l=1}^{N}i_{l}}e^{-\langle i,\Theta\rangle}\quad\text{and}\quad
X_{{t_{1}-i_{1}},\dots,{t_{N}-i_{N}}},$

respectively, where $i=(i_{1},\dots,i_{N})\in\\{0,1\\}^{N},i\neq\bf{0}.$ Then
the inner product of the vectors is

$\langle\hat{\Theta},\hat{X}^{-}_{t}\rangle=\sum_{\begin{subarray}{c}(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}\\\
i\neq\bf{0}\end{subarray}}(-1)^{1+\sum_{l=1}^{N}i_{l}}e^{-\langle
i,\Theta\rangle}X_{{t_{1}-i_{1}},\dots,{t_{N}-i_{N}}}.$

######  Example 2.13.

In the two-dimensional case we obtain that, for
$\Theta=(\theta_{1},\theta_{2})\in(0,\infty)^{2}$ and
$X=(X_{t})_{t\in\mathbb{Z}^{2}}$, the vectors $\hat{\Theta}$ and
$\hat{X}^{-}_{t}$ are given as

$\hat{\Theta}=(e^{-\theta_{1}},e^{-\theta_{2}},-e^{-\theta_{1}-\theta_{2}})\quad\text{and}\quad\hat{X}_{t}^{-}=\begin{pmatrix}X_{t_{1}-1,t_{2}}\\\
X_{t_{1},t_{2}-1}\\\ X_{t_{1}-1,t_{2}-1}\end{pmatrix}$

and we have

$\langle\hat{\Theta},\hat{X}^{-}_{t}\rangle=e^{-\theta_{1}}X_{t_{1}-1,t_{2}}+e^{-\theta_{2}}X_{t_{1},t_{2}-1}-e^{-\theta_{1}-\theta_{2}}X_{t_{1}-1,t_{2}-1}.$
(4)

######  Remark 2.14.

The vectors $\hat{\Theta}$ and $\hat{X}^{-}_{t}$ are unique up to permutations
of their elements.

######  Definition 2.15 (Class $\mathcal{G}_{{\Theta}}$).

Let $\Theta=(\theta_{1},\dots,\theta_{N})\in(0,\infty)^{N}$. Let
$G=(G_{t})_{t\in\mathbb{Z}^{N}}$ be a stationary increment field with
$G_{t}=0$ for all $t$ such that $\sum_{i=1}^{N}t_{i}\in\\{0,-1,\dots,-N+1\\}$.
If

$\lim_{M_{1}\to\infty}\dots\lim_{M_{N}\to\infty}\sum_{j_{1}=-M_{1}}^{t_{1}}\dots\sum_{j_{N}=-M_{N}}^{t_{N}}e^{\sum_{l=1}^{N}j_{l}\theta_{l}}\Delta_{(j_{1},\dots,j_{N})}G$
(5)

converges in probability defining an almost surely finite random variable for
every $t=(t_{1},\dots,t_{N})\in\mathbb{Z}^{N}$, then
$G\in\mathcal{G}_{{\Theta}}$.

######  Remark 2.16.

It turns out that the order of the limits in (5) is irrelevant, and any
permutation of the order leads to the convergence towards the same limiting
random variable that turns out to be the stationary field $X_{t}$, cf. proof
of Theorem 2.28.

######  Remark 2.17.

Condition $G_{t}=0$ for all $t$ such that
$\sum_{i=1}^{N}t_{i}\in\\{0,-1,\dots,-N+1\\}$ is rather peculiar and it
essentially means that $G$ has to vanish along discrete points of certain
$N-1$-dimensional planes. However, this condition is not required _a priori_
for the characterisation, but turns out to hold true and is also required to
obtain uniqueness of the representation, cf. Theorem 2.20 and Remark 2.29.

######  Remark 2.18.

Note that if $G\in L^{1}$ (or $\Delta G\in L^{1}$), then
$G\in\mathcal{G}_{{\Theta}}$ for every $\Theta$. Indeed, this can be seen from

$\begin{split}&\mathbb{E}|\sum_{j_{1}=-\infty}^{t_{1}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=1}^{N}j_{l}\theta_{l}}\Delta_{(j_{1},\dots,j_{N})}G|\\\
&\leq\sum_{j_{1}=-\infty}^{t_{1}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=1}^{N}j_{l}\theta_{l}}\mathbb{E}|\Delta_{(j_{1},\dots,j_{N})}G|\\\
&=\sum_{j_{1}=-\infty}^{t_{1}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=1}^{N}j_{l}\theta_{l}}\mathbb{E}|\Delta_{(1,\dots,1)}G|\\\
&\leq\sum_{j_{1}=-\infty}^{t_{1}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=1}^{N}j_{l}\theta_{l}}\sum_{(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}}\mathbb{E}|G_{1-i_{1},\dots,1-i_{N}}|.\end{split}$

######  Example 2.19.

In the two-dimensional case, we have $G\in\mathcal{G}_{{\Theta}}$ for
$\Theta=(\theta_{1},\theta_{2})\in(0,\infty)^{2}$ provided that $G$ has
stationary increments, $G_{t,-t}=G_{t,-t-1}=0$ and

$\lim_{M_{1}\to\infty}\lim_{M_{2}\to\infty}\sum_{j_{1}=-M_{1}}^{t_{1}}\sum_{j_{2}=-M_{2}}^{t_{2}}e^{j_{1}\theta_{1}}e^{j_{2}\theta_{2}}\Delta_{(j_{1},j_{2})}G$

exists as an almost surely finite random variable. That is, $G$ has stationary
increments and $G_{t_{1},t_{2}}$ is set to zero on lines $t_{1}=-t_{2}$ and
$t_{1}=-t_{2}-1$. The existence of such fields follow as a by-product of our
main results.

###  2.2 AR(1) type characterisation of stationary fields

Our main result is the following characterisation that is a natural extension
of the one-dimensional case presented in [15].

######  Theorem 2.20.

Let $\Theta=(\theta_{1},\dots,\theta_{N})\in(0,\infty)^{N}$. A field
$X=(X_{t})_{t\in\mathbb{Z}^{N}}$ is stationary if and only if the following
conditions are satisfied.

  1. (i)

$\lim_{m\to-\infty}e^{m\theta_{j}}X_{t_{1},\dots,t_{j-1},m,t_{j+1},\dots,t_{N}}\overset{\mathbb{P}}{\longrightarrow}0$

for every $j\in\\{1,\dots,N\\}$ and
$t_{1},\dots,t_{j-1},t_{j+1},\dots,t_{N}\in\mathbb{Z}$.

  2. (ii)

There exists $G=(G_{t})_{t\in\mathbb{Z}^{N}}\in\mathcal{G}_{{\Theta}}$ such
that

$X_{t}=\langle\hat{\Theta},\hat{X}_{t}^{-}\rangle+\Delta_{t}G\quad\text{for
every }t\in\mathbb{Z}^{N},$ (6)

where $\langle\hat{\Theta},\hat{X}_{t}^{-}\rangle$ is given by Definition
2.12.

Moreover, for a given $\Theta$, the stationary increment field
$G\in\mathcal{G}_{{\Theta}}$ in (6) is unique.

As a direct corollary we obtain the following version in a two-dimensional
case.

######  Corollary 2.21.

Let $\Theta=(\theta_{1},\theta_{2})\in(0,\infty)^{2}$. A field
$X=(X_{t})_{t\in\mathbb{Z}^{2}}$ is stationary if and only if the following
conditions are satisfied.

  1. (i)

$\lim_{m\to-\infty}e^{m\theta_{2}}X_{t_{1},m}\overset{\mathbb{P}}{\longrightarrow}0\quad\text{and}\quad\lim_{m\to-\infty}e^{m\theta_{1}}X_{m,t_{2}}\overset{\mathbb{P}}{\longrightarrow}0$

for every $t_{1}$ and $t_{2}$.

  2. (ii)

There exists $G=(G_{t})_{t\in\mathbb{Z}^{2}}\in\mathcal{G}_{{\Theta}}$ such
that

$X_{t}=\langle\hat{\Theta},\hat{X}_{t}^{-}\rangle+\Delta_{t}G\quad\text{for
every }t\in\mathbb{Z}^{2},$ (7)

where $\langle\hat{\Theta},\hat{X}_{t}^{-}\rangle$ is given by (4).

The stationary increment field $G\in\mathcal{G}_{{\Theta}}$ in (7) is unique.

The proof of Theorem 2.20 is split into a series of lemmas and auxiliary
theorems. We begin with the following result that is a version of Lamperti
theorem.

######  Theorem 2.22.

If $X=(X_{t})_{t\in\mathbb{Z}^{N}}$ is stationary, then
$(\mathcal{L}_{\Theta}X)_{e^{t}}$ is $\Theta$-self-similar. Conversely, if
$Y=(Y_{e^{t}})_{t\in\mathbb{Z}^{N}}$ is $\Theta$-self-similar, then
$(\mathcal{L}^{-1}_{\Theta}Y)_{t}$ is stationary.

###### Proof.

First, assume that $X$ is stationary. Set
$Y_{e^{t}}=(\mathcal{L}_{\Theta}X)_{e^{t}}$ and let $n\in\mathbb{N}$. Now

$\begin{split}(Y_{e^{t_{1}+s}},\dots,Y_{e^{t_{n}+s}})&=(e^{\langle
t_{1}+s,\Theta\rangle}X_{t_{1}+s},\dots e^{\langle
t_{n}+s,\Theta\rangle}X_{t_{n}+s})\\\ &\overset{\text{law}}{=}(e^{\langle
s,\Theta\rangle}e^{\langle t_{1},\Theta\rangle}X_{t_{1}},\dots,e^{\langle
s,\Theta\rangle}e^{\langle t_{n},\Theta\rangle}X_{t_{n}})\\\ &=(e^{\langle
s,\Theta\rangle}Y_{e^{t_{1}}},\dots,e^{\langle
s,\Theta\rangle}Y_{e^{t_{n}}}),\end{split}$

proving the first part of the claim. Next, assume that $Y$ is $\Theta$-self-
similar. Set $X_{t}=(\mathcal{L}^{-1}_{\Theta}Y)_{t}$ and let
$n\in\mathbb{N}$. Now

$\begin{split}(X_{t_{1}+s},\dots,X_{t_{n}+s})&=(e^{-\langle
t_{1}+s,\Theta\rangle}Y_{e^{t_{1}+s}},\dots,e^{-\langle
t_{n}+s,\Theta\rangle}Y_{e^{t_{n}+s}})\\\ &\overset{\text{law}}{=}(e^{-\langle
t_{1},\Theta\rangle}Y_{e^{t_{1}}},\dots,e^{-\langle
t_{n},\Theta\rangle}Y_{e^{t_{n}}})\\\
&=(X_{t_{1}},\dots,X_{t_{n}}),\end{split}$

completing the proof. ∎

The following lemma provides one of our key observations.

######  Lemma 2.23.

Let $(Y_{e^{t}})_{t\in\mathbb{Z}^{N}}$ be $\Theta$-self-similar. Set

$\Delta_{t}Y=\sum_{(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}}(-1)^{\sum_{l=1}^{N}i_{l}}Y_{e^{t_{1}-i_{1}},\dots,e^{t_{N}-i_{N}}}.$

For $\sum_{l=1}^{N}t_{l}\geq 1$, we set

$G_{t}=\sum_{k_{1}=1-t_{2}-\dots-
t_{N}}^{t_{1}}\sum_{k_{2}=1-k_{1}-t_{3}-\dots-
t_{N}}^{t_{2}}\dots\sum_{k_{N}=1-k_{1}-\dots-k_{N-1}}^{t_{N}}e^{-\langle
k,\Theta\rangle}\Delta_{k}Y,$

and, for $\sum_{l=1}^{N}t_{l}\leq 0$, we set

$G_{t}=(-1)^{N}\sum_{k_{1}=t_{1}+1}^{-t_{2}-\dots-
t_{N}-N+1}\sum_{k_{2}=t_{2}+1}^{-k_{1}-t_{3}-\dots-
t_{N}-N+2}\dots\sum_{k_{N}=t_{N}+1}^{-k_{1}-\dots-k_{N-1}}e^{-\langle
k,\Theta\rangle}\Delta_{k}Y.$

Here $\langle k,\Theta\rangle$ is the standard inner product and sums of the
type $\sum_{s_{2}}^{s_{1}}$ with $s_{1}<s_{2}$ are interpreted as empty sums.
Now

  1. (i)

$G_{t}=0$ for all $t$ such that $\sum_{l=1}^{N}t_{l}\in\\{0,-1,\dots,-N+1\\}$,

  2. (ii)

$\Delta_{t}G=e^{-\langle t,\Theta\rangle}\Delta_{t}Y$ for every
$t\in\mathbb{Z}^{N}$,

  3. (iii)

$G=(G_{t})_{t\in\mathbb{Z}^{N}}$ is a stationary increment field.

######  Remark 2.24.

It turns out that $G$ defined as above satisfies $G\in\mathcal{G}_{\Theta}$,
see also Lemma 2.31 below.

######  Example 2.25.

In the two-dimensional case, for $\Theta$-self-similar
$(Y_{e^{t}})_{t\in\mathbb{Z}^{2}}$, we denote

$\Delta_{t}Y=Y_{e^{t_{1}},e^{t_{2}}}-Y_{e^{t_{1}-1},e^{t_{2}}}-Y_{e^{t_{1}},e^{t_{2}-1}}+Y_{e^{t_{1}-1},e^{t_{2}-1}}.$

The field $G=(G_{t})_{t\in\mathbb{Z}^{2}}$ defined as

$G_{t_{1},t_{2}}=\begin{cases}\sum_{k_{1}=1-t_{2}}^{t_{1}}\sum_{k_{2}=1-k_{1}}^{t_{2}}e^{-\langle(k_{1},k_{2}),\Theta\rangle}\Delta_{k}Y,&\quad
t_{1}+t_{2}\geq 1\\\
\sum_{k_{1}=t_{1}+1}^{-t_{2}-1}\sum_{k_{2}=t_{2}+1}^{-k_{1}}e^{-\langle(k_{1},k_{2}),\Theta\rangle}\Delta_{k}Y,&\quad
t_{1}+t_{2}\leq 0,\end{cases}$

belongs to the class $\mathcal{G}_{\Theta}$. Here sums of the type
$\sum_{s+1}^{s}$ are interpreted as empty sums.

The proof of Lemma 2.23 is based on the following additional lemmas that we
prove first. The first one provides an auxiliary result on sums of binomial
coefficients. Although the result is quite elementary, we provide a proof for
the reader’s convenience.

######  Lemma 2.26.

We have the following identities:

$\sum_{m=0}^{\frac{M-1}{2}}\binom{M}{2m}=2^{M-1},\quad\sum_{m=0}^{\frac{M-1}{2}}\binom{M}{2m+1}=2^{M-1},\qquad\text{when
$M\geq 1$ is odd}.$
$\sum_{m=0}^{\frac{M}{2}}\binom{M}{2m}=2^{M-1},\quad\sum_{m=0}^{\frac{M}{2}-1}\binom{M}{2m+1}=2^{M-1},\qquad\text{when
$M\geq 2$ is even}.$

###### Proof.

In the odd case

$\sum_{m=0}^{\frac{M-1}{2}}\binom{M}{2m}=\binom{M}{0}+\ldots+\binom{M}{M-1}\quad\text{and}\quad\sum_{m=0}^{\frac{M-1}{2}}\binom{M}{2m+1}=\binom{M}{1}+\ldots+\binom{M}{M}.$

These sums are equal since $\binom{M}{k}=\binom{M}{M-k}$. In addition,

$\sum_{m=0}^{\frac{M-1}{2}}\binom{M}{2m}+\sum_{m=0}^{\frac{M-1}{2}}\binom{M}{2m+1}=2^{M}$

completing the proof of the first case. For the even case, we obtain

$\begin{split}\sum_{m=0}^{\frac{M}{2}-1}\binom{M}{2m+1}&=\binom{M}{1}+\ldots+\binom{M}{M-1}\\\
&=\binom{M-1}{0}+\binom{M-1}{1}+\dots+\binom{M-1}{M-2}+\binom{M-1}{M-1}=2^{M-1}.\end{split}$

Observing that

$\sum_{m=0}^{\frac{M}{2}}\binom{M}{2m}+\sum_{m=0}^{\frac{M}{2}-1}\binom{M}{2m+1}=2^{M}$

completes the proof. ∎

The following lemma sheds light on how the field $G$ of Lemma 2.23 is
constructed from a self-similar field $Y$.

######  Lemma 2.27.

  1. (i)

Let $\sum_{l=1}^{N}t_{l}\geq 1$. Then a term $e^{-\langle
j,\Theta\rangle}\Delta_{j}Y$ belongs to the sum defining $G_{t}$ in Lemma 2.23
if and only if

$j_{l}\leq t_{l}\quad\text{for every $l$ and}\quad\sum_{l=1}^{N}j_{l}\geq 1.$

  2. (ii)

Let $\sum_{l=1}^{N}t_{l}\leq-N$. Then a term $e^{-\langle
j,\Theta\rangle}\Delta_{j}Y$ belongs to the sum defining $G_{t}$ in Lemma 2.23
if and only if

$j_{l}\geq t_{l}+1\quad\text{for every $l$ and}\quad\sum_{l=1}^{N}j_{l}\leq
0.$

###### Proof.

Item $(i)$: By the upper bounds in the sum defining $G_{t}$, it is clear that
we have $j_{l}\leq t_{l}$ for all $l$. By the lower bound of the inmost
summation, we obtain $j_{N}\geq 1-j_{1}-\dots-j_{N-1}$. That is,
$\sum_{l=1}^{N}j_{l}\geq 1.$ We also observe that the other lower bounds of
the sum defining $G_{t}$ yield conditions

$j_{N-h}\geq 1-\sum_{l=1}^{N-h-1}j_{l}-\sum_{l=N-h+1}^{N}t_{l}\quad\text{for
every }h\in\\{0,\dots,N-1\\}.$

These conditions are satisfied since

$\sum_{l=1}^{N-h}j_{l}+\sum_{N-h+1}^{N}t_{l}\geq\sum_{l=1}^{N}j_{l}\geq 1.$

This completes the proof of the first item.

Item $(ii)$: By the lower bounds in the sum defining $G_{t}$, it is clear that
we have $j_{l}\geq t_{l}+1$ for all $l$. By the upper bound of the inmost
summation, we obtain also that $j_{N}\leq-j_{1}-\dots-j_{N-1}$. That is,
$\sum_{l=1}^{N}j_{l}\leq 0.$ We also observe that the other upper bounds of
the sum defining $G_{t}$ yield conditions

$j_{N-h}\leq-\sum_{l=1}^{N-h-1}j_{l}-\sum_{l=N-h+1}^{N}t_{l}-N+(N-h)\quad\text{for
every }h\in\\{0,\dots,N-1\\}.$

These conditions are satisfied since

$\begin{split}\sum_{l=1}^{N-h}j_{l}+\sum_{N-h+1}^{N}t_{l}&\leq\sum_{l=1}^{N-h}j_{l}+\sum_{N-h+1}^{N}(j_{l}-1)=\sum_{l=1}^{N}j_{l}-h\leq-h.\end{split}$

This completes the proof of the second item, and thus the whole proof is
completed. ∎

###### Proof of Lemma 2.23.

Item $(i)$: Let $\sum_{l=1}^{N}t_{l}\in\\{0,-1,\dots,-N+1\\}$. Then

$\sum_{l=2}^{N}t_{l}\in\\{-t_{1},-1-t_{1},\dots,-N+1-t_{1}\\}$

and for the upper bound of the first summation in the definition of $G_{t}$ it
holds that

$-\sum_{l=2}^{N}t_{l}-N+1\in\\{t_{1}-N+1,t_{1}-N+2,\dots,t_{1}\\}.$

Hence, $G_{t}$ is given by an empty sum.

Item $(ii)$: Recall that $G_{t}=0$ for all $t$ such that
$\sum_{i=1}^{N}t_{i}\in\\{0,-1,\dots,-N+1\\}$.

We treat the case $\sum_{l=1}^{N}t_{l}\geq 1$ first. Let $M$ be such that
$\sum_{l=1}^{N}t_{l}-M=1$. Then

$\Delta_{t}G=\sum_{\begin{subarray}{c}(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}\\\
\sum_{l=1}^{N}i_{l}\leq
M\end{subarray}}(-1)^{\sum_{l=1}^{N}i_{l}}G_{t_{1}-i_{1},\dots,t_{N}-i_{N}}.$
(8)

By Lemma 2.27, $\Delta_{t}G$ consists of terms $e^{-\langle
j,\Theta\rangle}\Delta_{j}Y$ with $j_{l}\leq t_{l}$ for every $l$ and
$\sum_{l=1}^{N}j_{l}\geq 1$. Let $m$ be the number of indices $l$ for which
$j_{l}=t_{l}$.
Assume that $m<N$. By Lemma 2.27, $e^{-\langle j,\Theta\rangle}\Delta_{j}Y$
belongs to summands of (8) that satisfy $j_{l}\leq t_{l}-i_{l}$ for every $l$.
That is, $m$ of the indices $i_{l}$ are zero while the remaining $N-m$ indices
may be zeros or ones. In addition,

$1\leq\sum_{l=1}^{N}j_{l}\leq\sum_{l=1}^{N}t_{l}-(N-m)$

giving

$N-m\leq\sum_{l=1}^{N}t_{l}-1=M.$ (9)

Now if $N-m$ is odd, then by Lemma 2.26 and (9), the number of terms
$e^{-\langle j,\Theta\rangle}\Delta_{j}Y$ in (8) with a positive sign is

$\displaystyle\binom{N-m}{0}+\ldots+\binom{N-m}{\min\\{N-m,M\\}-1}=\binom{N-m}{0}+\ldots+\binom{N-m}{N-m-1}=2^{N-m-1}.$

Thus, terms $e^{-\langle j,\Theta\rangle}\Delta_{j}Y$ cancel out in (8).
Similarly, if $N-m$ is even, then the number of terms
$e^{-j\Theta}\Delta_{j}Y$ with a positive sign is

$\binom{N-m}{0}+\dots+\binom{N-m}{\min\\{N-m,M\\}}=\binom{N-m}{0}+\dots+\binom{N-m}{N-m}.$

Again, by Lemma 2.26, terms $e^{-\langle j,\Theta\rangle}\Delta_{j}Y$ cancel
out in (8).
If $m=N$, we have that $e^{-\langle j,\Theta\rangle}\Delta_{j}Y=e^{-\langle
t,\Theta\rangle}\Delta_{t}Y$ belongs only to the summand of (8) with
$i=\bf{0}$. Hence we have shown that

$\Delta_{t}G=e^{-\langle t,\Theta\rangle}\Delta_{t}Y\quad\text{for every $t$
such that }\sum_{l=1}^{N}t_{l}\geq 1.$

This proves the claim for the case $\sum_{l=1}^{N}t_{l}\geq 1$.

Assume next that $\sum_{l=1}^{N}t_{l}\leq 0$ and let $M$ be such that
$\sum_{l=1}^{N}t_{l}-M=-N$. Then

$\Delta_{t}G=\sum_{\begin{subarray}{c}(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}\\\
\sum_{l=1}^{N}i_{l}\geq
M\end{subarray}}(-1)^{\sum_{l=1}^{N}i_{l}}G_{t_{1}-i_{1},\dots,t_{N}-i_{N}}.$
(10)

By Lemma 2.27, $\Delta_{t}G$ consists of terms $e^{-\langle
j,\Theta\rangle}\Delta_{j}Y$ with $j_{l}\geq t_{l}$ for every $l$ and
$\sum_{l=1}^{N}j_{l}\leq 0$. As before, let $m$ be the number of indices $l$
for which $j_{l}=t_{l}$. If $m<N$, then, by Lemma 2.27, $e^{-\langle
j,\Theta\rangle}\Delta_{j}Y$ belongs to summands of (10) that satisfy
$j_{l}\geq t_{l}-i_{l}+1$ for every $l$. That is, $m$ of the indices $i_{l}$
are equal to one while the remaining $N-m$ indices may be zeros or ones. In
addition,

$0\geq\sum_{l=1}^{N}j_{l}\geq\sum_{l=1}^{N}t_{l}+(N-m)=M-m$

giving $m\geq M$ and $N-m\leq N-M$. If $N-m$ is odd, by Lemma 2.26, the number
of terms $e^{-\langle j,\Theta\rangle}\Delta_{j}Y$ in (10) with the sign
$(-1)^{N+m}$ is equal to

$\displaystyle\binom{N-m}{0}+\dots+\binom{N-m}{N-m-1}=2^{N-m-1}.$

Thus, terms cancel out in (10). Similarly, terms cancel out when $N-m$ is
even. Finally, for the case $m=N$ we observe that $e^{-\langle
j,\Theta\rangle}\Delta_{j}Y=e^{-\langle t,\Theta\rangle}\Delta_{t}Y$ belongs
only to the summand of (10) with $i=\bf{1}$. The corresponding sign is
$(-1)^{N}(-1)^{N}=1$.
To conclude, we have shown that

$\Delta_{t}G=e^{-\langle t,\Theta\rangle}\Delta_{t}Y\quad\text{for every $t$
such that }\sum_{l=1}^{N}t_{l}\leq 0.$

This completes the proof of item $(ii)$.

Item $(iii)$: For $s=(s_{1},\dots,s_{N})\in\mathbb{Z}^{N}$,

$\begin{split}\Delta_{t+s}G&=e^{-\langle t+s,\Theta\rangle}\Delta_{t+s}Y\\\
&=e^{-\langle t,\Theta\rangle}e^{-\langle
s,\Theta\rangle}\sum_{(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}}(-1)^{\sum_{l=1}^{N}i_{l}}Y_{e^{t_{1}+s_{1}-i_{1}},\dots,e^{t_{N}+s_{N}-i_{N}}}\\\
&\overset{\text{law}}{=}e^{-\langle
t,\Theta\rangle}\sum_{(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}}(-1)^{\sum_{l=1}^{N}i_{l}}Y_{e^{t_{1}-i_{1}},\dots,e^{t_{N}-i_{N}}}\\\
&=e^{-\langle t,\Theta\rangle}\Delta_{t}Y=\Delta_{t}G.\end{split}$

Treating multidimensional distribution similarly completes the proof of item
$(iii)$. ∎

We are now ready to prove three results, Theorem 2.28, Theorem 2.30, and Lemma
2.32, that give us the main result of this article, Theorem 2.20.

######  Theorem 2.28.

Let $\Theta=(\theta_{1},\dots,\theta_{N})\in(0,\infty)^{N}$ and let
$X=(X_{t})_{t\in\mathbb{Z}^{N}}$ be a random field. If for some
$G=(G_{t})_{t\in\mathbb{Z}^{N}}\in\mathcal{G}_{{\Theta}}$ it holds that

$X_{t}=\langle\hat{\Theta},\hat{X}_{t}^{-}\rangle+\Delta_{t}G\quad\text{for
every }t\in\mathbb{Z}^{N},$ (11)

and

$\lim_{m\to-\infty}e^{m\theta_{j}}X_{t_{1},\dots,t_{j-1},m,t_{j+1},\dots,t_{N}}\overset{\mathbb{P}}{\longrightarrow}0$

for every $j\in\\{1,\dots,N\\}$ and
$t_{1},\dots,t_{j-1},t_{j+1},\dots,t_{N}\in\mathbb{Z}$, then
$X=(X_{t})_{t\in\mathbb{Z}^{N}}$ is stationary.

###### Proof.

Denote $Z_{t}=\Delta_{t}G$ and $i=(i_{1},\ldots,i_{N})$. From (11) and
Definition 2.12, we get

$X_{t}=\sum_{\begin{subarray}{c}(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}\\\
i\neq\bf{0}\end{subarray}}(-1)^{1+\sum_{l=1}^{N}i_{l}}e^{-\langle
i,\Theta\rangle}X_{{t_{1}-i_{1}},\dots,{t_{N}-i_{N}}}+Z_{t},$

which gives

$\begin{split}&X_{t}-\sum_{\begin{subarray}{c}(i_{1},\dots,i_{N-1})\in\\{0,1\\}^{N-1}\\\
i\neq\bf{0}\end{subarray}}(-1)^{1+\sum_{l=1}^{N-1}i_{l}}e^{-\langle(i_{1},\dots,i_{N-1},0),\Theta\rangle}X_{{t_{1}-i_{1}},\dots,t_{N-1}-i_{N-1},{t_{N}}}\\\
&=\sum_{(i_{1},\dots,i_{N-1})\in\\{0,1\\}^{N-1}}(-1)^{\sum_{l=1}^{N-1}i_{l}}e^{-\langle(i_{1},\dots,i_{N-1},0),\Theta\rangle}X_{{t_{1}-i_{1}},\dots,t_{N-1}-i_{N-1},{t_{N}}}\\\
&=\sum_{(i_{1},\dots,i_{N-1})\in\\{0,1\\}^{N-1}}(-1)^{\sum_{l=1}^{N-1}i_{l}}e^{-\langle(i_{1},\dots,i_{N-1},1),\Theta\rangle}X_{{t_{1}-i_{1}},\dots,t_{N-1}-i_{N-1},{t_{N}-1}}+Z_{t}.\end{split}$
(12)

Set

$Y_{t_{1},\dots,t_{N-1}}(t_{N})=\sum_{(i_{1},\dots,i_{N-1})\in\\{0,1\\}^{N-1}}(-1)^{\sum_{l=1}^{N-1}i_{l}}e^{-\langle(i_{1},\dots,i_{N-1},0),\Theta\rangle}X_{{t_{1}-i_{1}},\dots,t_{N-1}-i_{N-1},{t_{N}}}.$

Then, by iterating the recursive Equation (12), we get

$\begin{split}Y_{t_{1},\dots,t_{N-1}}(t_{N})=&e^{-\theta_{N}}Y_{t_{1},\dots,t_{N-1}}(t_{N}-1)+Z_{t}\\\
=&e^{-(n+1)\theta_{N}}Y_{t_{1},\dots,t_{N-1}}(t_{N}-n-1)+\sum_{j_{N}=0}^{n}e^{-j_{N}\theta_{N}}Z_{t_{1},\dots,t_{N-1},t_{N}-j_{N}}\\\
=&e^{-(n+1)\theta_{N}}Y_{t_{1},\dots,t_{N-1}}(t_{N}-n-1)+e^{-t_{N}\theta_{N}}\sum_{j_{N}=t_{N}-n}^{t_{N}}e^{j_{N}\theta_{N}}Z_{t_{1},\dots,t_{N-1},j_{N}},\end{split}$

for every $n\in\mathbb{N}$. Above

$\begin{split}e^{-(n+1)\theta_{N}}Y_{t_{1},\dots,t_{N-1}}(t_{N}-n-1)&=e^{-t_{N}\theta_{N}}e^{(t_{N}-n-1)\theta_{N}}Y_{t_{1},\dots,t_{N-1}}(t_{N}-n-1)\\\
&=e^{-t_{N}\theta_{N}}e^{m\theta_{N}}Y_{t_{1},\dots,t_{N-1}}(m)\end{split}$

by the change of variable $m=t_{N}-n-1$. Furthermore

$\begin{split}&e^{m\theta_{N}}Y_{t_{1},\dots,t_{N-1}}(m)\\\
=&e^{m\theta_{N}}\sum_{(i_{1},\dots,i_{N-1})\in\\{0,1\\}^{N-1}}(-1)^{\sum_{l=1}^{N-1}i_{l}}e^{-\langle(i_{1},\dots,i_{N-1},0),\Theta\rangle}X_{{t_{1}-i_{1}},\dots,t_{N-1}-i_{N-1},m},\end{split}$

which, by the assumptions, converges to zero in probability as $m\to-\infty$.
Hence

$\begin{split}&Y_{t_{1},\dots,t_{N-1}}(t_{N})\\\
&=\sum_{(i_{1},\dots,i_{N-1})\in\\{0,1\\}^{N-1}}(-1)^{\sum_{l=1}^{N-1}i_{l}}e^{-\langle(i_{1},\dots,i_{N-1},0),\Theta\rangle}X_{{t_{1}-i_{1}},\dots,t_{N-1}-i_{N-1},{t_{N}}}\\\
&=e^{-t_{N}\theta_{N}}\sum_{j_{N}=-\infty}^{t_{N}}e^{j_{N}\theta_{N}}Z_{t_{1},\dots,t_{N-1},j_{N}}\eqqcolon
Q_{t}^{(N)}.\end{split}$

In the following summations, let $(i_{1},\dots,i_{N-k-1})\in\\{0,1\\}^{N-k-1}$
and $(i_{1},\dots,i_{N-k-2})\in\\{0,1\\}^{N-k-2}$. We proceed by induction and
assume that for some $k\in\mathbb{N}\cup\\{0\\}$ it holds that

$\displaystyle\sum_{(i_{1},\dots,i_{N-k-1})}(-1)^{\sum_{l=1}^{N-k-1}i_{l}}e^{-\langle(i_{1},\dots,i_{N-k-1},0,\dots,0),\Theta\rangle}X_{{t_{1}-i_{1}},\dots,t_{N-k-1}-i_{N-k-1},t_{N-k},\dots,{t_{N}}}$
$\displaystyle=$ $\displaystyle
e^{-\sum_{l=0}^{k}t_{N-l}\theta_{N-l}}\sum_{j_{N-k}=-\infty}^{t_{N-k}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=0}^{k}j_{N-l}\theta_{N-l}}Z_{t_{1},\dots,t_{N-k-1},j_{N-k},\dots,j_{N}}\eqqcolon
Q_{t}^{(N-k)}.$ (13)

Now

$\displaystyle\sum_{(i_{1},\dots,i_{N-k-2})}\big{(}(-1)^{\sum_{l=1}^{N-k-2}i_{l}}e^{-\langle(i_{1},\dots,i_{N-k-2},0,\dots,0),\Theta\rangle}$
$\displaystyle\cdot
X_{{t_{1}-i_{1}},\dots,t_{N-k-2}-i_{N-k-2},t_{N-k-1},\dots,{t_{N}}}\big{)}$
$\displaystyle=$
$\displaystyle-\sum_{(i_{1},\dots,i_{N-k-2})}\big{(}(-1)^{1+\sum_{l=1}^{N-k-2}i_{l}}e^{-\langle(i_{1},\dots,i_{N-k-2},1,0,\dots,0),\Theta\rangle}$
$\displaystyle\cdot
X_{{t_{1}-i_{1}},\dots,t_{N-k-2}-i_{N-k-2},t_{N-k-1}-1,\dots,{t_{N}}}\big{)}+Q_{t}^{(N-k)}$
$\displaystyle=$
$\displaystyle\sum_{(i_{1},\dots,i_{N-k-2})}\big{(}(-1)^{\sum_{l=1}^{N-k-2}i_{l}}e^{-\langle(i_{1},\dots,i_{N-k-2},1,0,\dots,0),\Theta\rangle}$
$\displaystyle\cdot
X_{{t_{1}-i_{1}},\dots,t_{N-k-2}-i_{N-k-2},t_{N-k-1}-1,\dots,{t_{N}}}\big{)}+Q_{t}^{(N-k)}.$
(14)

Let $t^{*}=(t_{1},\dots,t_{N-k-2},t_{N-k},\dots,t_{N})$ and define

$\begin{split}&Y_{t^{*}}(t_{N-k-1})\\\
=&\sum_{(i_{1},\dots,i_{N-k-2})}(-1)^{\sum_{l=1}^{N-k-2}i_{l}}e^{-\langle(i_{1},\dots,i_{N-k-2},0,\dots,0),\Theta\rangle}X_{{t_{1}-i_{1}},\dots,t_{N-k-2}-i_{N-k-2},t_{N-k-1},\dots,{t_{N}}}.\end{split}$

Equation (14) gives

$\begin{split}Y_{t^{*}}(t_{N-k-1})=&e^{-\theta_{N-k-1}}Y_{t^{*}}(t_{N-k-1}-1)+Q_{t}^{(N-k)}\\\
=&e^{-(n+1)\theta_{N-k-1}}Y_{t^{*}}(t_{N-k-1}-n-1)\\\
&+\sum_{j_{N-k-1}=0}^{n}e^{-j_{N-k-1}\theta_{N-k-1}}Q^{(N-k)}_{t_{1},\dots,t_{N-k-2},t_{N-k-1}-j_{N-k-1},t_{N-k},\dots,t_{N}}\\\
=&e^{-(n+1)\theta_{N-k-1}}Y_{t^{*}}(t_{N-k-1}-n-1)\\\
&+e^{-t_{N-k-1}\theta_{N-k-1}}\sum_{\mathclap{j_{N-k-1}=t_{N-k-1}-n}}^{t_{N-k-1}}e^{j_{N-k-1}\theta_{N-k-1}}Q^{(N-k)}_{t_{1},\dots,t_{N-k-2},j_{N-k-1},t_{N-k},\dots,t_{N}}\end{split}$

for every $n\in\mathbb{N}$. Above

$\begin{split}e^{-(n+1)\theta_{N-k-1}}Y_{t^{*}}(t_{N-k-1}-n-1)=&e^{-t_{N-k-1}\theta_{N-k-1}}e^{(t_{N-k-1}-n-1)\theta_{N-k-1}}Y_{t^{*}}(t_{N-k-1}-n-1)\\\
=&e^{-t_{N-k-1}\theta_{N-k-1}}e^{m\theta_{N-k-1}}Y_{t^{*}}(m)\end{split}$

by the change of variable $m=t_{N-k-1}-n-1$. As before, the expression
converges to zero in probability as $m\to\\-\infty$. Hence, we obtain that

$\begin{split}&Y_{t^{*}}(t_{N-k-1})\\\
=&\sum_{(i_{1},\dots,i_{N-k-2})}\big{(}(-1)^{\sum_{l=1}^{N-k-2}i_{l}}e^{-\langle(i_{1},\dots,i_{N-k-2},0,\dots,0),\Theta\rangle}\\\
&\cdot
X_{{t_{1}-i_{1}},\dots,t_{N-k-2}-i_{N-k-2},t_{N-k-1},\dots,{t_{N}}}\big{)}\\\
=&e^{-t_{N-k-1}\theta_{N-k-1}}\sum_{\mathclap{j_{N-k-1}=-\infty}}^{t_{N-k-1}}e^{j_{N-k-1}\theta_{N-k-1}}Q^{(N-k)}_{t_{1},\dots,t_{N-k-2},j_{N-k-1},t_{N-k},\dots,t_{N}}\\\
=&e^{-t_{N-k-1}\theta_{N-k-1}}\Big{(}\sum_{\mathclap{j_{N-k-1}=-\infty}}^{t_{N-k-1}}e^{j_{N-k-1}\theta_{N-k-1}}e^{-\sum_{l=0}^{k}t_{N-l}\theta_{N-l}}\\\
&\cdot\sum_{j_{N-k}=-\infty}^{t_{N-k}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=0}^{k}j_{N-l}\theta_{N-l}}Z_{t_{1},\dots,t_{N-k-2},j_{N-k-1},\dots,j_{N}}\Big{)}\\\
=&e^{-\sum_{l=0}^{k+1}t_{N-l}\theta_{N-l}}\sum_{j_{N-k-1}=-\infty}^{t_{N-k-1}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=0}^{k+1}j_{N-l}\theta_{N-l}}Z_{t_{1},\dots,t_{N-k-2},j_{N-k-1},\dots,j_{N}},\end{split}$

which proves the induction step. Therefore choosing $k=N-2$ in (13) yields

$\begin{split}&\sum_{i_{1}\in\\{0,1\\}}(-1)^{i_{1}}e^{-\langle(i_{1},0,\dots,0),\Theta\rangle}X_{t_{1}-i_{1},t_{2},\dots,t_{N}}=X_{t}-e^{-\theta_{1}}X_{t_{1}-1,t_{2},\dots,t_{N}}\\\
&=e^{-\sum_{l=0}^{N-2}t_{N-l}\theta_{N-l}}\sum_{j_{2}=-\infty}^{t_{2}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=0}^{N-2}j_{N-l}\theta_{N-l}}Z_{t_{1},j_{2},\dots,j_{N}}\eqqcolon
Q_{t}^{(2)}\end{split}$

and we obtain a recursive equation

$X_{t}=e^{-\theta_{1}}X_{t_{1}-1,t_{1},\dots,t_{N}}+Q_{t}^{(2)}.$

By repeating the earlier procedure once more, we obtain that

$\begin{split}X_{t}=&e^{-t_{1}\theta_{1}}\sum_{j_{1}=-\infty}^{t_{1}}e^{j_{1}\theta_{1}}Q^{(2)}_{j_{1},t_{2},\dots,t_{N}}\\\
=&e^{-t_{1}\theta_{1}}\sum_{j_{1}=-\infty}^{t_{1}}e^{j_{1}\theta_{1}}e^{-\sum_{l=0}^{N-2}t_{N-l}\theta_{N-l}}\sum_{j_{2}=-\infty}^{t_{2}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=0}^{N-2}j_{N-l}\theta_{N-l}}Z_{j_{1},j_{2},\dots,j_{N}}\\\
=&e^{-\sum_{l=1}^{N}t_{l}\theta_{l}}\sum_{j_{1}=-\infty}^{t_{1}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=1}^{N}j_{l}\theta_{l}}Z_{j_{1},j_{2},\dots,j_{N}}\\\
=&e^{-\sum_{l=1}^{N}t_{l}\theta_{l}}\sum_{j_{1}=-\infty}^{t_{1}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=1}^{N}j_{l}\theta_{l}}\Delta_{j_{1},j_{2},\dots,j_{N}}G,\end{split}$

which, since $G\in\mathcal{G}_{{\Theta}}$, defines an almost surely finite
random variable. Hence it remains to prove that $X$ is stationary. To this
end, we show that the one dimensional distributions of $X$ are stationary. The
proof extends straightforwardly to multidimensional distributions. Let
$s=(s_{1},\dots,s_{N})\in\mathbb{Z}^{N}$. Note that

$X_{t}=\sum_{j_{1}=0}^{\infty}\dots\sum_{j_{N}=0}^{\infty}e^{-\sum_{l=1}^{N}j_{l}\theta_{l}}\Delta_{t_{1}-j_{1},\dots,t_{N}-j_{N}}G.$

Since $G$ is a stationary increment field, we have that

$\begin{split}&\sum_{j_{1}=0}^{M_{1}}\dots\sum_{j_{N}=0}^{M_{N}}e^{-\sum_{l=1}^{N}j_{l}\theta_{l}}\Delta_{t_{1}+s_{1}-j_{1},\dots,t_{N}+s_{N}-j_{N}}G\\\
\overset{\text{law}}{=}&\sum_{j_{1}=0}^{M_{1}}\dots\sum_{j_{N}=0}^{M_{N}}e^{-\sum_{l=1}^{N}j_{l}\theta_{l}}\Delta_{t_{1}-j_{1},\dots,t_{N}-j_{N}}G\end{split}$

for every $M_{1},\dots,M_{N}\in\mathbb{N}$. Moreover, since
$G\in\mathcal{G}_{{\Theta}}$, the iterated limits of both sides converge and
hence, the limits are equal in distribution. This gives

$\begin{split}X_{t+s}&=\sum_{j_{1}=0}^{\infty}\dots\sum_{j_{N}=0}^{\infty}e^{-\sum_{l=1}^{N}j_{l}\theta_{l}}\Delta_{t_{1}+s_{1}-j_{1},\dots,t_{N}+s_{N}-j_{N}}G\\\
&\overset{\text{law}}{=}\sum_{j_{1}=0}^{\infty}\dots\sum_{j_{N}=0}^{\infty}e^{-\sum_{l=1}^{N}j_{l}\theta_{l}}\Delta_{t_{1}-j_{1},\dots,t_{N}-j_{N}}G=X_{t}\end{split}$

and thus the proof is completed. ∎

######  Remark 2.29.

The property $G_{t}=0$ for all $t$ such that
$\sum_{i=1}^{N}t_{i}\in\\{0,-1,\dots,-N+1\\}$ of the class
$\mathcal{G}_{\Theta}$ is not utilized in the proof of Theorem 2.28. However,
in order to obtain uniqueness in the characterising Equation (11), we need to
pose this additional assumption, see Lemma 2.32.

######  Theorem 2.30.

Let $\Theta=(\theta_{1},\dots,\theta_{N})\in(0,\infty)^{N}$. Assume that the
field $X=(X_{t})_{t\in\mathbb{Z}^{N}}$ is stationary. Then there exists
$G=(G_{t})_{t\in\mathbb{Z}^{N}}\in\mathcal{G}_{{\Theta}}$ such that

$X_{t}=\langle\hat{\Theta},\hat{X}_{t}^{-}\rangle+\Delta_{t}G\quad\text{for
every }t\in\mathbb{Z}^{N}.$

###### Proof.

Applying Lamperti transformation gives

$\begin{split}\Delta_{t}X&=X_{t}-X_{t}^{-}=e^{-\langle
t,\Theta\rangle}Y_{e^{t}}-X^{-}_{t}=e^{-\langle
t,\Theta\rangle}(Y_{e^{t}}-\Delta_{t}Y)-X_{t}^{-}+e^{-\langle
t,\Theta\rangle}\Delta_{t}Y.\end{split}$

Hence

$X_{t}=e^{-\langle t,\Theta\rangle}(Y_{e^{t}}-\Delta_{t}Y)+e^{-\langle
t,\Theta\rangle}\Delta_{t}Y,$

where

$Y_{e^{t}}-\Delta_{t}Y=\sum_{(i_{1},\dots,i_{N})\in\\{0,1\\}^{N},i\neq\bf{0}}(-1)^{1+\sum_{l=1}^{N}i_{l}}Y_{e^{t_{1}-i_{1}},\dots,e^{t_{N}-i_{N}}}.$

Furthermore

$\begin{split}e^{-\langle
t,\Theta\rangle}(Y_{e^{t}}-\Delta_{t}Y)=&\sum_{(i_{1},\dots,i_{N})\in\\{0,1\\}^{N},i\neq\bf{0}}(-1)^{1+\sum_{l=1}^{N}i_{l}}e^{-\langle
t-i,\Theta\rangle}e^{-\langle
i,\Theta\rangle}Y_{e^{t_{1}-i_{1}},\dots,e^{t_{N}-i_{N}}}\\\
=&\sum_{(i_{1},\dots,i_{N})\in\\{0,1\\}^{N},i\neq\bf{0}}(-1)^{1+\sum_{l=1}^{N}i_{l}}e^{-\langle
i,\Theta\rangle}X_{{t_{1}-i_{1}},\dots,{t_{N}-i_{N}}}=\langle\hat{\Theta},\hat{X}_{t}^{-}\rangle,\end{split}$

where the last equality follows by Lamperti transformation. Let $G$ be the
field defined in Lemma 2.23. Then
$X_{t}=\langle\hat{\Theta},\hat{X}_{t}^{-}\rangle+\Delta_{t}G$. Moreover, the
proof of Theorem 2.28 and the property

$\lim_{m\to-\infty}e^{m\theta_{i}}X_{t_{1},\dots,t_{i-1},m,t_{i+1},\dots,t_{N}}\overset{\mathbb{P}}{\longrightarrow}0,$

for every $i\in\\{1,\dots,N\\}$ and
$t_{1},\dots,t_{i-1},t_{i+1},\dots,t_{N}\in\mathbb{Z}$, give us the
representation

$X_{t}=e^{-\sum_{l=1}^{N}t_{l}\theta_{l}}\sum_{j_{1}=-\infty}^{t_{1}}\dots\sum_{j_{N}=-\infty}^{t_{N}}e^{\sum_{l=1}^{N}j_{l}\theta_{l}}\Delta_{(j_{1},j_{2},\dots,j_{N})}G.$
(15)

Hence, $G\in\mathcal{G}_{{\Theta}}$. ∎

As a by-product, we observe the following result.

######  Lemma 2.31.

Let $G$ be defined as in Lemma 2.23. Then $G\in\mathcal{G}_{{\Theta}}$.

###### Proof.

By Lemma 2.23, it suffices to show that the limit (5) exists. Now since $Y$ is
$\Theta$-self-similar, there exists a stationary $X$ such that
$Y=\mathcal{L}_{\Theta}X$. Following the lines of the proof of Theorem 2.30,
we obtain Representation (15), where $G$ is defined as in Lemma 2.23. Hence
$G\in\mathcal{G}_{{\Theta}}$ by the proof of Theorem 2.30. ∎

######  Lemma 2.32.

Let $\Theta=(\theta_{1},\dots,\theta_{N})\in(0,\infty)^{N}$ be fixed. Then the
stationary increment field $G\in\mathcal{G}_{{\Theta}}$ in Theorems 2.28 and
2.30 is unique.

###### Proof.

Assume that $G,G^{\prime}\in\mathcal{G}_{{\Theta}}$ satisfy the recursive
Equation (11) of Theorems 2.28 and 2.30. Then,
$\Delta_{t}G=\Delta_{t}G^{\prime}$ for all $t\in\mathbb{Z}^{N}$. It remains to
show that this implies $G_{t}=G^{\prime}_{t}$ for all $t$, which we do by
induction.

Since $G_{t}=G^{\prime}_{t}=0$ for all $t$ such that
$\sum_{l=1}^{N}t_{l}\in\\{0,-1,\dots,-N+1\\}$, we may set an induction
assumption that there exists $s\in\mathbb{Z}$ such that $G_{t}=G^{\prime}_{t}$
for all $t$ such that $\sum_{l=1}^{N}t_{l}\in\\{s,s-1,\dots,s-N+1\\}$. In
order to complete the proof, it remains to show that $G_{t}=G^{\prime}_{t}$
under conditions $\sum_{l=1}^{N}t_{l}=s+1$ and $\sum_{l=1}^{N}t_{l}=s-N$. Let
$\sum_{l=1}^{N}t_{l}=s+1$. We now have

$\begin{split}\Delta_{t}G&=\sum_{(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}}(-1)^{\sum_{l=1}^{N}i_{l}}G_{t_{1}-i_{1},\dots,t_{N}-i_{N}}\\\
&=G_{t}+\sum_{\begin{subarray}{c}(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}\\\
i\neq\bf{0}\end{subarray}}(-1)^{\sum_{l=1}^{N}i_{l}}G_{t_{1}-i_{1},\dots,t_{N}-i_{N}},\end{split}$

where for the indices of $G$ it holds that

$\sum_{l=1}^{N}(t_{l}-i_{l})\in\\{s,s-1,\dots,s-N+1\\}.$

Therefore, by the induction assumption, $G_{t}=G_{t}^{\prime}$ for all $t$
such that $\sum_{l=1}^{N}t_{l}=s+1$. We proceed proving that
$G_{t}=G_{t}^{\prime}$ for all $t$ such that $\sum_{l=1}^{N}t_{l}=s-N$. Let
$\sum_{l=1}^{N}t_{l}=s-N$. Now

$\begin{split}&(-1)^{N}G_{t_{1},\dots,t_{N}}+\sum_{\begin{subarray}{c}(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}\\\
i\neq\bf{1}\end{subarray}}(-1)^{\sum_{l=1}^{N}i_{l}}G_{t_{1}+1-i_{1},\dots,t_{N}+1-i_{N}}\\\
&=(-1)^{N}G^{\prime}_{t_{1},\dots,t_{N}}+\sum_{\begin{subarray}{c}(i_{1},\dots,i_{N})\in\\{0,1\\}^{N}\\\
i\neq\bf{1}\end{subarray}}(-1)^{\sum_{l=1}^{N}i_{l}}G^{\prime}_{t_{1}+1-i_{1},\dots,t_{N}+1-i_{N}},\end{split}$

where

$\sum_{l=1}^{N}(t_{l}+1-i_{l})\in\\{s,s-1,\dots,s-N+1\\}.$

Hence, by the induction assumption, $G_{t}=G_{t}^{\prime}$ for all $t$ such
that $\sum_{l=1}^{N}t_{l}=s-N$. This completes the proof. ∎

The proof of our main theorem, Theorem 2.20, now follows by combining Theorems
2.28 and 2.30 together with Lemma 2.32.

###  2.3 Fractional Ornstein-Uhlenbeck fields on $\mathbb{Z}^{N}$

We illustrate our approach by constructing a stationary _fractional Ornstein-
Uhlenbeck field_ on $\mathbb{Z}^{2}$. For this purpose, let $G_{t_{1},t_{2}}$
be a discrete fractional Gaussian field on $(t_{1},t_{2})\in\mathbb{N}^{2}$. A
discrete fractional Gaussian field can be constructed from its continuous time
analogue by embedding to the space $\mathbb{N}^{2}$. That is, let $\tilde{X}$
be a self-similar Gaussian field having stationary rectangular increments such
that $\mathbb{E}\tilde{X}^{2}(1,1)=1$. Existence of such fields is studied in
[10]. We can now embed $\tilde{X}$ into $\mathbb{N}^{2}$ by considering values
of $\tilde{X}$ at points $(t_{1},t_{2})\in\mathbb{N}^{2}$. Now Definition 2.15
of $\mathcal{G}_{\Theta}$ allows us to extend $G$ into $\mathbb{Z}^{2}$ by
setting $G_{t,-t}=G_{t,-t-1}=0$ and requiring that $G$ has stationary
increments in the sense of Definition 2.8. As $G$ is Gaussian, in view of
Remark 2.18, it follows that our field $G=(G_{t})_{t\in\mathbb{Z}^{2}}$
belongs to $\mathcal{G}_{\Theta}$ for any $\Theta$. This allows us to define a
discrete time stationary _fractional Ornstein-Uhlenbeck field of the first
kind_ by Representation (6). On the other hand, using self-similarity of
$\tilde{X}$, we can define a stationary field through Lamperti transformation
and obtain, in view of Theorem 2.22, stationary _fractional Ornstein-Uhlenbeck
field of the second kind_ corresponding to a different noise
$G\in\mathcal{G}_{\Theta}$. Obviously, this approach can be used to define
fractional Ornstein-Uhlenbeck fields in arbitrary dimensions $N\geq 2$.

##  3 Conclusions

In this article we have provided a characterisation of discrete stationary
fields in terms of fields having stationary rectangular increments. Our
characterisation is analogous to the one-parameter case provided in [15]. As
an example, we have extended the notion of fractional Ornstein-Uhlenbeck
processes, introduced in [2] and [7], to the multi-parameter case.

Natural further prospects of our study are two folded. Firstly, it would be
interesting to generalise the characterisation into the continuous parameter
space $t\in\mathbb{R}^{N}$. In this case however, the increments $\Delta_{t}G$
are replaced with differentials, and the notion of differential $d_{t}G$ is
more subtle when $t\in\mathbb{R}^{N}$. This is a topic of a forthcoming
article. Secondly, parameter estimation is of paramount importance in
statistics, which in our setting would correspond to the estimation of the
parameter $\Theta$. This has been a topic of active research in the context of
fractional Ornstein-Uhlenbeck processes, see e.g. [6, 12] and the references
therein. For a more general setting in the context of multivariate time
series, parameter estimation is studied in [16]. We believe that ideas and
methods presented in [16] could be useful in our context as well.

## References

  * [1] Hermine Biermé, Mark M Meerschaert, and Hans-Peter Scheffler.  Operator scaling stable random fields.  Stoch. Process. their Appl., 117(3):312–332, 2007.
  * [2] P. Cheridito, H. Kawaguchi, and M. Maejima.  Fractional Ornstein-Uhlenbeck processes.  Electron. J. Probab., 8(3):1–14, 2003.
  * [3] Marianne Clausel.  Gaussian fields satisfying simultaneous operator scaling relations.  In Recent Developments in Fractals and Related Fields, pages 327–341. Springer, 2010.
  * [4] Paul Embrechts.  Selfsimilar processes.  In Princeton Series in Applied Mathematics. Princeton University Press, 2002.
  * [5] Marc G Genton, Olivier Perrin, and Murad S Taqqu.  Self-similarity and Lamperti transformation for random fields.  Stoch. Models, 23(3):397–411, 2007.
  * [6] Y. Hu and D. Nualart.  Parameter estimation for fractional Ornstein-Uhlenbeck processes.  Stat. Probab. Lett., 80(11–12):1030–1038, 2010.
  * [7] T. Kaarakka and P. Salminen.  Fractional Ornstein-Uhlenbeck processes.  Commun. Stoch. Anal., 5(1):121–133, 2011.
  * [8] J.W. Lamperti.  Semi-stable processes.  Trans. Amer. Math. Soc., 104:62–78, 1962.
  * [9] Vitalii Makogin and Yuliya Mishura.  Example of a Gaussian self-similar field with stationary rectangular increments that is not a fractional Brownian sheet.  Stoch. Anal. Appl., 33(3):413–428, 2015.
  * [10] Vitalii Makogin and Yuliya Mishura.  Gaussian multi-self-similar random fields with distinct stationary properties of their rectangular increments.  Stoch. Models, 35(4):391–428, 2019.
  * [11] Gennady Samorodnitsky, Murad S Taqqu, and RW Linde.  Stable non-Gaussian random processes: stochastic models with infinite variance.  Bull. London Math. Soc., 28(134):554–555, 1996.
  * [12] T. Sottinen and L. Viitasaari.  Parameter estimation for the Langevin equation with stationary-increment Gaussian noise.  Stat. Inference Stoch. Process., 21(3):569–601, 2018.
  * [13] Lauri Viitasaari.  Representation of stationary and stationary increment processes via Langevin equation and self-similar processes.  Stat. Prob. Lett., 115:45–53, 2016.
  * [14] Marko Voutilainen.  Modeling and estimation of multivariate discrete and continuous time stationary processes.  Front. Appl. Math. Stat., 6, 2020.
  * [15] Marko Voutilainen, Lauri Viitasaari, and Pauliina Ilmonen.  On model fitting and estimation of strictly stationary processes.  Mod. Stoch.: Theory Appl., 4(4):381–406, 2017.
  * [16] Marko Voutilainen, Lauri Viitasaari, Pauliina Ilmonen, Soledad Torres, and Ciprian Tudor.  Vector-valued generalized Ornstein–Uhlenbeck processes: Properties and parameter estimation.  Scand. J. Stat., 49(3):992–1022, 2022.