# Stochastic representation of processes with resetting

Marcin Magdziarz<EMAIL_ADDRESS>Faculty of Pure and Applied
Mathematics, Hugo Steinhaus Center, Wroclaw University of Science and
Technology, Wyspianskiego 27, 50-370 Wroclaw, Poland  Kacper Taźbierski
Faculty of Pure and Applied Mathematics, Hugo Steinhaus Center, Wroclaw
University of Science and Technology, Wyspianskiego 27, 50-370 Wroclaw, Poland

(September 7, 2021)

###### Abstract

In this paper we introduce a general stochastic representation for an
important class of processes with resetting. It allows to describe any
stochastic process intermittently terminated and restarted from a predefined
random or non-random point. Our approach is based on stochastic differential
equations called jump-diffusion models. It allows to analyze processes with
resetting both, analytically and using Monte Carlo simulation methods. To
depict the strength of our approach, we derive a number of fundamental
properties of Brownian motion with Poissonian resetting, such as: the Itô
lemma, the moment-generating function, the characteristic function, the
explicit form of the probability density function, moments of all orders,
various forms of the Fokker-Planck equation, infinitesimal generator of the
process and its adjoint operator. Additionally, we extend the above results to
the case of time-nonhomogeneous Poissonian resetting. This way we build a
general framework for the analysis of any stochastic process with intermittent
random resetting.

###### pacs:

05.40.Fb, 02.50-r

##  I Introduction

In recent years we have been observing a great interest in the field of
intermittent stochastic processes, in which the motion of a particle is
interrupted by random resetting to an initial state. Random motion with
resetting is typically observed in various searching strategies and foraging
patterns, where the living organism after an unsuccessful excursion returns to
the initial position and start the search again Campos  _et al._ (2014). The
justification of such strategy is that if one does not succeed in finding the
target within certain not very long time period, sometimes it is more
beneficial and safe to move back to the origin and start another excursion,
see Manrubia and Zanette (1999); Gelenbe (2010). On the other hand, in the
field of population dynamics resetting could be interpreted as an event
reducing population to its natural size according to the environmental
capacity Kyriakidis (1994). What makes the idea also interesting is the way
that random resetting completely changes the properties of the diffusion
process. In particular, mean first passage time of diffusion with resetting is
finite, which is in sharp contrast to the standard diffusion case (Evans and
Majumdar, 2011; Evans  _et al._ , 2020). This fact depicts the efficiency of
searching strategies with resetting mechanism Tong  _et al._ (2008). However,
if the resetting timer has an infinite mean then mean first passage time may
be infinite. Intermittent processes have found important application in other
fields, such as: phenotypic diversity, population growth and information in
fluctuating environments Kussell and Leibler (2005), ecology Viswanathan  _et
al._ (2011), enzymatic reactions Reuveni  _et al._ (2014) or computer science
and optimization problems Montanari and Zecchina (2002). The historical
background of processes with resetting can be found in Montero  _et al._
(2017). The interested reader is referred to a recent comprehensive review by
Evans et al. Evans  _et al._ (2020). Some latest results related to non-
Markovian resetting models can be found in Bodrova  _et al._ (2019); Zhou  _et
al._ (2020, 2021). We also refer the reader to the following relevant works
related to the restarted processes, published well before the terms ”restart”
and ”resetting” became popular Eliazar and Klafter (2004, 2005, 2006).

To characterize a process with resetting, one needs to deal with two sources
of randomness. The first one accounts for the particle motion between
resetting times, the second one is the point process determining the moments
of resetting. The usual description of motion with resetting is via CTRW with
additional resetting mechanism or via related deterministic partial
differential equations of the Fokker-Planck type (Evans  _et al._ , 2020).

Here we use another approach, which is founded on the theory of Lévy processes
and the corresponding stochastic differential equations (SDEs) – the so-called
jump-diffusion models. There is a rich history of research of Lévy processes.
The most well-known is obviously Brownian motion (or Wiener process). Other
important examples include $\alpha$-stable, Linnik, Mittag-Leffler, Gamma, or
Laplace processes Ken-Iti (1999). Their first applications arose quickly in
the physical Einstein (1905); Smoluchowski (1906) and financial Bachelier
(1900); Merton (1976); Black and Scholes (1973) sciences. They have naturally
spread to many other fields of science, e.g. as models in biology, chemistry,
data mining, statistics, etc. Barndorff-Nielsen  _et al._ (2012). In 1940’s,
thanks to Kiyoshi Itô, the theory of stochastic differential equations was
developed Itô (1944), Itô (1946). The stochastic Itô integral and the
celebrated Itô lemma allowed to integrate functions or processes with respect
to the Wiener process. Also, other types of stochastic integrals with
important role in physics, such as Stratanovich integral, were developed.
Jump-diffusion models are defined as SDEs, in which the Brownian noise term is
replaced or complemented by arbitrary Lévy noise Øksendal and Sulem (2019).
Since general Lévy processes have discontinuous trajectories Ken-Iti (1999),
therefore solutions of such Lévy driven SDEs display jumps. Jump-diffusion
models have found applications in various fields, such as: finance,
statistical physics, pattern theory and computational vision and many more,
see Øksendal and Sulem (2019); Eliazar and Klafter (2003) and references
therein. As we will show in the following, these processes are also tailor-
made to describe and analyze processes with resetting.

In this paper we introduce a general stochastic framework for processes with
resetting in the language of stochastic differential equations (SDEs) with
jumps (jump-diffusion models). Using the stochastic approach a version of the
Itô lemma for processes with resetting will be derived. Next, applying the
latter result we will derive the moment generating function (MGF) and the
Fourier transform for diffusion with resetting. Using the inverse Fourier
transform we will then determine the probability density function (PDF) of the
process in explicit form. The corresponding Fokker-Planck equation will be
also derived in various forms. We will discuss the equivalence of our Fokker-
Planck equation with the one given in Evans  _et al._ (2020). Moreover we will
find the infinitesimal generator and the adjoint operator for the Markovian
resetting process. In addition, we will extend the above results to the case
of time-nonhomogeneous Poissonian resetting. We will show that, depending on
the intensity function of the resetting events, we can observe: convergence to
stationary distribution, convergence to a point, subdiffusive or diffusive
behaviour, for long times.

It should be underlined that the presented here approach can be easily
generalized to the case of non-Markovian resetting times as well as other, not
necessarily Brownian, driving processes (Lévy processes, fractional processes,
etc.).

##  II Stochastic representation of processes with resetting

Let us consider a diffusing particle with initial position $x_{0}$ and
Poissonian resetting with rate $r$ to position $x_{R}$. In general the initial
position $x_{0}$ and the resetting position $x_{R}$ can be distinct, they can
be even random (in this paper we assume that $x_{0}$ and $x_{R}$ are constant,
however the methodology we use, can also be applied to the random case, which
will lead to different results and distributions). Position of such particle
at time $t$ is usually defined in the following way Evans  _et al._ (2020):

$\begin{gathered}x(t+\operatorname{\operatorname{d}\\!}t)=x_{R}\quad\text{with
probability }r\operatorname{\operatorname{d}\\!}t\\\
=x(t)+\xi(t)(\operatorname{\operatorname{d}\\!}t)^{\frac{1}{2}}\quad\text{with
probability }(1-r\operatorname{\operatorname{d}\\!}t),\end{gathered}$ (1)

where $\xi(t)$ is the standard Gaussian white noise with mean 0 and variance
$2D$. The typical trajectory of a diffusion process with resetting is depicted
in Fig. 1. The Fokker-Planck equation corresponding to (1) has the form Evans
_et al._ (2020)

$\partial_{t}p(x,t)=D\partial_{xx}p(x,t)+r\delta(x-x_{R})-rp(x,t),$ (2)

with initial condition $p(x,0)=\delta(x-x_{0})$. Here $p(x,t)$ is the PDF of
the diffusing particle with resetting and $\delta(\cdot)$ is the Dirac delta.
In what follows we will put for simplicity $D=1/2$. Using the above formalism
one can examine the interesting phenomenon of non-equilibrium steady states
Evans  _et al._ (2020). Even though the stationary state is independent of
time, there is a driving force in form of the resetting that creates the
probability flow. In a purely mathematical setting it means that the gradient
of stationary state is non-zero. This corresponds to the physical concept of
non-equilibrium steady state. A nice review on this topic can be found in Ge
_et al._ (2012).

Also in Evans  _et al._ (2020) an absorption of the process by a trap was
studied. For the standard diffusion, time to absorption follows Lévy
distribution, which has an infinite mean (being an $\alpha$-stable
distribution with $\alpha=1/2$). When the resetting is introduced, the mean
time to absorption becomes finite and we are able to find an optimal resetting
rate Evans  _et al._ (2020). This fact clearly shows that resetting can be
beneficial in the searching strategies.

In this paper we introduce a different approach to define and analyze
processes with resetting. Namely, we will use jump-diffusion processes to
define the stochastic dynamics with resetting. A jump-diffusion process is
defined as the solution of the following SDE Øksendal and Sulem (2019):

$\displaystyle
dX_{t}=\mu_{t}dt+\sigma_{t}dW_{t}+\nu_{t}dL_{t}\;,\;\;\;X_{0}=x_{0}.$ (3)

Here $W_{t}$ is the standard Brownian motion, $L_{t}$ the Lévy process, it is
usually the Poisson process and introduces jumps to the observed dynamics,
$\mu_{t}$, $\sigma_{t}$ and $\nu_{t}$ are the appropriate parameters of the
model, which can be in general space and time dependent. It is also assumed
that $W_{t}$ and $L_{t}$ are independent. Recall that Lévy process $L_{t}$ is
a stochastic process satisfying Øksendal and Sulem (2019):

  1. 1.

$L_{0}=0$,

  2. 2.

$L_{t}$ has independent increments,

  3. 3.

$L_{t}$ has stationary increments,

  4. 4.

$\displaystyle\mathop{\forall}_{\varepsilon,t>0}\lim_{h\to
0}P\left(|L_{t+h}-L_{t}|>\varepsilon\right)=0$, i.e. $L_{t}$ is continuous in
probability.

Since every Lévy process is a Markov process, solutions of SDE (3) are
Markovian. Typical examples of Lévy processes are: Brownian motion, Poisson
process, and $\alpha$-stable Lévy process. Brownian motion is the only Lévy
process (up to a constant drift) with continuous trajectories. All the other
Lévy processes have jumps. Poisson process is a non-decreasing Lévy process
with jumps of size $1$ and flat periods between jumps. Times between
consecutive jumps of the Poisson process are independent and drawn from the
exponential distribution with mean $r>0$. The constant $r$ is called intensity
of the Poisson process

Now let us consider the following particular case of jump-diffusion process:

$\operatorname{\operatorname{d}\\!}X_{t}=\operatorname{\operatorname{d}\\!}W_{t}+\left(x_{R}-X_{t}\right)\operatorname{\operatorname{d}\\!}N_{t},\quad
X_{0}=x_{0},$ (4)

where $N_{t}$ is the Poisson Process with intensity $r$, $x_{0}$ and $x_{R}$
are constants. We argue that (4) is the continuous-time analogue of diffusion
with Poissonian resetting defined in (1). Indeed, whenever $N_{t}$ has jump,
we get that $dN_{t}=1$ and therefore the term
$\left(x_{R}-X_{t}\right)\operatorname{\operatorname{d}\\!}N_{t}$ sends back
the process to the resetting position $x_{R}$. Between the resetting times,
i.e. when $N_{t}$ is constant, we get that $dN_{t}$=0. So in this case
$X_{t}=W_{t}$. Thus between consecutive resetting times the particle performs
Brownian motion. A typical trajectory obtained using (4) can be seen in Fig,
1. Clearly, we observe a diffusing particle with resetting to $x_{R}=2$.

Figure 1: An exemplary plot of a resetting process. The process starts at
$x_{0}=0$ and then diffuses normally until it moves instantaneously to
$x_{R}=2$ in every consecutive resetting event.

We can also write down the integrated form of (4):

$X_{t}=x_{0}+W_{t}+\int_{0}^{t}(x_{R}-X_{s})dN_{s}=x_{0}+W_{t}+\sum\displaylimits_{n=1}^{N_{t}}\left(x_{R}-X_{\tau_{n}}\right).$
(5)

Here $\tau_{n}=\max\\{t>0:N_{t}\leq n\\}$ are the consecutive resetting times.

In what follows we will show that the derived here stochastic representation
(4) constitutes a general framework for the analysis of processes with
resetting. It allows to study many key properties of intermittent processes
both analytically and using Monte Carlo simulation methods. Moreover, it
should be underlined that the introduced approach can be easily generalized to
the case of non-Poissonian resetting times as well as other arbitrary driving
processes. More precisely, suppose that we want to define a stochastic process
$\tilde{X}_{t}$ describing an arbitrary dynamics with resetting. Let
$T_{1},T_{2},...$ be an arbitrary sequence of positive random variables
describing times between consecutive resetting events. Moreover, let us assume
that $\tilde{W}_{t}$ is an arbitrary stochastic process describing the
particle motion between resetting events. Then $\tilde{X}_{t}$ can be defined
as the solution of the following SDE:

$\operatorname{\operatorname{d}\\!}\tilde{X}_{t}=\operatorname{\operatorname{d}\\!}\tilde{W}_{t}+\left(x_{R}-\tilde{X}_{t}\right)\operatorname{\operatorname{d}\\!}\tilde{N}_{t},\quad\tilde{X}_{0}={x}_{0},$
(6)

which is a straightforward generalization of (4). Here
$\tilde{N}_{t}=\max\\{n\in\mathbb{N}:\sum_{i=1}^{n}T_{i}\leq t\\}$. It counts
the number of resetting events up to time $t$. In particular, $\tilde{N}_{t}$
can be a renewal processes. Consequently, any process with resetting can be
written in the form of (6).

###  II.1 Itô lemma

Itô lemma is the main tool in the analysis of SDEs. It allows to find their
solutions as well as many of their key properties. Let us derive a version of
the Itô lemma corresponding to the resetting process (4). For completeness we
will derive it in a more general setting. Let $Z_{t}=f(X_{t})$, where $f$ is
appropriately smooth and $X_{t}$ is governed by the following SDE

$\displaystyle
dX_{t}=\mu_{t}dt+\sigma_{t}dW_{t}+(x_{R}-X_{t})dN_{t}\;,\;\;\;X_{0}=x_{0}.$
(7)

Note that for $\mu_{t}=0$ and $\sigma_{t}=1$ we recover the resetting process
(4). Using the Taylor expansion we get

$\displaystyle\operatorname{\operatorname{d}\\!}f(X_{t})=\sum\displaylimits_{i=1}^{\infty}\frac{f^{(i)}(X_{t})}{i!}\left(\operatorname{\operatorname{d}\\!}X_{t}\right)^{i}.$

According to classical Itô calculus we have that
$(\operatorname{\operatorname{d}\\!}W_{t})^{2}=\operatorname{\operatorname{d}\\!}t$
and $(\operatorname{\operatorname{d}\\!}W_{t})^{n}=0$ for $n>2$. Moreover
$(\operatorname{\operatorname{d}\\!}N_{t})^{n}=\operatorname{\operatorname{d}\\!}N_{t}$
for any $n\in\mathbb{N}$ since $N_{t}$ is the Poisson process. Therefore

$(\operatorname{\operatorname{d}\\!}X_{t})^{2}=\left(\sigma^{2}_{t}\operatorname{\operatorname{d}\\!}t+\left(x_{R}-X_{t}\right)^{2}\operatorname{\operatorname{d}\\!}N_{t}\right)$

and

$(\operatorname{\operatorname{d}\\!}X_{t})^{i}=\left(x_{R}-X_{t}\right)^{i}\operatorname{\operatorname{d}\\!}N_{t}$

for $i>2$. Consequently

$\displaystyle\operatorname{\operatorname{d}\\!}f(X_{t})$
$\displaystyle=f^{{}^{\prime}}(X_{t})\left(\mu_{t}\operatorname{\operatorname{d}\\!}t+\sigma_{t}\operatorname{\operatorname{d}\\!}W_{t}+\left(x_{R}-X_{t}\right)\operatorname{\operatorname{d}\\!}N_{t}\right)+$
$\displaystyle+\frac{1}{2}f^{{}^{\prime\prime}}(X_{t})\left(\sigma^{2}_{t}\operatorname{\operatorname{d}\\!}t+\left(x_{R}-X_{t}\right)^{2}\operatorname{\operatorname{d}\\!}N_{t})\right)+\sum\displaylimits_{i=3}^{\infty}\frac{f^{(i)}(X_{t})}{i!}\left(x_{R}-X_{t}\right)^{i}\operatorname{\operatorname{d}\\!}N_{t}$
$\displaystyle=f^{{}^{\prime}}(X_{t})\sigma_{t}\operatorname{\operatorname{d}\\!}W_{t}+\left(f^{{}^{\prime}}(X_{t})\mu_{t}+f^{{}^{\prime\prime}}(X_{t})\frac{\sigma_{t}^{2}}{2}\right)\operatorname{\operatorname{d}\\!}t+$
$\displaystyle+\left(f^{{}^{\prime}}(X_{t})\left(x_{R}-X_{t}\right)+\frac{1}{2}f^{{}^{\prime\prime}}(X_{t})\left(x_{R}-X_{t}\right)^{2}+\sum\displaylimits_{i=3}^{\infty}\frac{f^{(i)}(X_{t})}{i!}\left(x_{R}-X_{t}\right)^{i}\right)\operatorname{\operatorname{d}\\!}N_{t}$
$\displaystyle=f^{{}^{\prime}}(X_{t})\sigma_{t}\operatorname{\operatorname{d}\\!}W_{t}+\left(f^{{}^{\prime}}(X_{t})\mu_{t}+f^{{}^{\prime\prime}}(X_{t})\frac{\sigma_{t}^{2}}{2}\right)\operatorname{\operatorname{d}\\!}t+\left(\sum\displaylimits_{i=1}^{\infty}\frac{f^{(i)}(X_{t})}{i!}\left(x_{R}-X_{t}\right)^{i}+f(X_{t})-f(X_{t})\right)\operatorname{\operatorname{d}\\!}N_{t}$
$\displaystyle=f^{{}^{\prime}}(X_{t})\sigma_{t}\operatorname{\operatorname{d}\\!}W_{t}+\left(f^{{}^{\prime}}(X_{t})\mu_{t}+f^{{}^{\prime\prime}}(X_{t})\frac{\sigma_{t}^{2}}{2}\right)\operatorname{\operatorname{d}\\!}t+\left(\sum\displaylimits_{i=0}^{\infty}\frac{f^{(i)}(X_{t})}{i!}\left(x_{R}-X_{t}\right)^{i}-f(X_{t})\right)\operatorname{\operatorname{d}\\!}N_{t}.$

Finally, using the Taylor formula we get the generalised Itô formula

$\operatorname{\operatorname{d}\\!}Z_{t}=f^{{}^{\prime}}(X_{t})\sigma_{t}\operatorname{\operatorname{d}\\!}W_{t}+\left(f^{{}^{\prime}}(X_{t})\mu_{t}+f^{{}^{\prime\prime}}(X_{t})\frac{\sigma_{t}^{2}}{2}\right)\operatorname{\operatorname{d}\\!}t+\left(f(x_{R})-f(X_{t})\right)\operatorname{\operatorname{d}\\!}N_{t}.$

For $\mu_{t}=0$ and $\sigma_{t}=1$ we recover the Itô lemma for the resetting
process (4)

$\displaystyle\operatorname{\operatorname{d}\\!}f(X_{t})=\frac{1}{2}f^{{}^{\prime\prime}}(X_{t})\operatorname{\operatorname{d}\\!}t+f^{{}^{\prime}}(X_{t})\operatorname{\operatorname{d}\\!}W_{t}+\left(f(x_{R})-f(X_{t})\right)\operatorname{\operatorname{d}\\!}N_{t}.$
(8)

###  II.2 Moment generating function and Fourier transform

Now, using the above Itô lemma, let us calculate the MGF
$M_{t}(s)=\mathbb{E}\left(e^{sX_{t}}\right)$ of the resetting process $X_{t}$
given by (4). Here $\operatorname{\mathbb{E}}$ is the expected value. Applying
(8) for the function $f_{s}(x)=e^{sx}$ and putting $Z_{t}=f_{s}(X_{t})$ we
arrive at the following SDE:

$\operatorname{\operatorname{d}\\!}Z_{t}=\operatorname{\operatorname{d}\\!}f_{s}(X_{t})=sZ_{t}\operatorname{\operatorname{d}\\!}W_{t}+\frac{1}{2}s^{2}Z_{t}\operatorname{\operatorname{d}\\!}t+\left(f_{s}(x_{R})-Z_{t}\right)\operatorname{\operatorname{d}\\!}N_{t}.$
(9)

Note that $\mathbb{E}\left(Z_{t}\right)$ is equal to the MGF of $X_{t}$.
Taking the expected value of both sides of (9) we get

$\displaystyle\operatorname{\mathbb{E}}\left(\operatorname{\operatorname{d}\\!}Z_{t}\right)$
$\displaystyle=s{\operatorname{\mathbb{E}}\left(Z_{t}\operatorname{\operatorname{d}\\!}W_{t}\right)}+\frac{1}{2}s^{2}\operatorname{\operatorname{d}\\!}t\operatorname{\mathbb{E}}(Z_{t})+\operatorname{\mathbb{E}}\left(e^{sx_{R}}-f(X_{t})\right)\operatorname{\operatorname{d}\\!}N_{t}$
$\displaystyle=\frac{1}{2}s^{2}\operatorname{\operatorname{d}\\!}t\operatorname{\mathbb{E}}(Z_{t})+\sum\displaylimits_{i=0}^{1}\operatorname{\mathbb{E}}\left(\left(e^{sx_{R}}-f(X_{t})\right)\operatorname{\operatorname{d}\\!}N_{t}\Big{|}\operatorname{\operatorname{d}\\!}N_{t}=i\right)\operatorname{\mathbb{P}}\left(\operatorname{\operatorname{d}\\!}N_{t}=i\right)$
$\displaystyle=\frac{1}{2}s^{2}\operatorname{\operatorname{d}\\!}t\operatorname{\mathbb{E}}(Z_{t})+\operatorname{\mathbb{E}}\left(e^{sx_{R}}-f(X_{t})\Big{|}\operatorname{\operatorname{d}\\!}N_{t}=1\right)\operatorname{\mathbb{P}}\left(\operatorname{\operatorname{d}\\!}N_{t}=1\right)$
$\displaystyle=\frac{1}{2}s^{2}\operatorname{\operatorname{d}\\!}t\operatorname{\mathbb{E}}(Z_{t})+r\operatorname{\operatorname{d}\\!}te^{sx_{R}}-r\operatorname{\operatorname{d}\\!}t\operatorname{\mathbb{E}}\left(f(X_{t})\right).$

The sum in the second line of the above equalities is the consequence of the
total probability formula. Dividing both sides by
$\operatorname{\operatorname{d}\\!}t$, interchanging the order of integration
and differentiation and using the dominated convergence theorem (the process
has clearly lighter tails than the ordinary diffusion, and as such the theorem
can be applied), substituting $\operatorname{\mathbb{E}}(Z_{t})=M_{t}(s)$ and
rearranging the equation, we get a simple ordinary differential equation

$\displaystyle\partial_{t}M_{t}(s)+\left(r-\frac{1}{2}s^{2}\right)M_{t}(s)=re^{sx_{R}}.$

This yields the solution

$\displaystyle
M_{t}(s)=\frac{re^{sx_{R}}}{r-\frac{1}{2}s^{2}}+c(s)e^{-\left(r-\frac{1}{2}s^{2}\right)t}.$

Using conditions $M_{0}(s)=e^{sX_{0}}$ and $M_{t}(0)=1$ for all $t$ we get
that the unique constant $c(s)$ equals

$c(s)=\left(e^{sX_{0}}-\frac{r}{r-\frac{1}{2}s^{2}}e^{sx_{R}}\right).$

Figure 2: A plot of the derived mean of the resetting process (top left), MGF
(top right), the real part of the Fourier transform (bottom left) and the
imaginary part of the Fourier transform (bottom right) at a different time
points (0.1, 0.5, 1.5), as a blue, red and yellow lines respectively. The
process starts at 0 and resets with rate 1 to point 5. We can also see the
unusual behavior of the MGF at the point $\sqrt{2}$, due to the divergence of
the function.

Finally our MGF has the form

$\displaystyle
M_{t}(s)=\frac{re^{sx_{R}}}{r-\frac{1}{2}s^{2}}+\left(e^{sX_{0}}-\frac{re^{sx_{R}}}{r-\frac{1}{2}s^{2}}\right)e^{-\left(r-\frac{1}{2}s^{2}\right)t}.$
(10)

We can also notice that in case of the lack of resetting ($r=0$) the MGF
simplifies to the MGF of Brownian motion starting at $X_{0}$

$\frac{0e^{sx_{R}}}{0-\frac{1}{2}s^{2}}+\left(e^{sX_{0}}-\frac{0}{0-\frac{1}{2}s^{2}}e^{sx_{R}}\right)e^{-\left(0-\frac{1}{2}s^{2}\right)t}=e^{sX_{0}+t\frac{1}{2}s^{2}}.$

Additionally for large time $t\rightarrow\infty$ MGF of the process with
resetting converges to

$M(s)=\frac{e^{sx_{R}}}{1-\frac{1}{2r}s^{2}},$

which is a MGF of the Laplace distribution with mean $x_{R}$ and scale
parameter $(2r)^{-\frac{1}{2}}$.

The Fourier transform of the process with resetting
$\varphi_{t}(s)=\mathbb{E}\left(e^{isX_{t}}\right)$ can be derived analogously
or using the fact that $\varphi_{t}(s)=M_{t}(is)$, and equals

$\displaystyle\varphi_{t}(s)=\frac{re^{isx_{R}}}{r+\frac{1}{2}s^{2}}+\left(e^{isX_{0}}-\frac{re^{isx_{R}}}{r+\frac{1}{2}s^{2}}\right)e^{-\left(r+\frac{1}{2}s^{2}\right)t}.$
(11)

Plots of $M_{t}(s)$ and $\varphi_{t}(s)$ are presented in Fig. 2. We can
immediately notice that the MGF is well defined only for $|s|<\sqrt{2r}$. The
moment generating function and Fourier transform are closely related via
formula $\varphi_{t}(s)=M_{t}(is)$. However, the advantage of $\varphi_{t}(s)$
is that it is well defined for any $s\in\mathbb{R}$. The moment generating
function does not have this property. On the other hand, the advantage of
$M_{t}(s)$ is that it is real not complex. Both functions $M_{t}(s)$ (if well
defined) and $\varphi_{t}(s)$ determine uniquely the underlying distribution.

###  II.3 Explicit PDF

Now, using the inverse Fourier transform we are able to invert
$\varphi_{t}(s)$. This way we will arrive at the explicit formula for the PDF
$p(x,t)$ of the resetting process $X_{t}$. An exemplary plot of $p(x,t)$ is
presented in Fig. 3.

Figure 3: A histogram of $X_{t}$ obtained using Monte Carlo simulations
compared with the derived analytical PDF of the process at $t=0.1$. We observe
perfect agreement between both results. The process starts at 0 and resets to
point 3 with rate 1.

Let us split the Fourier transform (11) into three parts in the following way

$\varphi_{t}(s)=\varphi_{1}+\varphi_{2}+\varphi_{3},$

where $\varphi_{1}=\frac{r}{r+\frac{1}{2}s^{2}}e^{isx_{R}}$,
$\varphi_{2}=e^{isX_{0}}e^{-\left(r+\frac{1}{2}s^{2}\right)t}$,
$\varphi_{3}=-\frac{r}{r+\frac{1}{2}s^{2}}e^{isx_{R}}e^{-\left(r+\frac{1}{2}s^{2}\right)t}$.
Using the well known properties of the Fourier transform, we can already see
that the first term is just the Fourier transform of the Laplace distribution
with the corresponding PDF Billingsley (1986)

$\displaystyle\operatorname{\mathcal{F}}^{-1}\left\\{\varphi_{1}\right\\}=f_{\mathcal{L}}(x)=\sqrt{\frac{r}{2}}\exp\left(-\frac{|x-x_{R}|}{(2r)^{-1/2}}\right),$
(12)

see also Reed and Hughes (2002); Reed (2003) for the relationship between
Laplace distribution and stopping times of random motions.

We will calculate the second term directly from the definition of the inverse
transform

$\displaystyle\operatorname{\mathcal{F}}^{-1}\left\\{\varphi_{2}\right\\}$
$\displaystyle=\frac{1}{2\pi}\int_{\mathbb{R}}\displaylimits
e^{-isx}e^{isX_{0}-rt-\frac{t}{2}s^{2}}\operatorname{\operatorname{d}\\!}s=\frac{e^{-rt}}{2\pi}\int_{\mathbb{R}}\displaylimits
e^{i(X_{0}-x)s-\frac{t}{2}s^{2}}\operatorname{\operatorname{d}\\!}s$
$\displaystyle=\frac{e^{-rt-\frac{(x-X_{0})^{2}}{2t}}}{\sqrt{2\pi
t}}=e^{-rt}f_{W_{X_{0}}}(x),$

where $f_{W_{X_{0}}}(x)=\frac{e^{-\frac{(x-X_{0})^{2}}{2t}}}{\sqrt{2\pi t}}$
is the PDF of standard Brownian motion starting at $X_{0}$.

Let us now write $\varphi_{3}$ in the following way

$\varphi_{3}=-e^{-rt}\cdot\frac{r}{r+\frac{1}{2}s^{2}}e^{isx_{R}}\cdot
e^{-\frac{t}{2}s^{2}}.$

The first term of the product is independent of $s$, so it’s treated as a
scalar in the inverse Fourier transform. The second one is the Fourier
transform of the Laplace distribution and the third is simply the Fourier
transform of the normal distribution $\mathcal{N}(0,t)$. Having all this and
using the well-known fact that
$\operatorname{\mathcal{F}}\\{f\\}\operatorname{\mathcal{F}}\\{g\\}=\operatorname{\mathcal{F}}\\{f*g\\}$,
where $f*g$ is the convolution of functions, we finally get that the explicit
PDF $p(x,t)$ of the resetting process (4) equals

$p(x,t)=f_{\mathcal{L}}(x)+e^{-rt}\left(f_{W_{X_{0}}}(x)-f_{\mathcal{N}(0,t)}*f_{\mathcal{L}}(x)\right).$
(13)

Here $f_{\mathcal{N}(0,t)}(x)$ is the PDF of $\mathcal{N}(0,t)$. Note that if
$f$ and $g$ are two PDFs of some independent random variables, say V and W,
then $f*g$ is the PDF of $V+W$. We underline that the convolution of Laplace
and normal distributions can be calculated explicitly, see Geraci (2017).

###  II.4 Moments

Applying the above formula for PDF we are also able to calculate the mean of
the resetting process

$\operatorname{\mathbb{E}}(X_{t})=\int_{\mathbb{R}}xp(x,t)\operatorname{\operatorname{d}\\!}x=\int_{\mathbb{R}}xf_{\mathcal{L}}(x)\operatorname{\operatorname{d}\\!}x+e^{-rt}\left(\int_{\mathbb{R}}x\left(f_{W_{X_{0}}}(x)-f_{\mathcal{N}(0,t)}*f_{\mathcal{L}}(x)\right)\operatorname{\operatorname{d}\\!}x\right).$

Using the fact that the convolution of PDFs is the PDF of the sum of the
corresponding independent random variables, we get that

$\operatorname{\mathbb{E}}(X_{t})=x_{R}+e^{-rt}\left(X_{0}-x_{R}\right).$

Let us now derive the formula for the $n$-th moment of the resetting process.
For simplicity of the notation we will focus on the case $x_{R}=0$. We have

$\operatorname{\mathbb{E}}(X_{t}^{n})=\int_{\mathbb{R}}\displaylimits
x^{n}p(x,t)\operatorname{\operatorname{d}\\!}x=\int_{\mathbb{R}}\displaylimits
x^{n}f_{\mathcal{L}}(x)\operatorname{\operatorname{d}\\!}x+e^{-rt}\left(\int_{\mathbb{R}}\displaylimits
x^{n}f_{W_{X_{0}}}(x)\operatorname{\operatorname{d}\\!}x-\int_{\mathbb{R}}\displaylimits
x^{n}f_{\mathcal{N}(0,t)}*f_{\mathcal{L}}(x)\operatorname{\operatorname{d}\\!}x\right).$

The first integral is the $n$-th moment of the Laplace random variable $L$
with mean $0$ and scale parameter $(2r)^{-\nicefrac{{1}}{{2}}}$. Using the
scaling property it can be rewritten as $L=(2r)^{-\nicefrac{{1}}{{2}}}Y$,
where $Y\sim\mathcal{L}(0,1)$. It is well known, that the $n$-th moment of $Y$
equals $n!$ for even $n$ and otherwise zero. Therefore

$\operatorname{\mathbb{E}}(L^{n})=\operatorname{\mathbb{E}}\left(((2r)^{-1/2}Y^{n})\right)=(2r)^{-n/2}\operatorname{\mathbb{E}}Y^{n}=(2r)^{-n/2}n!$
(14)

for even $n$ and zero otherwise. The second integral being the $n$-th moment
of the normal distribution $\mathcal{N}(X_{0},{t})$ can be stated as follows
Winkelbauer (2014)

$\displaystyle\operatorname{\mathbb{E}}(W^{n}_{X_{0}})=$
$\displaystyle=\begin{cases}(\sqrt{2t})^{n}\frac{\Gamma((n+1)/2)}{\sqrt{\pi}}\Phi\left(-\frac{n}{2},\frac{1}{2};-\frac{X_{0}^{2}}{2t}\right)&\text{
for even }n\\\
X_{0}(\sqrt{t})^{n-1}2^{(n+1)/2}\frac{\Gamma(n/2+1)}{\sqrt{\pi}}\Phi\left(\frac{1-n}{2},\frac{3}{2};-\frac{X_{0}^{2}}{2t}\right)&\text{
for odd }n.\end{cases}$ (15)

Here $\Phi(a,b;c)$ is the Kummer’s confluent hypergeometric function
Winkelbauer (2014).

The case $X_{0}=0$ simplifies the results even more

$\operatorname{\mathbb{E}}(W^{n}_{0})=\operatorname{\mathbb{E}}(W^{n})=t^{\frac{n}{2}}(n-1)!!$
(16)

for even $n$ and zero otherwise.

The third integral is the most interesting. As we can see, it represents the
$n$-th moment of the sum of normal random variable $W\sim\mathcal{N}(0,{t})$
and Laplace random variable $L\sim\mathcal{L}(0,(2r)^{-\nicefrac{{1}}{{2}}})$.
Note that $W$ and $L$ are independent. Let us calculate the $n$-th moment of
their sum

$\int\displaylimits_{\operatorname{\mathbb{R}}}x^{n}f_{\mathcal{N}(0,t)}*f_{\mathcal{L}}\operatorname{\operatorname{d}\\!}x=\operatorname{\mathbb{E}}((W+L)^{n}).$

Utilizing the binomial theorem we get the following result

$\displaystyle\operatorname{\mathbb{E}}((W+L)^{n})$
$\displaystyle=\operatorname{\mathbb{E}}\left(\sum\displaylimits_{k=0}^{n}\binom{n}{k}W^{k}L^{n-k}\right)=\sum\displaylimits_{k=0}^{n}\binom{n}{k}\left(\operatorname{\mathbb{E}}(W^{k})\right)\left(\operatorname{\mathbb{E}}(L^{n-k})\right)$
$\displaystyle=\sum\displaylimits_{\begin{subarray}{c}k=0,\\\ k\text{
even}\end{subarray}}^{n}\binom{n}{k}\left(\operatorname{\mathbb{E}}(W^{k})\right)\left(\operatorname{\mathbb{E}}(L^{n-k})\right),$
(17)

for even $n$ and zero otherwise. This is due to the fact that that both $W$
and $L$ are independent and symmetric with mean $0$.

Gathering all the terms (14),(15) and (17) we get the final expression for the
$n-th$ moment of $X_{t}$

$\operatorname{\mathbb{E}}(X_{t}^{n})=\operatorname{\mathbb{E}}(L^{n})+e^{-rt}\operatorname{\mathbb{E}}(W_{X_{0}}^{n})-e^{-rt}\operatorname{\mathbb{E}}\left((W+L)^{n}\right).$

We note that another way of deriving moments of resetting process is by
calculating derivatives of the MGF $M_{t}(s)$ in eq. (10).

###  II.5 Fokker-Planck equations

Having the PDF (13) of the resetting process we can directly derive the
corresponding Fokker-Planck equation. Let us calculate the necessary
derivatives of (13):

$\displaystyle\partial_{t}p(x,t)$
$\displaystyle=\partial_{t}\left(f_{\mathcal{L}}+e^{-rt}\left(f_{W_{X_{0}}}-f_{\mathcal{N}(0,t)}*f_{\mathcal{L}}\right)\right)=\partial_{t}\left(e^{-rt}\left(f_{W_{X_{0}}}-f_{\mathcal{N}(0,t)}*f_{\mathcal{L}}\right)\right)$
$\displaystyle=-re^{-rt}\left(f_{W_{X_{0}}}-f_{\mathcal{N}(0,t)}*f_{\mathcal{L}}\right)+e^{-rt}\partial_{t}\left(f_{W_{X_{0}}}-f_{\mathcal{N}(0,t)}*f_{\mathcal{L}}\right)$
$\displaystyle=-re^{-rt}\left(f_{W_{X_{0}}}-f_{\mathcal{N}(0,t)}*f_{\mathcal{L}}\right)+e^{-rt}\left(\partial_{t}f_{W_{X_{0}}}-\left(\partial_{t}f_{\mathcal{N}(0,t)}\right)*f_{\mathcal{L}}\right)$
$\displaystyle=-rp(x,t)+rf_{\mathcal{L}}+e^{-rt}\left(\partial_{t}f_{W_{X_{0}}}-\left(\partial_{t}f_{\mathcal{N}(0,t)}\right)*f_{\mathcal{L}}\right)$

The second spatial derivative equals to

$\displaystyle\partial_{xx}p(x,t)$
$\displaystyle=\partial_{xx}\left(f_{\mathcal{L}}+e^{-rt}\left(f_{W_{X_{0}}}-f_{W}*f_{\mathcal{L}}\right)\right)=\partial_{xx}f_{\mathcal{L}}+e^{-rt}\left(\partial_{xx}f_{W_{X_{0}}}-\left(\partial_{xx}f_{W}\right)*f_{\mathcal{L}}\right).$

We can now note that the PDF of the Wiener process fulfills the standard
diffusion equation

$\partial_{t}f_{W_{X_{0}}}=\frac{1}{2}\partial_{xx}f_{W_{X_{0}}},$

so the time derivative can be rewritten as

$\displaystyle\partial_{t}p(x,t)$
$\displaystyle=-rp(x,t)+rf_{\mathcal{L}}+\frac{1}{2}\left(\partial_{xx}f_{\mathcal{L}}+e^{-rt}\left(\partial_{xx}f_{W_{X_{0}}}-\left(\partial_{xx}f_{\mathcal{N}(0,t)}\right)*f_{\mathcal{L}}\right)-\partial_{xx}f_{\mathcal{L}}\right)$
$\displaystyle=\frac{1}{2}\partial_{xx}p(x,t)+\left(r-\frac{1}{2}\partial_{xx}\right)f_{\mathcal{L}}-rp(x,t)$

The derivative (in the weak sense) of the density of Laplace distribution is

$\displaystyle\partial_{x}f_{\mathcal{L}}=\partial_{x}\sqrt{\frac{r}{2}}e^{-\sqrt{2r}|x-x_{R}|}=-re^{-\sqrt{2r}|x-x_{R}|}\partial_{x}\left(|x-x_{R}|\right)=-re^{-\sqrt{2r}|x-x_{R}|}\operatorname{sgn}(x-x_{R}),$

while the second spatial derivative is

$\displaystyle\partial_{xx}f_{\mathcal{L}}$
$\displaystyle=\sqrt{2r^{3}}e^{-\sqrt{2r}|x-x_{R}|}\operatorname{sgn}^{2}(x-x_{R})-re^{-\sqrt{2r}|x-x_{R}|}2\delta(x-x_{R})$
$\displaystyle=e^{-\sqrt{2r}|x-x_{R}|}\left(2r\delta(x-x_{R})-\sqrt{2r^{3}}\operatorname{sgn}^{2}(x-x_{R})\right)\stackrel{{\scriptstyle\text{a.e.}}}{{=}}f_{\mathcal{L}}\left(2\sqrt{2r}\delta(x-x_{R})-2r\right).$

Plugging this result to the expression for time derivative we arrive at the
final Fokker-Planck equation for the resetting process

$\partial_{t}p(x,t)=\frac{1}{2}\partial^{2}_{x}p(x,t)+\sqrt{2r}\delta(x-x_{R})f_{\mathcal{L}}-rp(x,t).$

The above equation is different from Evans  _et al._ (2020), see (2), however
we can check the equivalence of both equations by looking at them in the
frequency domain via Fourier transform. We only need to look at the component
$\sqrt{2r}\delta(x-x_{R})f_{\mathcal{L}}$. We have

$\displaystyle\operatorname{\mathcal{F}}\left\\{\sqrt{2r}\delta(x-x_{R})f_{\mathcal{L}}\right\\}$
$\displaystyle=\sqrt{2r}\int\displaylimits_{\operatorname{\mathbb{R}}}\sqrt{\frac{r}{2}}e^{-\sqrt{2r}|x-x_{R}|}\delta(x-x_{R})e^{isx}\operatorname{\operatorname{d}\\!}x$
$\displaystyle=re^{-\sqrt{2r}|x_{R}-x_{R}|}e^{isx_{R}}=re^{isx_{R}}=\operatorname{\mathcal{F}}\left\\{r\delta(x-x_{R})\right\\}.$

Both equations are equal in the frequency domain, implying that their
solutions coincide. One should also underline that the difference between both
equations is only in the Dirac delta term $\delta(x-x_{R})$, which affect the
solutions only on the set of points with zero Lebesgue measure. Since the
random variable $X_{t}$ is continuous, modification of its PDF on the set of
points with zero Lebesgue measure does not change the distribution.

Let us note that the above Fokker-Planck equation can be used to derive the
stationary PDF of the resetting process just by putting $\partial_{t}p(x,t)=0$
and calculating the corresponding $p$.

###  II.6 Infinitesimal generator

The infinitesimal generator is a key operator in the theory of Markov
processes Ethier and Kurtz (2009). It contains great deal of information about
the Markov process. In particular it can be used to find the corresponding
evolution equations. Applying the general formula for jump diffusion processes
Øksendal and Sulem (2019), we get that the infinitesimal generator
$\mathcal{A}$ of the resetting process (4) equals:

$\mathcal{A}g(x)=\frac{1}{2}\partial_{xx}g(x)+r\left(g(x_{R})-g(x)\right),$
(18)

for appropriately smooth function $g$. The first component
$\frac{1}{2}\partial_{xx}g(x)$ is the diffusive part of the generator. The
second component (jump part) of the above expression can be identified as a
generator of the compound Poisson process with jump sizes $x_{R}-X_{t}$, which
corresponds to the jumps of the process to the resetting point $x_{R}$ with
intensity $r$.

Let us now derive the adjoint generator of our resetting process. It follows
from the calculations below

$\displaystyle\int f(z)\mathcal{O}g(z)\operatorname{\operatorname{d}\\!}z$
$\displaystyle=\int g(x_{R})f(z)\operatorname{\operatorname{d}\\!}z=\int
g(u)\delta(u-x_{R})\operatorname{\operatorname{d}\\!}u\int
f(z)\operatorname{\operatorname{d}\\!}z$ $\displaystyle=\iint
g(u)\delta(u-x_{R})f(z)\operatorname{\operatorname{d}\\!}z\operatorname{\operatorname{d}\\!}u=\int
g(u)\delta(u-x_{R})\int
f(z)\operatorname{\operatorname{d}\\!}z\operatorname{\operatorname{d}\\!}u$
$\displaystyle=\int g(z)\delta(z-x_{R})\int
f(u)\operatorname{\operatorname{d}\\!}u\operatorname{\operatorname{d}\\!}z=\int
g(z)\mathcal{O}^{*}f(z)\operatorname{\operatorname{d}\\!}z,$

that the adjoint generator $\mathcal{A}^{*}$ equals

$\mathcal{A}^{*}g(x)=\frac{1}{2}\partial_{xx}g(x)+r\left(\delta(x-x_{R})\int\displaylimits_{\mathbb{R}}g(z)\operatorname{\operatorname{d}\\!}z-g(x)\right).$
(19)

From the general theory of Markov processes we know that the PDF of a Markov
process with adjoint generator $\mathcal{A}^{*}$ satisfies the following
Fokker-Planck equation Pavliotis (2014)

$\partial_{t}p=\mathcal{A}^{*}p.$

Therefore, applying (19) we obtain the following Fokker-Planck formula for the
resetting process

$\partial_{t}p(x,t)=\frac{1}{2}\partial_{xx}p(x,t)+r\left(\delta(x-x_{R})\int\displaylimits_{\mathbb{R}}p(z,t)\operatorname{\operatorname{d}\\!}z-p(x,t)\right).$

Since
$\int\displaylimits_{\mathbb{R}}p(z,t)\operatorname{\operatorname{d}\\!}z=1$
we obtain the same Fokker-Planck equation as the one derived in Evans  _et
al._ (2020) for resetting processes, cf. equation (2).

###  II.7 Nonhomogeneous in time Poissonian resetting

Non-homogeneous Poisson process (NPP) is a natural generalization of the
standard Poisson process to the case of time-dependent intensity. It is
defined as a counting process with the following properties Ross (1996):

  * •

$\tilde{N}_{0}=0$,

  * •

$\tilde{N}_{t}$ has independent increments,

  * •

the increments $\tilde{N}_{t+h}-\tilde{N}_{t}$ are Poisson distributed with
mean $\int_{t}^{t+h}r(s)ds$.

Here, the non-negative function $r(t)$ is called the intensity function of
$\tilde{N}_{t}$. For higher values of $r(t)$ we observe more jumps of
$\tilde{N}_{t}$, whereas small $r(t)$ gives less jumps on the average. For
$r(t)=r=const$ we recover the standard Poisson process.

Now, we introduce a modified SDE defining our resetting process

$\operatorname{\operatorname{d}\\!}X_{t}=\operatorname{\operatorname{d}\\!}W_{t}+(x_{R}-X_{t})\operatorname{\operatorname{d}\\!}\tilde{N}_{t}.$
(20)

We put $X_{0}=x_{R}=0$ for simplicity. $\tilde{N}_{t}$ here is the NPP with
intensity function $r(t)$ and $W(t)$ is the standard Brownian motion
independent of $\tilde{N}_{t}$.

In our further analysis we will focus on the power-law intensity function and
the impact it has on both the mean square displacement (MSD) and the overall
distribution. We will thus assume that the intensity function has the form
$r(t)=r\cdot(t+1)^{p}$. The previously studied homogeneous resetting is now a
special case with $p=0$, which makes it an interesting starting consideration
point. With $p=0$ the process displayed a nondegenerate stationary
distribution, so any increase in resetting frequency should have a profound
impact on the probability law, most likely creating a degenerate distribution.
On the other hand the decrease of the resetting intensity intuitively should
recover the classical diffusion at infinity. This will be verified in detail.
We will show that the diffusive, subdifussive, stationary and deterministic
behavior can be observed for long times depending on the parameter $p$.

Now let’s define
$R(t)=\int_{0}^{t}r(\omega)\operatorname{\operatorname{d}\\!}\omega$ as the
mean number of resets up to point $t$. Using exactly the same methods and
calculations as in the homogeneous case, based on Itô lemma, we conclude that
the Fourier transform of the process must be a solution of the following
differential equation

$\partial_{t}\varphi_{t}(s)+(r(t)+\frac{1}{2}s^{2})\varphi_{t}(s)=r(t).$

Integration yields the Fourier transform

$\varphi_{t}(s)=e^{-R(t)-\frac{1}{2}ts^{2}}\left(1+\int_{0}^{t}r(w)e^{R(w)+\frac{1}{2}ws^{2}}\operatorname{\operatorname{d}\\!}w\right)$
(21)

which after simple inspection can be inverted to obtain the density of the
process given in (20)

$f_{X_{t}}(x,t)=e^{-R(t)}f_{\mathcal{N}(0,t)}(x)+e^{-R(t)}\int\displaylimits_{0}^{t}r(\omega)e^{R(\omega)}f_{\mathcal{N}(0,t-\omega)}(x)\operatorname{\operatorname{d}\\!}\omega.$
(22)

We can easily extract the MSD simply by multiplying by $x^{2}$ and
integrating. We get

$MSD(X_{t})=e^{-R(t)}MSD(W_{t})+e^{-R(t)}\int\displaylimits_{0}^{t}r(\omega)e^{R(\omega)}MSD(W_{t-\omega})\operatorname{\operatorname{d}\\!}\omega,$

which after simple calculations turns into

$MSD(X_{t})=e^{-R(t)}\int\displaylimits_{0}^{t}e^{R(\omega)}\operatorname{\operatorname{d}\\!}\omega.$
(23)

Now, let us assume that the intensity function has the form
$r(t)=r\cdot(t+1)^{p}$, where $p$ can be any real number and $r$ is a positive
constant. Then we get
$R(t)=\int_{0}^{t}r(\omega)\operatorname{\operatorname{d}\\!}\omega=\frac{r}{p+1}\left(\left(t+1)^{p+1}-1\right)\right)$
if $p\neq-1$ and $R(t)=\ln{(t+1)^{r}}$ if $p=-1$.

We already know the asymptotic behavior of $X_{t}$ for $p=0$ from previous
sections. Let us now check the case $p=-1$. We have

$MSD(X_{t})=e^{-\ln{(t+1)^{r}}}\int\displaylimits_{0}^{t}e^{\ln{(\omega+1)^{r}}}\operatorname{\operatorname{d}\\!}\omega=\frac{1}{(t+1)^{r}}\int\displaylimits_{0}^{t}(\omega+1)^{r}\operatorname{\operatorname{d}\\!}\omega$

The MSD is now equal to $\frac{t+1}{r+1}-\frac{1}{(t+1)^{r}(r+1)}$, which for
large times gives us linear diffusive scaling.

Next, if $p\in(-1,0)$ then $r(t)$ is decreasing, but $R(t)$ is increasing,
giving us the most interesting case to investigate. We obtain

$MSD(X_{t})=e^{-\frac{r}{p+1}\left(\left(t+1)^{p+1}-1\right)\right)}\int\displaylimits_{0}^{t}e^{\frac{r}{p+1}\left(\left(\omega+1)^{p+1}-1\right)\right)}\operatorname{\operatorname{d}\\!}\omega=e^{-\frac{r}{p+1}(t+1)^{p+1}}\int\displaylimits_{0}^{t}e^{\frac{r}{p+1}(\omega+1)^{p+1}}\operatorname{\operatorname{d}\\!}\omega.$

Applying the de l’Hospital rule we get for large times

$MSD(X_{t})=\frac{\int\displaylimits_{0}^{t}e^{\frac{r}{p+1}(\omega+1)^{p+1}}\operatorname{\operatorname{d}\\!}\omega}{e^{\frac{r}{p+1}(t+1)^{p+1}}}\approx\frac{e^{\frac{r}{p+1}(t+1)^{p+1}}}{r(t+1)^{p}e^{\frac{r}{p+1}(t+1)^{p+1}}}\approx\frac{1}{r}(t+1)^{-p}.$

The process now displays subdiffusive behavior with exponent $-p\in(0,1)$.

Consequently, for $p<-1$ we obtain for large times

$MSD(X_{t})=e^{-\frac{r}{p+1}(t+1)^{p+1}}\int\displaylimits_{0}^{t}e^{\frac{r}{p+1}(\omega+1)^{p+1}}\operatorname{\operatorname{d}\\!}\omega\approx\int\displaylimits_{0}^{t}e^{\frac{r}{p+1}(\omega+1)^{p+1}}\operatorname{\operatorname{d}\\!}\omega$

Again, applying the de l’Hospital rule we get

$MSD(X_{t})\approx t$

for large times. As expected, the process displays normal diffusive scaling of
the MSD.

Examining the case $p>0$ in analogous way, we find that
$MSD(X_{t})\approx\frac{1}{r}(t+1)^{-p}$, which converges to 0 for large
times.

Let us now focus on the asymptotic distributions. Again, we expect different
results for different $p$. The case $p=0$ corresponds to the homogeneous
Poisson process and was analyzed previously, giving rise to the Laplace
stationary distribution.

For $p>0$ we get from (21) and from de l’Hospital rule that the Fourier
transform satisfies

$\varphi_{t}(s)\rightarrow 1$

for large times. This implies that $X_{t}\stackrel{{\scriptstyle
d}}{{\rightarrow}}x_{R}=0$ as $t\rightarrow\infty$. Here
$\stackrel{{\scriptstyle d}}{{\rightarrow}}$ means convergence in
distribution. This result agrees with our intuition. For $p>0$ the number of
resetting times increases with time, therefore $X_{t}$ keeps coming back to
$x_{R}$ more and more often.

For $p\in(-1,0)$ we get from (21) and from de l’Hospital rule that the Fourier
transform satisfies

$\varphi_{t}(s)\approx\frac{1}{1+\frac{s^{2}}{2r(t+1)^{p}}}$

for large times. This means that the distribution of $X_{t}$ is asymptotically
Laplace with subdiffusive $MSD(X_{t})\approx t^{-p}$ for large times.

Ananlogous result is obtained for $p=-1$:

$\varphi_{t}(s)\approx\frac{1}{1+\left(\frac{t+1}{2r}\right)s^{2}}$

for large times, with linear in time asymptotic MSD.

Figure 4: A Monte-Carlo simulation results (solid lines) compared with the
theoretical asymptotics (dotted lines) of the MSD of Brownian motion under
nonhomogeneous Poisson resetting with $r(t)=r\cdot(t+1)^{p}$, calculated for
different $p\leq 0$. The parameter $\mu$ here is the corresponding power-law
exponent of the theoretical asymptotic MSD. The time horizon of the
simulations was of order $10^{2}$, while the number of samples was equal to
$10^{4}$. We can observe the convergence of the anomalous diffusion exponents
to the theoretical ones. Figure 5: A Monte-Carlo simulation (solid lines) of
the MSD of the process with nonhomogeneous Poisson resetting with intensity
function $r(t)=r\cdot(t+1)^{p}$, calculated for different $p>0$. Asterixes are
the corresponding theoretical values of the MSD obtained via numerical
evaluation of (23). We can see the convergence of the MSD to 0, indicating the
convergence to the degenerate distribution concentrated in $x_{R}=0$. The
convergence rate increases with the increase of the $p$ parameter.

Looking at $p<-1$ we get that for large times the corresponding characteristic
function satisfies

$\varphi_{t}(s)\approx
e^{-R(t)}e^{-ts^{2}/2}+\frac{r(t)(t+1)}{p}\frac{1}{1+\frac{t+1}{2p}s^{2}},$

which is a combination of Gaussian and Laplace distributions. The
corresponding MSD is asymptotically linear in time. Note that for
$p\rightarrow-\infty$, i.e. when the resetting events vanish, only the
diffusive Brownian part is left in the above formula.

Summarizing the above findings for large times:

  * •

for $p>0$ we observe degenerate distribution concentrated in $x_{R}$ with
$MSD(X_{t})\approx const$;

  * •

for $p=0$ a stationary Laplace distribution with $MSD(X_{t})\approx const$ is
observed;

  * •

for $p\in(-1,0)$ we get a non-stationary Laplace distribution with
subdiffusive $MSD(X_{t})\approx t^{-p}$;

  * •

for $p=-1$ we get a non-stationary Laplace distribution with diffusive
$MSD(X_{t})\approx t$;

  * •

for $p<-1$ we observe a combination of normal and Laplace distributions with
diffusive $MSD(X_{t})\approx t$. As $p\rightarrow-\infty$, the standard
Brownian diffusion is recovered.

In Figs. 4–5 we can observe the approximated MSD with a fitted (anomalous)
diffusion exponent. We observe perfect agreement between simulation and
theory.

##  III Conclusions

Summarizing, we have introduce a general stochastic representation for
processes with resetting. It allows to analyze such processes using stochastic
differential equations, Itô lemma and related tools. We have shown robustness
of our approach for the case of homogeneous, and nonhomogeneous Poissonian
resetting. We would like to underline that the derived in this paper results
build a general framework for the analysis of stochastic processes with
intermittent random resetting. It allows to analyze processes with resetting
both, analytically and using Monte Carlo simulation methods. Moreover, the
presented approach can be easily generalized to the case of non-Markovian
resetting times as well as other, not necessarily Brownian, driving processes
(Lévy flights, Lévy processes, fractional processes, etc.).

###### Acknowledgements.

This research was partially supported by NCN Sonata Bis-9 grant nr
2019/34/E/ST1/00360.

## References

  * Campos  _et al._ (2014) D. Campos, F. Bartumeus, V. Méndez, and X. Espadaler, Journal of The Royal Society Interface 11, 20130859 (2014).
  * Manrubia and Zanette (1999) S. C. Manrubia and D. H. Zanette, Physical Review E 59, 4945 (1999).
  * Gelenbe (2010) E. Gelenbe, Physical Review E 82, 061112 (2010).
  * Kyriakidis (1994) E. Kyriakidis, Statistics & Probability Letters 20, 239 (1994).
  * Evans and Majumdar (2011) M. R. Evans and S. N. Majumdar, Journal of Physics A: Mathematical and Theoretical 44, 435001 (2011).
  * Evans  _et al._ (2020) M. R. Evans, S. N. Majumdar, and G. Schehr, Journal of Physics A: Mathematical and Theoretical 53, 193001 (2020).
  * Tong  _et al._ (2008) H. Tong, C. Faloutsos, and J. Pan, Knowledge and Information Systems 14, 327 (2008).
  * Kussell and Leibler (2005) E. Kussell and S. Leibler, Science 309, 2075 (2005).
  * Viswanathan  _et al._ (2011) G. M. Viswanathan, M. G. Da Luz, E. P. Raposo, and H. E. Stanley,  _The physics of foraging: an introduction to random searches and biological encounters_ (Cambridge University Press, 2011).
  * Reuveni  _et al._ (2014) S. Reuveni, M. Urbakh, and J. Klafter, Proceedings of the National Academy of Sciences 111, 4391 (2014).
  * Montanari and Zecchina (2002) A. Montanari and R. Zecchina, Physical Review Letters 88, 178701 (2002).
  * Montero  _et al._ (2017) M. Montero, A. Masó-Puigdellosas, and J. Villarroel, The European Physical Journal B 90, 1 (2017).
  * Bodrova  _et al._ (2019) A. S. Bodrova, A. V. Chechkin, and I. M. Sokolov, Physical Review E 100, 012119 (2019).
  * Zhou  _et al._ (2020) T. Zhou, P. Xu, and W. Deng, Physical Review Research 2, 013103 (2020).
  * Zhou  _et al._ (2021) T. Zhou, P. Xu, and W. Deng, arXiv preprint arXiv:2105.12571 (2021).
  * Eliazar and Klafter (2004) I. Eliazar and J. Klafter, Physica A: Statistical Mechanics and its Applications 334, 1 (2004).
  * Eliazar and Klafter (2005) I. Eliazar and J. Klafter, Journal of statistical physics 118, 177 (2005).
  * Eliazar and Klafter (2006) I. Eliazar and J. Klafter, Physica A: Statistical Mechanics and its Applications 367, 106 (2006).
  * Ken-Iti (1999) S. Ken-Iti,  _Lévy Processes and Infinitely Divisible Distributions_ (Cambridge university press, 1999).
  * Einstein (1905) A. Einstein, Annalen der Physik 322, 549 (1905).
  * Smoluchowski (1906) M. Smoluchowski, Annalen der Physik 326, 756 (1906).
  * Bachelier (1900) L. Bachelier, Annales scientifiques de l’École Normale Supérieure 3e série, 17, 21 (1900).
  * Merton (1976) R. Merton, Journal of Financial Economics 3, 125 (1976).
  * Black and Scholes (1973) F. Black and M. Scholes, Journal of Political Economy 81, 637 (1973).
  * Barndorff-Nielsen  _et al._ (2012) O. E. Barndorff-Nielsen, T. Mikosch, and S. I. Resnick,  _Lévy processes: Ttheory and Applications_ (Springer Science & Business Media, 2012).
  * Itô (1944) K. Itô, Proceedings of the Imperial Academy 20, 519 (1944).
  * Itô (1946) K. Itô, Proceedings of the Japan Academy 22, 32 (1946).
  * Øksendal and Sulem (2019) B. Øksendal and A. Sulem,  _Applied Stochastic Control of Jump Diffusions; 3rd ed._, Universitext (Springer, Cham, 2019).
  * Eliazar and Klafter (2003) I. Eliazar and J. Klafter, Journal of Statistical Physics 111, 739 (2003).
  * Ge  _et al._ (2012) H. Ge, M. Qian, H. Qian, C. Hartnack, H. Oeschler, Y. Leifels, E. Bratkovskaya, and J. Aichelin, Physics Reports 510, 201 (2012).
  * Billingsley (1986) P. Billingsley,  _Probability and Measure_ , 2nd ed. (John Wiley and Sons, 1986).
  * Reed and Hughes (2002) W. J. Reed and B. D. Hughes, Physical Review E 66, 067103 (2002).
  * Reed (2003) W. J. Reed, Physica A: Statistical Mechanics and its Applications 319, 469 (2003).
  * Geraci (2017) M. Geraci, “Mixed-effects models using the normal and the laplace distributions: A $\mathbf{2\times 2}$ convolution scheme for applied research,” (2017), arXiv:1712.07216 [stat.ME] .
  * Winkelbauer (2014) A. Winkelbauer, ArXiv abs/1209.4340 (2014).
  * Ethier and Kurtz (2009) S. N. Ethier and T. G. Kurtz,  _Markov Processes: Characterization and Convergence_ , Vol. 282 (John Wiley & Sons, 2009).
  * Pavliotis (2014) G. A. Pavliotis,  _Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations_ , Vol. 60 (Springer, 2014).
  * Ross (1996) S. M. Ross,  _Stochastic Processes_ (Wiley, 1996).


# Parametric triggering of vortices in toroidally trapped rotating Bose-
Einstein Condensates

M. Arivazhagan  P. Muruganandam<EMAIL_ADDRESS>N. Athavan
<EMAIL_ADDRESS>Department of Physics, H. H. The Rajah’s College
(Autonomous), Pudukkottai 622001, (Affiliated to Bharathidasan University,
Tiruchirapalli), Tamilnadu, India  Department of Physics, Bharathidasan
University, Tiruchirappalli 620024, Tamilnadu, India

###### Abstract

We study the creation of vortices by triggering the rotating Bose-Einstein
condensates in a toroidal trap with trap parameters such as laser beam waist
and Gaussian potential depth. By numerically solving the time-dependent Gross-
Pitaevskii equation in two dimensions, we observe a change in vortex structure
and a considerable increase in the number of vortices when the waist of the
irradiated laser beam is in consonance with the area of the condensate as we
vary the Gaussian potential depth. By computing the root mean square radius of
the condensate, we confirm the variation in the number of vortices generated
as a function of the ratio between the root-mean-square radius of the
condensate and the laser beam waist. In particular, the number of hidden
vortices reaches the maximum value when the above ratio is close to the value
0.7. We find the variation in the number of vortices is rapid for deeper
Gaussian potentials, and we conclude that the larger beam waist and deeper
Gaussian potentials generate more vortices. Further, we calculate the number
of vortices using the Feynman rule with Thomas Fermi approximation and compare
them with the numerical results. We also observe that the critical rotation
frequency decreases with an increase in depth of Gaussian potential.

###### keywords:

Bose-Einstein condensates, vortices, toroidal trap, Gross-Pitaevskii equation

##  1 Introduction

A most intriguing occurrence in the realm of rotating Bose-Einstein
condensates (BECs) is the creation of quantized vortices like superfluidity of
liquid helium [1]. Quantum vortices are nonlinear physical entities that arise
due to the excitation of BECs [2, 3, 4, 5, 6, 7]. Experimentally vortices in
BEC have been generated by either applying rotating magnetic traps or laser
stirring, and they are usually formed above a critical rotation frequency. In
a BEC, the vortices have quantized circulation, and their interaction on
inhomogeneous density background have been studied extensively [8, 9, 10]
complementing the theoretical investigations [11, 12] over the past few
decades.

On the other hand, the observation of BEC in a toroidal geometry has opened up
a new window for remarkable studies about their fundamental properties
including persistent currents [13, 14], ring dark solitons [15, 16, 17], weak
links [18, 19, 20], self interference [21], symmetry breaking [22], atomic-
phase interference [23], decay of persistent current via phase slips [24],
creation of a vortex-antivortex pair, multiple quantized vortices [25, 26, 27,
28], hidden vortices [29, 30], transition from vortex excitation to solid-body
like motion [31], transition from quantum anti-resonance to dynamical
localization [32], half-quantum vortex and complex vortex lattice [33].
Toroidally trapped Bose-Einstein condensates have been experimentally realized
by superimposing a laser beam on a harmonically trapped condensate, thus
creating a harmonic plus Gaussian potential [34, 35, 36, 37], which also
favours the occurrence of vortices under suitable conditions.

By using a toroidal trap geometry, one can show that the formed condensate can
either be like a disc or an annulus depending upon the parameters of the
system, such as rotational frequency ($\Omega$), intensity ($V_{0}$) and waist
of Gaussian potential depth of the laser beam [38]. To be precise, the
condensate bears a disc-like structure when the rotational frequency
($\Omega$) is notably smaller than radial trap frequency ($\omega_{\perp}$),
and on the other hand, rapid rotation ($\Omega/\omega_{\perp}\sim 1$) makes
the condensate to expand resulting in a large annulus with a vortex lattice
inside the condensate, and a large circulation within the central hole [39,
40]. Moreover, the vortex dynamics of rotational BECs trapped by harmonic plus
Gaussian potential for varied rotational velocity has been recently analysed
in reference [38]. However, they have investigated the dynamics of the
vortices for different rotational frequencies by keeping the beam waist
constant. In contrast to these above investigations, we intend to study the
vortex dynamics of BECs in such a toroidal trap as a function of trap
parameters, namely the potential depth and waist of Gaussian potential under
constant rotation. In particular, by altering the trap parameters suitably, we
explore the possibility of creating more number of new vortices in the
rotating BECs and the maximum extent of vortex number for a particular choice
of beam waist. Additionally, we analyze the relationship between the critical
rotational frequency and the depth of the Gaussian potential.

By numerically solving the quasi-two-dimensional Gross-Pitaevskii equation, we
calculate the number of vortices as a function of the potential (Gaussian)
depth for different waist lengths of Gaussian laser, and compare them with
that of calculated values by adhering to the Feynmann rule with the Thomas-
Fermi (TF) approximation. We have found that the increase in Gaussian
potential depth hastens the vortex formation in BECs. Also, we have observed
that the critical rotation frequency decreases when the depth of the Gaussian
potential increases for the same beam waist. We have calculated other
associated quantities, such as the expectation value of angular momentum, root
mean square (RMS) radius, and chemical potential of the BEC as a function of
Gaussian potential depth ($V_{0}$) separately for three different beam waist
values.

This paper is organised as follows. In section 2, we provide the theoretical
background by describing the dynamics of rotating BEC confined in an axially
symmetric toroidal trap with the aid of the mean-field Gross-Pitaevskii
equation. Then, in section 3, we numerically illustrate the formation of
vortices and their dynamical behaviour with respect to the change in the
quantities which shape the Gaussian potential. We also investigate the role
played by the waist of the Gaussian laser beam on the enhancement of the
number of vortices. Finally, in section 4 we provide a summary and conclusion.

##  2 Theoretical background

###  2.1 Evolution Equation

According to the mean-field approximation, the dynamics of a rotating BEC with
$N$ atoms, each of mass $m$ at an absolute zero temperature in a rotating
frame can be described by the dimensionless form of the 2D Gross-Pitaevskii
(GP) equation, which reads as [41, 42]

$\displaystyle i\frac{\partial\phi(r,t)}{\partial
t}=\left[-\frac{1}{2}\nabla^{2}+V(r)+g|\phi(r,t)|^{2}-\Omega
L_{z}\right]\phi(r,t),$ (1)

with $i=\sqrt{(-1)}$. Here $V(r)$, $r\equiv(x,y)$, is the external trapping
potential, $g=4\pi aN/\sqrt{2\pi}d_{z}$ is a two-dimensional coupling
parameter, $\phi(r,t)$ is the wave function at time $t$, and $a$ is atomic
scattering length which can be tuned via Feshbach resonance.
$L_{z}=-i(x\partial_{y}-y\partial_{x})$ corresponds to the $z$-component of
the angular momentum due to the rotation of the BEC about $z$ axis with
angular velocity $\Omega$. The radial and axial trap frequencies of the
harmonic potential, $V_{ha}$, are $\omega_{\perp}$ and $\omega_{z}$,
respectively, and are related with the trap aspect ratio,
$\lambda\equiv\omega_{z}/\omega_{\perp}$ and $d_{z}=1/\sqrt{\lambda}$. It is
to be noted that in Eq. (1) length is measured in units of characteristic
harmonic oscillator length $\sqrt{\hbar/m\omega_{\perp}}$, frequency in units
of $\omega_{\perp}$ and time in units of $\omega_{\perp}^{-1}$.

The normalization condition of the wave function is

$\displaystyle\int|\phi(r,t)|^{2}d\tau=1,$ (2)

where $d\tau$ is the spatial integral element and the chemical potential in
this case for the state $\phi(r,t)$ can be written as

$\displaystyle\mu=\int\left[\frac{1}{2}|\nabla\phi|^{2}+V(r)|\phi|^{2}+g|\phi|^{4}-{\phi^{*}}\Omega
L_{z}\phi\right]d\tau,$ (3)

where ${\phi^{*}}$ is the complex conjugate of the wave function $\phi$ and
$L_{z}$ is the $z$ component of the angular momentum.

The dimensionless 2D toroidal trap potential $V(r)$ is given by

$\displaystyle V(r)=\frac{r^{2}}{2}+V_{0}\exp(-2r^{2}/w_{0}^{2}),$ (4)

with Gaussian potential depth $V_{0}$ with units of $\hbar\omega_{\perp}$, and
$w_{0}$ is the waist of laser beam in units of $\sqrt{\hbar/m\omega_{\perp}}$.

###  2.2 Root mean square radius of the condensate

We determine the rms radius ($r_{rms}$) with respect to $V_{0}$ in the mean-
field TF regime. When the interaction energy is larger than the kinetic
energy, the kinetic energy can be neglected and in doing so we enter into the
TF regime in Eq. (1) and therefore the approximation yields the expression for
the TF density as [38]

$\displaystyle
g|\phi_{TF}|^{2}=\mu_{TF}-\frac{1}{2}(1-\Omega^{2})r^{2}-V_{0}e^{-\frac{2r^{2}}{w_{0}^{2}}},$
(5)

where $\mu_{TF}$ is the chemical potential in the TF regime.

As a next step we deduce $r_{rms}$ for the two limiting values of laser waist
$w_{0}$. When $w_{0}\gg r$, or in other words, the laser beam waist is
comparatively larger than the outer TF radius of the condensate,
$2r^{2}/w_{0}^{2}$ is very small, that is $(\ll 1)$, and hence the expansion
of the exponential term in the expression for the TF density in Eq. (5) yields

$\displaystyle
g|\phi_{TF}|^{2}=\mu_{TF}^{{}^{\prime}}-\frac{1}{2}\left(1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}\right)r^{2},$
(6)

where $\mu_{TF}^{{}^{\prime}}=\mu_{TF}-V_{0}$.

When $4V_{0}/w_{0}^{2}$ becomes larger than $1-\Omega^{2}$, we consider the
next higher order term of the exponential function in Eq. (5).

It may be noted that the density of the condensate $|\phi_{TF}|^{2}$ vanishes
at radius $R_{TF}$, which in turn leads to the following chemical potential
relation with frequency of rotation ($\Omega$) from Eq. (6),

$\displaystyle\mu_{TF}^{{}^{\prime}}(\Omega)=\frac{1}{2}\left(1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}\right)R_{TF}^{2}(\Omega).$
(7)

By using the above Eq. (7), we get the ratio of square of the TF radius with
and without rotation as

$\displaystyle\frac{R^{2}_{TF}(\Omega)}{R^{2}_{TF}(0)}=\frac{\mu_{TF}^{{}^{\prime}}(\Omega)}{\mu_{TF}^{{}^{\prime}}(0)}\left(\frac{1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}}{1-\frac{4V_{0}}{w_{0}^{2}}}\right)^{-1}.$
(8)

After substituting the expression for $|\phi_{TF}|^{2}$ in Eq. (6) into the
normalization condition in Eq. (2) and integrating it over the domain between
$0$ and $R_{TF}$, we get

$\displaystyle\mu_{TF}^{{}^{\prime}}(\Omega)-\frac{1}{4}\left(1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}\right)R_{TF}^{2}(\Omega)=\frac{g}{\pi
R_{TF}^{2}(\Omega)}.$ (9)

We now substitute for $R_{TF}^{2}$ from Eq. (7) in Eq. (9) and get

$\displaystyle\mu_{TF}^{{}^{\prime}}{}^{2}(\Omega)=\frac{g}{\pi}\left(1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}\right).$
(10)

Then the ratio of chemical potentials with and without rotation can be written
as

$\displaystyle\frac{\mu_{TF}^{{}^{\prime}}(\Omega)}{\mu_{TF}^{{}^{\prime}}(0)}=\left(\frac{1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}}{1-\frac{4V_{0}}{w_{0}^{2}}}\right)^{\frac{1}{2}}.$
(11)

Now putting Eq. (11) in Eq. (8), we get the TF radius as a function of
rotation ($\Omega$) as

$\displaystyle
R_{TF}(\Omega)=R_{TF}(0)\left(\frac{1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}}{1-\frac{4V_{0}}{w_{0}^{2}}}\right)^{-\frac{1}{4}}.$
(12)

In order to get an expression for $R_{TF}(0)$, we equate Eqs. (7) and (10)
after substituting $\Omega=0$ and obtain
$R_{TF}(0)=\left[4g/\pi(1-\frac{4V_{0}}{w_{0}^{2}})\right]^{\frac{1}{4}}$. By
using the same analysis, one can obtain the expression for the outer radius of
the condensate $(R_{2})$ and subsequently the inner radius $(R_{1})$ after
eliminating the higher order terms in Eq. (18). If $w_{0}\ll r$, that is the
laser beam waist is relatively small compared to the outer TF radius of the
condensate, $2r^{2}/w_{0}^{2}$ becomes very large, that is $(\gg 1)$, and the
exponential term in Eq. (5) can be neglected. Finally, we get radius of the
condensate as

$\displaystyle R_{TF}(\Omega)=R_{TF}(0)(1-\Omega^{2})^{-1/4},$ (13)

where $R_{TF}(0)=(4g/\pi)^{1/4}$. One can easily obtain the expressions for
outer and inner radii of the condensate using the expressions given in Eqs.
(21) and (22) of Appendix A, respectively. $R_{TF}(0)$ is related to root mean
square radius $r_{rms}$ by $R_{TF}(0)=\sqrt{3}r_{rms}$ [43, 44].

###  2.3 Vortex Number

The dimensionless form of Feynman rule under equilibrium for the number of
vortices in a superfluid rotating in a rigid container is given by
$N_{v}=\Omega R_{\perp}^{2}(\Omega)$, where $R_{\perp}(\Omega)$ is the radius
of the rotating fluid [45]. As one needs a radius value for the numerical
estimation of the number of vortices, root-mean-square radius value may be
conveniently chosen for obtaining the number of vortices closer to the values
suggested by Feynman rule as in ref.[44].

Substituting the expression for $R_{TF}^{2}(\Omega)$ from Eqs. (12) and (13)
the Feynman’s relation (number of vortices) for the case of rotating BEC
trapped under harmonic and Gaussian potential with large laser waist and with
small laser waist are

$\displaystyle N_{v}=\Omega
R_{\perp}^{2}(0)\left(\frac{1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}}{1-\frac{4V_{0}}{w_{0}^{2}}}\right)^{-1/2},$
(14)

and

$\displaystyle N_{v}=\Omega R_{\perp}^{2}(0)(1-\Omega^{2})^{-1/2},$ (15)

respectively.

###  2.4 Expectation value of Angular momentum

The birth of the new vortices can be identified by calculating the increase in
the value of average angular momentum per atom ($\langle L_{z}\rangle$) as we
increase the Gaussian potential depth ($V_{0}$). $\langle L_{z}\rangle$ is
calculated by using the following expression

$\displaystyle\left\langle
L_{z}\right\rangle=i\int\phi^{*}(r,t)\left(y\partial_{x}-x\partial_{y}\right)\phi(r,t)\,d\tau.$
(16)

##  3 Numerical illustration of vortex dynamics

To numerically illustrate the vortex dynamics of BECs under toroidal trap
which mimics the real situation, we consider the respective dimensionless
equivalent numerical values of the parameters used by the authors in reference
[34] during their experimental investigations on the rotating properties of
Bose-Einstein condensate of ${}^{87}Rb$ atoms confined in toroidal geometry
with a weak laser beam. The characteristic harmonic oscillator length of
${}^{87}Rb$ atoms is $1.24\mu m$ and the trap aspect ratio of the condensate
used in their experiment was equivalent to $\lambda=0.146$. Precisely, they
have used a laser power of $1.2$mW with waist $w_{0}=25\mu m$ from which one
can obtain the dimensionless values for the Gaussian potential depth $V_{0}$
and the laser waist $w_{0}$ as $24.8$ and $20.2$, respectively. Further, we
consider $\Omega=0.5$ to match the experimental value of 37.75Hz and the two-
dimensional coupling parameter $g\sim 500$ for $60000$ ${}^{87}Rb$ atoms. All
our numerical simulations have been carried out with $dx=dy=0.2$ (space step)
and $dt=0.005$ (time step).

###  3.1 Effect of varying the depth and waist of Gaussian potential on the
shape of the condensate and number of vortices

To start with, we obtain the ground state wavefunction by solving Eq. (1)
using imaginary time propagation in the presence of harmonic and Gaussian
potential without rotation $(\Omega=0)$. Then, the obtained ground state wave
function has been allowed to evolve in real time with constant rotation
(frequency of rotation $\Omega=0.5$). The dissipation is introduced by
replacing ‘$i$’ with ‘($i-\gamma$)’ in the time dependent equation (1) to
achieve the equilibrium state, where $\gamma$ ($\sim 10^{-5}$) accounts for
the dissipation strength. We have fixed the coupling parameter $g=500$. The
shape of the condensate can be tailor made by choosing the appropriate values
for the intensity and the waist of the Gaussian beam such that one can get
either a disc or an annular shaped condensate. At first instance, we have
investigated the effect on the number of vortices $N_{v}$ for various laser
waist values $w_{0}=20.2$, 9.43, 4.71 and 2.36 with Gaussian potential depth
fixed at $V_{0}=50$. When we take the largest laser waist $w_{0}=20.2$, we
obtain the corresponding rms radius of the condensate ($r_{rms}$) as 3.83 and
the density of the condensate is maximum at the centre (disk) with eight
vortices, as shown in Fig. 1(a), and the corresponding phase plot has been
shown in Fig. 1(e).

Figure 1: Contour plots of the density distributions $|\phi|^{2}$ showing the
steady vortex state in a rotating BEC under toroidal trap with $g=500$ and
$\Omega=0.5$. For $V_{0}=50$, (a) $w_{0}=20.2$, (b) $w_{0}=9.43$, (c)
$w_{0}=4.71$ and (d) $w_{0}=2.36$. (e)-(h) are the corresponding phase
distributions of (a)-(d).

Next, when we reduce $w_{0}$ to $9.43$, the rms radius ($r_{rms}$) increases
to $6.96$ and the condensate becomes annular with no visible vortices
($N_{vv}$) as shown in Fig. 1(b). However, with the aid of the phase plot, we
identify the presence of $20$ hidden vortices ($N_{hv}$) which can be seen
from Fig. 1(f). Further reduction in $w_{0}$ to $4.71$ and then to $2.36$
resulted in the decline of the rms radius ($r_{rms}$) to $5.48$ and $3.98$,
respectively. By using the density plot Figs. 1(c)-1(d), we show that the size
of the annulus got decreased, and the corresponding phase plot Figs. 1(g)-1(h)
show that the number of hidden vortices reduces to $12$ and $6$, respectively.

Figure 2: Plot of the number of vortices as a function of the ratio of rms
radius of condensate and waist of Gaussian beam ($r_{rms}/w_{0}$) for
$V_{0}=50$ and 150. Red colour curves indicate the visible vortices ($N_{vv}$)
and green colour curves represent the hidden vortices ($N_{hv}$).

The plot, for the number of vortices $(N_{v}=N_{vv}+N_{hv})$ with respect to
the ratio between the rms radius of condensate ($r_{rms}$) and the waist of
Gaussian beam ($w_{0}$), in Fig. 2 clearly shows the number of hidden vortices
($N_{hv}$) is maximum when $r_{rms}/w_{0}$ is close to value of 0.7 and
$N_{v}$ decreases as $r_{rms}/w_{0}$ deviates from that value. Interestingly,
from Fig. 2 one can observe an abrupt variation in the number of vortices for
deeper Gaussian potential.

Next, we investigate the effect on the number of vortices $N_{v}$ for a select
set of Gaussian potential depths in multiples of 30, say, $V_{0}=30$, $60$,
$90$ and $120$ with different Gaussian beam waists, for instance,
$w_{0}=20.2$, $4.71$ and $0.47$.

Figure 3: Contour plots of the density distributions $|\phi|^{2}$ showing the
steady vortex state in a rotating BEC under toroidal trap with $g=500$ and
$\Omega=0.5$. For $w_{0}=20.2$, (a) $V_{0}=30$, (b) $V_{0}=60$, (c) $V_{0}=90$
and (d) $V_{0}=120$. (e)-(h) are the corresponding phase distributions of
(a)-(d)

.

In the case of the largest laser waist $w_{0}=20.2$, we have seen that the
density of the condensate is maximum at the centre for Gaussian depths
$V_{0}=30$ and $V_{0}=60$ and the corresponding number of vortices in the
condensate are eight and twelve as given in Figs. 3(a) and 3(b), respectively.
Figs. 3(e) and 3(f) show the respective phase plots of Figs. 3(a) and 3(b).
When we increase the Gaussian depth to $V_{0}=90$, the number of vortices
increases to $22$ with local minimum density at the centre of condensate as
shown in Fig. 3(c). Fig. 3(g) shows the respective phase pattern of Fig. 3(c).
Finally, when we increase the Gaussian depth to $V_{0}=120$, the condensate
takes up the annular shape with the emergence of 14 visible vortices along the
ring as shown in Fig. 3(d). In fact, we have observed, by using the phase plot
Fig. 3(h), $24$ hidden vortices in the annulus which are clustered together,
forming a large hole at the centre.

Now, by considering the intermediate value of laser waist, that is,
$w_{0}=4.71$, we note that, the condensate becomes annular and vortices are
hard to see in the condensate as exposed in Figs. 4(a)-4(d).

Figure 4: Contour plots of the density distributions $|\phi|^{2}$ showing the
steady vortex state in a rotating BEC under toroidal trap with $g=500$ and
$\Omega=0.5$. For $w_{0}=4.71$, (a) $V_{0}=30$, (b) $V_{0}=60$, (c) $V_{0}=90$
and (d) $V_{0}=120$. (e)-(h) are the corresponding phase distributions of
(a)-(d).

For $V_{0}=30$ and $60$, we identify $12$ hidden vortices from the phase plot
Figs. 4(e)-4(f). When $V_{0}=90$ and $120$, we witnessed more number of hidden
vortices (about $16$ to $20$ ) by analysing the corresponding phase plots
shown in Figs. 4(g)-4(h).

By considering the smallest laser waist of Gaussian laser beam $w_{0}=0.47$ to
suit the conditions for condensate formation, we plot the density plots
5(a)-5(c), which show zero density at the centre forming thin annular having
four visible vortices that appear in the condensate for $V_{0}=30$, $60$ and
$90$, two visible and two hidden vorices are found in 5(d) for $V_{0}=120$.
However the total number of vortices in the condensate are invariable in all
four values. Figs. 5(e)-5(h) show the respective phase plots of Figs.
5(a)-5(d).

Figure 5: Contour plots of the density distributions $|\phi|^{2}$ showing the
steady vortex state in a rotating BEC under toroidal trap with $g=500$ and
$\Omega=0.5$. For $w_{0}=0.47$, (a) $V_{0}=30$, (b) $V_{0}=60$, (c) $V_{0}=90$
and (d) $V_{0}=120$. (e)-(h) are the respective phase distributions of (a)-(d)

.

From the three sets of the laser waist values ($w_{0}=20.2$, $4.71$ and
$0.47$) along with the corresponding four sets of Gaussian potential depths
($V_{0}=30$, $60$, $90$ and $120$), we conclude that, the number of vortices
is larger for the intermediate laser waist, that is for $w_{0}=4.71$. However,
when the Gaussian potential is hiked from the minimum value, say $V_{0}=30$,
the growth rate of the number of vortices is larger for the larger laser
waist, that is for $w_{0}=20.2$. This is because when Gaussian potential depth
($V_{0}$) increases, the rms radius of the condensates approaches the value of
the waist of the laser beam corresponding to $w_{0}=20.2$. For $w_{0}=4.71$,
the rms radius value exceeds the value of the waist resulted in poor
overlapping. Finally, when we take $w_{0}=0.47$, rms radius is still larger
than the waist which means that a pinpoint like-overlap of Gaussian beam on
the condensate takes place.

Table 1: The comparison of vortex numbers with respect to Gaussian potential depth ($V_{0}$) for different laser waists ($w_{0}$) $w_{0}$ | $N_{v}$
---|---
$V_{0}=30$ | $V_{0}=60$ | $V_{0}=90$ | $V_{0}=120$ | $V_{0}=150$
20.20 | 8 | 12 | 22 | 38 | 48
4.71 | 12 | 12 | 16 | 20 | 24
0.47 | 4 | 4 | 4 | 4 | 4

Hence, even for a larger Gaussian potential depth, we could not find any
elevation vortex number. The comparison between the number of vortices for
different choices of waist and potential depth are presented in table 1.

###  3.2 Root mean square radius and chemical potential

To understand the variation in the number of vortices with Gaussian potential
depth ($V_{0}$), we calculate the rms radius ($r_{rms}$) for various values of
$V_{0}$. We plot the calculated values of rms radius, number of vortices and
chemical potential as a function of Gaussian potential depth for the two
extreme values of $w_{0}$, say, $20.2$ and $0.47$ in Figs. (6), (7) and (8),
respectively.

Figure 6: Plot of rms radius ($r_{rms}$) as a function of Gaussian potential
depth ($V_{0}$) for $w_{0}=20.2$ and $w_{0}=0.47$. Red colour curves represent
the numerically calculated values and blue colour curves denote the TF
approximated values. Figure 7: Plot of the number of vortices ($N_{v}$) as a
function of Gaussian potential depth ($V_{0}$) for $w_{0}=20.2$ and
$w_{0}=0.47$. Red colour curves represent the numerically calculated values
and blue colour curves denote the TF approximated values. Figure 8: Plot of
the chemical potential ($\mu$) versus Gaussian potential depth ($V_{0}$) for
$w_{0}=20.2$ and $w_{0}=0.47$. Red colour curves represent the numerically
calculated values and blue colour curves denote the TF approximated values.

In Fig. 6, we plot $r_{rms}$ as a function of $V_{0}$ for largest and smallest
$w_{0}$ values, that is, $w_{0}=20.2$ and $0.47$. The curves drawn in dashed
line with red open circles are for the numerically calculated values, while
the solid line with blue filled circles correspond to the TF approximated
values. Fig. 6 shows that the rms radius, ($r_{rms}$), increases exponentially
for $w_{0}=20.2$ while it remains constant for $w_{0}=0.47$. Consequently, the
constancy of rms radius for $w_{0}=0.47$ irrespective of the values of
Gaussian potential depth reveals the fact that the number of vortices within
the condensate remains constant for all values of $V_{0}$ which can be seen
from Fig. 7.

As far as the chemical potential is concerned, it linearly increases with
Gaussian potential depth for $w_{0}=20.2$ as shown in Fig. 8 due to the
heavier overlap of the condensate and the beam. While the curve stays flat for
$w_{0}=0.47$ as depicted in Fig. 8 because of a minimal overlap between the
condensate and the beam. In all three figures, namely Figs. 6, 7 and 8, it can
be seen that the numerically calculated values corroborated very well with
that of the TF approximation.

###  3.3 Expectation value of angular momentum

To ascertain the birth of new vortices as a function of waist and potential
depth of the laser beam, we calculate the expectation values of angular
momentum as a function of $V_{0}$ from the numerical solution. Fig. 9, depicts
the variation of $\langle L_{z}\rangle$ with respect to $V_{0}$ for the two
values of $w_{0}$, namely, $w_{0}=20.2$ and $0.47$.

Figure 9: Plot of the expectation value of angular momentum ($\langle
L_{z}\rangle$) versus the Gaussian potential depth ($V_{0}$) for $w_{0}=20.2$
and $w_{0}=0.47$.

In Fig. 9, the expectation values of angular momentum ($\langle L_{z}\rangle$)
seem to increase monotonically for $w_{0}=20.2$, which signifies the birth of
new vortices within the condensate. On the other hand, understandably, for the
thin beam due to the poor overlap of the beam with the condensate, these
$\langle L_{z}\rangle$ values for $w_{0}=0.47$ remains almost constant,
denoting no more new vortices generated within the condensate as shown in Fig.
9.

###  3.4 Free energy

Another significant quantity that demonstrates the creation of new vortices as
the potential depth increases is the free energy $F$, which can be estimated
as [46]

$\displaystyle F=E-\mu-\Omega\langle L_{z}\rangle,$ (17)

where $E$ and $\mu$ correspond to the energy and chemical potential in the
absence of rotation. Formation of vortices can be seen in terms of the free
energy by computing the characteristic vortex nucleation time. Therefore, we
plot a graph between the change in free energy ($\Delta F=F(t)-F(0)$) with
respect to time

Figure 10: The time evolution of the free energy for three different values of
$V_{0}$ plotted as the difference between the instant free energy and the free
energy at $t=0$ with $w_{0}=20.2$

and observe that during the formation of a vortex, the free energy of the
condensate rapidly decreases for some time and it stabilizes once the vortex
reaches the equilibrium state. Also, from Fig. 10, we observe that as the
Gaussian potential depth increases, the characteristic nucleation time
decreases. Moreover, the free energy fall-off is steeper for a higher depth of
the Gaussian potential ($V_{0}$).

###  3.5 Critical Rotation frequency

Even though we keep the rotation frequency constant for our analysis, we are
interested to know whether the considered value of rotational frequency is
fair enough to suit the conditions for vortex formation. In addition to the
above, we intend to study the relationship between the critical rotational
frequency and the Gaussian potential depth for the formation of a single
vortex. We calculate the critical rotational frequencies for different
Gaussian potentials for two different waist values, viz $w_{0}=20.2$ and
$0.47$, and also for the entire range of potential depth of interest.

Figure 11: Plot of critical rotational frequency ($\Omega_{c}$) as a function
of Gaussian potential depth ($V_{0}$) for $w_{0}=20.2$ and $w_{0}=0.47$

In Fig. 11, we show the variation of critical rotational frequency as a
function of $V_{0}$ for $w_{0}=20.2$ and $0.47$, respectively, with
sufficiently long time evolution to ensure the equilibrium. It is clearly seen
from Fig. 11 that our choice of rotational frequency is well above the
critical rotational frequency for the entire range, which bears a downtrend
linear relationship with $V_{0}$ for $w_{0}=20.2$ and remains almost constant
with $V_{0}$ when $w_{0}=0.47$.

##  4 Summary and Conclusion

In this paper, we have studied the creation of vortices by triggering the
rotating Bose-Einstein condensates formed inside a toroidal trap with the trap
parameters such as laser beam waist and Gaussian potential depth. By
numerically solving the time-dependent Gross-Pitaevskii equation in two
dimensions, we have shown that the vortex structure has undergone change when
we vary the Gaussian potential depth. Also, we noticed a considerable increase
in the number of vortices when the waist of the irradiated laser beam is in
good agreement with the root mean square of the condensate. By computing the
root mean square radius of the condensate, we have confirmed the variation in
the number of vortices created as a function of the ratio between the root
mean square radius of the condensate and the laser beam waist. The number of
hidden vortices reaches the maximum value when the above ratio is close to
$0.7$. In addition, we have found the variation in the number of vortices is
abrupt for a deep Gaussian potential. In a nutshell, we conclude that the
larger beam waist and a deep Gaussian potential have generated a lot of
vortices. Also, we calculated the number of vortices using the Feynman rule
with the Thomas Fermi approximation and compared them with the numerical
results. Finally, we have observed that the critical rotation frequency
decreases with an increase in the depth of Gaussian potential. Our results may
throw some light on the applicability of rotating BECs towards the fabrication
of quantum computers. Precisely, our work may be useful for the creation of
phase vortex qubits, as discussed in ref. [47].

### Acknowledgements

The authors thank Dr. R. Kishor Kumar for helpful discussions. The work of
P.M. is supported by CSIR under Grant No. 03(1422)/18/EMR-II, DST-SERB under
Grant No. CRG/2019/004059, DST-FIST under Grant No. SR/FST/PSI-204/2015(C),
MoE RUSA 2.0 (Physical Sciences) and DST-PURSE (Phase-II) Programmes.

## References

  * [1] L. Onsager, Statistical hydrodynamics, Nuovo Cim 6 (S2) (1949) 279–287.  doi:10.1007/bf02780991.
  * [2] K. W. Madison, F. Chevy, W. Wohlleben, J. Dalibard, Vortex formation in a stirred bose-einstein condensate, Phys. Rev. Lett. 84 (5) (2000) 806–809.  doi:10.1103/physrevlett.84.806.
  * [3] J. R. Abo-Shaeer, C. Raman, J. M. Vogels, W. Ketterle, Observation of vortex lattices in bose-einstein condensates, Science 292 (5516) (2001) 476–479.  doi:10.1126/science.1060182.
  * [4] M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S. Hall, C. E. Wieman, E. A. Cornell, Vortices in a bose-einstein condensate, Phys. Rev. Lett. 83 (13) (1999) 2498–2501.  doi:10.1103/physrevlett.83.2498.
  * [5] C. Raman, J. R. Abo-Shaeer, J. M. Vogels, K. Xu, W. Ketterle, Vortex nucleation in a stirred bose-einstein condensate, Phys. Rev. Lett. 87 (21) (nov 2001).  doi:10.1103/physrevlett.87.210402.
  * [6] Z. Dutton, M. Budde, C. Slowe, L. V. Hau, Observation of quantum shock waves created with ultra- compressed slow light pulses in a bose-einstein condensate, Science 293 (5530) (2001) 663–668.  doi:10.1126/science.1062527.
  * [7] B. P. Anderson, P. C. Haljan, C. A. Regal, D. L. Feder, L. A. Collins, C. W. Clark, E. A. Cornell, Watching dark solitons decay into vortex rings in a bose-einstein condensate, Phys. Rev. Lett. 86 (14) (2001) 2926–2929.  doi:10.1103/physrevlett.86.2926.
  * [8] P. Mason, N. G. Berloff, Motion of quantum vortices on inhomogeneous backgrounds, Phys. Rev. A 77 (3) (mar 2008).  doi:10.1103/physreva.77.032107.
  * [9] H. M. Nilsen, G. Baym, C. J. Pethick, Velocity of vortices in inhomogeneous bose-einstein condensates, Proc. Natl. Acad. Sci. 103 (21) (2006) 7978–7981.  doi:10.1073/pnas.0602541103.
  * [10] W. Li, M. Haque, S. Komineas, Vortex dipole in a trapped two-dimensional bose-einstein condensate, Phys. Rev. A 77 (5) (may 2008).  doi:10.1103/physreva.77.053610.
  * [11] L. M. Pismen, Vortices in Nonlinear Fields: From Liquid Crystals to Superfluids, from Non-Equilibrium Patterns to Cosmic Strings, Oxford University Press, 1999.
  * [12] A. L. Fetter, A. A. Svidzinsky, Vortices in a trapped dilute bose-einstein condensate, J. Phys.: Condens. Matter 13 (12) (2001) R135–R194.  doi:10.1088/0953-8984/13/12/201.
  * [13] J. Javanainen, S. M. Paik, S. M. Yoo, Persistent currents in a toroidal trap, Phys. Rev. A 58 (1) (1998) 580–583.  doi:10.1103/physreva.58.580.
  * [14] R. Dubessy, T. Liennard, P. Pedri, H. Perrin, Critical rotation of an annular superfluid bose-einstein condensate, Phys. Rev. A 86 (1) (2012) 011602.  doi:10.1103/physreva.86.011602.
  * [15] L. A. Toikka, J. Hietarinta, K.-A. Suominen, Exact soliton-like solutions of the radial gross-pitaevskii equation, J. Phys. A: Math. Theor. 45 (48) (2012) 485203\.  doi:10.1088/1751-8113/45/48/485203.
  * [16] L. A. Toikka, K.-A. Suominen, Snake instability of ring dark solitons in toroidally trapped bose-einstein condensates, Phys. Rev. A 87 (4) (2013) 043601\.  doi:10.1103/physreva.87.043601.
  * [17] L. A. Toikka, O. Kärki, K.-A. Suominen, Creation and revival of ring dark solitons in an annular bose-einstein condensate, J. Phys. B: At. Mol. Opt. Phys. 47 (2) (2014) 021002.  doi:10.1088/0953-4075/47/2/021002.
  * [18] A. Ramanathan, K. C. Wright, S. R. Muniz, M. Zelan, W. T. Hill, C. J. Lobb, K. Helmerson, W. D. Phillips, G. K. Campbell, Superflow in a toroidal bose-einstein condensate: An atom circuit with a tunable weak link, Phys. Rev. Lett. 106 (13) (2011) 130401.  doi:10.1103/physrevlett.106.130401.
  * [19] K. C. Wright, R. B. Blakestad, C. J. Lobb, W. D. Phillips, G. K. Campbell, Driving phase slips in a superfluid atom circuit with a rotating weak link, Phys. Rev. Lett. 110 (2) (2013) 025302.  doi:10.1103/physrevlett.110.025302.
  * [20] F. Piazza, L. A. Collins, A. Smerzi, Critical velocity for a toroidal bose-einstein condensate flowing through a barrier, J. Phys. B: At. Mol. Opt. Phys. 46 (9) (2013) 095302.  doi:10.1088/0953-4075/46/9/095302.
  * [21] L. A. Toikka, Self-interference of a toroidal bose-einstein condensate, New J. Phys. 16 (4) (2014) 043011.  doi:10.1088/1367-2630/16/4/043011.
  * [22] A. Oliinyk, I. Yatsuta, B. Malomed, A. Yakimenko, Symmetry breaking in interacting ring-shaped superflows of bose-einstein condensates, Symmetry 11 (10) (2019) 1312.  doi:10.3390/sym11101312.
  * [23] B. P. Anderson, K. Dholakia, E. M. Wright, Atomic-phase interference devices based on ring-shaped bose-einstein condensates: Two-ring case, Phys. Rev. A 67 (3) (2003) 033601.  doi:10.1103/physreva.67.033601.
  * [24] S. Moulder, S. Beattie, R. P. Smith, N. Tammuz, Z. Hadzibabic, Quantized supercurrent decay in an annular bose-einstein condensate, Phys. Rev. A 86 (1) (2012) 013629.  doi:10.1103/physreva.86.013629.
  * [25] M. Cozzini, B. Jackson, S. Stringari, Vortex signatures in annular bose-einstein condensates, Phys. Rev. A 73 (1) (2006) 013603.  doi:10.1103/physreva.73.013603.
  * [26] U. R. Fischer, G. Baym, Vortex states of rapidly rotating dilute bose-einstein condensates, Phys. Rev. Lett. 90 (14) (2003) 140402.  doi:10.1103/physrevlett.90.140402.
  * [27] A. D. Jackson, G. M. Kavoulakis, E. Lundh, Phase diagram of a rotating bose-einstein condensate with anharmonic confinement, Phys. Rev. A 69 (5) (2004) 053619.  doi:10.1103/physreva.69.053619.
  * [28] A. Aftalion, I. Danaila, Giant vortices in combined harmonic and quartic traps, Phys. Rev. A 69 (3) (2004) 033608.  doi:10.1103/physreva.69.033608.
  * [29] L. Wen, H. Xiong, B. Wu, Hidden vortices in a Bose-Einstein condensate in a rotating double-well potential, Phys. Rev. A 82 (5) (2010) 053627.  doi:10.1103/PhysRevA.82.053627
  * [30] L. H. Wen, X. B. Luo, Formation and structure of vortex lattices in a rotating double-well Bose-Einstein condensate, Laser Phys. Lett. 9 (8) (2012) 618–624.  doi:10.7452/lapl.201210044
  * [31] M. Ögren, G. Drougakis, G. Vasilakis, W. von Klitzing, G. M. Kavoulakis, Stationary states of bose-einstein condensed atoms rotating in an asymmetric ring potential, J. Phys. B: At. Mol. Opt. Phys. 54 (14) (2021) 145303.  doi:10.1088/1361-6455/ac1647.
  * [32] Q. Wang, W. Zhao, L. Wen, Dynamics of kicked spin-orbit-coupled bose-einstein condensates, Results Phys. 25 (2021) 104317.  doi:10.1016/j.rinp.2021.104317.
  * [33] H. Yang, Q. Wang, L. Wen, Ground states of dipolar spin-orbit-coupled bose-einstein condensates in a toroidal trap, J. Phys. Soc. Jpn. 88 (6) (2019) 064001.  doi:10.7566/jpsj.88.064001.
  * [34] V. Bretin, S. Stock, Y. Seurin, J. Dalibard, Fast rotation of a bose-einstein condensate, Phys. Rev. Lett. 92 (2004) 050403.  doi:10.1103/PhysRevLett.92.050403.
  * [35] S. Stock, V. Bretin, F. Chevy, J. Dalibard, Shape oscillation of a rotating bose-einstein condensate, Europhys. Lett. 65 (5) (2004) 594–600.  doi:10.1209/epl/i2003-10179-x.
  * [36] C. Ryu, M. Andersen, P. Cladé, V. Natarajan, K. Helmerson, W. Phillips, Observation of persistent flow of a bose-einstein condensate in a toroidal trap, Phys. Rev. Lett. 99 (26) (2007) 260401.  doi:10.1103/physrevlett.99.260401.
  * [37] C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis, B. P. Anderson, Spontaneous vortices in the formation of bose-einstein condensates, Nature 455 (7215) (2008) 948–951.  doi:10.1038/nature07334.
  * [38] A. Aftalion, P. Mason, Rotation of a bose-einstein condensate held under a toroidal trap, Phys. Rev. A 81 (2) (2010) 023607.  doi:10.1103/physreva.81.023607.
  * [39] P. Mason, N. G. Berloff, Dynamics of quantum vortices in a toroidal trap, Phys. Rev. A 79 (4) (2009) 043620.  doi:10.1103/physreva.79.043620.
  * [40] A. L. Fetter, Rotating trapped bose-einstein condensates, Rev. Mod. Phys. 81 (2) (2009) 647–691.  doi:10.1103/revmodphys.81.647.
  * [41] S. K. Adhikari, P. Muruganandam, Effect of an impulsive force on vortices in a rotating bose-einstein condensate, Phys. Lett. A 301 (3-4) (2002) 333–339.  doi:10.1016/s0375-9601(02)00980-5.
  * [42] W. Bao, H. Wang, An efficient and spectrally accurate numerical method for computing dynamics of rotating bose-einstein condensates, J. Comput. Phys. 217 (2) (2006) 612–626.  doi:10.1016/j.jcp.2006.01.020.
  * [43] M. Abad, M. Guilleumas, R. Mayol, M. Pi, D. M. Jezek, Energy barriers for vortex nucleation in dipolar condensates, Laser Phys. 20 (5) (2010) 1190–1196.  doi:10.1134/s1054660x1009001x.
  * [44] R. Kishor Kumar, P. Muruganandam, Vortex dynamics of rotating dipolar Bose–Einstein condensates, J. Phys. B: At. Mol. Opt. Phys. 45 (2012) 215301 (9pp).  doi:10.1088/0953-4075/45/21/215301
  * [45] R. Feynman, Chapter II application of quantum mechanics to liquid helium, in: Prog. Low. Temp. Phys., Elsevier, 1955, pp. 17–53.  doi:10.1016/s0079-6417(08)60077-3.
  * [46] A. Kato, Y. Nakano, K. Kasamatsu, T. Matsui, Vortex formation of a bose-einstein condensate in a rotating deep optical lattice, Phys. Rev. A 84 (5) (2011) 053623.  doi:10.1103/physreva.84.053623.
  * [47] K. T Kapale, J. P Dowling, Vortex Phase Qubit: Generating Arbitrary, Counterrotating, Coherent Superpositions in Bose-Einstein Condensates via Optical Angular Momentum Beams, Phys. Rev. Lett. 95 (17) (2005) 173601.  doi:10.1103/PhysRevLett.95.173601

##  Appendix A Expressions for inner and outer radii of an annular condensate

For $w_{0}\gg r$, one can obtain the outer radius of the condensate ($R_{2})$
by the following equation [38]

$\displaystyle\frac{8V_{0}}{15w_{0}^{8}}R_{2}^{10}-\frac{V_{0}}{w_{0}^{6}}R_{2}^{8}+\frac{4V_{0}}{3w_{0}^{4}}R_{2}^{6}+\frac{1}{4}\left(1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}\right)R_{2}^{4}-\frac{g}{\pi}=0$
(18)

and the inner radius of the condensate $(R_{1})$ can be obtained by the
relation

$\displaystyle
R_{1}\sim\frac{\xi}{\sqrt{\left[\left(1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}\right)R_{2}^{2}+\frac{4V_{0}}{w_{0}^{4}}R_{2}^{4}-\frac{8V_{0}}{3w_{0}^{6}}R_{2}^{6}+\frac{4V_{0}}{3w_{0}^{8}}R_{2}^{8}-\frac{8V_{0}}{15w_{0}^{1}0}R_{2}^{10}\right]}}$
(19)

where $\xi$ is the quantum circulation of the condensate which can be obtained
from the following equation

$\displaystyle\frac{g\Omega}{\pi\xi}\sim\frac{1}{2}\left(1-\Omega^{2}-\frac{4V_{0}}{w_{0}^{2}}\right)R_{2}^{2}\left[\ln(R_{2}^{2})-2\right]+\frac{V_{0}}{w_{0}^{4}}R_{2}^{4}\left[2\ln(R_{2}^{2})-3\right]-\frac{4V_{0}}{8w_{0}^{6}}R_{2}^{6}\left[3\ln(R_{2}^{2})-4\right]$
(20)

For $w_{0}\ll r$, one can get the outer radius of the condensate by the
following relation

$\displaystyle
R_{2}\sim\left[\left(\frac{4g}{\pi}+2V_{0}w_{0}^{2}\right)\frac{1}{(1-\Omega^{2})}\right]^{1/4}$
(21)

and hence the inner radius of the condensate can be expressed as a quadratic
equation in $(R_{1}^{2})$ as follows:

$\displaystyle\frac{4V_{0}}{w_{0}^{2}}R_{1}^{4}+[(1-\Omega^{2})R_{2}^{2}-2V_{0}]R_{1}^{2}-\xi^{2}=0$
(22)

with the quantum circulation

$\displaystyle\xi\sim\frac{2g\Omega}{\pi(1-\Omega^{2})R_{2}^{2}[\ln(R_{2}^{2})-2]+8.59V_{0}}.$
(23)