# Noise-induced instabilities in a stochastic Brusselator

Maximilian Engel  Freie Universität Berlin<EMAIL_ADDRESS>and
Guillermo Olicón-Méndez  Freie Universität Berlin  maximilian.engel@fu-
berlin.de

###### Abstract.

We consider a stochastic version of the so-called Brusselator - a mathematical
model for a two-dimensional chemical reaction network - in which one of its
parameters is assumed to vary randomly. It has been suggested via numerical
explorations that the system exhibits noise-induced synchronization when time
goes to infinity. Complementing this perspective, in this work we explore some
of its finite-time features from a random dynamical systems perspective. In
particular, we focus on the deviations that orbits of neighboring initial
conditions exhibit under the influence of the same noise realization. For
this, we explore its local instabilities via finite-time Lyapunov exponents.
Furthermore, we present the stochastic Brusselator as a fast-slow system in
the case that one of the parameters is much larger than the other one. In this
framework, an apparent mechanism for generating the stochastic instabilities
is revealed, being associated to the transition between the slow and fast
regimes.

###### Key words and phrases:

Finite-time Lyapunov exponents, stochastic oscillator, Brusselator

###### 2020 Mathematics Subject Classification:

34C26, 34E13, 34E15, 37G10, 37N25

##  1\. Introduction

In the past few decades there has been an increased interest in the effects of
noisy perturbations on deterministic dynamical systems and their bifurcations,
i.e. critical parameter-dependent changes of qualitative behaviour in the
system. In many cases, such a change of behaviour mirroring the deterministic
bifurcation can be observed in the density of the respective stationary
distribution (if it exists). The change of its shape, typically observed in
the local maxima of such density functions, has been called a P-bifurcation,
standing for phenomenological bifurcation [1].

Another direction has been addressing the qualitative change in the structure
of the objects which remain invariant under the random dynamics. In this
context, a seminal insight was obtained at the hand of the pitchfork
bifurcation with additive noise [7], stating that for any such perturbation of
the pitchfork normal form the system admits a unique random equilibrium
independently from the parameter value, whereas the deterministic system
exhibits a bifurcation from one to two locally attracting equilibria. Even in
the bistable scenario, when the system is perturbed with additive noise the
unique random equilibrium is globally attracting, from which it follows that
any two trajectories of the system approach each other as time grows. This
phenomenon is commonly known as _noise-induced synchronization_. However, it
has been argued that the random bifurcation persists in a more subtle manner,
being detected by changes of finite-time stability and the dichotomy spectrum
[4]. This point of view, measuring a random bifurcation in particular by means
of _finite-time Lyapunov exponents_ (FTLEs), has been extended to Hopf
bifurcations in normal form [10] and the infinite-dimensional situation
recently [3]. This perspective puts emphasis on the transient dynamics rather
than on the asymptotic behaviour of the orbits, which becomes of particular
relevance in the applied sciences.

Despite the apparent asymptotic stabilizing effect that noisy perturbations
often have on deterministic systems, there has been recent progress on showing
noise-induced chaos for particular classes of stochastically perturbed
oscillators [11], implying a stochastic bifurcation from synchronization to
chaotic attractors. Building up on numerical investigations [24] and
_Furstenberg-Khasminskii_ formulas for Lyapunov exponents [17, 18], it was
shown lately [5] that a positive Lyapunov exponent can be obtained for
appropriate parameter combinations in a model of Hopf bifurcations with shear
terms and additive noise.

In this work we turn our attention to a stochastic version of the Brusselator.
This system is an example of a chemical reaction network whose reaction rate
equation exhibits a Hopf bifurcation, i.e. the emergence of sustained
oscillations after passing a certain critical parameter value. The model has
proven to be dynamically very rich. For instance, when considering the space-
time dynamics of the Brusselator as a reaction-diffusion equation the system
exhibits successive bifurcations which induce dissipative structures [28].
Furthermore, under the influence of a periodic forcing, it may exhibit chaotic
behaviour [25]. Here, we explore a complementary route by assuming that one of
the system’s parameters is perturbed by Brownian noise, yielding a stochastic
differential equation (SDE) as considered in [2].

We employ the perspective of _random dynamical systems_ (RDS) in order to
understand the trajectory-wise behaviour of the solutions of this SDE on the
state space $X=\mathbb{R}^{2}$. Taking the time set $\mathbb{R}$, the RDS
$(\theta,\varphi)$ consists of the noise model
$(\Omega,\mathcal{B},\mathbb{P},\theta)$, where $(\Omega,\mathcal{B})$ is the
canonical path space with Wiener measure $\mathbb{P}$, and a
$\mathbb{P}$-invariant dynamical system
$\theta=(\theta_{t})_{t\in\mathbb{R}}$, and the family of maps
$\varphi^{t}_{\omega}:X\rightarrow X$ which determines the state of the system
at time $t\geq 0$ for a fixed noise realization $\omega\in\Omega$. The RDS
approach allows to compare two trajectories with the same noise realizations
$\omega\in\Omega$. In other words, given that the map $\varphi^{t}_{\omega}$
determines the evolution of the system, we want to compare
$\varphi^{t}_{\omega}(x_{0},y_{0})$ and $\varphi^{t}_{\omega}(x_{1},y_{1})$
for distinct elements $(x_{0},y_{0})$ and $(x_{1},y_{1})$ of the state space
$X$. Specifically, the RDS approach analyses whether these two trajectories
synchronize or diverge asymptotically, determining the form of a _random
attractor_ [8]. We refer the reader to [1] for a comprehensive exposition of
RDS theory.

While it has been suggested numerically in [2] that the Brusselator with
parametric noise exhibits noise-induced synchronization as
$t\rightarrow\infty$, we adopt a different point of view by quantifying
finite-time instabilities via the FTLEs. We observe numerically that the Hopf
bifurcation of the stochastic Brusselator persists in a similar manner as in
the cases previously mentioned – the system exhibits instabilities in finite-
time windows.

We perform our numerical explorations in two complementary directions: we
calculate the FTLE for a fixed noise realization $\omega\in\Omega$, first by
fixing a time length $T>0$ and varying the intial condition, and secondly by
fixing an initial condition $(x,y)\in\mathbb{R}^{2}_{+}$ and seeing the
evolution of the FTLE through time windows of increasing length. In the former
scenario, it is shown that the state space has regions of instability (i.e.
where the FTLE is positive) which are not only localized around the
deterministic unstable equilibrium. In the latter, we show that the FTLEs may
attain positive values for arbitrarily large time windows whenever one of the
parameters is large enough. It is precisely this last setting where the
Brusselator can be formulated as a slow-fast system. We sketch the different
time scales involved in the deterministic setting, and how they extend to the
stochastic picture.

### Structure of the paper

The rest of the paper is structured as follows. In Section 2 we introduce the
notion of a chemical reaction network, and its different modeling approaches,
depending on the volume scaling under consideration. In particular, we present
the Brusselator as a chemical reaction network, and the model that we analyse
throughout the rest of this paper. In Section 3 we introduce the notion of
finite-time Lyapunov exponents, as an indicator of instabilities in bounded
time windows. In Section 4, we perform a change of coordinates which turns the
stochastic Brusselator into a slow-fast system in standard form. We further
sketch a possible mechanism which generates the finite-time instabilities
observed via the finite-time Lyapunov exponents. Last, in Section 5 we provide
some concluding remarks with an outlook for future work.

##  2\. The Brusselator: a chemical oscillator

Consider the general model of $L$ chemical species $S_{1},\dots,S_{L}$ such
that the state of the whole reaction system at time $t\geq 0$ is given by
$Z(t)=(Z_{1},\dots,Z_{L}(t))\in\mathbb{N}_{0}^{L}$ with $Z_{l}(t)$ denoting
the number of molecules of type $S_{l}$ at time $t$. Reactions $R_{k}$,
changing the number of molecules of a subset of species, occur at rates
$\alpha_{k}(x)$ which are defined as propensity funtions involving reaction
rate constants $\gamma_{k}$. Assuming that these functions depend only on the
current value of $Z(t)$, we obtain the _time-change representation of the
Markov jump process_

(2.1)
$Z(t)=Z(0)+\sum_{k=1}^{K}\mathcal{P}_{k}\left(\int_{0}^{t}\alpha_{k}(Z(s))\mathrm{d}s\right)\nu_{k},$

where $\nu_{k}$ is the state change vector of reaction $R_{k}$ and
$\mathcal{P}_{k}$, $k=1,\dots,K$, denote independent Poisson processes with
unit rate. For high accumulations of reaction events, the Poisson variables
can be approximated by normal random variables, yielding the _chemical
Langevin equation_ in Itô form

(2.2)
$\mathrm{d}Z(t)=\sum_{k}\alpha_{k}(Z(t))\nu_{k}\mathrm{d}t+\frac{1}{\sqrt{V}}\sum_{k}\sqrt{\alpha_{k}(Z(t))}\nu_{k}\mathrm{d}W_{k}(t),$

where each $W_{k}(t)$ denotes a one-dimensional Wiener process and the state
space for solutions of equation (2.2) is now extended to $\mathbb{R}^{L}$. In
the thermodynamic limit where volume and volume-dependent molecule numbers
scale in a fixed proportion approaching infinity, we may average over the
noise terms to obtain the purely deterministic version of equation (2.2), i.e.
the ODE

(2.3) $\mathrm{d}Z(t)=\sum_{k}\alpha_{k}(Z(t))\cdot\nu_{k}\mathrm{d}t,$

also called _reaction rate equation_. For chemical reaction networks with
different volume-dependent scalings, the three description levels (2.1), (2.2)
and (2.3) can also be combined into hybrid models [30]. For an introductory
exposition of chemical reaction networks from these different perspectives see
[31].

In this work we focus on the chemical reaction network

(2.4) $\begin{array}[]{lrcl}R_{1}:&A&\rightarrow&X,\\\
R_{2}:&B+X&\rightarrow&Y+D,\\\ R_{3}:&2X+Y&\rightarrow&3X,\\\
R_{4}:&X&\rightarrow&E.\end{array}$

We refer to (2.4) as the Brusselator reaction network. As classically
considered, we suppose that $A$ and $B$ satisfy the infinite pool assumption,
for which an inexhaustible supply of both reactants is injected into the
system so that their number of molecules (or their concentration when properly
rescaled with the volume of the system) are assumed to be constant parameters.
Under this assumption, system (2.4) becomes an open system. It is relevant to
note here that both $E$ and $D$ are products of the system which do not
interact back with the rest of the reactants. For this reason, we focus only
on the variables $X$ and $Y$ in the rest of the paper.

The Brusselator was originally proposed in [27] as an example of a chemical
system involving only two species which under suitable conditions is driven
out of the thermodynamic equilibrium. The exploration of such type of reaction
networks becomes relevant since many biochemical systems depend essentially on
two main species. For such systems which involve only uni- and bimolecular
reactions, it is known that if the rate functions $\alpha_{k}$ obey the law of
mass action the dynamics given by their reaction rate equations (cf. (2.3))
does not admit attracting limit cycles, see [28, Chapter 7]. The intermediate
reaction channel $R_{3}$ thus provides a simple trimolecular step for which
interesting dynamics may emerge. Notice, however, that the reaction $R_{3}$ is
physically unrealistic, making the Brusselator only a prototype of chemical
instabilities.

###  2.1. The reaction jump process

Following [31], for each of the reaction channels $R_{k}$ ($k=1,2,3,4$) we
define the vectors $\vec{v}_{k}$ as

$\begin{array}[]{cc}\begin{array}[]{rcl}\nu_{1}&=&(1,0,0,0),\\\
\nu_{2}&=&(-1,1,1,0),\end{array}&\begin{array}[]{rcl}\nu_{3}&=&(1,-1,0,0),\\\
\nu_{4}&=&(-1,0,0,1).\\\ \end{array}\end{array}$

According to the law of mass action, the corresponding reaction rates
$\alpha_{k}$ for each reaction channel $R_{k}$ are given as [9]

$\begin{array}[]{cc}\begin{array}[]{rcl}\alpha_{1}(Z)&=&\gamma_{1}A\\\
\alpha_{2}(Z)&=&\gamma_{2}V^{-1}B\cdot
X\end{array}&\begin{array}[]{rcl}\alpha_{3}(Z)&=&\gamma_{3}V^{-2}\cdot
X(X-1)Y\\\ \alpha_{4}(Z)&=&\gamma_{4}\cdot X,\\\ \end{array}\end{array}$

where $Z=(X,Y,D,E)$, and, abusing the notation, the variables $Z_{i}(t)$
denote the number of molecules of species $Z_{i}\in\\{X,Y,D,E\\}$ at time $t$.
The constants $\gamma_{i}$ for $i=1,2,3,4$ are positive real numbers,
$A,B\in\mathbb{N}$, and $V>0$ denotes the volume of the system. The evolution
of the system is thus given by the Markov jump process (2.1).

The Brusselator reaction network, modeled as a Markov jump process, has been
numerically studied in [15], whose oscillatory behaviour was compared with the
Lotka-Volterra system. The numerical scheme employed is called the direct
method, and splits the jump process into two stochastic processes
$(Z_{n}),(T_{n})$, where $T_{n}$ is the time at which a reaction $R_{k}$
occurs and $Z_{n}=Z(T_{n})$ is the state of the system which remains constant
until another reaction occurs. Both quantities evolve according to a random
difference equation, and thus define an RDS, see [14] for general results in
this direction.

###  2.2. Reaction rate equations

By considering the thermodynamic limit when $V\rightarrow\infty$, the reaction
rate equation (2.3) for the Brusselator reaction network reads as

(2.5) $\begin{array}[]{rcl}\dot{X}&=&a-(1+b)X+X^{2}Y,\\\
\dot{Y}&=&bX-X^{2}Y.\end{array}$

System (2.5) has been thoroughly studied [28], exhibiting the following key
properties: it has a unique equilibrium point $(a,b/a)$ which is stable
whenever $b<b_{crit}:=1+a^{2}$ and unstable when $b>b_{crit}$. An attracting
limit cycle exists precisely in the latter case due to a Hopf bifurcation, see
Figure 1.

\begin{overpic}[width=433.62pt]{TimeSeries_LargeB.png}
\put(52.0,-3.0){\scriptsize$t$} \end{overpic} (a)

\begin{overpic}[width=433.62pt]{Cycles.png} \put(-2.0,40.0){\scriptsize$y$ }
\put(50.0,-2.0){\scriptsize$x$} \end{overpic} (b)

Figure 1. In (A), the time series of two different orbits of the system for
two different parameter values respectively. Oscillatory behaviour arises when
$b>b_{crit}$, and exhibits a clear time scale separation when $b\gg b_{crit}$.
In (B), different profiles of the attracting limit cycle are depicted
depending on the value of $b$. The limit cycles grow as $b$ becomes larger.

###  2.3. The chemical Langevin equation

Consider now the rescaled variable $(X,Y)\mapsto\frac{1}{V}(X,Y)$. For large
values of $V$ the Markov jump process can be approximated by the chemical
Langevin equation (2.2), which for the Brusselator reads as the Itô SDE

(2.6)
$\begin{array}[]{lll}dX_{t}&=&\left(\tilde{\gamma}_{1}\tilde{A}-\tilde{\gamma}_{2}\tilde{B}X_{t}+\tilde{\gamma}_{3}X_{t}(X_{t}-\frac{1}{V})Y_{t}-\tilde{\gamma}_{4}X_{t}\right)dt+\sqrt{\frac{\tilde{\gamma}_{1}\tilde{A}}{V}}d\tilde{W}_{t}^{1}\\\
&&-\sqrt{\frac{\tilde{\gamma}_{2}\tilde{B}}{V}X_{t}}d\tilde{W}_{t}^{2}+\sqrt{\frac{\tilde{\gamma}_{3}}{V}X_{t}(X_{t}-\frac{1}{V})Y_{t}}d\tilde{W}_{t}^{3}-\sqrt{\frac{\tilde{\gamma}_{4}}{V}X_{t}}d\tilde{W}_{t}^{4},\\\
\\\
dY_{t}&=&\left(\tilde{\gamma}_{2}\tilde{B}X_{t}-\tilde{\gamma}_{3}X_{t}(X_{t}-\frac{1}{V})Y_{t}\right)dt+\sqrt{\frac{\tilde{\gamma}_{2}\tilde{B}}{V}X_{t}}d\tilde{W}_{t}^{2}\\\
&&-\sqrt{\frac{\tilde{\gamma}_{3}}{V}X_{t}(X_{t}-\frac{1}{V})Y_{t}}d\tilde{W}_{t}^{3},\end{array}$

where the rescaled variables $\tilde{A}=\frac{1}{V}A$ and
$\tilde{B}=\frac{1}{V}B$ are assumed to be constant parameters independent of
the volume $V$, while for $i=1,2,3,4$ the quantities $\tilde{\gamma}_{i}$ are
positive constants, and $\tilde{W}_{t}^{i}$ are mutually independent one-
dimensional Wiener processes.

While solutions of the Markov jump process and the reaction rate equations are
typically well defined for all $t\geq 0$, this is very often not the case for
the chemical Langevin equation. In particular, solutions of (2.6) may not
exist for all $t>0$, since the noisy perturbations push the trajectories
through the boundary of the domain $\mathbb{R}_{+}^{2}$ with positive
probability. This is clearly an inconsistency with the chemical interpretation
of $X$ and $Y$ as concentrations. There have been attempts to resolve this
issue. For example, in [29] a modified Langevin equation has been proposed
mainly for uni- and bimolecular reactions which respect the positivity of
solutions. However, note that the Brusselator involves a trimolecular
reaction, namely $R_{3}$. As a different strategy, a method with reflecting
boundary conditions has been considered in [22]. This last proposal does not
come from first principles, and some interpretation of the original problem
may be lost.

###  2.4. Brusselator with parametric noise

In the rest of this work, we focus on a stochastic version of equation (2.5),
where the parameter $b$ is replaced by a perturbation with Brownian noise of
Stratonovich type $b+\sigma\circ dW_{t}$. More precisely, in this work we
study the system [2]

(2.7)
$\begin{array}[]{rcl}\mathrm{d}x_{t}&=&\left(a-(1+b)x_{t}+x_{t}^{2}y_{t}\right)\mathrm{d}t-\sigma
x_{t}\circ\mathrm{d}W_{t},\\\
\mathrm{d}y_{t}&=&\left(bx_{t}-x_{t}^{2}y_{t}\right)\mathrm{d}t+\sigma
x_{t}\circ\mathrm{d}W_{t}.\\\ \end{array}$

We consider the stochastic Brusselator (2.7) to be an intermediate model
between the deterministic system (2.5) and the chemical Langevin equation
(2.6), guaranteeing the existence of solutions that remain in the positive
quadrant $\mathbb{R}_{+}^{2}$ for all positive times. The stochastic
perturbations of this form still induce intriguing dynamical features that we
will explore in the rest of this work.

The solutions of (2.7) induce a RDS as indicated in the Introduction. Formally
speaking, we consider the noise space
$\Omega=C_{0}(\mathbb{R},\mathbb{R}^{2})$ endowed with is canonical
$\sigma$-algebra $\mathcal{B}$, the Wiener measure $\mathbb{P}$, and the shift
maps $\theta=(\theta_{t})_{t\in\mathbb{R}}$ given by

$\theta_{t}\omega(\cdot)=\omega(\cdot+t)-\omega(t),$

which are $\mathbb{P}$-invariant and ergodic. The state of the system is given
by a measurable function
$\varphi:\mathbb{R}_{+}\times\Omega\times\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$,
called the cocycle map over $\theta$, such that for all $t,s\in\mathbb{R}_{+}$
and $\omega\in\Omega$,

$\varphi^{0}_{\omega}=\operatorname{id},\qquad\varphi^{t+s}_{\omega}=\varphi^{t}_{\theta_{s}\omega}\circ\varphi^{s}_{\omega},$

where we adopt the notation $\varphi^{t}_{\omega}=\varphi(t,\omega,\cdot)$.
The pair $(\theta,\varphi)$ is a smooth RDS if the cocycle map
$\varphi^{t}_{\omega}$ is $C^{1}$ for all $t\in\mathbb{R}_{+}$ and
$\omega\in\Omega$.

We recall the following result.

######  Theorem 2.1.

The SDE (2.7) generates a smooth RDS in $\mathbb{R}^{2}_{+}$. Moreover, system
(2.7) admits at most one stationary distribution.

###### Proof.

See [2]. ∎

Theorem 2.1 is a consequence of standard existence results for SDEs. In fact,
it follows that solutions exist locally backwards in time, so that for each
$(\omega,x,y)\in\Omega\times\mathbb{R}_{+}^{2}$ there is
$\tau_{-}(\omega,x,y)<0$ for which the RDS can be extended up to the interval
$(\tau_{-}(\omega,x,y),\infty)$. Nevertheless, we emphasize that the shift
maps $\theta_{t}$ on the base noise space $\Omega$ are defined for all
$t\in\mathbb{R}$, and it is sufficient for the cocycle map to be defined in
$\mathbb{R}_{+}$.

######  Remark 2.2.

In [2] it is conjectured that in fact, for any $a,b,\sigma>0$ there is a
unique stationary distribution $\rho$ with corresponding invariant measure
$\mu$ for the RDS, whose sample measures $\mu_{\omega}$ are given by delta
distributions $\mu_{\omega}=\delta_{A(\omega)}$. In particular, the random
equilibrium point $A(\omega)$ would be such that for all $t\geq 0$

$\varphi^{t}_{\omega}(A(\omega))=A(\theta^{t}\omega),\quad\mathbb{P}-a.s.$

It has been further conjectured that the top Lyapunov exponent is negative,
implying that this random equilibrium is exponentially attracting [26], and
thus every pair of trajectories synchronizes exponentially fast as
$t\rightarrow\infty$. While it has been shown in [2] that there are initial
conditions for which solutions do not exist for all negative times with
positive probability, numerical explorations suggest that, in fact, the random
equilibrium $A(\omega)$ constitutes the only trajectory whose dynamics can be
extended for all $t<0$.

In order to visualise the trajectories of the system we employ the Euler-
Maruyama method [19] which approximates the solution of a general Itô SDE

(2.8) $dX_{t}=f(X_{t},t)+\sum_{i=1}^{m}{g_{i}(X_{t},t)\;\mathrm{d}W_{t}^{i}},$

with vector fields $g_{i}=(g_{1i},g_{2i},\ldots,g_{ni})^{\top}$ and a fixed
initial condition $X_{0}$ at $t=0$, via a recursive random map

(2.9) $X_{n+1}=X_{n}+hf(X_{n},t_{n})+\sum_{i=1}^{m}g_{i}(X_{n},t_{n})\;\Delta
W_{n}^{i},$

where $\left(\Delta W_{n}^{i}\right)_{n=0}^{\infty}$ are independently and
identically distributed normal random variables with mean $0$ and variance $h$
for step size $h>0$ sufficiently small, and $t_{n}=nh$.

Hence, for numerical solutions we convert (2.7) to its Itô version using the
Wong-Zakai corrections. More precisely, for an Itô SDE in $\mathbb{R}^{n}$ of
the form (2.8) where
$f,g_{i}:\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}^{n}$ for
$i=1,2,\ldots,m$ are smooth, let
$g:\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}^{n\times m}$ be the
matrix whose columns consist of the vector fields $g_{i}$ for $i=1,\ldots,m$,
and let $W_{t}:=(W_{t}^{1},\ldots,W_{t}^{m})^{\top}$. The Itô and Stratonovich
integrals are component-wise related for $l=1,2,\ldots,n$ by

(2.10) $\left(\int_{0}^{t}{g(X_{s},s)\circ
dW_{s}}\right)_{l}=\left(\int_{0}^{t}{g(X_{s},s)dW_{s}}\right)_{l}+\frac{1}{2}\int_{0}^{t}\sum_{k=1}^{m}\sum_{j=1}^{n}\frac{\partial
g_{lk}}{\partial X_{j}}(X_{t},t)\cdot g_{jk}(X_{t},t)$

By using (2.10) for (2.7) we obtain the Itô SDE

(2.11)
$\begin{array}[]{rcl}dx&=&[a-(1+b)x+x^{2}y+\frac{\sigma^{2}}{2}x]dt-\sigma
xdW_{t},\\\ dy&=&[bx-x^{2}y-\frac{\sigma^{2}}{2}x]dt+\sigma
xdW_{t},\end{array}$

which we can use for numerical simulations of the induced random dynamical
system. We start by exploring the _two-point motion_ , i.e. comparing two
trajectories $\varphi_{\omega}^{t}(x_{0},y_{0})$ and
$\varphi_{\omega}^{t}(x_{1},y_{1})$ for a fixed noise path $\omega\in\Omega$.
Figure 2 depicts such two trajectories for parameter values before and after
the Hopf bifurcation.

\begin{overpic}[width=433.62pt]{StocTraj_Eq.png}
\put(52.0,-2.0){\scriptsize$x$} \put(-2.0,40.0){\scriptsize$y$} \end{overpic}
(a)

\begin{overpic}[width=433.62pt]{StocTraj_Oscil.png}
\put(52.0,-2.0){\scriptsize$x$} \put(-2.0,40.0){\scriptsize$y$} \end{overpic}
(b)

Figure 2.  Here we portray two different orbits of (2.11) with neighbouring
initial conditions, which are subject to the same noise realization. We take a
step size $h=10^{-3}$ and iterate the Euler-Maruyama scheme from $t_{0}=0$ to
$T=150$. We take $\sigma=0.1$ and $a=1$. In (A), we consider the parameter
$b=b_{crit}-1$ and the orbits tend to stay around the deterministic
equilibrium point. In (B), we have $b=b_{crit}+1$, and the orbits exhibit
fluctuations around the deterministic limit cycle.

The qualitative picture in Figure 2 for $b>b_{crit}$ has some interesting
features. On the one hand, given an initial condition
$(x,y)\in\mathbb{R}_{+}^{2}$ the trajectory $\varphi^{t}_{\omega}(x,y)$
clearly resembles the triangle-like shape of the deterministic limit cycle,
cf. Figure 1 – starting from the bottom right, moves to the left in a slow
time scale until it takes a fast transition towards a more vertical-like
movement. Once the trajectory is closer to the vertical axis, it heads towards
its uppermost part in an ultraslow fashion. Once on its apex, the system
“slides down” almost in a straight diagonal line back to the bottom right in a
fast time scale, where the trajectory repeats the process. However, due to the
random perturbation of the system, the trajectory seems to choose different
highest points, and times to reach them on each repetition. A similar
behaviour for the Brusselator has been observed previously in [15] at the
level of the Markov jump process. A more detailed description of the dynamics
seen in Figure 2 in terms of a fast-slow system will be given in Section 4,
see also Figure 7.

In order to explore in more detail dynamical features in terms of noise-
induced instabilities or synchronization, we calculate the distance between
trajectories

$d_{t}\left((x_{0},y_{0}),(x_{1},y_{1})\right)=\|\varphi_{\omega}^{t}(x_{0},y_{0})-\varphi_{\omega}^{t}(x_{1},y_{1})\|,$

as shown in Figure 3. We observe that, in both cases before and after the
bifurcation, any two arbitrary orbits approach each other eventually as
$t\rightarrow\infty$. However, the fashion in which this is done is
qualitatively different. Indeed, when $b<b_{crit}$ the convergence between
orbits is uniform, while in the case when $b>b_{crit}$ the orbits deviate
repeatedly from each other, but they do so with less strength as time grows.
In spite of that, the amplitude of their distance may become larger again at
intermediate time windows.

\begin{overpic}[width=433.62pt]{Dist_Eq.png} \put(52.0,-2.0){\scriptsize$t$}
\put(-4.0,35.0){\scriptsize$d_{t}$} \end{overpic} (a)

\begin{overpic}[width=433.62pt]{Dist_Oscil.png}
\put(52.0,-2.0){\scriptsize$t$} \put(-4.0,35.0){\scriptsize$d_{t}$}
\end{overpic} (b)

Figure 3. For the cases above, we depict here the distance between the two
different orbits as a function of time. We consider $b=b_{crit}-1$ and
$b=b_{crit}+1$ in (A) and (B), respectively. The rest of the parameters were
chosen to be the same as in Figure 2. Note the big difference between the
scales in both axes.

##  3\. Finite-time Lyapunov exponents

As mentioned in Remark 2.2 and sketched in Figure 3, it was conjectured that
any two trajectories of the system, under the influence of the same noise
realization, synchronize exponentially fast for any parameter combination of
$a,b$ and $\sigma$. However, Figure 3 still shows a clear qualitative
difference between the case $b<b_{crit}$ and $b>b_{crit}$. Specifically, when
$b>b_{crit}$ there seems to be a large time window where the two trajectories
deviate from each other recurrently before their distance finally starts
approaching zero.

A useful quantity to study local finite-time stability in RDS are the finite-
time Lyapunov exponents (FTLEs). For their computation, we consider the
linearization along an orbit $\varphi^{t}_{\omega}(x,y)$ of (2.7). It is
straightforward to see that

$v_{t}=\text{D}\varphi^{t}_{\omega}(x,y)\cdot v_{0},\quad
v_{0}\in\mathbb{R}^{2},$

satisfies the variational equation for each $(x,y)\in\mathbb{R}_{+}^{2}$ and
$\mathbb{P}$-a.e $\omega\in\Omega$ which, for system (2.7), reads

(3.1)
$\mathrm{d}v=J\left(\varphi_{\omega}^{t}(x,y)\right)v\mathrm{d}t+\sigma\left[\begin{array}[]{cc}-1&0\\\
1&0\end{array}\right]v\circ\mathrm{d}W_{t},$

where $J$ is the Jacobian matrix of the drift given by

$J(x,y)=\left[\begin{array}[]{cc}-(1+b)+2xy&x^{2}\\\
b-2xy&-x^{2}\end{array}\right].$

We can also consider the matrix form of (3.1) as the SDE

(3.2)
$\mathrm{d}\Phi=J\left(\varphi^{t}_{\omega}(x,y)\right)\Phi\mathrm{d}t+\sigma\left[\begin{array}[]{cc}-1&0\\\
1&0\end{array}\right]\Phi\circ\mathrm{d}W_{t},$

where $\Phi\in\mathbb{R}^{2\times 2}$.

######  Definition 3.1.

Let $\varphi^{t}_{\omega}(x,y)$ be a trajectory of (2.7), and let $\Phi_{t}$
be the solution of (3.2) at time $t$, with initial condition
$\Phi_{0}=\operatorname{id}$. The (maximal) finite-time Lyapunov exponent
(FTLE) of (2.7) for $(\omega,x,y)$ at time $T>0$ is given by

(3.3) $\lambda^{T}(\omega,x,y):=\frac{1}{T}\ln\left\|\Phi_{T}\right\|,$

where $\|\cdot\|$ is the operator norm or, equivalently, its largest singular
value.

Note that the classical maximal Lyapunov exponent is given by

$\lambda(\omega,x,y)=\limsup_{T\rightarrow\infty}\lambda^{T}(\omega,x,y).$

Again, in order to solve numerically (3.2) above one must convert the system
into its Itô version using (2.10). We thus get the system

$\begin{array}[]{rcl}\mathrm{d}x&=&[a-(1+b)x+x^{2}y+\frac{\sigma^{2}}{2}x]\mathrm{d}t-\sigma
x\mathrm{d}W_{t}\\\
\mathrm{d}y&=&[bx-x^{2}y-\frac{\sigma^{2}}{2}x]\mathrm{d}t+\sigma
x\mathrm{d}W_{t}.\\\
\mathrm{d}v&=&\left\\{J\left(\varphi_{\omega}^{t}(x_{0},y_{0})\right)+\frac{\sigma^{2}}{2}\left[\begin{array}[]{cc}1&0\\\
-1&0\end{array}\right]\right\\}v\mathrm{d}t++\sigma\left[\begin{array}[]{cc}-1&0\\\
1&0\end{array}\right]v\cdot\mathrm{d}W_{t}.\end{array}$

By fixing a noise realization, we compute numerically the FTLEs throughout the
state space by implementing the Euler-Maruyama method, see (2.9). In Figure 4,
the FTLE (3.3) is numerically calculated over a mesh of fine grid of initial
conditions in the state space. We use the period of the deterministic limit
cycle $\tau$ as a reference, so that $T\approx\tau/2$. For $b>b_{crit}$
(Figure 4(A)), we observe an area of instability around the limit cycle within
a certain region whose role is explored in more detail later, when exploring
the multiscale dynamics. For $b<b_{crit}$ (Figure 4(B)), no such instabilities
can be found. Note that the change of behaviour in terms of the FTLEs is
really measured around the critical attracting object, i.e. the equilibrium
before the bifurcation and the limit cycle after the bifurcation. This is in
contrast to where the change of stability is only measured around the
equilibrium (see e.g. [10]).

In Figure 5 we explore the variation of the FTLE with respect to $T$, by
fixing both $\omega\in\Omega$ and $(x,y)\in\mathbb{R}^{2}_{+}$. Here, an
interesting phenomenon is revealed: as the parameter $b$ increases, it seems
to allow the FTLEs to attain positive values at much larger values of $T$. In
other words, as $b$ grows the system exhibits finite-time instabilities in
larger time windows. It remains an open question whether this can be achieved
in the limit as $T\rightarrow\infty$ for an adequate rescaling of both $b$ and
$\sigma$, showing that for some values of the parameters the classical
Lyapunov exponent is positive. In the following, we explore the case in which
$b\rightarrow\infty$.

\begin{overpic}[width=433.62pt]{FTLEmap_final2.png}
\put(50.0,-2.0){\scriptsize$x$} \put(-2.0,40.0){\scriptsize$y$} \end{overpic}
(a)

\begin{overpic}[width=433.62pt]{FTLEmapeq_2.png}
\put(50.0,-2.0){\scriptsize$x$} \put(-2.0,40.0){\scriptsize$y$} \end{overpic}
(b)

Figure 4. We depict a landscape of the maximal FTLE as given in Definition
3.1. Here we took $a=1$ and $\sigma=0.1$. The time window length considered
was $T\approx\tau/2\equiv\frac{1}{2\Omega}$, where $\Omega$ is the dominant
frequency obtained when applying the fast Fourier transform to a time series
of the $x$ coordinate. We compute the FTLE throughout the state space in a
neighbourhood of the deterministic invariant sets over a mesh of $100\times
100$ points. In (A), we took $b=b_{crit}+2$. For comparison, we portray in (B)
the FTLE landscape for $b=b_{crit}-1$, for which the deterministic dynamics
admits a globally attracting equilibrium.

\begin{overpic}[width=433.62pt]{FTLE_series3.png}
\put(50.0,-3.0){\scriptsize$T$} \put(-3.0,35.0){\scriptsize$\lambda^{T}$}
\end{overpic} (a)

\begin{overpic}[width=433.62pt]{FTLE_series2.png}
\put(50.0,-3.0){\scriptsize$T$} \put(-2.0,35.0){\scriptsize$\lambda^{T}$}
\end{overpic} (b)

\begin{overpic}[width=433.62pt]{FTLEseries4.png}
\put(50.0,-3.0){\scriptsize$T$} \put(-2.0,35.0){\scriptsize$\lambda^{T}$}
\end{overpic} (c)

Figure 5. Here the FTLE is plotted as a function of $T$ for a fixed initial
condition, considering $\sigma=0.1$, $a=1$, and a time step $h=10^{-3}$. In
(A) we took $b=b_{crit}+1$, while in (B) we considered $b=b_{crit}+7$. The
time windows in which the FTLE explores the positive real line is considerably
larger in the latter case. As a comparison, in (C) we took $b=b_{crit}-1$, for
which the system admits a globally attracting equilibrium.

##  4\. A slow-fast structure for the Brusselator

Consider again the deterministic system (2.7), with a change of coordinates
(see also [23]) given by

(4.1) $\left(\begin{array}[]{c}x\\\
y\end{array}\right)\mapsto\left(\begin{array}[]{c}y\\\ x+y\end{array}\right)$

We thus obtain the SDE

(4.2)
$\begin{array}[]{rcl}\mathrm{d}x&=&\left[b(y-x)-x(y-x)^{2}\right]\mathrm{d}t+\sigma(y-x)\circ\mathrm{d}W_{t},\\\
\mathrm{d}y&=&\left[a-(y-x)\right]\mathrm{d}t,\end{array}$

whose state space is now given by the set $\\{(x,y):x\geq 0\ y\geq x\\}$.
Based on the observations made in the previous section, we are interested in
the case $b\gg b_{crit}$. Hence, we write $b=\frac{a}{\epsilon}$ with
$\epsilon\ll 1$, so that system (4.2) reads as

(4.3) $\begin{array}[]{rcl}\epsilon\mathrm{d}x&=&\left[a(y-x)-\epsilon
x(y-x)^{2}\right]\mathrm{d}t+\epsilon\sigma(y-x)\circ\mathrm{d}W_{t}\\\
\mathrm{d}y&=&\left[a-(y-x)\right]\mathrm{d}t.\end{array}$

Observe that (4.3) has the structure of a singularly perturbed system in
standard form, representing the dynamics on a slow time scale. By considering
the fast time scale $\tau=t/\epsilon$, and recalling that
$\tilde{W}_{\tau}=\epsilon^{-1/2}W_{\epsilon\tau}$ is again a Wiener process
for all $\epsilon$, we have the equivalent system (in the distributional
sense)

(4.4) $\begin{array}[]{rcl}\mathrm{d}x&=&\left[a(y-x)-\epsilon
x(y-x)^{2}\right]\mathrm{d}\tau+\sqrt{\epsilon}\sigma(y-x)\circ\mathrm{d}\tilde{W}_{\tau}\\\
\mathrm{d}y&=&\epsilon\left[a-(y-x)\right]\mathrm{d}\tau.\end{array}$

By taking $\epsilon=0$ in (4.3), we obtain the reduced problem

(4.5) $\begin{array}[]{rcl}0&=&a(y-x),\\\
\mathrm{d}y&=&\left[a-(y-x)\right]\mathrm{d}t,\end{array}$

which is, in fact, a deterministic algebraic-differential equation. This
yields the _critical manifold_

$\mathcal{S}_{0}=\\{(x,y)\;|\;x=y\\},$

see also Figure 6.

###  4.1. Deterministic dynamics for $\epsilon\ll 1$

As a first approximation, we sketch the dynamics in the deterministic system
by considering $\sigma=0$ in (4.3)–(4.4), that is we consider the
deterministic slow-fast system in standard form

(4.6)
$\begin{array}[]{cc}\begin{array}[]{rcl}\epsilon\dot{x}&=&a(y-x)-\epsilon
x(y-x)^{2}\\\
\dot{y}&=&a-(y-x)\end{array}&\begin{array}[]{rcl}x^{\prime}&=&a(y-x)-\epsilon
x(y-x)^{2}\\\ y^{\prime}&=&\epsilon[a-(y-x)],\end{array}\end{array}$

where $\dot{z}:=dz/dt$ and $z^{\prime}:=dz/d\tau$. In the next proposition we
show that $\mathcal{S}_{0}$ is partially attracting.

######  Proposition 4.1.

The critical manifold $\mathcal{S}_{0}$ is normally hyperbolic. Its stable
foliation has leaves

$W_{s}(z)=span\\{(1,0)\\},$

for each $z\in\mathcal{S}_{0}$. Furthermore, the dynamics of (4.5) is given by

$x(t)=y(t)=x_{0}+at$

###### Proof.

Let $h(x,y)=(a(y-x),0)$. Its Jacobian matrix is calculated as

$Dh(x,y)=\left[\begin{array}[]{cc}-a&a\\\ 0&0\end{array}\right],$

for which the nontrivial eigenvalue is simply $\lambda=-a<0$, with eigenvector
$(1,0)$. By considering $x=y$, then $\mathrm{d}y/\mathrm{d}t=a$ and the result
follows. ∎

Fenichel’s theorem (see for instance [21, Theorem 2.4.2]) implies that for any
compact submanifold $\hat{\mathcal{S}}_{0}$ of $\mathcal{S}_{0}$ there is
$\epsilon_{0}$ sufficiently small such that for each
$\epsilon\in(0,\epsilon_{0})$ there is a normally hyperbolic locally invariant
manifold $\hat{\mathcal{S}}_{\epsilon}$ which is an
$\mathcal{O}(\epsilon)$-perturbation of $\mathcal{S}_{0}$. Since
$\mathcal{S}_{0}$ is normally hyperbolic, we can take $\hat{\mathcal{S}}_{0}$
arbitrarily long. Yet, the flow makes a turn and traverses the nullcline
$\dot{x}=0$, namely

$\mathcal{N}=\left\\{(x,y):y=x+\frac{a}{\epsilon x}\right\\}.$

Since the nullcline has $\mathcal{S}_{0}$ as its asymptote as
$x\rightarrow\infty$, whenever $\epsilon\ll 1$ the point in which the flow
traverses the nullcline is still very close to $\mathcal{S}_{0}$, but its
$x$-value is large. Hence, the dominating term for $x^{\prime}$ is the cubic
term $-\epsilon x(y-x)^{2}$ (see (4.6)), and the flow heads to the left until
it reaches $\mathcal{N}$ again. After this, the flow follows $\mathcal{N}$
closely with a slow time scale, up to a turning point where it takes a fast
transition back to a neighborhood of $\mathcal{S}_{0}$, and the cycle is
reapeated.

It is worth noticing that due to Proposition 4.1, no loss of hyperbolicity
anticipates the turning from $\mathcal{S}_{0}$. The fact that the nullcline
$\mathcal{N}$ approaches $\mathcal{S}_{0}$ as $x\rightarrow\infty$ suggests
that the analysis should be performed near “infinity”, by compactifying the
state space. To the best of our knowledge, no rigorous singular perturbation
analysis of the Brusselator has been done before, for which we refer to [13],
where by compactifying the state space, the “line at infinity” becomes a
nonhyperbolic invariant manifold, which can be desingularized using an
appropriate rescaling of the time variables. Similar arguments have been
implemented in [20, 16].

###  4.2. Stochastic Brusselator for $\epsilon\ll 1$.

Let us sketch the stochastic dynamics of (4.3)–(4.4) for $\sigma>0$. First,
notice that the noisy perturbation in (4.3) vanishes precisely on
$\mathcal{S}_{0}$, which is an attracting normally hyperbolic invariant
manifold for the deterministic system. It is intuitively clear that the
trajectories of (4.3) starting close to $\mathcal{S}_{0}$, will remain so for
some long time, whenever both $\epsilon$ and $\sigma$ are sufficiently small.
Similarly to the deterministic case, our numerical explorations suggest that
if a trajectory which is following closely the critical manifold
$\mathcal{S}_{0}$ traverses the nullcline $\mathcal{N}$, it enters the fast
regime and moves practically horizontally to the left until reaching
$\mathcal{N}$ again. The process then fluctuates around $\mathcal{N}$ until
being sufficiently close to $\mathcal{S}_{0}$ again, see Figure 6.

Observe that this description matches completely with the one given in Section
2 for Figure 2 by reverting the change of coordinates (4.1). A matching
between the corresponding time scales in systems (2.7) and (4.3) is portrayed
in Figure 6.

\begin{overpic}[width=433.62pt]{Orbits.png} \put(50.0,-3.0){\scriptsize$x$}
\put(-3.0,35.0){\scriptsize$y$} \end{overpic} (a)

\begin{overpic}[width=433.62pt]{Orbits_zoomright.png}
\put(50.0,-3.0){\scriptsize$x$} \put(-2.0,35.0){\scriptsize$y$} \end{overpic}
(b)

Figure 6. Two different orbits of system (4.3) with initial conditions $(2,2)$
and $(2.1,2.1)$ are obtained numerically and presented in (A) in blue and
yellow, respectively. The critical manifold $\mathcal{S}_{0}$ is portrayed in
red, the vertical nullcline in blue, and the deterministic limit cycle in
black. A closer look at the turning points is shown in (B). The parameters
taken here are $a=1$, $b=b_{crit}+1$, and $\sigma=0.2$

#### A mechanism for noise-induced (finite-time) instabilities

At last, let us sketch the mechanism by which the Brusselator exhibits finite-
time instabilities. We refer the reader to Figure 6 for the time scale regimes
of the state space.

Let us consider $(x_{0},y_{0})$ and $(x_{1},y_{1})$ two initial conditions for
the system (4.3), which are close to each other and also sufficiently close to
$\mathcal{S}_{0}$. Let $\tau_{\mathcal{N}}(\omega,x_{i},y_{i})$ be the first
hitting time of the trajectory $\varphi_{\omega}^{t}(x_{i},y_{i})$ ($i=0,1$)
to the nullcline $\mathcal{N}$, so that for
$t<\min\\{\tau_{\mathcal{N}}(\omega,x_{0},y_{0}),\tau_{\mathcal{N}}(\omega,x_{1},y_{1})\\}$,
it yields

$d_{t}((x_{0},y_{0}),(x_{1},y_{1}))\approx\|(x_{0},y_{0})-(x_{1},y_{1})\|^{2}.$

In other words, the ultraslow regime (I) close to $\mathcal{S}_{0}$ is of a
neutral nature, and no finite-time instabilities are observed in such regime.

It is worth noticing that for $\epsilon$ sufficiently small, in Figure 6 the
transitions from the ultraslow regime (I) to the subsequent fast regime (II)
may occur at different altitudes and at different times, generating an almost
instantaneous deviation between the trajectories if $\epsilon\ll 1$. Indeed,
the trajectories $\varphi_{\omega}^{t}(x_{i},y_{i})$ ($i=0,1$) may traverse
the deterministic nullcline $\mathcal{N}$ in points which are sufficiently
away from each other, and therefore at different times. When one of them, say
$\varphi^{t}_{\omega}(x_{0},y_{0})$, enters the fast regime first, it moves
quickly towards $\mathcal{N}$ while $\varphi^{t}_{\omega}(x_{1},y_{1})$ may
remain still close to $\mathcal{S}_{0}$.

The slow regime (III) and the fast regime (IV) are the only parts of the state
space where the system might present noise-induced stabilization. Not only a
rigorous analysis of the FTLEs on the regimes (I)-(IV) would give an
indication of the finite-time instabilities that the system exhibits, but it
may be of use in order to prove or disprove the conjecture that the classical
Lyapunov exponent is negative.

\begin{overpic}[width=433.62pt]{scales2.png} \put(50.0,-3.0){\scriptsize$x$}
\put(-3.0,35.0){\scriptsize$y$} \put(35.0,65.0){\tiny(II) fast}
\put(9.0,24.0){\tiny\begin{turn}{90.0} (III) slow \end{turn}}
\put(20.0,11.0){\tiny(IV) fast} \put(50.0,21.0){\tiny\begin{turn}{45.0} (I)
ultra slow \end{turn}} \par\end{overpic} (a)

\begin{overpic}[width=433.62pt]{scales_original2.png}
\put(50.0,-3.0){\scriptsize$x$} \put(-2.0,35.0){\scriptsize$y$}
\put(9.0,34.0){\tiny\begin{turn}{90.0} (I) ultra slow\end{turn}}
\put(46.0,56.0){\tiny\begin{turn}{-41.0}(II) fast\end{turn}}
\put(52.0,10.0){\tiny(III) slow} \put(14.0,22.0){\tiny\begin{turn}{-18.0}(IV)
fast\end{turn}} \end{overpic} (b)

Figure 7. The time scale separation in the stochastic Brusselator. In (A), an
example of an orbit of the system (4.3) is portrayed, following the regimes
(I), (II), (III), and (IV) in that order (counterclockwise). In (B), we
present a prototypical trajectory in the original coordinates (clockwise), see
(2.7). Each regime (I)-(IV) is identified via the change of coordinates (4.1)

##  5\. Concluding remarks and outlook

In this work we have studied finite-time instabilities for a stochastic
Brusselator via FTLEs. In [2] it was conjectured that parametric noise
destroys the Hopf bifurcation, in the sense that for any $\sigma>0$ the system
(2.7) admits a unique random equilibrium which is exponentially attracting.
While this conjecture remains open, our numerical explorations show that
finite-time instabilities occur in arbitrarily large time intervals as the
parameter $b$ grows.

In order to understand better the role of large $b$ already for the
deterministic system, we have sketched a suitable singular perturbation
analysis in Section 4, introducing the parameter of time scale separation
$\varepsilon\approx 1/b$. In [13], under suitable rescalings of $x,y$ and $t$,
we study the existence of a critical cycle for the deterministic system from
which the attracting periodic orbits are
$\mathcal{O}(\epsilon)$-perturbations.

We consider it possible that an appropriate scaling of both the parameters $b$
and $\sigma$ may lead to noise-induced chaos (cf. [5]), giving a negative
answer to the conjecture posed in [2]. Along these lines, there exist
extensions of singular perturbation analysis in the stochastic setting, see
for instance [21, Chapter 15]. This analysis for the stochastic Brusselator is
left for future work.

While we have focused on the stochastic Brusselator (2.7) defining a global
RDS, it is still possible to consider the chemical Langevin equation (2.6) and
perform a similar exploration, despite solutions not being defined for all
$t>0$. A suitable treatment of this situation may involve the notion of an
absorbed Markov process so that the process stops when one of the reactants
becomes zero. In such a context, stationary distributions cannot exist and
have to be replaced by the concept of quasi-stationary distributions, for an
introduction to the topic see [6]. On this basis, one can study changes of
stability via conditioned Lyapunov exponents [12].

At the level of the Markov jump process (2.1), the trajectories portrayed in
[15] resemble those obtained in Figure 2. As mentioned in [15], these noisy
oscillations have “…the apparent inability of the system to consistently
retrace its previous path on the diagonal leg of the limit cycle.” In our
setting, this property emerges as a consequence of the transition from a slow
regime into a fast regime of the system, as exhibited in Figure 6. To the best
of our knowledge, this point of view has not been studied for the Markov jump
process of the Brusselator reaction network, yielding another intriguing
direction for future work.

### Acknowledgements

We acknowledge the support of Deutsche Forschungsgemeinschaft (DFG) through
CRC 1114 and under Germany’s Excellence Strategy – The Berlin Mathematics
Research Center MATH+ (EXC-2046/1, project 390685689, in particular subproject
AA1-8). M. E. additionally thanks the DFG-funded SPP 2298 for supporting his
research. G. O.-M. also thanks FU Berlin for a 3-month Forschungsstipendium.
The authors gratefully acknowledge Peter Szmolyan and Bernd Krauskopf, for
insightful discussions.

## References

  * [1] L. Arnold, _Random dynamical systems_ , Springer, Berlin, 1998.
  * [2] L. Arnold, G. Bleckert, and K. R. Schenk-Hoppé, _The stochastic Brusselator: parametric noise destroys Hopf bifurcation_ , Springer, New York, NY, 2015.
  * [3] A. Blumenthal, M. Engel, and A. Neamtu, _On the pitchfork bifurcation for the Chafee-Infante equation with additive noise_ , arXiv:2108.11073 (2021).
  * [4] M. Callaway, T. S. Doan, J. S. W. Lamb, and M. Rasmussen, _The dichotomy spectrum for random dynamical systems and pitchfork bifurcations with additive noise_ , Ann. Inst. Henri Poincaré Probab. Stat. 53 (2017), no. 4, 1548–1574. MR 3729628
  * [5] D. Chemnitz and M. Engel, _Positive Lyapunov exponent in the Hopf normal form with additive noise_ , arXiv:2212.06547 (2022), To appear in Communications in Mathematical Physics.
  * [6] P. Collet, S. Martínez, and J. San Martín, _Quasi-stationary distributions. Markov chains, diffusions and dynamical systems_ , Springer Berlin, Heidelberg, 2012.
  * [7] H. Crauel and F. Flandoli, _Additive noise destroys a pitchfork bifurcation_ , J. Dyn. Differ. Equ. 10 (1998), 259––274.
  * [8] H. Crauel and P. E. Kloeden, _Nonautonomous and random attractors_ , Jahresbericht der Deutschen Mathematiker-Vereinigung 117 (2015), 173–206.
  * [9] S. Dey, D. Das, and P. Parmananda, _Intrinsic noise induced resonance in presence of sub-threshold signal in Brusselator_ , Chaos: An Interdisciplinary Journal of Nonlinear Science 21 (2011), no. 3.
  * [10] T.S. Doan, M. Engel, J.S.W. Lamb, and M. Rasmussen, _Hopf bifurcation with additive noise_ , Nonlinearity 31 (2018), no. 10, 4567–4601.
  * [11] M. Engel, J.S.W. Lamb, and M. Rasmussen, _Bifurcation analysis of a stochastically driven limit cycle_ , Comm. Math. Phys. 365 (2019), no. 3, 935–942.
  * [12] by same author, _Conditioned Lyapunov exponents for random dynamical systems_ , Trans. Am. Math. Soc. 372 (2019), no. 9, 6343–6370.
  * [13] M. Engel and G. Olicón-Méndez, _A singular perturbation analysis of the Brusselator_ , in preparation (2023).
  * [14] M. Engel, G. Olicón-Méndez, N. Unger, and S. Winkelmann, _Synchronization and random attractors for reaction jump processes_ , arXiv:2207.00602 (2022).
  * [15] D. T. Gillespie, _Exact stochastic simulation of coupled chemical reactions_ , J. Phys. Chem. 81 (1977), no. 25, 2340–2361.
  * [16] I. Gucwa and P. Szmolyan, _Geometric singular perturbation analysis of an autocatalator_ , Discrete Contin. Dyn. Syst.–S 2 (2009), 783–806.
  * [17] P. Imkeller and C. Lederer, _An explicit description of the Lyapunov exponents of the noisy damped harmonic oscillator_ , Dynam. Stability Systems 14 (1999), no. 4, 385–405. MR 1746113
  * [18] by same author, _Some formulas for Lyapunov exponents and rotation numbers in two dimensions and the stability of the harmonic oscillator and the inverted pendulum_ , Dyn. Syst. 16 (2001), no. 1, 29–61. MR 1835906
  * [19] P. Kloeden and E. Platen, _Numerical solution of stochastic differential equations_ , Springer Berlin, Heidelberg, 1992.
  * [20] C. Kuehn, _Normal hyperbolicity and unbounded critical manifolds_ , Nonlinearity 27 (2014), 1351––1366.
  * [21] by same author, _Multiple time scale dynamics_ , Springer Cham, 2015.
  * [22] S. C. Leite and R.J. Williams, _A constrained Langevin approximation for chemical reaction networks_ , Ann. Appl. Probab. 29 (2019), no. 3, 1541–1608.
  * [23] X. Li, J. Hou, and Y. Shen, _Slow-fast effect and generation mechanism of Brusselator based on coordinate transformation_ , Open Phys. 14 (2016), 261––26.
  * [24] K. K. Lin and L. S. Young, _Shear-induced chaos_ , Nonlinearity 21 (2008), no. 5, 899–922. MR 2412320
  * [25] K. Lu, Q. Wang, and L.S. Young, _Strange attractors for periodically forced parabolic equations_ , vol. 224, Mem. Am. Math. Soc., 2013.
  * [26] Julian Newman, _Synchronisation of almost all trajectories of a random dynamical system_ , Discrete Contin. Dyn. Syst. 40 (2020), no. 7, 4163–4177. MR 4097539
  * [27] I. Prigogine and R. Lefever, _Symmetry breaking instabilities in dissipative systems. II._ , J. Chem. Phys. 48 (1968), 1665–1700.
  * [28] I. Prigogine and G. Nicolis, _Self-organisation in nonequilibrium systems_ , Wiley, 1977.
  * [29] J. Wilkie and Y. M. Wong, _Positivity preserving chemical Langevin equations_ , Chem. Phys. 353 (2008), 132–138.
  * [30] S. Winkelmann and C. Schütte, _Hybrid models for chemical reaction networks: Multiscale theory and application to gene regulatory systems_ , J. Chem. Phys. 147 (2017), no. 11, 114115.
  * [31] by same author, _Stochastic dynamics in computational biology_ , Springer, 2020.


# Effect of gallium termination on InGaAs wetting layer properties in droplet
epitaxy InGaAs quantum dots

D Fricker1,5, P Atkinson2, M Lepsa1,3, Z Zeng1,5, A Kovács4,
L Kibkalo4, RE Dunin-Borkowski4,5, and BE Kardynał1,5 1 Peter Grünberg
Institute 9, Forschungszentrum Jülich, 52425 Jülich, Germany  2 Institut des
Nano Sciences de Paris, CNRS UMR 7588, Sorbonne Université, 75005 Paris,
France  3 Peter Grünberg Institute 10, Forschungszentrum Jülich, 52425 Jülich,
Germany  4 Ernst Ruska-Centre for Microscopy and Spectroscopy with Electrons,
Peter Grünberg Institute 5, Forschungszentrum Jülich, 52428 Jülich, Germany  5
Department of Physics, RWTH Aachen University, 52074 Aachen, Germany
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Self-assembled quantum dots based on III-V semiconductors have excellent
properties for applications in quantum optics. However, the presence of a 2D
wetting layer which forms during the Stranski-Krastanov growth of quantum dots
can limit their performance. Here, we investigate wetting layer formation
during quantum dot growth by the alternative droplet epitaxy technique. We use
a combination of photoluminescence excitation spectroscopy, lifetime
measurements, and transmission electron microscopy to identify the presence of
an InGaAs wetting layer in these droplet epitaxy quantum dots, even in the
absence of distinguishable wetting layer photoluminescence. We observe that
increasing the amount of Ga deposited on a GaAs (100) surface prior to the
growth of InGaAs quantum dots leads to a significant reduction in the emission
wavelength of the wetting layer to the point where it can no longer be
distinguished from the GaAs acceptor peak emission in photoluminescence
measurements.

Keywords: droplet epitaxy, InAs/GaAs, wetting layer properties, wetting layer
characterization

##  1 Introduction

The field of semiconductors as a platform for non-classical light sources has
been rapidly evolving over the years as a result of growing interest in
quantum communication. Photoluminescence from a single InAs/GaAs quantum dot
(QD) grown by molecular beam epitaxy (MBE) was first reported in 1994 [1].
These QDs were grown in the Stranski-Krastanov (SK) growth mode, in which
quantum dot formation is driven by strain relaxation of the InAs layer when
the layer thickness exceeds a critical thickness [2]. These SK QDs have been a
workhorse of semiconductor quantum optics and used extensively for the
fabrication of light sources for quantum communication [3]. Since then,
numerous remarkable optical phenomena have been observed in semiconductor-
based QDs as their fabrication has been continually improved to achieve better
control over QD formation [4]. However, one disadvantage of the SK growth mode
is the presence of a wetting layer (WL), which leads to decoherence processes
in the QDs that are highly undesired in some applications [5, 6, 7]. One
method to decrease the influence of this wetting layer on the QDs is to use an
AlGaAs capping layer which eliminates confined electronic states in the
conduction band in the wetting layer [7].
In recent years, however, a different growth mode, the droplet epitaxy (DE)
technique, has been developed, yielding QDs with comparable optical properties
as SK-grown QDs [8, 9]. This technique has also been reported to result in
wetting-layer free QDs [10].
In this report, we show that the emission wavelength of the WL can be
controlled by the amount of Ga used to prepare a Ga-terminated surface prior
to the growth of InGaAs dots by droplet epitaxy. We study wetting layer
electronic states using photoluminescence excitation (PLE) and time-resolved
photoluminescence (TR PL). This has allowed the presence of a WL to be
identified even when the WL emission cannot be clearly identified in
photoluminescence (PL) spectroscopy. The In distribution in the wetting layers
grown on different Ga-terminated surfaces is studied using transmission
electron microscopy (TEM).

##  2  Epitaxy

The investigated samples were grown on undoped (100) GaAs substrates in a
molecular beam epitaxy system from MBE Komponenten, using an arsenic valve
cracker cell for $As_{4}$ flux control and a pyrometer to measure the
substrate temperature. After thermal deoxidation of the substrate, a 500 nm
thick GaAs buffer layer was grown at 630 ∘C with a growth rate of 1 monolayer
(ML)/s. A 120 nm thick AlAs/GaAs superlattice was incorporated in the middle
of the GaAs buffer layer in order to smooth the surface. Following buffer
layer growth, the surface was prepared for QD growth. Firstly, the substrate
temperature was reduced to 575 ∘C and then the $As_{4}$ valve was closed for 5
minutes. Afterwards, the substrate was slowly cooled to 350 ∘C over 20 minutes
to allow the residual arsenic in the chamber to be pumped away. When the
chamber background pressure reached 8$\times$10-9 mbar, 1.8 ML of Ga was
deposited on the GaAs surface with a low flux of 0.06 ML/s, which should
result in a Ga-terminated GaAs surface [11, 12, 13]. On this surface 1.4 ML of
In was then deposited with a flux of 0.045 ML/s to form In droplets. The Ga
and In growth rates and layer thicknesses are given as equivalent for GaAs and
InAs, respectively. The droplets were immediately crystallized with an
$As_{4}$ beam flux of 2$\times$10-5 mbar to form InGaAs QDs. Simultaneously,
slow heating of the wafer was started. Then, after 30 minutes, when the
substrate temperature reached 530 ∘C the QDs were capped with a 2 nm thick
GaAs capping layer, followed by the growth of a 98 nm thick GaAs layer at 580
∘C. Both capping layers were grown with a growth rate of 1 ML/s.
The deposition amounts given here (1.8 ML Ga, 1.4 ML In) correspond to the
deposition amount at the wafer centre. However, the wafers were grown without
substrate rotation during the Ga and In deposition steps, resulting in an
expected $\pm 25\,\%$ variation in deposition amount across a 4-inch wafer
[14]. After Ga deposition, the substrate was rotated to ensure that the In
flux gradient was at 90∘ to the Ga flux gradient. This technique allows the
effect of different Ga amounts during the Ga-termination step, without
changing other growth parameters, to be studied on a single wafer.
To check that the observed variation across the wafer was not due to a
temperature gradient, a second wafer with identical growth conditions but a
nominal Ga amount of 2.16 ML was grown. Later, it will be shown that the
emission near the centre of the 2.16 ML Ga wafer (wafer 2) overlaps well with
the emission near the edge of the 1.8 ML Ga wafer (wafer 1) where the Ga
amount was $\sim 2.2\,ML$ Ga. This indicates that a temperature gradient from
the centre to the edge of the wafer cannot account for the observed shift in
emission wavelength.

##  3 Sample characterization

###  3.1 Optical characterization

Mapping of the PL across the wafer was carried out with meV resolution, using
a cold finger He cryostat mounted on a manual stage allowing large movements
of the sample.
More detailed $\mu$-PL, PLE and TR PL measurements were carried out in a
closed-cycle cryostat, where the samples were mounted on a computer-controlled
xyz-stage offering sub-$\mu m$ resolution. An achromatic lens with a numerical
aperture (NA) of 0.81 focused the laser beam on the sample and collected the
resulting PL signal. A continuous wave (CW) Ti:sapphire laser with a tunable
wavelength was used to perform the $\mu$-PL and PLE measurements.
To perform the TR PL measurements, a pulsed diode laser with a tuneable pulse
repetition rate up to 80 MHz and an emission wavelength of 660 nm was used.
The laser pulse width was $\sim 200\,ps$. The emitted PL was spectrally
selected with wavelength-tunable filters and guided to a single-photon
avalanche photodiode (SPAD), which has 30 ps time resolution.

###  3.2 Structural characterization

Structural characterization of the samples was carried out using high-angle
annular dark-field (HAADF) scanning TEM (STEM) imaging. Cross-sectional
specimens were prepared using focused Ga ion beam (FIB) sputtering in a dual
beam scanning electron microscope. The surface damage was reduced by low-
energy ($<1$ keV) Ar ion beam sputtering [15]. The images were recorded using
an electron probe aberration corrected microscope operated at 200 kV [16]. The
inner annular dark-field detector semi-angle used was 69 mrad, resulting in
HAADF imaging, which is sensitive to chemical composition as the image
intensity scales approximately as $I\sim Z^{2}$, where Z is the atomic number
($Z_{In}=49,Z_{Ga}=31$).

Figure 1: (a) Normalize$\mu$-PL spectra from three different positions on
wafer 1 along the Ga gradient. (b) Normalized wetting layer signal measured on
wafer 1 along the Ga gradient in 1 mm steps. (c) Central wavelength of the
Gaussian function fitted to the signal assigned to wetting layer emission for
wafers 1 and 2. Wafer 1 had a deposition of 1.8 ML Ga at the wafer centre,
while wafer 2 had a deposition of 2.16 ML Ga at the centre. Both wafers had
the same deposition amount of In (1.4 ML) at the centre.

##  4 Results and discussion

Figure 1(a) shows normalized $\mu$-PL spectra obtained at different positions
on the wafer along the Ga gradient. The positions correspond to Ga amounts of
1.4 ML, 1.8 ML and 2.2 ML. Thanks to the use of the gradient technique during
QD layer formation, these three spectra were taken from different positions
along the Ga gradient from wafer 1. The deposited InAs amount was 1.4 ML. The
two PL peaks at 818 nm and 830 nm originate from radiative recombination of
the GaAs free exciton and the free electron-neutral acceptor ($(A^{0},e)$),
respectively [17, 18]. The positions of these GaAs PL lines are, as expected,
independent of the Ga amount deposited during the Ga-termination step.
A decrease in emission wavelength of the WL peak can be observed when the Ga
amount increases from 1.4 ML to 1.8 ML. When the Ga amount is 2.2 ML, the
wetting layer peak can no longer be resolved in the PL spectrum. PL
measurements in fine steps along the Ga gradient track the development of the
WL emission wavelength with Ga amount in more detail. The spectra shown in
Figure 1(b) were recorded in 1 mm steps along the Ga-gradient on wafer 1, with
the deposited In amount fixed at 1.4 ML InAs. Each spectrum in Figure 1(b) was
fitted with Gaussian functions. The central wavelengths of the fitted
functions are plotted as black points in Figure 1(c). The position on the
wafer is converted into the expected amount of Ga. The WL emission wavelength
continuously decreases from 847 to 834 nm, with the Ga amount increasing from
1.35 ML to 1.85 ML. At around 1.87 ML the WL emission wavelength changes
abruptly from 834 nm to 836 nm, and remains in the 836$\pm$ nm range until no
WL is observed for larger amounts of Ga of 2.2 ML. The central wavelength of
WL emission on wafer 2 is shown in red in Figure 1(c) and shows similar
behaviour but with no abrupt change in WL emission wavelength at any point.
There is a good overlap of the WL emission for wafers 1 and 2 for the Ga
deposition amount between 1.9-2.1 ML, despite the fact that these Ga
deposition amounts occur at different spatial positions on the two wafers.
This means that any effect of a difference in substrate temperature with
position on the wafer can be neglected, so that the effect of the Ga
deposition amount on WL formation is considered in the following discussion
only.
There are three possible explanations for the absence of a WL PL signal.
Firstly, it is possible that the formation of the WL has been suppressed by
the Ga deposition amount. It has previously been demonstrated that, on an
arsenic-terminated surface, the first 0.75-1.75 ML of group III deposition
bonds to excess arsenic on the surface, forming a wetting layer and droplet
formation only starts to occur once the surface is metal-rich [12]. However,
under our growth conditions and with the residual arsenic overpressure in our
MBE chamber, it is possible that the transition to a fully Ga-terminated
surface was achieved [13] only once 2.2 ML Ga was deposited. It could be
expected that all of the deposited In on a Ga-terminated surface forms
droplets and no InAs WL is formed.
A second possible explanation is that the wetting layer emission persists but
overlaps with the acceptor signal, so that differentiation in PL measurements
is not possible for samples with more than 2.2 ML of deposited Ga. The two PL
signals can however be distinguished by time-resolved PL, since WL and
acceptor-mediated PL occur with different decay times.
The third possible explanation for the absence of the WL PL signal is that
there is rapid relaxation of photo-excited carriers from the wetting layer
into the QDs, so that PL signal from the wetting layer is not observable. The
presence of a wetting layer in the latter case can be probed using PLE.

###  4.1 PLE

PLE measurements were performed with a wavelength-tunable Ti:Sapphire laser in
the $\mu$-PL set-up. The PLE data are shown in Figure 2 in black. The data
represent the normalized QD emission intensity as a function of excitation
wavelength on three different positions on the wafer. At each position, the
PLE spectrum of several QDs, with QD emission wavelength ranging from 855 nm
to 890 nm, were studied. No significant differences were observed in their PLE
spectra, so the data are shown for a single representative QD. The PL of the
sensing QD is shown in the inset. PLE spectra are shown together with PL
spectra from the same position and spectral range, measured using an
excitation wavelength of 780 nm.
Figure 2(a) shows the data for 1.4 ML of deposited Ga amount from wafer 1. The
PLE data demonstrate a clear increase in QD intensity when the laser light is
resonant with the GaAs bandgap at 818 nm and with the wetting layer states.
Among the wetting layer states, both the light hole $WL_{LH}$ and the heavy
hole $WL_{HH}$ states are observed [19]. The $WL_{HH}$ state is also
observable in PL measurements and there is a Stokes shift of 2 nm (3.5 meV)
between the PLE and PL peak positions. Figures 2(b) and (c) show the data for
2.2 ML and 2.5 ML of Ga from wafer 2, respectively. For 2.2 ML, the PL from WL
emission at a wavelength of 833 nm partially overlaps with the $(A^{0},e)$
signal, which can be seen from fitting with two Gaussian functions. The 2.5 ML
sample does not show a detectable wetting layer in the PL spectrum. The PLE
data nevertheless show a clear WL peak at 833 nm. Thus, it can be shown that a
WL exists even though there is no distinguishable WL PL signal. Both the
Stokes shift between the PL and PLE signal of the $WL_{HH}$ states and the
wavelength difference between heavy and light hole WL states decreases with
the wavelength of the PLE $WL_{HH}$ signal, probably as a result of weaker
confinement of the different states within the wetting layer. Other PLE
measurements on wafers 1 and 2 at positions with different WL wavelengths show
the same effect as that presented here.

Figure 2: Normalized PL and PLE intensities for (a) 1.4 ML, (b) 2.2 ML and (c)
2.5 ML of deposited Ga amounts for Ga termination. The identified states for
the PL data are the GaAs free exciton, $(A^{0},e)$ and WL recombinations,
whereas the PLE data identifies GaAs, $WL_{LH}$ and $WL_{HH}$ states. The
inset of each figure shows the PL spectrum of the QD used for the PLE
measurements.

###  4.2 TR PL

Time-resolved photoluminescence measurements were performed in the $\mu$-PL
setup using a pulsed diode laser. Figure 3 shows lifetime measurements on two
samples with different Ga amounts: a 1.4 ML Ga sample from wafer 1, where the
acceptor (in black) and WL (in red) signals are spectrally well separated (as
seen in Figure 2(a)) and a 2.2 ML Ga sample from wafer 2 (in blue), where the
WL cannot be spectrally well separated by PL measurements (shown in Figure
2(b)). Exponential decay fitting to the 1.4 ML sample shows that the decay
time of $(A^{0},e)$ is $\tau_{(A^{0},e)}=13$ ns, which is significantly
greater than the wetting layer lifetimes. The latter is found by fitting
double exponential decay functions to be $\tau_{WL_{1}}=0.3$ ns and
$\tau_{WL_{2}}=1.2$ ns. The fast WL decay time may due to exciton relaxation
to the QD or other losses, while the slower WL decay time is likely to
correspond to the radiative emission from the WL, which is visible in the PL
measurement. Similar measurements (not shown here) at different positions on
wafers 1 and 2, where the WL can be resolved in the PL spectrum, result in
lifetimes in the range $\tau_{(A^{0},e)}=9-13$ ns, $\tau_{WL_{1}}=0.2-0.6$ ns
and $\tau_{WL_{1}}=1.2-1.8$ ns. It is therefore clear that the wetting layer
emission can be distinguished from the acceptor signal due to the significant
difference in decay time.
The fitting of a double-exponential decay function to the mixed $(A^{0},e)$+WL
emission of the 2.2 ML sample gives $\tau_{Mix_{1}}=0.3$ ns and
$\tau_{Mix_{2}}=1.3$ ns. Both lifetimes extracted from the decay of the single
emission line seen in the 2.2 ML sample ,$\tau_{Mix}$ (0.3 ns and 1.3 ns ),
are similar to the lifetimes of the WL radiative emission seen in the 1.4 ML
sample, $\tau_{WL}$ (0.3 ns and 1.2 ns). In the mixed state the long lifetime
of the acceptor excitation $\tau_{(A^{0},e)}$ (13 ns) could not be found.
Further, measurements (not shown here) result in lifetimes of mixed states in
the range $\tau_{Mix_{1}}=0.2-0.4$ ns and $\tau_{Mix_{2}}=1.2-1.6$ ns, which
is in the range of decay times for resolved WL emission for Ga deposition
below 2.2 ML and is one order of magnitude smaller than that for $(A^{0},e)$
emission. These observations indicate that, for Ga deposition larger than 2.2
ML, where only the acceptor peak appears to be visible in the PL, this peak is
in fact composed of emission from $(A^{0},e)$ and the WL.

Figure 3:  Normalized lifetime measurements for a 1.4 ML Ga sample with
spectrally separable $(A^{0},e)$ (black) and WL (red) PL signals with
corresponding exponential fits and for the 2.2 ML Ga with spectrally
overlapping $(A^{0},e)$ and WL signals (blue), respectively. The signal that
contains WL emission is fitted with a double exponential decay and shows
significantly shorter lifetimes than pure $(A^{0},e)$ emission.

###  4.3 TEM

Three samples, with deposited Ga amounts of 1.8 ML Ga, 2.16 ML Ga and 2.5 ML
Ga, were investigated with TEM. Both samples with higher deposited Ga amounts
are from wafer 2. The STEM imaging conditions were identical in studies of the
specimens. HAADF STEM images of the three samples are shown in Figure 4(a-c).
In all three samples, a continuous bright line is observed, indicating a
continuous InGaAs wetting layer. The sample with the lowest amount of
deposited Ga shows the weakest contrast of the wetting layer relative to the
surrounding GaAs. This may be due to a small difference in In deposition
amount between wafer 1 and 2 despite having nominally the same amount or
because of an artefact from the measurements and preparation. Figure 4(d)
shows the normalized line scan intensity (averaged over $\sim$ 70 nm along the
wetting layer) along the growth direction for all three samples. All
distributions show a sharp increase in In signal at the start of the In
deposition. The In segregates during overgrowth with GaAs, giving a more
gradual drop in In concentration along the growth direction. This profile can
be observed for all three samples. The thicknesses of these wetting layers are
the same for all Ga deposition amounts. Note that the intensity fluctuation in
the linescan of the sample with 1.8 ML Ga is above the noise in the TEM
measurements.
The data presented here provide strong evidence for the presence of a WL in
droplet-epitaxy-grown InAs dots on GaAs substrates, even in the case of Ga-
terminated surfaces. This may be due to In exchange with Ga atoms on the Ga-
terminated surface, resulting in a partially In-terminated surface during the
initial stages of In deposition prior to droplet formation. During the
recrystallisation step, this In-Ga-terminated surface would form an InGaAs
wetting layer. Alternatively, outdiffusion from the In droplet during the
recrystallisation step may result in the formation of a thin WL between the
dots. This kind of outdiffusion has been observed with GaAs droplets on AlGaAs
surfaces, leading to nano-disks with a diameter of up to several 100 nm
forming around the dots [20]. We have not seen evidence for this disk
formation when carrying out AFM measurements on surface droplet epitaxy dots
grown under the conditions investigated here. However, given the greater
diffusion length of In compared to Ga atoms on GaAs [21, 22], we cannot
exclude that this has occurred but that the ring is so large and small in
height that it has not been observed in AFM measurements.

Figure 4:  HAADF STEM images for (a) 1.8 ML, (b) 2.16 ML and (c) 2.5 ML of
deposited Ga amount. The bright horizontal line indicates in all samples a
continuous InGaAs layer, corresponding to the wetting layer. d) Normalized
line scans for the three samples recorded perpendicular to the WL in HAADF
STEM images.

##  5 Conclusions

We have demonstrated that the Ga deposition amount prior to In droplet epitaxy
is an effective tool to tune the WL emission wavelength. This emission
wavelength was reduced until the WL signal overlapped with the GaAs free-
electron-acceptor signal in PL measurements, such that these signals could not
be spectrally distinguished. However, the presence of WL absorption peaks
observed in PLE unambiguously revealed the WL present for all studied Ga
amounts. The presence of the WL was also confirmed by time-resolved
measurements of the emission decay lifetime, since the WL decay lifetime was
found to be an order of magnitude shorter than the free-electron acceptor
decay lifetime. These measurement techniques show a clear presence of a WL
signal for all investigated samples, indicating that the depositions amounts
investigated, 1.4 ML to 2.5 ML Ga deposition followed by 1.4 ML In deposited
at 350 ∘C did not result in the suppression of wetting layer formation during
In droplet epitaxy. This is in contrast with a previous report of WL-free InAs
droplet epitaxy dots [10] which may indicate the effect of subtle differences
in growth parameters (e.g., chamber background pressure during droplet
deposition) due to the different MBE equipment used. In addition, TEM
measurements demonstrated no significant change in the thickness of the WL for
different deposited Ga amounts. This observation shows that the shift in the
WL emission wavelength observed is mainly an effect of changing composition of
InGaAs within the WL along the Ga-gradient. This work demonstrates the
sensitivity of complementary PLE, TR and TEM measurements to identify the
presence of a WL in the absence of a WL emission signature in
photoluminescence measurements.

We thank members of the ML4Q Cluster for stimulating discussions and Christoph
Krause and Benjamin Bennemann for technical support in growing the sample.
Access to the Helmholtz Nano Facility (HNF) and the Ernst Ruska-Centre (ER-C)
at Forschungszentrum Jülich is acknowledged. We acknowledge funding from
Germany’s Excellence Strategy-Cluster of Excellence Matter and Light for
Quantum Computing (ML4Q) EXC 2004/1-390534769.

## References

## References

  * [1] J. Y. Marzin, J. M. Gérard, A. Izraël, D. Barrier, and G. Bastard.  Photoluminescence of single InAs quantum dots obtained by self-organized growth on GaAs.  Physical review letters, 73(5):716–719, August 1994.
  * [2] D. Leonard, K. Pond, and and P. M. Petroff.  Critical layer thickness for self-assembled InAs islands on GaAs.  Physical review B, 50(16), 1994.
  * [3] P Intallura, M Ward, Otabek Karimov, Zhiliang Yuan, P See, P Atkinson, David Ritchie, and A Shields.  Quantum communication using single photons from a semiconductor quantum dot emitting at a telecommunication wavelength.  Journal of Optics A: Pure and Applied Optics, 11:054005, 03 2009\.
  * [4] Peter Lodahl, Sahand Mahmoodian, and Soren Stobbe.  Interfacing single photons and single quantum dots with photonic nanostructures.  Reviews of Modern Physics, 87(2):347–400, 2015.
  * [5] L. Seravalli, G. Trevisi, P. Frigeri, S. Franchi, M. Geddo, and G. Guizzetti.  The role of wetting layer states on the emission efficiency of InAs/InGaAs metamorphic quantum dot nanostructures.  Nanotechnology, 20(27), 2009.
  * [6] Q. Q. Wang, A. Muller, P. Bianucci, E. Rossi, Q. K. Xue, T. Takagahara, C. Piermarocchi, A. H. MacDonald, and C. K. Shih.  Decoherence processes during optical manipulation of excitonic qubits in semiconductor quantum dots.  Physical Review B, 72(3), 2005.
  * [7] Matthias C. Löbl, Sven Scholz, Immo Söllner, Julian Ritzmann, Thibaud Denneulin, András Kovács, Beata E. Kardynał, Andreas D. Wieck, Arne Ludwig, and Richard J. Warburton.  Excitons in InGaAs quantum dots without electron wetting layer states.  Communications Physics, 2(1), 2019.
  * [8] Massimo Gurioli, Zhiming Wang, Armando Rastelli, Takashi Kuroda, and Stefano Sanguinetti.  Droplet epitaxy of semiconductor nanostructures for quantum photonic devices.  Nature materials, 18(8):799–810, 2019.
  * [9] M. Anderson, T. Müller, J. Huwer, J. Skiba-Szymanska, A. B. Krysa, R. M. Stevenson, J. Heffernan, D. A. Ritchie, and A. J. Shields.  Quantum teleportation using highly coherent emission from telecom C-band quantum dots.  npj Quantum Information, 6(1), 2020.
  * [10] Joanna Skiba-Szymanska, R. Mark Stevenson, Christiana Varnava, Martin Felle, Jan Huwer, Tina Müller, Anthony J. Bennett, James P. Lee, Ian Farrer, Andrey B. Krysa, Peter Spencer, Lucy E. Goff, David A. Ritchie, Jon Heffernan, and Andrew J. Shields.  Universal Growth Scheme for Quantum Dots with Low Fine-Structure Splitting at Various Emission Wavelengths.  Physical Review Applied, 8(1):1–8, 2017.
  * [11] Takaaki Mano, Katsuyuki Watanabe, Shiro Tsukamoto, Hiroshi Fujioka, Masaharu Oshima, and Nobuyuki Koguchi.  New Self-Organized Growth Method for InGaAs Quantum Dots on GaAs(001) Using Droplet Epitaxy.  Japanese Journal of Applied Physics, 38(Part 2, No. 9A/B):L1009–L1011, 1999.
  * [12] S. Sanguinetti, K. Watanabe, T. Tateno, M. Gurioli, P. Werner, M. Wakaki, and N. Koguchi.  Modified droplet epitaxy GaAs/AlGaAs quantum dots grown on a variable thickness wetting layer.  Journal of Crystal Growth, 253(1-4):71–76, 2003.
  * [13] Masafumi Jo, Takaaki Mano, and Kazuaki Sakoda.  Two-Step Formation of Gallium Droplets with High Controllability of Size and Density.  Crystal Growth $\&$ Design, 11(10):4647–4651, 2011.
  * [14] Flux distribution simulations carried out by MBE Komponenten.
  * [15] Paul E. Fischione, Robert E.A. Williams, Arda Genç, Hamish L. Fraser, Rafal E. Dunin-Borkowski, Martina Luysberg, Cecile S. Bonifacio, and András Kovács.  A Small Spot, Inert Gas, Ion Milling Process as a Complementary Technique to Focused Ion Beam Specimen Preparation.  Microscopy and Microanalysis, 23(4):782–793, 2017.
  * [16] András Kovács, Roland Schierholz, and Karsten Tillmann.  FEI Titan G2 80-200 CREWLEY.  Journal of large-scale research facilities JLSRF, 2, 2016.
  * [17] Kazuhiro Kudo, Yunosuke Makita, Ichiro Takayasu, Toshio Nomura, Toshihiko Kobayashi, Tomio Izumi, and Tokue Matsumori.  Photoluminescence spectra of undoped GaAs grown by molecular-beam epitaxy at very high and low substrate temperatures.  Journal of Applied Physics, 59(3):888–891, 1986.
  * [18] D. Bimberg, H. Munzel, A. Steckenborn, and and J. Christen.  Kinetics of relaxation and recombination of nonequilibrium carriers in GaAs: Carrier capture by impurities.  Physical Review B, 31(12), 1985.
  * [19] E. S. Moskalenko, K. F. Karlsson, P. O. Holtz, B. Monemar, W. V. Schoenfeld, J. M. Garcia, and P. M. Petroff.  Acceptor-induced threshold energy for the optical charging of InAs single quantum dots.  Physical Review B, 66(19), 2002.
  * [20] Sergio Bietti, Claudio Somaschini, Luca Esposito, Alexey Fedorov, and Stefano Sanguinetti.  Gallium surface diffusion on GaAs (001) surfaces measured by crystallization dynamics of Ga droplets.  Journal of Applied Physics, 116(11):114311, 2014.
  * [21] A. Wakejima.  Lateral variation of indium content in InGaAs grown on GaAs channeled substrates by molecular beam epitaxy.  Journal of Vacuum Science $\&$ Technology B: Microelectronics and Nanometer Structures, 12(2):1102, 1994.
  * [22] S. Koshiba, Y. Nakamura, M. Tsuchiya, H. Noge, H. Kano, Y. Nagamune, T. Noda, and H. Sakaki.  Surface diffusion processes in molecular beam epitaxial growth of GaAs and AlAs as studied on GaAs (001)–(111) B facet structures.  Journal of Applied Physics, 76(7):4138–4144, 1994.


Victoria Toptun, Igor Chilingarian, Kirill Grishin, Ivan Katkov, Ivan
Zolotukhin, Vladimir Goradzhanov, Mariia Demianenko, and Ivan Kuzmun
11institutetext: Sternberg Astronomical Institute, Lomonosov Moscow State
University, Universitetsky pr. 13, Moscow 119234, Russia
22institutetext: Center for Astrophysics - Harvard and Smithsonian, 60 Garden
Street MS09, Cambridge, MA 02138, USA
33institutetext: Department of Physics, M.V. Lomonosov Moscow State
University, 1 Vorobyovy Gory, Moscow 119991, Russia
44institutetext: Center for Astro, Particle, and Planetary Physics, NYU Abu
Dhabi, PO Box 129188 Abu Dhabi, UAE
55institutetext: Moscow Institute of Physics and Technology, Institutskiy per.
9, Dolgoprudny 141701, Russia

# Confirmation of intermediate-mass black holes candidates with X-ray
observations

Victoria Toptun 1133 Igor Chilingarian 1122 Kirill Grishin 11 Ivan Katkov 1144
Ivan Zolotukhin 11 Vladimir Goradzhanov 1133 Mariia Demianenko 1155
Ivan Kuzmun 33

###### Abstract

The origin of supermassive black holes (SMBH) in galaxy centers still remains
uncertain. There are two possible ways of their formation - from massive
($10^{5}-10^{6}$ $M_{\sun}$ ) and low-mass (100 $M_{\sun}$ ) BH nuclei. The
latter scenario should leave behind a large number of intermediate mass black
holes (IMBH, $10^{2}-10^{5}$ $M_{\sun}$ ). The largest published sample of
bona-fide IMBH-powered AGN contains 10 objects confirmed in X-ray. Here we
present a new sample of 15 bona-fide IMBHs, obtained by confirming the
optically selected IMBH candidates by the presence of radiation from the
galactic nucleus in the X-ray range, which increases the number of confirmed
IMBHs at the centers of galaxies by 2.5 times. In the same way, 99 black holes
with masses of $2\cdot 10^{5}-10^{6}$ $M_{\sun}$ were confirmed. The sources
of X-ray data were publicly available catalogs, archives of data, and our own
observations on XMM-Newton, Chandra and Swift. The Eddington coefficients for
30% of the objects from both samples turned out to be close to critical, from
0.5 to 1, which is an unusually high fraction. Also for the first time for
light-weight SMBH the correlations between the luminosity in the [OIII]
emission line or the broad component of the H$\alpha$ line and the luminosity
in the X-ray range were plotted.

###### keywords:

supermassive black holes, active galaxies, x-ray observations

##  1 Introduction

Massive black holes in galaxy centres are believed to co-evolve with the
spheroids of their hosts (Kormendy & Ho 2013), growing via coalescences during
galaxy mergers (Merritt & Milosavljević 2005) and by gas accretion (Volonteri
2012). The discovery of quasars in the early Universe (z $>6.3$, only 750-900
Myr after the Big Bang) hosting super-massive black holes (SMBHs) as heavy as
$10^{10}$$M_{\sun}$ (Mortlock et al. 2011; Wu et al. 2015) cannot be explained
by gas accretion onto stellar mass black hole seeds ($<100$$M_{\sun}$ ) alone.
Alternatively, the rapid inflow and subsequent direct collapse of gas clouds
(Loeb & Rasio 1994; Begelman et al. 2006) can form massive seeds (M
$>10^{5}-10^{6}$$M_{\sun}$ ). The latter scenario solves the SMBH early
formation puzzle butleads to a gap in the present-day black hole mass function
in the IMBH regime ($100<M_{BH}<2\cdot 10^{5}$$M_{\sun}$ ), whereas stellar
mass seeds should leave behind a large number of IMBHs. Therefore, the elusive
IMBH population holds a clue to the understanding of SMBH formation. Now the
largest published sample of bona-fide IMBH-powered AGN contains 10 objects
confirmed in X-ray.

##  2 Sample description

The work is based on a sample of 1928 galaxies with active nuclei with virial
masses of black holes up to $10^{6}$$M_{\sun}$ , 305 of which are IMBHs with
masses less than $2\cdot 10^{5}$ $M_{\sun}$ . The sample was obtained as part
of the work on the article (Chilingarian et al. 2018) by analyzing the spectra
of about a million galaxies from the SDSS DR7 (Abazajian et al. 2009) survey.
A broad component of the $H_{\alpha}$ emission line was identified in the
spectra, and the virial masses of the central black holes were calculated from
it’s parameters. A method for selecting BHs of intermediate masses and light-
weight SMBH was proposed in the work (Chilingarian et al. 2018). Also only
galaxies with well-measured central BH masses and reliable broad component in
lines were selected, which ratio of the BH mass to the mass error is more than
3, important lines do not overlay on the sky lines due to the redshift, and
falls into the AGN or transition region on the BPT (Baldwin et al. 1981)
diagnostic diagram.

##  3 Data analysis methods

###  3.1 X-ray data collection, calibration and analysis

To confirm the AGN nature of an object it is necessary to detect the X-ray
emission, produced by the accretion disk of the black hole. X-ray data were
obtained from the Chandra Source Catalog 2, 4XMM-DR10, Second Swift XRT Point
Source Catalog, Second ROSAT all-sky survey, from the x-ray data archives as
Chandra Data Archive and XMM-Newton science archive and our own observations
at XMM-Newton, Chandra and Swift. For calibration of raw observational data
and extraction of spectra was used special software: Science Analysis System
(SAS) for XMM-Newton data, Chandra Interactive Analysis of Observations for
Chandra and XSelect for Swift.

The next important step was extracting the spectrum of the object and fitting
its parameters. Using the parameters and form of the spectrum we can once
again confirm its AGN nature. Spectral data was fitted with powerlaw and
photoelectric absorption models with XSPEC and Sherpa packages in Python 3. It
turned out that most of our objects have a photon index over 2, that is
typical for AGN.

Figure 1: Examples of spectra. The model is shown with a red line

###  3.2 X-ray luminosity of the stellar population

It is worth remembering that the object’s X-ray luminosity should be higher
than expected from the stellar population of the host galaxy. X-ray stellar
population mainly consists of low-mass and high-mass X-ray binaries. Their
contribution to luminosity can be calculated using the following relations

$L_{HX}^{gal}=\alpha M_{*}+\beta SFR$ $\alpha=(9.05\pm 0.37)\times
10^{28}\,ergs\,s^{-1}\,M_{\odot}^{-1}$ $\beta=(1.62\pm 0.22)\times
10^{39}\,ergs\,s^{-1}\,(M_{\odot}\,yr^{-1})^{-1}$

###  3.3 Eddington luminosities

The luminosity of the accretion disk directly depends on this accretion rate.
Therefore, in order to understand how high the rate of accretion onto a black
hole is, it is necessary to calculate the Eddington ratio. The masses are
known from optical observations and the Eddington luminosity is easily
calculated

$L_{edd}=\frac{4\pi GMm_{p}c}{\sigma_{T}}\approx
10^{38}\frac{M}{M_{\odot}}\>ergs\cdot s^{-1}$

The bolometric luminosity values are determined from the luminosity in the
X-ray range through the bolometric correction coefficient defined as

$K(L_{X})=\frac{L_{bol}}{L_{X}}=15.33\left[1+\left(\frac{\log(L_{X}/L_{\odot}}{11.48}\right)^{16.20}\,\right]$

##  4 Results and conclusions

###  4.1 New IMBHs with X-ray confirmation

From the sample of 1928 candidates for low-mass black holes, 124 black holes
with masses less than $10^{6}\,$$M_{\sun}$ were confirmed in the X-ray range,
of which 24 are intermediate masses black holes with masses $M_{BH}<2\cdot
10^{5}\,$$M_{\sun}$ , and their x-ray luminosities were calculated. For 55
objects from this sample, in addition to the luminosity, the parameters of
their X-ray spectra were determined. Thus, confirmed intermediate-mass black
holes sample was increased by a factor of 2.5, and this number of IMBHs is a
strong argument for the formation of supermassive black holes from low-mass
progenitors. Most of the sources have ”soft” X-ray spectra – for 39 of 55
objects the photon index is $\Gamma>2$. For a large number of objects the
Eddington ratio turned out to be very high. Thus, for 39 of 123 BHs it is more
than 30%. This points on very high accretion rates on many IMBH.

###  4.2 Correlation between X-ray luminosity and luminosity in emission lines

For supermassive black holes have found a correlation between the luminosity
in the X-ray range and the optical emission lines [OIII] and $H_{\alpha}$
(Heckman et al. 2005; Ueda et al. 2015; Jin et al. 2012). We verified this
correlation for light-weight SMBH and IMBH for the first time. The fluxes in
lines were obtained together with the black holes masses by processing optical
spectra. Correlations were approximated by linear regression

$\log(L_{emis.\>Line})=a\cdot\log(L_{X})+b$

The results of this approximation are shown in the table 1.

Table 1: Aproximation results Line | a | b | error
---|---|---|---
[OIII] | $0.339\pm 0.045$ | $25.950\pm 1.898$ | 0.475
$H_{\alpha}$ | $0.344\pm 0.039$ | $25.680\pm 1.642$ | 0.411

It turned out that for low-mass black holes the slope of the regression is
much smaller than for the sample of supermassive black holes. It could be
caused by a change in the geometry of the region of X-ray radiation - this and
other possible hypotheses requires further research. However, it is already
clear that these correlations are extremely useful - it can be used to
estimate the X-ray luminosity for calculating the exposures necessary to
accumulate a sufficient number of photons for spectral analysis.

Figure 2: Correlations between the luminosity in the [OIII] emission line or
the broad component of the H$\alpha$ line and the luminosity in the X-ray
range

###  4.3 Conclusions about the ways of black holes formation

The key result of this study is a significant increase in sample of IMBH, and
this confirms that SMBH grew from their lower mass progenitors. Also we show
that at low mass regime the number of objects emitting close to the Eddington
limit is unexpectedly large. This indicates that BH can grow rapidly,
including due to accretion of matter, not only due to merging and host
galaxies bulges growth.

Acknowledgements. This project is supported by the RScF Grant 17-72-20119.

## References

  * Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Begelman et al. (2006) Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, MNRAS, 370, 289
  * Chilingarian et al. (2018) Chilingarian, I. V., Katkov, I. Y., Zolotukhin, I. Y., et al. 2018, ApJ, 863, 1
  * Heckman et al. (2005) Heckman, T. M., Ptak, A., Hornschemeier, A., & Kauffmann, G. 2005, ApJ, 634, 161
  * Jin et al. (2012) Jin, C., Ward, M., & Done, C. 2012, MNRAS, 422, 3268
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511
  * Loeb & Rasio (1994) Loeb, A. & Rasio, F. A. 1994, ApJ, 432, 52
  * Merritt & Milosavljević (2005) Merritt, D. & Milosavljević, M. 2005, Living Reviews in Relativity, 8
  * Mortlock et al. (2011) Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616
  * Ueda et al. (2015) Ueda, Y., Hashimoto, Y., Ichikawa, K., et al. 2015, ApJ, 815, 1
  * Volonteri (2012) Volonteri, M. 2012, Science, 337, 544
  * Wu et al. (2015) Wu, X.-B., Wang, F., Fan, X., et al. 2015, Nature, 518, 512