#  Quasi-Reliable Estimates of Effective Sample Size

Youhan Fang
Purdue University
<EMAIL_ADDRESS>Yudong Cao
Harvard University
<EMAIL_ADDRESS>Robert D. Skeel
Purdue University / Arizona State University

###### Abstract

The efficiency of a Markov chain Monte Carlo algorithm might be measured by
the cost of generating one independent sample, or equivalently, the total cost
divided by the effective sample size, defined in terms of the integrated
autocorrelation time. To ensure the reliability of such an estimate, it is
suggested that there be an adequate sampling of state space—to the extent that
this can be determined from the available samples. A possible method for doing
this is derived and evaluated.

##  1 Introduction

Markov chain Monte Carlo (MCMC) methods are very widely used to estimate
expectations of random variables, since they are often the only methods
applicable to complicated and high dimensional distributions. Their accuracy
is low relative to computational effort, so it is particularly important to
have reliable error estimates. Due to the correlated nature of MCMC samples, a
variance estimate for $N$ samples uses the effective sample size
$\mathit{ESS}$ rather than $N$ in the denominator. The effective sample size
for a random variable can be computed from its integrated autocorrelation
time. Proposed in this paper is a more reliable estimate of $\mathit{ESS}$,
particularly for multimodal distributions.

This study is motivated by an interest in designing more efficient MCMC
samplers. For evaluating performance of MCMC methods, there is wide support
among both statisticians, e.g., Charles Geyer (1992, page 11) [4], and
physicists, e.g., Alan Sokal (1997, page 13) [15] for using as a metric the
cost of an independent sample:

$\displaystyle\frac{N}{\mathit{ESS}}\times\mbox{cost per
step}=\tau\times\mbox{cost per step},$

where $\tau$ is the integrated autocorrelation time.

Similarly, for evaluating accuracy of MCMC estimates, it is deemed essential,
e.g., Alan Sokal (page 11) [15] again, to know the autocorrelation time.
Nonetheless, MCMC is often used for problems for which $\mathit{ESS}$ is very
small, perhaps even less than 1. How should such problems be handled? One
option is to avoid claims of sampling and instead designate the simulation as
“exploration”. Another option is to change the problem: For example, use a
simpler model. A coarsened model with error bars is better than a refined
model without error bars (because then the limitations are more explicit).
Alternatively, artificially limit the extent of state space, and adjust
conclusions accordingly.

A computed estimate of the integrated autocorrelation time might be marred by
incomplete sampling—at least in the case of multimodal distributions. Without
additional information, there is no foolproof way of detecting convergence, so
reliability is impractical in general. So the goal is instead quasi-
reliability, which is defined to mean apparent good coverage of state
space—ensuring thorough sampling of that part of state space that has already
been explored, to minimize the risk of missing an opening to an unexplored
part of state space. More concretely, it might mean thorough sampling of those
modes that have already been visited, to minimize the risk of missing an
opening to yet another mode.

This article makes a couple of contributions towards more reliable estimates
of effective sample size:

  1. 1.

A quantity is derived for measuring thoroughness of sampling that is simple to
describe and relatively easy to estimate, namely, the longest autocorrelation
time over all possible functions defined on state space, denoted by
$\tau_{\max}$. This quantity has two advantages over estimates of $\tau$ only
for functions of interest. (i) It is more likely to indicate the existence of
a transition path to a mode that has not been visited. (ii) It is more
reliably estimated with finite sampling than is $\tau$ for just any function:
even though a function of interest might have a small value of $\tau$, a
reliable estimate of its autocorrelation function and integrated
autocorrelation time $\tau$ may well require a number of samples dictated by
$\tau_{\max}$.

  2. 2.

An algorithm is constructed for estimating $\tau_{\max}$, assuming the
availability of a method for estimating the integrated autocorrelation time
$\tau$ for an arbitrary function, and evidence is presented for the utility of
this algorithm.

Section 2 of this article establishes notation by reviewing needed concepts
and methods. It presents the definition and role of the integrated
autocorrelation $\tau$. It defines a forward transfer operator that comprises
transition probabilites of the MCMC process. And it presents the notion of a
lag window for estimating $\tau$ for a finite number of samples $N$.

Section 3 is the heart of this article. It begins with an example illustrating
how it is possible to have an apparently reliable esimate of $\tau$ that is
far too small. It then defines thorough sampling to mean that (with, say, 95%
confidence) the fraction of samples from any given subset of state space
should differ from the actual probability by no more than some specified
tolerance $\mathit{tol}$. It is shown that this is achieved if
$N\geq\tau_{\max}/\mathit{tol}^{2}$. An approximate maximum is obtained by
considering arbitrary linear combinations of a finite set of basis functions.
This results in a generalized eigenvalue problem. The matrices are shown to be
symmetric not only for reversible MCMC methods but also for a class of
“modified” reversible MCMC methods that includes generalized hybrid Monte
Carlo samplers and reversible (inertial) Langevin integrators.

Section 4 tests the methodology on 4 examples. A 1-dimensional Gaussian
distribution is used to confirm the correctness of the theory. A sum of 2
Gaussians is used as a simple example of a multimodal problem, which
demonstrates the advantage of seeking the longest autocorrelation time. A
2-layer single-node neural network is used to capture the flavor of machine
learning applications. It yields a large enough estimate of $\tau_{\max}$ to
serve as a warning that an apparently large sample size is inadequate, unlike
the $\tau$ estimate for the quantity of interest. Finally a logistic
regression example confirms again the wisdom of seeking the longest
autocorrelation time $\tau_{\max}$.

Section 5 is a brief discussion and conclusion.

Initial attempts to estimate the integrated autocorrelation time employed a
small computer program called acor [5]. Due to concerns about its reliability,
Appendix A derives an alternative method for defining the lag window under a
couple of simplyfying assumptions, one of which becomes more credible for
observables having a large value of $\tau$. Experiments indicate slightly more
reliable estimates compared to acor.

##  2 Preliminaries and Background

Given a probability density $\rho_{\mathrm{q}}(\mathbf{q})$,
$\mathbf{q}\in\mathbb{R}^{d}$, known only up to a multiplicative factor, the
aim is to compute observables $\mathbb{E}[u(\mathbf{Q})]$, which are
expectations for specified functions $u(\mathbf{q})$. Note here the use of
upper case for random variables.

Consider a Markov chain
$\mathbf{Q}_{0}\rightarrow\mathbf{Q}_{1}\rightarrow\cdots\rightarrow\mathbf{Q}_{N-1}$
that is ergodic with stationary density $\rho_{\mathrm{q}}(\mathbf{q})$. Also,
assume $\mathbf{Q}_{0}\sim\rho_{\mathrm{q}}(\mathbf{q})$. To estimate an
observable, use

$\displaystyle\mathbb{E}[u(\mathbf{Q})]\approx\overline{U}=\frac{1}{N}\sum_{n=0}^{N-1}u(\mathbf{Q}_{n}),$

but use just one realization
$\mathbf{q}_{0}\rightarrow\mathbf{q}_{1}\rightarrow\cdots\rightarrow\mathbf{q}_{N-1}$,
assuming samples are renumbered to omit those prior to equilibration/burnin.

###  2.1 Variance of the estimated mean

There are a variety of methods for estimating the variance of the estimated
mean. Most appealing [4] are those that expoit the fact that the sample is
generated by a Markov chain. The variance of the estimated mean is given by

$\displaystyle\mathrm{Var}[\overline{U}]=\frac{1}{N}\mathrm{Var}[u(\mathbf{Q})]\left(1+2\sum_{k=1}^{N-1}\left(1-\frac{k}{N}\right)\frac{C(k)}{C(0)}\right)$

where the autocovariances

$\displaystyle
C(k)=\mathbb{E}[(u(\mathbf{Q}_{0})-\mu)(u(\mathbf{Q}_{k})-\mu)],$

with $\mu=\mathbb{E}[u(\mathbf{Q})]$. In the limit $N\rightarrow\infty$,

$\mathrm{Var}[\overline{U}]=\frac{1}{N}\mathrm{Var}[u(\mathbf{Q})]\tau+o(\frac{1}{N})$

where $\tau$ is the integrated autocorrelation time.

$\tau=1+2\sum_{k=1}^{+\infty}\frac{C(k)}{C(0)}.$ (1)

An estimate of covariances is provided by Ref. [11, pp. 323–324]:

$\displaystyle
C_{N}(k)=\frac{1}{N}\sum_{n=0}^{N-k-1}(u(\mathbf{q}_{n})-\overline{u})(u(\mathbf{q}_{n+k})-\overline{u}).$

###  2.2 Forward transfer operator

Some MCMC samplers augment state variables $\mathbf{q}$ with auxiliary
variables $\mathbf{p}$, e.g., momenta, extend the p.d.f. to
$\rho(\mathbf{q},\mathbf{p})$ so that

$\displaystyle\int\rho(\mathbf{q},\mathbf{p})\mathrm{d}\mathbf{p}=\rho_{\mathrm{q}}(\mathbf{q}),$

and make moves in extended state space
$\mathbf{z}_{0}\rightarrow\mathbf{z}_{1}\rightarrow\cdots\rightarrow\mathbf{z}_{N-1}$
where $\mathbf{z}=(\mathbf{q},\mathbf{p})$.

Associated with the MCMC propagator is an operator $\mathcal{F}$, which maps a
relative density $u_{n-1}=\rho_{n-1}/\rho$ for $\mathbf{Z}_{n-1}$ to a
relative density $u_{n}=\rho_{n}/\rho$ for $\mathbf{Z}_{n}$:

$\displaystyle u_{n}=\mathcal{F}u_{n-1}$

where

$\displaystyle\mathcal{F}u_{n-1}(\mathbf{z})=\frac{1}{\rho(\mathbf{z})}\int\rho(\mathbf{z}|\mathbf{z}^{\prime})u_{n-1}(\mathbf{z}^{\prime})\rho(\mathbf{z}^{\prime})\mathrm{d}\mathbf{z}^{\prime}$

and $\rho(\mathbf{z}|\mathbf{z}^{\prime})$ is the transition probability
density for the chain [14]. Note that $u\equiv 1$ is an eigenfunction of
$\mathcal{F}$ for eigenvalue $\lambda_{1}=1$.

The error in the probability density, starting from
$\rho_{0}(\mathbf{q})=\delta(\mathbf{q}-\mathbf{q}_{0})$, is

$\displaystyle\frac{1}{N}\sum_{n=0}^{N-1}\rho_{n}-\rho=\rho\frac{1}{N}(1-\mathcal{F}_{0})^{-1}(1-\mathcal{F}_{0}^{N})\frac{\rho_{0}-\rho}{\rho}$

where $\mathcal{F}_{0}$ is $\mathcal{F}$ with the eigenvalue 1 removed, i.e.,
$\mathcal{F}_{0}u=\mathcal{F}u-\mathbb{E}[u(\mathbf{Q})]$. From this, one sees
that the error is proportional to the reciprocal of the spectral gap
$|1-\lambda_{2}|$, where $\lambda_{2}$ is the nonunit eigenvalue of
$\mathcal{F}$ nearest 1. This, however, is not the relevant quantity—in
general (see 3.1).

###  2.3 Error in estimates of integrated autocorrelation time

The standard deviation of the statistical error in the estimate of $\tau$ from
Eq. (1) grows with the number of terms taken; more specifically, it is
approximately $\sqrt{M/N}\tau$ where $M$ is the number of terms [15, Eqs.
(3.19)]. Therefore, one uses instead

$\displaystyle\tau\approx 1+2\sum_{k=1}^{N-1}w(k)\frac{C_{N}(k)}{C_{N}(0)}$

where $w(k)$ is a decreasing weight function known as a lag window.

The tiny program acor uses a lag window that is $1$ for $k$ from $0$ to $M-1$
and $0$ for the rest where $M$ is the smallest number that exceeds the
estimated $\tau$ by a factor of 10. It requires the number of samples $N$ to
exceed $100\tau$.

##  3 Thorough Sampling

Depending on estimates of the ESS for just the observables of interest can be
deceptive.

Consider a mixture of 2 Gaussians,

$\displaystyle-\log\rho(q_{1},q_{2})=$
$\displaystyle\frac{1}{2}(36(q_{1}+1)^{2}+(q_{2}-3)^{2})$ (2)
$\displaystyle+\frac{1}{2}((q_{1}-2)^{2}+36q_{2}^{2})+\mbox{const},$

whose combined basin has an L shape with a barrier at the corner of the L. For
a sampler, consider the use of Euler-Maruyama Brownian dynamics (w/o
rejections), defined in Sec. 3.3. A realization of a sampling trajectory is
given in Figure 1. For the realization shown there, the trajectory makes the
first transition around $N=1000$.

(a) 5000 samples (b) Trajectories of $q_{1}$ and $q_{2}$

Figure 1:  Example trajectory for the L-shaped mixture of Gaussians
distribution.

Figure 2 shows the estimated integrated autocorrelation time and the effective
sample size as a function of the sample size. To consider the question of
thorough sampling, suppose that $\mathit{ESS}=100$ is regarded adequate for
estimating observables. Just before the $q_{1}$ curve jumps, the ESS is much
greater than 100 indicating thorough sampling. This is misleading since after
the jump the ESS drops dramatically. On the other hand, the $q_{2}$ curve more
reliably detects thorough sampling. Note that if the trajectory starts from
the vertical leg of the ”L”, the behavior of $q_{1}$ and $q_{2}$ are swapped.
This implies that a general observable that is independent of the starting
point and can capture the “hardest” move is needed for better detecting
convergence.

(a) Autocorrelation time (b) Effective sample size

Figure 2:  The autocorrelation time and the effective sample size vs. the
number of samples in the L-shaped mixture of Gaussians problem.

###  3.1 Definition of thorough sampling

Good coverage of state space might be defined as follows: for any subset $A$
of (unextended) state space, an estimate $\overline{1_{A}}$ of
$\mathbb{E}[1_{A}(\mathbf{Q})]=\Pr(\mathbf{Q}\in A)$ satisfies

$\mathrm{Var}[\overline{1_{A}}]\leq\frac{1}{4}\mathit{tol}^{2}.$

This would ensure 95% confidence that the proportion of samples for any subset
$A$ is not off by more than $\mathit{tol}$. Since

$\mathrm{Var}[\overline{1_{A}}]\approx\tau_{A}\frac{1}{N}\mathrm{Var}[1_{A}(\mathbf{Q})]\leq\frac{1}{4N}\tau_{A},$

it is enough to have

$\frac{1}{4N}\tau_{A}\leq\frac{1}{4}\mathit{tol}^{2}\quad\mbox{for all }A.$
(3)

There may be permutations $P$ such that
$\rho_{\mathrm{q}}(P\mathbf{q})=\rho_{\mathrm{q}}(\mathbf{q})$ and
$u(P\mathbf{q})=u(\mathbf{q})$ for all interesting $u$. Then, it is
appropriate to weaken the definition of good coverage by including only those
$A$ for which $1_{A}(P\mathbf{q})=1_{A}(\mathbf{q})$.

To facilitate analysis, introduce the inner product

$\langle
v,u\rangle=\int\overline{v(\mathbf{z})}u(\mathbf{z})\rho(\mathbf{z})\,\mathrm{d}\mathbf{z}.$

Note that $\langle 1,u\rangle$ is the expectation $\mathbb{E}[u(\mathbf{Z})]$
of $u(\mathbf{Z})$. One can show by induction that the cross covariance

$\mathbb{E}[v(\mathbf{Z}_{0})u(\mathbf{Z}_{k})]=\langle\mathcal{F}^{k}v,u\rangle.$

and, in particular, $C(k)=\langle\mathcal{F}^{k}u,u\rangle$.

For simplicity, instead of just indicator functions $1_{A}(\mathbf{q})$,
consider arbitrary functions in

$\displaystyle W=$
$\displaystyle\\{u=u(\mathbf{q})~{}|~{}\mathbb{E}[u(\mathbf{Q})]=0,$
$\displaystyle u(P\mathbf{q})=u(\mathbf{q})\mbox{ for symmetries }P\\},$

and define thorough sampling using

$\tau_{\max}=\max_{u\in
W}\left(1+2\sum_{k=1}^{+\infty}\frac{C(k)}{C(0)}\right)$

in Eq. (3). The maximum autocorrelation time can be rewritten as

$\displaystyle\tau_{\max}=$ $\displaystyle\max_{u\in
W}\left(1+2\sum_{k=1}^{+\infty}\frac{\langle\mathcal{F}^{k}u,u\rangle}{\langle
u,u\rangle}\right)$ $\displaystyle=$ $\displaystyle\max_{u\in W}\frac{\langle
g(\mathcal{F})u,u\rangle}{\langle u,u\rangle}$

where

$g(\lambda)=1+2\lambda/(1-\lambda).$

For reversible samplers, the spectral gap is intimately related to
$\tau_{\max}$. If $\tau_{\max}$ were the maximum over all $u(\mathbf{z})$,
then

$\tau_{\max}=\frac{1+\lambda_{2}}{1-\lambda_{2}}.$

where $\lambda_{2}$ is the second largest eigenvalue of the transfer operator.

###  3.2 Spatial discretization

For computational purposes, the idea is to seek a function that comes close to
maximizing the autocorrelation time by considering a linear combination
$u(\mathbf{q})=\mathbf{a}^{\mathsf{T}}\mathbf{u}(\mathbf{q})$ of given basis
functions $u_{i}\in W$ where the $a_{i}$ are to be determined. (In practice,
only the _estimated_ mean of the function $u(\mathbf{q})$ vanishes.) Then its
autocovariance

$C(k)=\langle\mathcal{F}^{k}u,u\rangle=\mathbf{a}^{\mathsf{T}}C_{k}\mathbf{a}$

where

$C_{k}=\langle\mathcal{F}^{k}\mathbf{u},\mathbf{u}^{\mathsf{T}}\rangle=\mathbb{E}[\mathbf{u}(\mathbf{Q}_{0})\mathbf{u}(\mathbf{Q}_{k})^{\mathsf{T}}]$

and

$\tau_{\max}\approx\max_{\mathbf{a}}\frac{\mathbf{a}^{\mathsf{T}}K\mathbf{a}}{\mathbf{a}^{\mathsf{T}}C_{0}\mathbf{a}}\quad\mbox{where
}K=C_{0}+2\sum_{k=1}^{+\infty}C_{k}$

—a generalized eigenvalue problem (arising from maximizing a Rayleigh
quotient).

Linear functions $u_{i}(\mathbf{q})=q_{i}$ are suggested as a general choice.

###  3.3 Symmetry of cross covariance matrices

The matrix $C_{0}$ is symmetric positive definite, if the components of
$\mathbf{u}$ are linearly independent. The eigenvalue problem is well
conditioned if $K$ is a symmetric matrix. Symmetry of $K$ does not hold for
every sampler. Even when $K$ is symmetric, its finitely sampled estimate is
not, so using the symmetric part can reduce sampling error.

A reversible MCMC samplers is defined to be one that satisfies detailed
balance:

$\rho(\mathbf{z}^{\prime}|\mathbf{z})\,\rho(\mathbf{z})=\rho(\mathbf{z}|\mathbf{z}^{\prime})\,\rho(\mathbf{z}^{\prime})$

Examples of reversible samplers include

  1. 1.

a Brownian sampler, the Euler-Maruyama discretization of Brownian dynamics
coupled with a Metropolis-Rosenbluth-Rosenbluth-Teller-Teller [10] acceptance
test (known as MALA in the statistics literature [13]),

  2. 2.

hybrid Monte Carlo [2],

  3. 3.

generalized hybrid Monte Carlo [6], with each step followed by a momentum
flip,

  4. 4.

a reversible Langevin integrator, with every fixed number of steps followed by
a momentum flip (see below). Note that reversible has a different meaning for
integrators than it does for samplers.

These last two samplers augment state space with momentum variables.

The momentum flip of the last two enumerated samplers defeats the purpose of
augmenting state space with momenta. Indeed, reversing the momentum flip leads
to more rapid mixing. More generally, a modified reversible propagator
effectively couples the reversible substep
$\mathbf{z}_{n-1}\rightarrow\bar{\mathbf{z}}_{n}$ with a substep,

$\mathbf{z}_{n}=R(\bar{\mathbf{z}}_{n}),$

where $R$ satisfies $R\circ R=\mathrm{id}$ and $\rho\circ R=\rho$. For a
momentum flip, $R(\mathbf{q},\mathbf{p})=(\mathbf{q},-\mathbf{p})$.

Following is an example of a reversible Langevin integrator. Let
$F(\mathbf{q})=-\nabla_{\mathbf{q}}\log\rho(\mathbf{q},\mathbf{p})$. Then each
step of the integrator consists of the following five stages:

  * B:

$\mathbf{p}:=\mathbf{p}+\frac{1}{2}\Delta tF(\mathbf{q})$

  * A:

$\mathbf{q}:=\mathbf{q}+\frac{1}{2}\Delta t\mathbf{p}$

  * O:

$\mathbf{p}:=\sqrt{1-2\gamma\Delta t}\mathbf{p}+\sqrt{2\gamma\Delta
t}\mathcal{N}(0,I)$

  * A:

$\mathbf{q}:=\mathbf{q}+\frac{1}{2}\Delta t\mathbf{p}$

  * B:

$\mathbf{p}:=\mathbf{p}+\frac{1}{2}\Delta tF(\mathbf{q})$

The dynamics of this integrator has a precise stationary density. It differs
from the desired density [8] by $\mathcal{O}(\Delta t^{2})$. Such error is
much smaller than statistical error. The Euler-Leimkuhler-Matthews integrator
[7] is the special case $\gamma=1/(2\Delta t)$ of this, and it has the
remarkable property of retaining second order accuracy. Since applying a
momentum flip has no effect in this case, it is a reversible MCMC propagator.
(The Euler-Leimkuhler-Matthews integrator can be expressed as a Markov chain
in $\mathbf{q}$-space if post-processing is employed [16].)

For a modified reversible propagator, the forward transfer operator is a
product

$\mathcal{F}=\mathcal{R}\bar{\mathcal{F}},$

where $\bar{\mathcal{F}}$ is the forward transfer operator of the reversible
substep and $\mathcal{R}$ is the operator for the added substep
$\mathbf{z}_{n}=R(\bar{\mathbf{z}}_{n})$:

$u_{n}=\mathcal{R}\bar{u}_{n-1}$

where

$\mathcal{R}\bar{u}_{n-1}(\mathbf{z})=\frac{1}{\rho(\mathbf{z})}\int\delta(\mathbf{z}-R(\mathbf{z}^{\prime}))\bar{u}_{n-1}(\mathbf{z}^{\prime})\rho(\mathbf{z}^{\prime})\mathrm{d}\mathbf{z}^{\prime}.$

It can be shown that the operators $\bar{\mathcal{F}}$ and $\mathcal{R}$ are
self adjoint:

$\langle\bar{\mathcal{F}}g,f\rangle=\langle
g,\bar{\mathcal{F}}f\rangle\quad\mbox{and}\quad\langle\mathcal{R}g,f\rangle=\langle
g,\mathcal{R}f\rangle.$

Additionally, assume
$\mathcal{R}\mathbf{u}(\mathbf{q})=\mathbf{u}(\mathbf{q})$. Then it is
straightforward to show that matrices $C_{k}$ are symmetric.

###  3.4 A method for estimating maximum autocorrelation time

Begin with a guess for $\mathbf{a}$, e.g., choose $\mathbf{a}$ to be a vector
of zeros except for a one corresponding to the observable with the the longest
autocorrelation time. The following algorithm might converge:

  1. 1.

With $c(k)=\mathbf{a}^{\mathsf{T}}C(k)\mathbf{a}$, use the method of Sec. A.1
to select $w(k)$, $k=1,2,\ldots,M-1$, where $M$ is heuristically chosen to be
large but not so large to impact computational efficiency,

  2. 2.

Set

$K=C_{0}+2\sum_{k=1}^{M-1}w(k)C_{k}.$

  3. 3.

Choose $\mathbf{a}$ to maximize
$\mathbf{a}^{\mathsf{T}}K\mathbf{a}/(\mathbf{a}^{\mathsf{T}}C_{0}\mathbf{a})$.

  4. 4.

Repeat.

The number of samples $N$ needed for an estimate of $\tau_{\max}$ depends on
$\tau_{\max}$ itself; for this one can use a formula suggested by others,
e.g., acor requires $N\geq 100\tau$, roughly.

To improve the convergence of the algorithm, one might prohibit $\tau_{\max}$
from decreasing from one iteration to the next.

It is important to keep the method simple, i.e., the time complexity must not
exceed the complexity for sampling. And keeping this modest is aided by use of
the FFT to estimate autocovariances.

##  4 Numerical Experiments

###  4.1 A Gaussian distribution

To confirm the correctness of the theory, consider a simple test problem for
sampling where the target distribution is the standard Gaussian. The sampler
is Brownian dynamics with a discrete time step of length $\Delta t=0.02$. This
step size is much smaller than needed in practice. It is chosen in this way to
make the discretization error negligible, thereby permitting the use of
analytical results derived for exact Brownian dynamics. The Markov chain is
obtained by subsampling the original trajectory at intervals of $0.1$ (every
5th point). This gives the true eigenvalues of the transfer operator [12] as
$1,\exp(-0.1),\exp(-0.2),\exp(-0.3),\ldots$. The corresponding eigenfunctions
are the Hermite polynomials $H_{i}(q)=2qH_{i-1}(q)-H^{\prime}_{i-1}(q)$ with
$H_{0}(q)=1$.

Consider three observables

$H_{3}+H_{2}+H_{1},\quad H_{3}-H_{2}+H_{1},\quad-H_{3}+H_{2}+H_{1}.$

The goal of this test is to recover the theoretical value

$\tau_{\max}=(1+\exp(-0.1))/(1-\exp(-0.1))\approx 20.0,$

and its corresponding observable, which is $H_{1}(q)$, by using the linear
combination of given observables. Twelve independent runs are performed to
evaluate the reliability of the method.

Figure 3 shows the estimated $\tau_{\max}$ for all 12 runs for increasing
sample size. It can be seen that the estimated value converges to the true
value. And the variance is comparatively small when $N>3\times 10^{3}.$ Table
1 shows the coefficients of the linear combinations of given observables. It
can be seen that the theoretical maximizing observable $H_{1}(q)$ is
successfully recovered.

Figure 3: The estimated maximum autocorrelation time in the 1D Gaussian problem. Table 1: The weights of the linear combination of the three observables with $a_{1}$ normalized to 1. | $N=10^{3}$ | $N=10^{4}$ | $N=10^{5}$ | $N=10^{6}$
---|---|---|---|---
$a_{1}$ | -0.036 | -0.008 | 0.003 | 0.010
$a_{2}$ | 0.972 | 1.024 | 1.008 | 0.994
$a_{3}$ | 1.000 | 1.000 | 1.000 | 1.000

###  4.2 An L-shaped mixture of two Gaussian distributions

Consider again the mixture of two Gaussians with target density given by eqn.
(2) sampled with discrete Brownian dynamics with time step size $\Delta
t=0.02$ and subsampled at time interval $0.1$ (every 5 points). The given
observables are $q_{1}$ and $q_{2}$. Theoretically, $q_{1}$ should have larger
autocorrelation time than $q_{2}$ for small $N$ while the chain remains in the
horizontal leg of the ”L”, since the frequency is smaller in the direction of
$q_{1}$ than it is in the direction of $q_{2}$. By the same reasoning, the
weight for $q_{1}$ in the linear combination forming the maximizing observable
should be greater than the weight for $q_{2}$, when $N$ is small. Due to
symmetry, on the other hand, both auto-correlation times and the magnitudes of
the weight for $q_{1}$ and $q_{2}$ should be equal when $N$ large.

Figure 4 shows ESS and $\tau$ of single observables and the maximizing
observable. It can be seen that $\tau_{\max}$ is always greater than the
larger $\tau$ of individual observables. As is discussed in previous section,
this tells us that obtaining $\tau_{\max}$ is helpful for better convergence
detection. Table 2 shows the weights of linear combinations for different
sample sizes. The weighting obtained meets expectations, ultimately giving
equal weight to $q_{1}$ and $q_{2}$.

(a) Auto-correlation time (b) Effective sample size

Figure 4:  Effective sample size and autocorrelation time vs. number of samples in the L-shape mixture of Gaussians.  Table 2: The weights of the linear combination of the $q_{1}$ and $q_{2}$ with $a_{1}$ normalized to $1$. | $N=10^{3}$ | $N=10^{4}$ | $N=10^{5}$ | $N=10^{6}$
---|---|---|---|---
$a_{1}$ | 1.000 | 1.000 | 1.000 | 1.000
$a_{2}$ | 0.089 | 0.085 | -1.207 | -1.051

###  4.3 A one-node neural network

Consider the simplest Bayesian neural network regression model [1], having a
single node in the hidden layer:

$y\approx u(\mathbf{q};x)=q_{3}\tanh(q_{1}x+q_{2})+q_{4},$

where $(x,y)$ represents a data point. Suppose a total of $100$ data points
are given as shown by the large dots of Figure 5 (a). The posterior
distribution is

$-\log\rho(\mathbf{q})=\frac{1}{2}\beta\sum_{i=1}^{100}(y_{i}-u(\mathbf{q};x_{i}))^{2}+\frac{1}{2}\alpha\|\mathbf{q}\|^{2}+\mbox{const},$

where $\beta=2.5$ and $\alpha=0.8$. Suppose the sampler is discretized
Brownian dynamics with $\Delta t=0.01$, and let the Markov chain be obtained
by subsampling the original trajectory at time interval $0.1$ (every 10
points). The given observables for finding the maximizing observable are the
parameters of the model $q_{1}$, $q_{2}$, $q_{3}$ and $q_{4}$. The prediction
for the first data point

$\bar{u}(x_{1}):=\frac{1}{N}\sum_{i=1}^{N}u(\mathbf{q}_{i};x_{1})$

is also examined as an example of an observable of interest.

Figure 5 (a) shows the overall prediction $\bar{u}(x)$ resulting from the
regression model, and Figures 5 (b) and (c) show the trajectories of
observables. It can be seen that when the sampling is insufficient, the
prediction is inaccurate (the asymmetry implies only one mode is explored). It
can also be seen that the maximizing observable makes many fewer transitions
than $q_{1}$. This is because the maximizing observable tends to make the
“hardest” moves, hence has longest autocorrelation time.

Figure 6 (a) and (b) show the ESS and $\tau$ for single observables and the
maximizing observable. It can be seen that the ESS of the prediction
$\bar{u}(x_{1})$ is obviously misleading for convergence detection, since it
indicates sufficient sampling when $N\approx 100$. It can also be seen that
the maximum $\tau$ is significantly larger than the longest $\tau$ of
individual observables. This is a result of combining 4 observables linearly.

Figure 6 (c) shows the mean squared error of the regression converging when
$N>2.0\times 10^{4}$. This coincides with the convergence of $\tau_{\max}$. On
the other hand, the the $\tau$ for the prediction might suggest stopping far
too early and the $\tau$ for $q_{1}$ might also be too optimistic.

(a) The prediction (b) Trajectory of $q_{1}$ (c) Trajectory of $u_{\max}$

Figure 5: The resulting prediction and the trajectories of observables in the
one-node neural network problem.

(a) $\tau$ (b) ESS (c) MSE

Figure 6: Autocorrelation times, effective sample sizes, and the mean squared
error for the one-node neural network model.

###  4.4 Logistic regression

Consider a Bayesian logistic regression model [3]. The logistic function maps
a linear combination of features $\mathbf{x}$ to a probability:

$\sigma(\mathbf{q};\mathbf{x})=1/(1+\exp(-\mathbf{q}^{\mathsf{T}}\mathbf{x})).$

The posterior distribution is

$\displaystyle-\log\rho(\mathbf{q})=$
$\displaystyle\beta\sum_{i=1}^{n}(y_{i}\log(\sigma_{i})+(1-y_{i})\log(1-\sigma_{i}))$
$\displaystyle+\frac{1}{2}\alpha\|\mathbf{q}\|^{2}+\mbox{const.},$

where $y_{i}$ is the class label of data example $i$, $\beta=1.0$,
$\alpha=0.1$, and $n$ is the total number of data examples. When predicting
the label of a given data example $\mathbf{x}$, use the average value

$\bar{\sigma}(\mathbf{x})=\frac{1}{N}\sum_{i=1}^{N}\sigma(\mathbf{q}_{i};\mathbf{x})$

over all samples. The label $t$ of a data example $\mathbf{x}$ is predicted to
be $1$ if $\bar{\sigma}(\mathbf{x})\geq 0.5$ and $0$ if
$\bar{\sigma}(\mathbf{x})<0.5$. For data, use the Australian Credit Approval
dataset from the UCI machine learning data repository [9]. It contains 690
data examples and provides 15 features for each data example. Half of the
examples in the dataset are extracted to form the training set, and the rest
of them are in the testing set. The sampler is Brownian dynamics with $\Delta
t=0.05$. This step size best balances computation cost and sampling accuracy
(measured by the training error). As for the one-node neural network problem,
the model parameters $q_{i}$ are used as observables for estimating
$\tau_{\max}$ and the prediction for the first data point in the training set
as an observable of interest.

Figure 7 shows $\tau$, the ESS, and the training/testing error in the logistic
regression problem. The training error is defined as

$err_{\mathrm{train}}=\frac{1}{n}\sum_{i=1}^{n}1_{t(\mathbf{x}_{i})\neq
y_{i}}$

The testing error is defined similarly except the data examples are from the
testing set.

In practice, one may monitor convergence by observing the convergence of the
training error. This approach, however, can be very expensive for large data
sets. One has to evaluate the prediction function for all data examples.
Instead, one might observe the convergence of observables such as the
prediction function of a single data example or a single weight. But as is
shown in the figure, it is risky to only observe one of them. Both the
training and testing error do not converge until the number of samples exceeds
$1000$, but both the prediction and the weight$q_{1}$ have an almost constant
autocorrelation time equal to 1 (implying independent samples). This is
obviously misleading. On the other hand, the maximizing observable does not
show convergence until $N>1000$. This demonstrates that the maximzing
observable is a more reliable indicator of convergence.

(a) $\tau$ (b) ESS (c) Training error

Figure 7: The autocorrelation times, the effective sample sizes, and the
training/testing error in the logistic regression problem.

##  5 Discussion and conclusions

There is considerable agreement on the importance of estimating the integrated
autocorrelation time $\tau$. It is illustrated here and acknowledged by others
that better estimates of $\tau$ are needed. Derived and presented here are
computational methods based on estimating the greatest possible $\tau$, which
for some problems give more reliable estimates at a cost that is only a small
fraction of that required to generate the samples. The approach suggested here
is quite general, and there are opportunities for refinement, for example,
using better basis functions.

## References

  * [1] C. M. Bishop.  Neural Network for Pattern Recognition.  Clarendon Press, 1996.
  * [2] S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth.  Hybrid Monte Carlo.  Phys. Rev. B, 195:216–222, 1987.
  * [3] A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin.  Bayesian Data Analysis.  New York: Chapman and Hall, 2004.
  * [4] C. J. Geyer.  Practical Markov chain Monte Carlo.  Stat. Sci., 7:473–511, 1992.
  * [5] J. Goodman.  Acor, statistical analysis of a time series, Spring 2009.
  * [6] A. M. Horowitz.  A generalized guided Monte-Carlo algorithm.  Phys. Lett. B, 268:247–252, 1991.
  * [7] B. Leimkuhler and C. Matthews.  Rational construction of stochastic numerical methods for molecular sampling.  Appl. Math. Res. Express, pages 34–56, 2013.
  * [8] B. Leimkuhler, C. Matthews, and G. Stoltz.  The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics.  IMA J. Numer. Anal., 36:13–79, 2016.
  * [9] M. Lichman.  UCI machine learning repository, 2013.
  * [10] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller.  Equation of state calculations by fast computing machines.  J. Chem. Phys., 21:1087–1092, 1953.
  * [11] M. B. Priestly.  Spectral Analysis and Time Series.  Academic Press, 1981.
  * [12] H. Risken.  The Fokker-Planck Equation.  Springer, 1996.
  * [13] G. O. Roberts and R. L. Tweedie.  Exponential convergence of Langevin diffusions and their discrete approximations.  Bernoulli, 2:341–363, 1996.
  * [14] C. Schütte and W. Huisinga.  Biomolecular conformations as metastable sets of Markov chain.  Proceedings of the 38th Annual Allerton Conference on Communication, Control and Computing,, pages 1106–1115, 2000.
  * [15] A. Sokal.  Monte Carlo methods in statistical mechanics: Foundations and new algorithms.  In C. DeWitt-Morette, P. Cartier, and A. Folacci, editors, Functional Integration: Basics and Applications, pages 131–192. Springer US, Boston, MA, 1997.
  * [16] G. Vilmart.  Postprocessed integrators for the high order integration of ergodic SDEs.  SIAM J. Sci. Comput., 37:A201–A220, 2015.

##  Appendix A Marginally Better Estimates of Integrated Autocorrelation Time

###  A.1 An optimal lag window

The autocorrelation function can be highly oscillatory, which would invalidate
the method proposed here. To obtain a much more slowly varying function, use
values based on doubling: $v_{i}=u(\mathbf{q}_{2i})+u(\mathbf{q}_{2i+1})$. In
the case of a reversible sampler, this can be shown to yield a positive
decreasing convex autocorrelation function [4, Sec. 3.3]. It is
straightforward to show that the autocorrelation time of $u(\mathbf{q})$ can
be recovered from that of $v$:

$\tau_{u}=\frac{1}{2}\frac{C_{v}(0)}{C_{u}(0)}\tau_{v}.$

####  A.1.1 Form of lag window

Assume the use of a lag window with $0\leq w(k)\leq 1$ and $w(k)$
nonincreasing. Assume that covariance estimates $C_{N}(k)$ satisfy

$C_{N}(k)=C(k)+\sigma C(0)\eta_{k}$

for some nonnegative value $\sigma$ where the $\eta_{k}$ are independent
standard normally distributed random values. The estimated autocorrelation
time is

$\tau^{\mathrm{est}}=1+2\sum_{k=1}^{+\infty}w(k)\frac{c(k)+\sigma\eta_{k}}{1+\sigma\eta_{0}}$

where $c(k)=C(k)/C(0)$ is the ACF and $w(k)=0$ for $k\geq M$. Assuming
$\sigma\ll 1$, the error in the integrated autocorrelation time, to a first
approximation, is

$-2R+2\sigma\sum_{k=1}^{+\infty}w(k)\eta_{k}-2\sigma\eta_{0}\sum_{k=1}^{+\infty}w(k)c(k),$

where

$R=\sum_{k=1}^{+\infty}(1-w(k))c(k).$

The expectation of the square of the error is

$4R^{2}+4\sigma^{2}\sum_{k=1}^{+\infty}w(k)^{2}+4\sigma^{2}\left(\sum_{k=1}^{+\infty}w(k)c(k)\right)^{2},$

which is minimized for a given value of $m$ if

$w(k)=\frac{c(k)}{\sigma^{2}}\left(R-\sigma^{2}\sum_{k=1}^{+\infty}w(k)c(k)\right).$
(4)

Note that, for small enough $\sigma$, the right-hand side is positive.

The ACF $c(k)$ is unknown, so, for the purpose of defining the lag window, it
is suggested to use the specific model

$C_{N}(k)=c_{0}\lambda^{k}+\sigma c_{0}\eta_{k},$

where $\lambda$, $c_{0}$, and $\sigma$ are determined as in Section A.1.2 that
follows. The goal is to capture the behavior of the tail of the ACF, which
will be achieved by using enough points $M$ in the fitting. Then
$c(k)=\lambda^{k}$, so $w(k)$ is proportional to $\lambda^{k}$, suggesting the
choice

$w(k)=\min\\{1,\lambda^{k-m}\\}$ (5)

where $m$ is to be determined optimally.

Using the specific model for $C_{N}(k)$, the expectation of the error squared
becomes

$\frac{4c_{0}\lambda^{2}}{(1-\lambda^{2})^{2}}\left(\mu^{2}+\sigma^{2}(1+\lambda-\mu)^{2}\right)+\frac{4\sigma^{2}\log\mu}{\log\lambda}-\frac{4\sigma^{2}}{1-\lambda^{2}}$

where $\mu=\lambda^{m}$. Differentiating with respect to $\mu$ and setting the
result to zero leads to a quadratic equation for $\mu$, which can be shown to
have a positive and negative root. If the positive root exceeds 1, this
indicates that the sample size $N$ is not large enough.

####  A.1.2 Fitting the autocorrelation function

With residuals

$r_{k}=c_{0}\lambda^{k}-C_{N}(k),$

the parameters $c_{0}$ and $0<\lambda<1$ minimize $\sum_{k=0}^{M-1}r_{k}^{2}$.
The value of $c_{0}$ that minimizes this must satisfy

$\sum_{k=0}^{M-1}\lambda^{k}r_{k}=0,$

whence $c_{0}={Q(\lambda)}/{P(\lambda)}$, where

$P(\lambda)=\sum_{k=0}^{M-1}\lambda^{2k}\mbox{ and
}Q(\lambda)=\sum_{k=0}^{M-1}C_{N}(k)\lambda^{k}.$

Substituting this into the objective function gives

$\sum_{k=0}^{M-1}r_{k}^{2}=-\frac{Q(\lambda)^{2}}{P(\lambda)}+\sum_{k=0}^{M-1}c(k)^{2}.$
(6)

The optimal value of $\lambda$ is obtained by minimizing this rational
function. (One can use golden section search beginning with the interval
$[0,1]$. We used bisection applied to the gradient.) The ML estimate of the
standard deviation $\sigma$ of the noise is obtained from the RMS of the
resulting residual.

One can simply use the exponential model for the tail
$\sum_{k=M}^{N-1}w(k)C_{N}(k)$ of $\tau^{\mathrm{est}}$.

###  A.2 A comparison with the lag window of acor

Figure 8 shows the estimated $\tau$’s of $H_{1}$ and $H_{1}+H_{2}+H_{3}$ as
well as $\tau_{\max}$. It can be seen that when $N$ is small, both methods
underestimate $\tau$, but the new one gives larger estimates. When $N$ is
large, both methods converge to the true value, and the new method has a
smaller variance. Overall, the proposed lag window is more reliable than acor.

(a) $H_{1}$ (b) $H_{1}+H_{2}+H_{3}$ (c) $\tau_{\max}$

Figure 8: Estimated $\tau$ and $\tau_{\max}$ by the new lag window and acor
for the 1D standard Gaussian.

###  A.3 The lag window for extremely small ESS

It might be difficult to get a good lag window when the number of samples is
extremely insufficient. Figure 9 shows the lag window obtained by the new
method and acor in the one-node neural network problem. For acor, the lag
window obviously contains not enough data. And for the proposed method, it
contains perhaps too much data since some of the data is already in the
negative value range due to statistical error. The noise for the auto-
correlation function estimator is highly correlated when $\tau$ is large. In
this example, the effective sample size is less than 10. In such cases, $\tau$
may be underestimated by both methods.

Figure 9: An example of ACF and lag window for the one-node neural network
problem.


# Self-accelerating solitons

Boris A. Malomed1,2 1Department of Physical Electronics, School of Electrical
Engineering, Faculty of Engineering, and Center for Light-Matter Interaction,
Tel Aviv University, Tel Aviv 69978, Israel
2Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 7D, Arica,
Chile

###### Abstract

Basic models which give rise to one- and two-dimensional (1D and 2D) solitons,
such as the Gross-Pitaevskii (GP) equations for BEC, feature the Galilean
invariance, which makes it possible to generate families of moving solitons
from quiescent ones. A challenging problem is to find models admitting stable
self-accelerating (SA) motion of solitons. SA modes are known in linear
systems in the form of Airy waves, but they are poorly localized states. This
brief review presents two-component BEC models which make it possible to
predict SA solitons. In one system, a pair of interacting 1D solitons with
opposite signs of the effective mass is created in a binary BEC trapped in an
optical-lattice potential. In that case, opposite interaction forces, acting
on the solitons with positive and negative masses, produce equal
accelerations, while the total momentum is conserved. The second model is
based on a system of GP equations for two atomic components, which are
resonantly coupled by a microwave field. The latter model produces an exact
transformation to an accelerating references frame, thus predicting 1D and 2D
stable SA solitons, including vortex rings.

Introduction. – A basic property of one- and two-dimensional (1D and 2D)
equations which produce solitons is the Galilean invariance, which generates
solitons moving with an arbitrary velocity from a quiescent one. A challenging
issue is to construct localized states moving at a constant acceleration,
rather than constant velocity Belanger ; Parker . A well-known fact is that
the linear Schrödinger equation, $iu_{z}+(1/2)u_{xx}=0$ (it is written in the
form of the paraxial propagation equation in optics, with propagation distance
$z$ and transverse coordinate $x$) admits self-accelerating (SA) solutions in
the form of Airy waves Berry . Later, it was predicted Siviloglu-
Christodoulides  and experimentally demonstrated, originally in optics
Siviloglu et al 2007 , and then in electron beams electron waves , plasmonics
Airy plasmonics , Bose-Einstein condensates (BECs) kli , acoustics acc , gas
discharge gas-discharge , and water waves Airy water  that truncated Airy
waves (TAWs) can be created in these media, the respective exact solution of
the linear Schrödinger equation being

$\displaystyle u_{\mathrm{TAW}}\left(x,z\right)=u_{0}\mathrm{Ai}\left(\alpha
x-\frac{\alpha^{4}}{4}z^{2}+i\aleph\alpha^{2}z\right)\exp\left(\aleph\alpha
x\right)$
$\displaystyle\times\exp\left(-\frac{i\alpha^{6}}{12}z^{3}+\frac{i\alpha^{3}}{2}xz-\frac{\aleph\alpha^{4}}{2}z^{2}+\frac{i\aleph^{2}\alpha^{2}}{2}z\right),$
(1)

where $\mathrm{Ai}$ is the Airy function, and constants $u_{0}$, $\alpha$, and
$\aleph>0$ define, respectively, the amplitude, internal scale, and truncation
of the Airy wave. Accordingly, the input which generates solution (1) is
$u\left(x;z=0\right)=u_{0}\mathrm{Ai}\left(\alpha
x\right)\exp\left(\aleph\alpha x\right)$. The truncation factor is necessary
to make the total norm (power, in terms of optics) of the input finite:

$N_{\mathrm{TAW}}\equiv\int_{-\infty}^{+\infty}|u(x)|^{2}dx=u_{0}^{2}\left(\sqrt{8\pi\aleph}\alpha\right)^{-1}\exp\left(2\aleph^{3}/3\right),$
(2)

while the TAW’s momentum,

$P=i\int_{-\infty}^{+\infty}u_{x}^{\ast}udx,$ (3)

is zero, in spite of the self-acceleration featured by this wave (both $N$ and
$P$ are dynamical invariants of the Schrödinger equation).

In fact, the truncation leads to gradual degradation of the TAW, as shown, in
particular, by factor $\exp\left(-\aleph\alpha^{4}z^{2}/2\right)$ in solution
(1). Another source of degradation is the action of nonlinearity, as TAW is
the eigenmode of the linear medium. Effects of nonlinearity on the Airy waves
were considered in many works Ellenbogen2 -Thawatchai , chiefly demonstrating
decay into solitons.

A possibility to create well-localized (unlike the Airy waves) SA two-
component pulses in optics was elaborated in terms of a system of coupled
nonlinear Schrödinger (NLS) equations with opposite signs of the group-
velocity dispersion (GVD) in them Peschel theory , aiming to create two pulse
components with _opposite signs of their effective masses_. In this case, the
opposite interaction forces with which the coupled components act on each
other give rise to identical signs of the acceleration. A solution for such an
SA bound state was constructed approximately, taking an unperturbed NLS
soliton in the anomalous-GVD component, and applying the Thomas-Fermi
approximation to the normal-GVD one. The so predicted optical SA pulses were
demonstrated experimentally in a roughly similar temporal-domain form, using a
pair of fiber loops with different lengths, coupled to each other at one point
Peschel-experiment .

It is relevant to mention that a soliton moving with a constant acceleration
(although not of the SA type) can be readily predicted taking the single NLS
equation with the self-attractive cubic term and an attractive or repulsive
local defect which moves with acceleration $a$:

$iu_{z}+(1/2)u_{xx}+|u|^{2}u=-\varepsilon\delta\left(x-az^{2}/2\right)u,$ (4)

where $\delta$ is the delta-function, $\varepsilon>0$ or $\varepsilon<0$
corresponding to the attractive or attractive defect, respectively. In terms
of the spatial-domain NLS equation (4) in optics, the “accelerating” defect
represents a narrow parabolic ($x=az^{2}/2$) stripe in the $\left(x,z\right)$
plane, with the locally increased ($\varepsilon>0$) or decreased
($\varepsilon<0$) value of the refractive index. It is convenient to rewrite
Eq. (4) in the co-moving reference frame, applying the corresponding boost
transformation Belanger :

$\displaystyle\xi\equiv x-az^{2}/2,u(x,z)\equiv
v\left(\xi,z\right)\exp\left(ia\xi z+ia^{2}z^{3}/6\right),$ (5) $\displaystyle
iv_{z}+(1/2)v_{\xi\xi}+|v|^{2}v=a\xi v-\varepsilon\delta(\xi)v.$ (6)

In the moving frame, the soliton may be approximated by the simple stationary
solution to Eq. (6),

$v_{\mathrm{sol}}(\xi,z)=(N/2)\mathrm{sech}\left((N/2)\left(\xi-\xi_{0}\right)\right)\exp\left(iN^{2}z/8\right),$
(7)

where $N$ is the soliton’s norm (see Eq. (2)), and $\xi_{0}$ is a shift of the
soliton from the position of the defect. Treating the terms on the right-hand
side of Eq. (6) by perturbations pert , the soliton is considered as a quasi-
particle under the action of an effective potential,

$U\left(\xi_{0}\right)=Na\xi_{0}-\left(\varepsilon
N^{2}/4\right)\mathrm{sech}^{2}\left(N\xi_{0}/2\right).$ (8)

This potential has an equilibrium position with sign
$\mathrm{sgn}\left(\xi_{0}\right)=-\mathrm{sgn}\left(a\varepsilon\right)$,
which exists, for given strength $\varepsilon$ of the defect, if the
acceleration does not exceed a critical value,

$\left|a_{\max}\right|=|\varepsilon|N^{2}/\left(6\sqrt{3}\right).$ (9)

The corresponding largest equilibrium value of the shift is
$\left(\left|\xi_{0}\right|\right)_{\max}=N^{-1}\ln\left(\left(\sqrt{3}+1\right)^{2}/2\right)\approx
1.32/N$. Note that the largest acceleration with which the soliton can be
dragged by the moving defect, as given by Eq. (9), does not depend on the sign
of $\varepsilon$.

A different result is produced for dragging solitons by the local defect
moving with constant acceleration in the framework of the NLS equation with
the quintic, rather than cubic, nonlinearity (which may also be realized in
optical media Cid ),

$iu_{z}+(1/2)u_{xx}+|u|^{4}u=-\varepsilon\delta\left(x-az^{2}/2\right)u.$ (10)

In the uniform space ($\varepsilon=0$), Eq. (10) gives rise to commonly known
1D Townes solitons Salerno ,
$u=\left(3k\right)^{1/4}\sqrt{\mathrm{sech}\left(2\sqrt{2k}x\right)}\exp\left(ikz\right)$,
with arbitrary propagation constant $k>0$. This soliton family is degenerate,
as its norm takes a single value, which does not depend on $k$,
$N=\sqrt{3/2}\pi/2$, and the family is completely unstable against the onset
of the critical collapse Sulem-Sulem . However, all the solitons are
_stabilized_ by the interaction with the quiescent attractive defect ($a=0$,
$\varepsilon>0$) Wang . Then, the above consideration can be developed for the
solitons of Eq. (10) pulled by the defect with constant acceleration. In
particular, the largest acceleration which can be supported by the defect with
given $\varepsilon>0$ is $\left|a_{\max}\right|=(4/\pi)\varepsilon k$, cf. Eq.
(9).

The objective of this perspective is to produce a brief summary of results
which predict possibilities of true SA motion of 1D and 2D solitons in
specific BEC models, one based on the spatially-periodic optical-lattice (OL)
potential, and another one making use of a binary BEC whose components are
resonantly coupled by a microwave (MW) field.

Co-accelerating bound states of solitons with positive and negative masses. –
A two-component model which allows one to predict stable SA bound states of
solitons with positive and negative effective masses is represented by a
system of Gross-Pitaevskii (GP) equations for wave functions $\phi$ and $\psi$
of the binary BEC, including the spatially-periodic potential of the OL type
in each equation, with strengths $U_{1}$ and $U_{2}$ HS :

$\displaystyle i\phi_{t}$ $\displaystyle=$
$\displaystyle-(1/2)\phi_{xx}-\left[g_{1}|\phi|^{2}+\gamma|\psi|^{2}+U_{1}\cos\left(2\pi
x\right)+fx\right]\phi,$ $\displaystyle i\psi_{t}$ $\displaystyle=$
$\displaystyle-(1/2)\psi_{xx}-\left[\gamma|\phi|^{2}-g_{2}|\psi|^{2}+U_{2}\cos\left(2\pi
x\right)+fx\right]\psi.$ (11)

Here, $g_{1}>0$ and $-g_{2}<0$ are coefficients of the self-interaction of the
components, implying that their signs are made opposite by means of the
Feshbach resonance applied to one of the components FR , and $\gamma>0$ is the
coefficient of the cross-attraction. The OL period in Eq. (11) is set equal to
$1$ by means of rescaling. The system also includes a possibility to consider
the action of gravity, with strength $f$, on both components.

Analytical considerations. To provide opposite signs of the effective mass for
solitons in components $\phi$ and $\psi$, it is natural to consider the case
when quasi-wavenumbers of wave functions $\phi$ and $\psi$ are set to be
close, respectively, to the center and edge of the first OL’s Brillouin zone.
The respective effective masses, calculated by means of the known methods for
the linear GP equation Pu ; we ; gapsol , are, respectively,

$M_{1}=\frac{2\pi^{3}+U_{1}^{2}+\pi^{2}\sqrt{4\pi^{4}+2U_{1}^{2}}}{10\pi^{4}+U_{1}^{2}-3\pi^{2}\sqrt{4\pi^{4}+2U_{1}^{2}}},$
(12) $-M_{2}=U_{2}/\left(U_{2}-2\pi^{2}\right)$ (13)

($M_{2}$ is defined with sign minus, to focus on the relevant case of the
negative mass, which means $M_{2}>0$), and the wave functions themselves are
approximated by

$\phi(x)=\Phi(x)\frac{1+2a\cos\left(2\pi
x\right)}{\sqrt{1+2a^{2}}},\psi(x)=\sqrt{2}\Psi(x)\cos\left(\pi x\right),$
(14)

with $a\equiv\sqrt{\left(\pi^{2}/U_{1}\right)^{2}+1/2}-\pi^{2}/U_{1}$, where
$\Phi(x)$ and $\Psi(x)$ are slowly varying (in comparison with $\cos(2\pi x)$)
envelopes. Actually, the negative mass is a characteristic feature of gap
solitons, generated by the interplay of the OL potential and self-repulsion we
.

The substitution of expressions (14) in Eqs. (11) leads, by means of the
averaging procedure, to equations governing the slow evolution of the envelope
amplitudes, which do not include the OL potential,

$\displaystyle i\frac{\partial\Phi}{\partial t}$ $\displaystyle=$
$\displaystyle-\frac{1}{2M_{1}}\frac{\partial^{2}\Phi}{\partial
x^{2}}-\left(G_{1}|\Phi|^{2}+\Gamma|\Psi|^{2}+fx\right)\Phi,$ (15)
$\displaystyle i\frac{\partial\Psi}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{1}{2M_{2}}\frac{\partial^{2}\Psi}{\partial
x^{2}}-\left(\Gamma|\Phi|^{2}-G_{2}|\Psi|^{2}+fx\right)\Psi,$ (16)

with effective nonlinearity coefficients,

$\displaystyle
G_{1}=g_{1}\frac{1+12a^{2}+6a^{4}}{(1+2a^{2})^{2}},G_{2}=\frac{3}{2}g_{2},$
$\displaystyle\Gamma=\frac{1+2a^{2}+2a}{1+2a^{2}}>0.$ (17)

Sign minus is eliminated in front of the second derivative in Eq. (16)
according to the definition of the respective effective mass in Eq. (13). Note
that Eqs. (15) and (16) conserve the total momentum of wave functions $\Phi$
and $\Psi$ (cf. Eq. (3)), in spite of the possibility to produce SA states, as
shown below.

It is obvious that Eqs. (15) and (16) may indeed feature opposite signs of the
effective masses, if $M_{1}$ and $M_{2}$ are both positive (or both negative),
and opposite signs of the self-interaction in the two components, if $G_{1}$
and $G_{2}$ are both positive (or both negative) too. These sign combinations
open the way to the creation of coupled pairs of solitons with opposite signs
of their dynamical masses.

Using the Lagrangian structure of Eqs. (15) and (16), the variational
approximation (VA) can be applied to the system, using the following ansatz
for two-component solitons:

$\left\\{\Phi,\Psi\right\\}=\sqrt{N_{1,2}}\left(2\alpha_{1,2}/\pi\right)^{1/4}\exp\left[i\varphi_{1,2}-\alpha_{1,2}\left(x-\xi_{1,2}\right)^{2}+ik_{1,2}\left(x-\xi_{1,2}\right)\right],$
(18)

with norms $N_{1,2}$, widths $\alpha_{1,2}^{-1/2}$, central coordinates
$\xi_{1,2}$, momenta $k_{1,2}$, and phases $\varphi_{1,2}$. The VA procedure
progress  leads to equations of motion for the coordinates in the absence of
gravity ($f=0$) HS :

$\displaystyle\frac{d^{2}\xi_{1,2}}{dt^{2}}=\frac{N_{2,1}}{M_{1,2}}\alpha\exp\left[-\frac{2\alpha_{1}\alpha_{2}}{\alpha_{1}+\alpha_{2}}(\xi_{1}-\xi_{2})^{2}\right](\xi_{2}-\xi_{1}),$
(19)
$\displaystyle\alpha\equiv\left(2\Gamma/\sqrt{\pi}\right)\left(2\alpha_{1}\alpha_{2}/(\alpha_{1}+\alpha_{2})\right)^{3/2},\rho\equiv
N_{2}/M_{1}-N_{1}/M_{2}.$ (20)

Further, in the linear approximation with respect to
$\Delta\xi\equiv\xi_{2}-\xi_{1}$ Eq. (19) gives rise to the following
equations for $\Delta\xi$ and the mean coordinate,
$\Xi\equiv\left(\xi_{1}+\xi_{2}\right)/2$:

$\frac{d^{2}\Delta\xi}{dt^{2}}=-\alpha\rho\Delta\xi,\frac{d^{2}\Xi}{dt^{2}}=n\Delta\xi,n\equiv\frac{\alpha}{2}\left(\frac{N_{2}}{M_{1}}+\frac{N_{1}}{M_{2}}\right).~{}$
(21)

In the case of $\rho>0$ (see Eq. (20)), the solution of Eqs. (21) predicts
shuttle motion of the bound state of the solitons with positive and negative
masses, combined with its inner oscillations:

$\Delta\xi(t)=\Delta\xi_{0}\cos\left(\sqrt{\alpha\rho}t\right),~{}\Xi(t)=\frac{n}{\alpha\rho}\Delta\xi_{0}\left[1-\cos\left(\sqrt{\alpha\rho}t\right)\right],$
(22)

where $\Delta\xi_{0}$ is an arbitrary amplitude of the oscillations. Examples
of the shuttle motion produced by numerical simulations of Eqs. (15) and (16)
are displayed below in Figs. 1(b) and 2(a).

The SA motion of the two-soliton bound state is predicted at $\rho=0$, i.e.,
as it follows from Eq. (20), for relation $N_{1}/N_{2}=M_{2}/M_{1}$ between
the norms and effective masses of the two components. This condition provides
cancellation of the force driving the evolution of $\Delta\xi$. In this case,
Eqs. (22) demonstrate that the separation $\Delta\xi_{0}$ between centers of
the components remains constant, and the solitons move with identical
accelerations, so that $\Xi(t)=(n/2)\Delta\xi_{0}t^{2}$. Examples of the SA
motion produced by simulations of Eqs. (15) and (16) are presented below in
Figs. 1(a) and 2(a). Finally, in the case of $\rho<0$ Eqs. (21) predict the
exponential growth of $\Delta\xi(t)$, i.e., splitting of the two-soliton pair,
which is also confirmed by direct numerical solutions.

Further, the VA was extended to include the gravity terms in Eqs. (15) and
(16). A straightforward analysis predicts the shift of the point at which the
SA state appears from $\rho=0$ at $f=0$ to
$\rho=-\left(M_{1}M_{2}\alpha\Delta\xi\right)^{-1}\left(M_{1}+M_{2}\right)\exp\left[2\alpha_{1}\alpha_{2}\left(\alpha_{1}+\alpha_{2}\right)^{-1}\left(\Delta\xi\right)^{2}\right]f$.
Thus, the gravity may be used to control the occurrence of the SA regime. The
appropriate value of $f$ can be tuned, in turn, by varying the angle,
$\Theta$, between the vertical axis and direction of the quasi-1D waveguide
into which the BEC is loaded: $f=f_{\max}\cos\Theta$.

Numerical results. Typical examples of the dynamics of the bound states of
solitons with effective positive and negative masses, produced by simulations
of Eqs. (15) and (16) with parameters and inputs

$M_{1}=1,G_{1}=0.9,\Gamma=0.1,G_{2}=0.1+M_{2}^{-1},$ (23)
$\Phi_{0}=\mathrm{sech~{}}x,\Psi_{0}=\mathrm{sech}\left(x-\Delta\xi\right),\Delta\xi=0.1$
(24)

(i.e., the soliton norms are $N_{1}=N_{2}=2$) are displayed in Fig. 1. In
particular, $M_{2}^{-1}=1$ and $M_{2}^{-1}=0.96$ in Figs. 1(a) and (b)
correspond, respectively to $\rho=0$ and $\rho=0.08$ (see Eq. (20)), which, in
agreement with the VA prediction, produce the SA and shuttle motion,
respectively. For values of $M_{2}^{-1}>1$, i.e., $\rho<0$, the simulations
demonstrate splitting of the two-soliton state, also as predicted by the VA
(not shown here).

Figure 1: The evolution of $\left|\Phi\left(x,t\right)\right|$ and
$\left|\Psi\left(x,t\right)\right|$ (nearly overlapping solid and dashed
lines, respectively) produced by simulations of Eqs. (16) and (15) with input
(24) and parameters (23), where $M_{2}^{-1}=1$ (a) and $M_{2}^{-1}=0.96$ (b)
(as per Ref. HS ).

Additional numerical results are presented in Fig. 2(a), which shows the law
of motion of central coordinates for $M_{2}^{-1}=1$, $0.98$, and $0.92$, i.e.,
$\rho_{1;0.98;0.92}=0;0.04;0.16$, respectively. Note that ratio
$\sqrt{\rho_{0.92}/\rho_{0.98}}=2$ exactly corresponds to the ratio of periods
$T_{0.98}/T_{0.92}=2$, in agreement with Eq. (22). Further, Fig. 2(b) compares
the numerically measured acceleration of the SA pair and its VA-predicted
counterpart, as given by Eq. (19). In the presence of the gravity in Eqs. (16)
and (15), the direct simulations also accurately corroborate the VA
predictions HS .

Figure 2: (a) Trajectories of centers of the $\Phi$ and $\Psi$ components
(solid and dashed lines, respectively), produced by the the same simulations
as in Fig. 1, with $M_{2}^{-1}=1$, $0.92$, and $0.98$, as indicated. (b) The
acceleration of the SA bound state, $a_{0}$, vs. the constant separation,
$\Delta\xi\equiv x_{0}$, between the bound solitons with the positive and
negative masses, in the case of $\rho=0$ (see Eq. (20). Rhombuses: results of
the simulations of Eqs. (16) and (15) with input (24) and parameters (23),
where $M_{2}^{-1}=1$. The dashed curve: the VA prediction produced by Eq. (19)
(as per Ref. HS ).

The existence of robustly moving SA bound states suggest a possibility to
consider collisions between such modes moving with opposite accelerations HS .
A typical example, viz., collision between the same SA state which is
displayed above in Figs. 1(a) and 2(a) and its mirror image is displayed in
Fig. 3(a). It shows that the colliding soliton pairs pass through each other,
causing increase of the separation between the bound solitons in each pair,
from $\Delta\xi=0.10$ to $\Delta\xi\approx 0.17$. This, in turn, leads to the
increase of the co-acceleration, in agreement with Eq. (19) and Fig. 2(b).

Figure 3: (a) The collision of the SA pair, shown in Figs. 1(a) and 2(a), and
its mirror image moving with the opposite self-acceleration. (b) Trajectories
of motion of centers of components $\phi$ and $\psi$ (the continuous and
dashed lines, respectively) produced by simulations of the underlying system
(11), see details in the text (as per Ref. HS ).

As concerns direct simulations, the results presented above were produced by
Eqs. (16) and (15), which were derived from the full system (11) as equations
for the slow evolution of amplitudes $\Phi$ and $\Psi$ in ansatz (14). Direct
simulations of the full equations produce similar results HS . As an example,
Fig. 3(b) displays trajectories of centers of the two components, as obtained
from simulations of Eqs. (11) with $U_{1}=0$, $U_{2}=8$, $f=0$ (no gravity)
and input given by Eq. (14), with $\Phi(x)$ and $\Psi(x)$ replaced by
amplitudes $A=0.150$, $0.094$, $0.075$ and $B=0.150$. As seen in the figure,
the stable SA bound state is generated by $A=0.094$, while the VA predicts,
for the same input, the SA state at $A=0.102$. The corresponding relative
error $\approx$ $0.08$ demonstrates the accuracy of the VA, in comparison with
the simulations of the full underlying system (11). The shuttle motion and
splitting are observed at $A=0.075$ and $0.150$, respectively, in agreement
with the VA. In Ref. Yulik , similar simulations were performed for the system
including, in addition to the OL, the quadratic potential in each component.
It produces opposite forces acting on the positive- and negative-mass
components, leading to their splitting.

If the solitons are considered in a more general form than defined by Eq.
(14), the formation of the SA bound states of solitons with positive and
negative masses is still possible (in particular, because gap solitons are
mobile in the general case we ). However, this possibility was not explored in
detail.

Self-acceleration of vortex rings (VRs) in microwave-coupled binary BEC. – The
existence of stable 2D SA solitons, including VRs, in a binary BEC formed by
two atomic states $\phi_{\uparrow,\downarrow}$, resonantly coupled by the
magnetic component $H$ of the MW field, was predicted in Ref. Qin2 . The
respective system of scaled GP equations for $\phi_{\uparrow,\downarrow}$ and
the Poisson equation for $H$ are Qin1

$\displaystyle i\frac{\partial\phi_{\downarrow}}{\partial t}$ $\displaystyle=$
$\displaystyle\left(-\frac{1}{2}\nabla^{2}+\eta\mathbf{-}\beta_{c}\left|\phi_{\uparrow}\right|^{2}-\beta_{s}\left|\phi_{\downarrow}\right|^{2}\right)\phi_{\downarrow}-\pi
H^{\ast}\phi_{\uparrow},$ $\displaystyle
i\frac{\partial\phi_{\uparrow}}{\partial t}$ $\displaystyle=$
$\displaystyle\left(-\frac{1}{2}\nabla^{2}-\eta\mathbf{-}\beta_{c}\left|\phi_{\downarrow}\right|^{2}-\beta_{s}\left|\phi_{\uparrow}\right|^{2}\right)\phi_{\uparrow}-\pi
H\phi_{\downarrow},$ (25)
$\nabla^{2}H=-\phi_{\downarrow}^{\ast}\phi_{\uparrow},$ (26)

where $\ast$ stands for the complex conjugate, $2\eta$ is detuning of the MW
from the transition between the atomic states $\left|\uparrow\right\rangle$
and $\left|\downarrow\right\rangle$, $\beta_{c}$ and $\beta_{s}$ being,
respectively, strengths of the cross- and self-interaction of the two
components. Equations (25) are supplemented by the normalization condition,
$\int\left(\left|\phi_{\uparrow}\right|^{2}+\left|\phi_{\downarrow}\right|^{2}\right)d\mathbf{r}=1$.
In the symmetric system with $\eta=0$, Eqs. (25) coalesce into a single one
for $\phi_{\downarrow}=\phi_{\uparrow}\equiv\phi$, supplement by the
accordingly simplified Poisson equation for real $H$,

$i\frac{\partial\phi}{\partial
t}=\left[-\frac{1}{2}\nabla^{2}-\beta\left|\phi\right|^{2}-\pi
H\right]\phi,\nabla^{2}H=-|\phi|^{2},$ (27)

with $\beta\equiv\beta_{c}+\beta_{s}$ and normalization
$\int\left|\phi(\mathbf{r})\right|^{2}d\mathbf{r}=1/2$.

Numerical solutions of Eq. (27) demonstrate that it gives rise to VR solutions
with chemical potential $\mu$ and vorticity $S=0,1,2,...$, in the form of
$\phi=\exp\left(-i\mu t+iS\theta\right)\Phi_{S}\left(r\right)$, where $r$ and
$\theta$ are the polar coordinates. The solutions exist in the region of
$\beta<\beta_{\max}(S)$, where $\beta_{\max}(S=0)\approx 11.7$ corresponds to
the commonly known value for Townes solitons Sulem ; Fibich , and
$\beta_{\max}(S\geq 1)$ are the respective values for the Townes solitons with
embedded vorticity Volkov , which are well approximated by
$\beta_{\max}(S)\approx 8\sqrt{3}\pi S$ Qin1 . In these intervals, the
solitons with $S=0$ are completely stable, while the VRs are stable in
narrower regions, $\beta<\beta_{\mathrm{st}}(S)$. An analytical approximation
yields $\beta_{\mathrm{st}}(S)=2\sqrt{6}\pi S\approx\allowbreak 15.4S$, while
an empirical formula for the numerical results is
$\beta_{\mathrm{st}}^{\mathrm{(num)}}(S)=15S-4$ Qin1 . Examples of stable VRs
with $S=1$ and $5$ are presented in Fig. 4. In intervals
$\beta_{\mathrm{st}}(S)<\beta\leq\beta_{\max}(S)$ the VRs are unstable against
splitting into necklace-shaped arrays of fragments. It is relevant to stress
that the growth of $\beta_{\mathrm{st}}^{\mathrm{(num)}}(S)$ $\char
126\relax~{}S$ for the “giant” VRs (ones with large values of $S$), produced
by Eq. (27), makes them much more robust modes than their counterparts with
smaller $S$. This feature is opposite to what was previously found in those
models which are able to produce stable VRs with $S>1$ Quiroga ; Pego ;
Borovkova ; Driben ; Sudharshan ; Reyna ; Zhang ; book .

Figure 4: Examples of the distribution of $\left|\phi(r)\right|$ in stable VRs
(vortex rings) with indicated values of $S$ and $\beta$. The initial shapes of
the VR ($\tau=0$) are compared to the outputs produced by simulations of Eq.
(27) at time $\tau=100$ (as per Ref. Qin1 ).

The finding reported in Ref. Qin2  is that Eqs. (27) are _exactly invariant_
with respect to a boost transformation, from the quiescent reference frame to
one which moves, in the 2D plane $\left(x,y\right)$, with vectorial
acceleration $\mathbf{a}=\left(a_{x},a_{y}\right)$, combined with a constant
velocity, $\mathbf{V}=\left(V_{x},V_{y}\right)$. The coordinates, wave
functions, and magnetic field in the accelerating frame are

$\left\\{x^{\prime},y^{\prime}\right\\}=\left\\{x.y\right\\}-V_{x,y}t-\left(a_{x,y}/2\right)t^{2},$
(28)
$\phi^{\prime}\left(x^{\prime},y^{\prime},t\right)=\phi\left(x,y,t\right)\exp\left[-i\left(a_{x}x+a_{y}y\right)t-i\left(V_{x}x+V_{y}y\right)+i\chi(t)\right],$
(29)
$\chi(t)=(1/6)\left[a_{x}^{-1}\left(V_{x}+a_{x}t\right)^{3}+a_{y}^{-1}\left(V_{y}+a_{y}t\right)^{3}\right],$
(30)
$H^{\prime}\left(x^{\prime},y^{\prime},t\right)=H\left(x,y,t\right)+\pi^{-1}\left(a_{x}x+a_{y}y\right).$
(31)

Actually, Eqs. (28)-(31) are a generalization of the usual Galilean boost for
the accelerating reference frame. Note that the solution of the 2D Poisson
equation in system (27), with the source represented by a quiescent soliton,
has the standard asymptotic form far from the region where the soliton is
located (it is determined by the Green’s function for the 2D Laplacian):

$H(\mathbf{r})\approx-(1/2\pi)\left(\int|\phi(\mathbf{r}^{\prime})|^{2}d\mathbf{r}^{\prime}\right)\ln
r.$ (32)

The difference of the magnetic-field component (31) of the SA soliton from its
quiescent counterpart (32) is the presence of the terms linear in $x$ and $y$,
which implies that the SA motion can be maintained by the properly constructed
background magnetic field. This field provides a reservoir of the momentum
which makes the self-acceleration possible.

According to Eq. (28), coordinates $\left(x_{c},y_{c}\right)$ of the center of
the stable 2D soliton moves as $x_{c}=V_{x}t+(1/2)a_{x}t^{2}$,
$y_{c}=V_{y}t+(1/2)a_{y}t^{2}$. This is a curvilinear trajectory in the 2D
plane: at small $t$, it is close to a straight line with slope
$x/y=V_{x}/V_{y}$, while at $t\rightarrow\infty$ it is close to a line with a
different slope, $x/y=a_{x}/a_{y}$. In particular, in the case of
$a_{x}=V_{y}=0$, the trajectory is a parabola:

$y_{c}=\left(a_{y}/2V_{x}^{2}\right)x_{c}^{2}.$ (33)

The analytical results are corroborated by Fig. 5, which displays stably
moving VRs produced by simulations of Eq. (27). The numerical solutions
demonstrate exactly the same SA motion of the VRs as predicted by Eq. (33).

Figure 5: (a) The plot of $\left|\phi\right|^{2}$ for a self-accelerating VR
with $a_{x}=10$, $a_{y}=0$, $V_{x.y}=0$ (see Eq. (28)), produced by the
numerical solution of Eq. (27) with $\beta=0$. (b) The same for $\beta=10$.
Right panels in (a) and (b) show the time dependence of coordinate
$\left\langle x\right\rangle$ of the VR’s center, cf. Eq. (33) (as per Qin2 ).

Essentially the same results were obtained for the SA motion of solitons
produced by the 1D version of Eq. (27), as well as for the 2D system based on
Eqs. (25) and (26) including the spin-orbit coupling between components
$\phi_{\downarrow}$ and $\phi_{\uparrow}$ Qin2 .

Conclusion. – This perspective represents basic models which make it possible
to predict counterintuitive regimes of motion of stable 1D and 2D solitons
(including vortex rings) with SA (self-acceleration). The corresponding models
represent two-component BECs. In one case, a pair of interacting 1D solitons
with opposite signs of the effective masses can be created in the binary BEC
loaded in the OL potential HS . In that case, the opposite interaction forces,
applied to the solitons with opposite signs of the mass, produce the SA
motion, similar to what was previously predicted Peschel theory  and
experimentally realized Peschel-experiment  in a nonsoliton form. In the
second case, the system of GP equations for the matter-wave components,
resonantly coupled by the magnetic component of the MW (microwave) field,
admits an exact transition to the accelerating references frame, thus
predicting stable 2D and 1D SA solitons, including VRs Qin2 .

A challenging issue is a possibility of experimental realization of the
predicted SA states. The theory also has a potential for development. In
particular, as concerns the SA pairs of solitons with positive and negative
masses, it is relevant to extend the analysis for the 2D system.

Acknowledgments – I thank Prof. Xiaojie Chen, EPL Deputy Editor, for the
invitation to write a perspective article. This work was supported, in part,
by the Israel Science Foundation through grant No. 1286/17.

## References

  * (1) Gagnon L. and Be1anger P. A., Opt. Lett. 15 (1990) 466-468.
  * (2) Parker D. F., Sophocleous C. and Radha C., J. Phys. A: Math. Gen. 35 (2002) 1283.
  * (3) Berry M. V. and Balazs N. L., Am. J. Phys. 47 (1979) 264-267.
  * (4) Siviloglou G. A. and Christodoulides D. N., Opt. Lett. 32 (2007) 979-981.
  * (5) Siviloglou G. A. et al., Phys. Rev. Lett. 99 (2007) 213901.
  * (6) Voloch-Bloch N. et al., Nature 494 (2013) 331-335.
  * (7) Minovich A. E. et al., Laser & Photonics Reviews 8 (2013) 221-232.
  * (8) Efremidis N. K., Paltoglou V. and von Klitzing W., Phys. Rev. A 87 (2013) 043637.
  * (9) Zhang P. et al., Nature Commun. 5 (2014) 4316.
  * (10) Clerici M. et al., Science Advances 1 (2015) 1400111.
  * (11) Fu S. et al., Phys. Rev. Lett. 115 (2015) 034501.
  * (12) Ellenbogen T. et al., Nature Phot. 3 (2009) 395-398.
  * (13) Hu Y. et al., Opt. Lett. 35, 3952-3954 (2010).
  * (14) Jia S. et al., Phys. Rev. Lett. 104, 253904 (2010).
  * (15) Kaminer I., Segev M., and Christodoulides D. N., Phys. Rev. Lett. 106 (2011) 213903.
  * (16) Lotti A. et al., Phys. Rev. A 84 (2011) 021807.
  * (17) Fattal Y., Rudnick A., and Marom D. M., Opt. Exp. 18 (2011) 17298-17307.
  * (18) Mayteevarunyoo T. and Malomed B. A., Opt. Lett. 40 (2015) 4947-4950.
  * (19) Batz S. and Peschel U., Phys. Rev. Lett. 110 (2013) 193901.
  * (20) Wimmer M. et al., Nature Phys. 9 (2013) 780-784.
  * (21) Kivshar Yu. S. and Malomed B. A., Rev. Mod. Phys. 61 (1989) 763-915.
  * (22) Reyna A. S. and C. B. de Araújo, Adv. Opt. Phot. 9 (2017) 720-774.
  * (23) Abdullaev F. Kh. and Salerno M., Phys. Rev. A 72 (2005) 033617.
  * (24) Sulem C. and Sulem P.-L., The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse (Springer, New York, 1999).
  * (25) Wang L., Malomed B. A. and Yan Z., Phys. Rev. E 99, 052206 (2019).
  * (26) Sakaguchi H. and Malomed B. A., Phys. Rev. E 99 (2019) 022216\.
  * (27) Chin C. et al., Rev. Mod. Phys. 82 (2010) 1225.
  * (28) Pu H. et al., Phys. Rev. A 67 (2003) 043605.
  * (29) H. Sakaguchi and Malomed B. A., J. Phys. B 37, 1443-1459 (2004).
  * (30) Brazhnyi V. A. and Konotop V. V., Mod. Phys. Lett. B 18 (2004) 627.
  * (31) Malomed B. A, Progr. Optics 43 (2002) 71-193.
  * (32) Bludov Yu. V. and García-Ñustes M. A., J. Phys. B: At. Mol. Opt. Phys. 50 (2017) 135004.
  * (33) Qin J. et al., Phys. Rev. A 99 (2019) 023610.
  * (34) Qin J., Dong G., and Malomed B. A., Phys. Rev. A 94 (2016) 053611.
  * (35) Sulem C. and Sulem P.-L., The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse (Springer, New York, 1999).
  * (36) Fibich G., The Nonlinear Schrödinger Equation: Singular Solutions and Optical Collapse (Springer, Heidelberg, 2015).
  * (37) Kruglov V. I. et al., J. Phys. A: Math. Gen. 21 (1988) 4381-4395.
  * (38) Quiroga-Teixeiro M. and Michinel H., J. Opt. Soc. Am. B 14 (1997) 2004-2009 (1997).
  * (39) Pego R. L. and Warchall H. A., J. Nonlinear Sci. 12 (2002) 347-394 (2002).
  * (40) Borovkova O. V. et al., Phys. Rev. E 84 (2011) 035602(R).
  * (41) Driben R. et al., Phys. Rev. Lett. 112 (2014) 020404.
  * (42) Sudharsan J. B. et al., Phys. Rev. A 92 (2015) 053601.
  * (43) Reyna A. S. et al., Phys. Rev. A 93 (2016) 013840.
  * (44) Zhang H. et al., Opt. Lett. 44 (2019) 3098-3101.
  * (45) Malomed B. A., Multidimensional Solitons (AIP Publishing, Melville, 2022).