# Multilevel Monte Carlo and its Applications in Financial Engineering

Devang Sinha  Department of Mathematics, Indian Institute of Technology
Guwahati, Guwahati-781039, Assam, India, e-mail<EMAIL_ADDRESS>Siddhartha
P. Chakrabarty  Department of Mathematics, Indian Institute of Technology
Guwahati, Guwahati-781039, Assam, India, e-mail<EMAIL_ADDRESS>Phone:
+91-361-2582606, Fax: +91-361-2582649

###### Abstract

In this article, we present a review of the recent developments on the topic
of Multilevel Monte Carlo (MLMC) algorithm, in the paradigm of applications in
financial engineering. We specifically focus on the recent studies conducted
in two subareas, namely, option pricing and financial risk management. For the
former, the discussion involves incorporation of the importance sampling
algorithm, in conjunction with the MLMC estimator, thereby constructing a
hybrid algorithm in order to achieve reduction for the overall variance of the
estimator. In case of the latter, we discuss the studies carried out in order
to construct an efficient algorithm in order to estimate the risk measures of
Value-at-Risk (VaR) and Conditional Var (CVaR), in an efficient manner. In
this regard, we briefly discuss the motivation and the construction of an
adaptive sampling algorithm with an aim to efficiently estimate the nested
expectation, which, in general is computationally expensive.

Keywords: MLMC; Option Pricing; Importance Sampling; Risk Management; Adaptive
Sampling

  * •

Brief overview of the Multilevel Monte Carlo estimator.

  * •

Review of the importance sampling algorithm to reduce the overall variance of
the MLMC estimator associated with the option pricing problems in financial
engineering.

  * •

Review of extension of adaptive sampling algorithm to multilevel paradigm,
with the aim of improving the computational complexity while estimating the
Value-at-Risk (VaR) and Conditional VaR (CVaR).

Subject area Mathematics and Statistics More specific subject area
Computational Finance Name of the reviewed methodology Importance Sampling for
Option Pricing and Adaptive Sampling in Efficient Risk Estimation. Keywords
MLMC; Option Pricing; Importance Sampling; Risk Management; Adaptive sampling
Resource availability N/A Review question 1\. What are recent development in
the field of Multilevel Monte Carlo? 2\. How have these developments led to
the improved efficiency of the existing estimator under various financial
applications? 3\. What are the shortcomings of the presented studies? 4\. What
are the recent developments catering to these shortcoming?

Table 1: Specifications Table

##  1 Introduction

In the broader area of computational finance, the mere establishment of the
existence of solution to a problem is not sufficient towards achieving the
tangible financial goals, for the problem that has been posed. Accordingly, as
is the case for many applications, we seek a solution that (in practice)
happens to be an approximation to the actual solution being sought. To this
end, we begin by observing that for problems in quantitative finance, one can
arrive at either an analytical or (possibly) a semi-analytical solution, in
only a handful of cases. Therefore, for the most part, one needs to devise
efficient methods to arrive at the desired and appropriate solution to the
posed problem which, in turn, necessitates the resorting to computational
techniques. At the heart of this paper, lies a specific computational
technique, widely used in the finance industry, namely, the Monte Carlo
simulation approach. Accordingly, we begin our presentation with a brief
narrative about this approach, with the main focus of the discussion being
steered towards the recent research and development in the area of Multilevel
Monet Carlo (MLMC). The interested readers may refer to [1, 2, 3, 4, 5, 6, 7]
for greater clarity on the approach of MLMC and the key developments with
respect to algorithms, as well as the applications in financial engineering
problems. In this article we primarily focus on the importance sampling
approach developed and studied in [8, 9] and also on how MLMC has led to the
development of algorithms for efficient risk estimation in the field of
financial risk management, discussed by the authors in [10]. However, we give
a brief overview of Monte Carlo and MLMC before directing our discussion
towards the aforesaid specific topics.

Monte Carlo methods have become one of the driving computing tools in the
finance industry. The necessity of simulating high-dimensional stochastic
models, which in turn may be attributed to the linear development in the
complexity corresponding to the size of the problem itself, is one of the
primary reasons that this approach is becoming the critical computational
strategy in the industry. The main objective of this method, in case of
computational finance is to reach the necessary degree of accuracy, which is
coupled with a high computational cost. More specifically, we intend to
approximate $\mathbf{E}[Y]$ where, $Y=G(X)$ is functional of the random
variable $X$. The traditional Monte Carlo approach requires a computational
complexity of an order of $O(\epsilon^{-3})$ to attain the root mean square
(RMS) error of $O(\epsilon)$ in a biased context [1]. This limitation led to
the introduction of the multilevel framework in [1] to address this issue and
achieve $O(\epsilon^{-2})$ computational complexity in the biased framework.

The idea behind the multilevel architecture is to employ independent standard
Monte Carlo on various resolution levels and use the differences as the
control variate for the Monte Carlo simulation at its most granular level,
which in mathematical terms is given by,

$\mathbf{E}[Y_{L}]=\mathbf{E}[Y_{1}]+\sum\limits_{l=2}^{L}\mathbf{E}[Y_{l}-Y_{l-1}],\leavevmode\nobreak\
\text{where}\leavevmode\nobreak\ Y_{l}=G(X_{l}).$ (1)

Using the standard Monte Carlo as the estimator to approximate the expectation
on the right hand side of (1), we obtain,

$\widehat{Y}_{L}=\frac{1}{N_{1}}\sum\limits_{k=1}^{N_{1}}Y_{1}^{k}+\sum\limits_{l=2}^{L}\frac{1}{N_{l}}\sum\limits_{k=1}^{N_{l}}\left(Y_{l}^{k}-Y_{l-1}^{k}\right),$
(2)

and therefore, $\mathbf{E}[Y_{L}]\approx\widehat{Y}_{L}$. Here $X_{l}$ is the
approximation of the random variable $X$ on level $l$ and this approximation
is contingent on the application under consideration. For example, if the
underlying stochastic process is driven by a stochastic differential equation
(SDE), then $X_{l}$ is the approximation of $X$, with some time discretization
parameter $h_{l}$. With all the preludes being presented in the preceding
discussion, we are now in a position to examine the following the seminal
result due to Giles [1].

######  Theorem 1.

Let $G$ denote a functional of the random variable X, and let $Y_{l}=G(X_{l})$
denote the corresponding level $l$ numerical approximation. If there exist
independent estimators $Z_{l}$, based on $N_{l}$ Monte Carlo samples, and
positive constants $\alpha,\beta,\gamma,c_{1},c_{2},c_{3},c_{4}$ such that
$\displaystyle{\alpha\geq\frac{1}{2}\min\left(\alpha,\beta\right)}$ and

  1. 1.

$\displaystyle{\left|\mathbf{E}\left[Y_{l}-Y_{l-1}\right]\right|\leq
c_{1}2^{-\alpha l}}$.

  2. 2.

$\displaystyle{\mathbf{E}\left[Z_{l}\right]=\begin{cases}\mathbf{E}\left[Y_{1}\right],\leavevmode\nobreak\
l=0,\\\ \mathbf{E}\left[Y_{l}-Y_{l-1}\right],\leavevmode\nobreak\
l>0.\end{cases}}$.

  3. 3.

$\displaystyle{\mathbf{V}[Z_{l}]\leq c_{2}N_{l}^{-1}2^{-\beta l}}$.

  4. 4.

$\displaystyle{C_{l}\leq c_{3}N_{l}2^{\gamma l}}$, where $C_{l}$ is the
computational complexity of $Z_{l}$,

then there exists a positive constant $c_{4}$ such that for any
$\epsilon<e^{-1}$, there are values $L$ and $N_{l}$ for which the multilevel
estimator $\displaystyle{Z=\sum\limits_{l=1}^{L}Z_{l}}$,has a MSE with bound,

$MSE\equiv\mathbf{E}\left[\left(Z-\mathbf{E}[Y]\right)^{2}\right]<\epsilon^{2},$

with a computational complexity $C$, having the bound,

$\mathbf{E}[C]\leq\begin{cases}c_{4}\epsilon^{-2},\leavevmode\nobreak\
\beta>\gamma,\\\
c_{4}\epsilon^{-2}\left(\log\epsilon\right)^{2},\leavevmode\nobreak\
\beta=\gamma,\\\
c_{4}\epsilon^{-2-\frac{(\gamma-\beta)}{\alpha}},\leavevmode\nobreak\
0<\beta<\gamma.\end{cases}$

It is quite evident from the above theorem that the computational complexity
is driven by the strong convergence of the estimator i.e.,
$\mathbf{V}[Z_{l}]$. Therefore, one of the main challenges while developing a
MLMC based estimator is to study the order of strong convergence of the
underlying approximation. With this brief introduction of MLMC, we now direct
our discussion towards its recent developments, pertaining to algorithm and
financial applications.

##  2 Importance Sampling Multilevel Algorithm

Since the advent of MLMC in literature, one of the directions of its
progression has been through various attempts to combine this algorithm, with
the already existing variance reduction techniques. For instance, Giles, in
[3, 11] studied and analyzed the combination of antithetic variates and MLMC
in order to bypass the Levy area simulation, encountered while using Milstein
discretization scheme, in order to simulate higher dimensional SDEs. However,
in our discussion we primarily focus on the combination of importance sampling
algorithm and Multilevel estimator.

The idea of incorporating importance sampling with multilevel estimators is
derived from the seminal paper by Arouna [12]. Arouna’s idea relied upon the
parametric change of measure and using a search algorithm to approximate the
optimal change of the measure parameter, in order to minimize the variance of
the standard Monte Carlo estimator. Before we discuss the research undertaken
in the area of multilevel pertaining to importance sampling algorithm, we give
a brief overview of the parametric importance sampling approach.

Consider a general problem of estimating $\mathbf{E}[G(X)]$, where $X$ is a
$d$-dimensional random variable. If $f(x)$ is the multivariate density
function, then,

$\mathbf{E}[G(X)]=\int G(x)f(x)dx=\int G(x+\theta)f(x+\theta)dx=\int
h(\theta,x)f(x)dx,$

where, $\displaystyle{h(\theta,x)=\frac{G(x+\theta)f(x+\theta)}{f(x)}}$. This
implies that, $\mathbf{E}[G(X)]=\mathbf{E}[h(\theta,X)]$. Therefore, we need
to determine the optimal value of $\theta$ such that $\text{Var}[h(\theta,X)]$
is minimum. Mathematically this is represented as,

$\theta^{*}=\arg\min_{\theta\in\mathbf{R}^{d}}\text{Var}[h(\theta,X)].$ (3)

In order to solve the above problem, one can resort to the usage of the
Robbins-Monro algorithm that deals with a sequence of random variable
$\left(\theta_{i}\right)_{i\in\mathbb{N}}$, which approximate $\theta^{*}$
accurately. However, the convergence of this algorithm requires certain
restrictive conditions, which are known as the non explosion condition (given
in [13]),

$\mathbf{E}[h^{2}\left(\theta,X\right)]\leq
C\left(1+|\theta|^{2}\right)\leavevmode\nobreak\ \text{for
all}\leavevmode\nobreak\ \theta\in\mathbf{R}^{d}.$

In order to deal with this restrictive condition, the authors in [14, 15]
introduced a truncation based procedure which was furthered in [16, 17]. An
unconstrained procedure to approximate $\theta^{*}$, by using the regularity
of the density function in an extensive manner, was introduced in [18] along
with the proof of convergence of the algorithm. Beside the stochastic
approximation algorithm, one can also use deterministic algorithm such as
sample average approximation, which, though being computationally expensive,
provides for a better approximation to $\theta^{*}$.

In problems dealing with the pricing of the options, for the most part the
underlying stochastic process $(X_{t})_{0\leq t\leq T}$ with $T>0$ being a
finite time horizon, is governed by some SDEs. The general form of these SDEs
is given as follows:

$dX_{t}=b(X_{t})dt+\sum\limits_{j=1}^{q}\sigma_{j}(X_{t})dW_{t}^{j},\leavevmode\nobreak\
X_{0}=x\in\mathbb{R}^{d},$ (4)

where, $\displaystyle{W:=\begin{pmatrix}W_{1}&W_{2}&\dots&W_{q}\end{pmatrix}}$
is a $q$-dimensional Brownian motion on a filtered probability space
$\left(\Omega,(\mathcal{F}_{t}\right)_{0\leq t\leq T},\mathbb{P})$, with
$b:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and
$\sigma_{j}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ being the functions
satisfying the following condition:

$\forall x,y\in\mathbb{R}^{d},\leavevmode\nobreak\ \lvert
b(x)-b(y)\rvert+\sum\limits_{j=1}^{q}\lvert\sigma_{j}(x)-\sigma_{j}(y)\rvert<K_{b,\sigma}\lvert
x-y\rvert,\leavevmode\nobreak\ \text{where}\leavevmode\nobreak\
K_{b,\sigma}>0.$ ($\text{C}_{1}$)

Assumption ($\text{C}_{1}$) ensures the existence and the uniqueness of
solution to (4). For the most part, constructing an analytical or semi-
analytical solution to (4) is not possible, and therefore we need to rely on
discretization schemes such as Euler or Milstein in order to simulate the
SDEs. For detailed discussion on these discretization schemes, the interested
readers may refer to [19]. Further, following the idea of [12], we consider a
family of stochastic process $\left(X_{t}(\theta)\right)_{0\leq t\leq T}$,
with $\theta\in\mathbb{R}^{d}$, being governed by the following SDE:

$dX_{t}(\theta)=(b(X_{t}(\theta))+\sigma(X_{t}(\theta))\theta)dt+\sum\limits_{j=1}^{q}\sigma_{j}(X_{t}(\theta))dW_{t}^{j},\leavevmode\nobreak\
\sigma(x)=\begin{pmatrix}\sigma_{1}(x)&\dots&\sigma_{q}(x)\end{pmatrix}.$ (5)

As a consequence of the Girsanov’s Theorem, we know that there exists a risk-
neutral probability measure $\mathbf{P}_{\theta}$, which is equivalent to
$\mathbf{P}$ such that,

$\frac{d\mathbf{P}_{\theta}}{{d\mathbf{P}|}_{\mathcal{F}_{t}}}=\exp{\left(-\langle\theta,W_{t}\rangle-\frac{1}{2}\lvert\theta\rvert^{2}t\right)},$
(6)

under which the process $\displaystyle{\left(\theta t+W_{t}\right)_{0\leq
t\leq T}}$ is a Brownian motion. Therefore,

$\mathbf{E}_{\mathbf{P}}\left[G(X_{T})\right]=\mathbf{E}_{\mathbf{P}_{\theta}}\left[G(X_{T}(\theta))\right]=\mathbf{E}_{\mathbf{P}}\left[G(X_{T}(\theta))e^{-\langle\theta,W_{T}\rangle-\frac{1}{2}\lvert\theta\rvert^{2}T}\right]$
(7)

Therefore, following the discussion above we have,

$\mathbf{E}\left[G(X_{T})\right]=\mathbf{E}\left[h(\theta,X_{T})\right].$

here,
$h(\theta,X_{T})=G(X_{T}(\theta))e^{-\langle\theta,W_{T}\rangle-\frac{1}{2}\lvert\theta\rvert^{2}T}$.
Now the idea of importance sampling Monte Carlo method is to estimate
$\mathbf{E}\left(G(X_{T})\right)$, where $\theta$ is given by,

$\theta^{*}=\arg\min_{\theta\in\mathbb{R}^{d}}\text{Var}\leavevmode\nobreak\
\left(G(X_{T}(\theta))e^{-\langle\theta,W_{T}\rangle-\frac{1}{2}\lvert\theta\rvert^{2}T}\right).$
(8)

In the context of Multilevel estimator, we present two approaches studied in
[8, 9], adapting the ideas studied by authors in [12, 13] and extending it to
multilevel scenarios. Under the parametric change of measure, the general
multilevel estimator is given defines as,

$\mathbf{E}[Y_{L}]=\mathbf{E}[Y_{1}^{\theta_{1}}]+\sum\limits_{l=2}^{L}\mathbf{E}[Y_{l}^{\theta_{l}}-Y_{l-1}^{\theta_{l}}],\leavevmode\nobreak\
\text{where}\leavevmode\nobreak\
Y_{l}^{\theta}=G(X_{l}^{\theta})e^{-\langle\theta_{l},W_{T}^{l}\rangle-\frac{1}{2}\lvert\theta_{l}\rvert^{2}T}.$
(9)

Under the framework of multilevel estimator, the parametric importance
sampling estimator looks like,

$\displaystyle\widehat{Y}_{L}^{\theta}=\frac{1}{N_{1}}\sum\limits_{k=1}^{N_{1}}Y_{1}^{k,\theta_{1}}+\sum\limits_{l=2}^{L}\frac{1}{N_{l}}\sum\limits_{k=1}^{N_{l}}(Y_{l}^{k,\theta_{l}}-Y_{l-1}^{k,\theta_{l}}).$
(10)

Considering the variance of the above estimator, we have [8],

$\displaystyle\text{Var}[\widehat{Y}_{L}^{\theta}]=\frac{1}{N_{1}}\mathbf{v}_{1}(\theta_{1})+\sum\limits_{l=2}^{L}\frac{1}{N_{l}}\sum\limits_{k=1}^{N_{l}}\frac{(M-1)T}{M^{l}}\mathbf{v}_{l}(\theta_{l}),$
(11)

where,

$\mathbf{v}_{1}(\theta_{1})=\text{Var}[Y_{1}^{\theta_{1}}]\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\
\mathbf{v}_{l}(\theta_{l})=\text{Var}[Y_{l}^{\theta_{1}}-Y_{l-1}^{\theta_{1}}].$

Therefore, as discussed, in order to solve the problem of minimizing the
overall variance of the estimator described above, we intend to minimize the
variance at each level of resolution, i.e., we aim at approximating
$\theta_{l}^{*}$ for $l=1,\dots,L$, such that,

$\theta_{1}^{*}=\arg\min_{\theta\in\mathbf{R}^{d}}v_{1}(\theta_{1})\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\
\theta_{l}^{*}=\arg\min_{\theta\in\mathbf{R}^{d}}v_{l}(\theta_{l}).$ (12)

Further, pertinent to the discussion carried out in [13] and another
application of the Girsanov’s Theorem, the above problem can be reformulated
as,

$\displaystyle\theta_{1}^{*}$ $\displaystyle=$
$\displaystyle\arg\min_{\theta_{1}\in\mathbf{R}^{d}}\mathbf{E}\left[G(X_{1})^{2}e^{-\langle\theta_{1},W_{T}^{1}\rangle+\frac{1}{2}\lvert\theta_{1}\rvert^{2}T}\right]$
$\displaystyle\theta_{l}^{*}$ $\displaystyle=$
$\displaystyle\arg\min_{\theta\in\mathbf{R}^{d}}\mathbf{E}\left[\frac{M^{l}}{(M-1)T}\left(G(X_{l})-G(X_{l-1})\right)^{2}e^{-\langle\theta_{l},W_{T}^{l}\rangle+\frac{1}{2}\lvert\theta_{l}\rvert^{2}T}\right].$
(13)

We present below the two algorithm namely, the sample average approximation
and stochastic approximation, in order to approximate the $\theta_{l}$’s as
the solution to (2).

###  2.1 Sample Average Approximation

The sample average approximation deals with approximating the above
expectations using $\widetilde{N}_{l}$ sample paths. More specifically,

$\mathbf{E}\left[G(X_{1})^{2}e^{-\langle\theta_{1},W_{T}^{1}\rangle+\frac{1}{2}\lvert\theta_{1}\rvert^{2}T}\right]\approx\frac{1}{\widetilde{N}_{1}}\sum\limits_{j=1}^{\widetilde{N}_{1}}G(X_{1}^{k})^{2}e^{-\langle\theta_{1},W_{T}^{1,k}\rangle+\frac{1}{2}\lvert\theta_{1}\rvert^{2}T}\equiv\mathcal{V}_{1},$
(14)

and,

$\displaystyle\mathbf{E}\left[\left(G(X_{l})-G(X_{l-1})\right)^{2}e^{-\langle\theta_{l},W_{T}^{l}\rangle+\frac{1}{2}\lvert\theta_{l}\rvert^{2}T}\right]$
(15) $\displaystyle\approx$
$\displaystyle\frac{1}{\widetilde{N}_{l}}\sum\limits_{j=1}^{\widetilde{N}_{l}}\frac{M^{l}}{(M-1)T}\left(G(X_{l}^{k})-G(X_{l-1}^{k})\right)^{2}e^{-\langle\theta_{l},W_{T}^{l,k}\rangle+\frac{1}{2}\lvert\theta_{l}\rvert^{2}T}\equiv\mathcal{V}_{l}.$

Having approximated the expectation in the minimization problem, the authors
used the standard Newton-Raphson algorithm on the functions $\mathcal{V}_{1}$
and $\mathcal{V}_{l}$ in order to approximate $\theta_{l}^{*}$ for
$l=1,\dots,L$. In [13] it is proved that if the functional $G(X)$ satisfies
the non-degeneracy conditions i.e., $\mathbf{P}((G(X^{1}_{T})\neq 0)>0$ and
$\mathbf{P}\left(((G(X^{l}_{T})-G(X^{l-1}_{T}))\neq 0\right)>0$ and further
assuming they have finite second moment, then by Lemma 2.1 in [13],
$\mathcal{V}_{1}$ and $\mathcal{V}_{l}$ are infinitely continuously
differentiable. Moreover, both $\mathcal{V}_{1}$ and $\mathcal{V}_{l}$ are
both strongly convex, thus implying the existence of the unique minimum
$\theta_{1}^{*}$ and $\theta_{l}^{*}$ as the solution to equation (2).

###  2.2 Adaptive Stochastic Approximation

Under the stochastic approximation, studied in [9] the aim of determining the
optimal change of parameter $\theta_{l}^{*}$ for $l=1,\dots,L$ is carried out
using the Robbins-Monro algorithm. Here, we briefly describe the algorithm.
Consider a compact convex set $\Theta\subset\mathbb{R}^{q}$ such that
$0\in\text{int}(\Theta)$. Then the recursive algorithm with projection is
defined as follows,

$\theta_{l}^{n+1}=\textbf{Proj}_{\Theta}\left[\theta_{l}^{n}-\gamma_{n+1}H_{l}(\theta_{l}^{n},Y_{l},W_{T}^{l})\right],$
(16)

where,
$\text{Proj}_{\Theta}(\theta)=\min_{\theta\in\Theta}\lvert\theta-\theta_{0}\rvert$.
The sequence $(\gamma_{n})_{n\geq 1}$ is a decreasing sequence of positive
real numbers satisfying,

$\sum\limits_{n=1}^{\infty}\gamma_{n}=\infty\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\
\sum\limits_{i=1}^{\infty}\gamma_{n}^{2}<\infty.$ (17)

Also,

$H_{l}(\theta_{l}^{n},Y_{l},W_{T}^{l})=\left\\{\begin{array}[]{ll}\left(\theta_{1}T-W_{T}^{1}\right)\left(G(X_{1})^{2}e^{-\langle\theta_{1},W_{T}^{1}\rangle+\frac{1}{2}\lvert\theta_{1}\rvert^{2}T)}\right),&l=1,\\\
\left(\theta_{l}T-W_{T}^{l}\right)\left[\frac{M^{l}}{(M-1)T}\left(G(X_{l})-G(X_{l-1})\right)^{2}e^{-\langle\theta_{l},W_{T}^{l}\rangle+\frac{1}{2}\lvert\theta_{l}\rvert^{2}T}\right],&l=2,\dots,L.\end{array}\right\\}$
(18)

The algorithm described above is the constrained version of the Robbins-Monro
algorithm. The inclusion of the projection operator in the recursive algorithm
is to satisfy the non-explosion condition described above. Similar to the
discussion carried out in the previous section, if the non-degeneracy
conditions are satisfied i.e., $\mathbf{P}\left(G(X^{1}_{T})\neq 0\right)>0$
and $\mathbf{P}\left(\left(G(X^{l}_{T})-G(X^{l-1}_{T})\right)\neq 0\right)>0$,
further assuming the finite second moment of $G(X_{1})$ and
$G(X_{l})-G(X_{l-1})$, we can conclude the convergence of the
$\theta_{l}^{*}$, constructed recursively using equation (16), for various
level of resolutions.

The term adaptive is used in the sense that, the estimation of the optimal
importance sampling parameter and the multilevel Monte Carlo run
simultaneously. The multilevel estimator in this case is given as follows,

$\widehat{Y}_{L}^{\theta}=\frac{1}{N_{1}}\sum\limits_{k=1}^{N_{1}}Y_{1}^{k,\theta_{1}^{k-1}}+\sum\limits_{l=2}^{L}\frac{1}{N_{l}}\sum\limits_{k=1}^{N_{l}}\left(Y_{l}^{k,\theta_{l}^{k-1}}-Y_{l-1}^{k,\theta_{l}^{k-1}}\right).$
(19)

However, for the purpose of the practical implementation, one needs to stop
the approximations procedure after finite number of iterations.

Having approximated the $\theta_{l}^{*}$ for $l=1,\dots,L$, we use the
multilevel algorithm described by equation (16) to estimate our expectation.
It is quite evident from the way the algorithms have been described that the
importance sampling algorithm combined with a multilevel estimator is more
computationally complex than the standard multilevel algorithm. However, the
variance reduction achieved by these combinations compensates for the high
computational complexity. That is, the hybrid algorithm achieves the desired
RMS error much faster than the MLMC estimator. The studies carried out in [8,
13] demonstrate the accuracy of the hybrid importance sampling multilevel
algorithm over standard multilevel algorithm, through a series of numerical
examples, where the underlying SDEs are multi-dimensional. The slight drawback
of the sample average approximation method, though more stable than the
adaptive stochastic algorithm, is the slow convergence rate to the optimal
value. As for the stochastic approximation, the algorithm is sensitive to the
learning parameter $\gamma_{n}$ and therefore is unstable. It may be pointed
out that the study performed above only deals with the Euler Multilevel Monte
Carlo, restricted to the use of Euler discretization to simulate the
underlying SDEs. More recently, a study carried out by authors in [20]
generalize this approach, undertaking higher order discretization schemes such
as Milstein to simulate the underlying SDEs. The interested reader can refer
to the references mentioned therein to get a more rigorous understanding of
this hybrid algorithm.

##  3 MLMC and Efficient Risk Estimation.

Risk measurement and consequent management is one of the essential components
of financial engineering. The computation of the former (risk measures) for a
financial portfolio is both challenging and computationally intensive, which
may be ascribed to computations involving nested expectation, which entails
multiple evaluations of the loss to the portfolio, for distinct risk
scenarios. Further, the cost of computing loss of portfolio entailing
thousands of derivatives becomes progressively expensive with an increase in
the size of the portfolio [21]. Value-at-Risk (VaR), Conditional VaR (CVaR),
and the likelihood of a large loss are the necessary risk metrics used to
estimate the risk of a financial portfolio. At the core of these estimation,
is the necessity of evaluating the nested expectation, given by,

$\eta=\mathbf{E}\left[H\left(\mathbf{E}[X\lvert Y]\right)\right]$ (20)

where, $H$ is the Heaviside function. More specifically, suppose we need to
compute the probability of the expected loss being greater than
$L_{\eta}\in\mathbf{R}$, i.e., we are interested in the following computation:

$\eta=\mathbf{E}\left[H(\mathbf{E}[\Delta\lvert R_{\tau}]-L_{\eta})\right],$
(21)

where $\mathbf{E}[\Delta\lvert R_{\tau}]$ is the expected loss in a risk-
neutral world, with $R_{\tau}$ being a possible risk scenario at some short
risk (time) horizon $\tau$. Also, $\Delta$ is the average loss of many losses
incurred from different financial derivatives, depending upon similar
underlying assets [21], that is,

$\Delta=\frac{1}{K}\sum\limits_{i=1}^{K}\Delta_{i},$ (22)

where $K$ is the total number of derivatives and $\Delta_{i}$ is the loss from
the $i$-th derivative. The average is considered to ensure the boundedness of
$\Delta$, when the portfolio size of $K$ increases. A standard and straight
forward way to estimate the nested expectation (20) is the usage of Monte
Carlo method. This involves, simulating $M$ independent scenarios of the risk
parameter $R_{\tau}$, and for each risk scenario, $N$ total loss samples,
which are independent. This method was explored in [22] and an extended
analysis was carried out in [23]. The total computational cost to perform the
above simulation is $O(\max(K\epsilon^{-2},\epsilon^{-3}))$ in order to
achieve the root-mean-squared (RMS) error of $\epsilon$ [21]. In order to
handle this issue we present the ideas studied in [10] under the realm of
MLMC.

###  3.1 Adaptive Sampling Multilevel estimator

As mentioned in the previous section, the cost of the standard Monte Carlo to
achieve the root-mean-squared error of $\epsilon$ is $O(\epsilon^{-3})$. To
improve the computational complexity, the authors in [24] developed an
efficient through the adaptation of the sample size required in the inner
sampler of Monte Carlo, to the particular outer sampler random variable. Under
certain conditions, the authors were able to achieve the $O(\epsilon^{-5/2})$
computational complexity to achieve the RMS of $\epsilon$. Giles in [10]
extended this approach to the multilevel framework and was able to achieve
$O\left(\epsilon^{-2}\lvert\log\epsilon\rvert^{2}\right)$ computational cost
for a RMS error tolerance $\epsilon$. Before presenting the work initiated by
Giles, we put forth a brief review of the studies carried out in [22] and
[24].

The authors in [22], estimated the inner expectation of the equation (20),
i.e., $\mathbf{E}[X\lvert Y=y]$, for a given $y$, using the unbiased Monte
Carlo estimator, with $N$ sample paths, as given by,

$\widehat{Z}_{N}(y)=\frac{1}{N}\sum\limits_{n=1}^{N}x_{n}(y),$ (23)

where, $\\{x_{n}(y)\\}_{n}$ are the mutually independent samples from the
random variable $X$, conditioned on $Y=y$. Again using the Monte Carlo for the
outer expectation, we have,

$\eta\approx\frac{1}{M}\sum\limits_{m=1}^{M}H\left(\widehat{Z}_{N}(y_{m})\right),$
(24)

where $\\{y_{m}\\}_{m}$ are the mutually independent samples from the random
variable $Y$. Further, they proved that if the two random variables
$\mathbf{E}[X\lvert Y]$ and $\widehat{Z}_{N}$ have the joint density
$d_{N}(y,z)$ and assuming that for $i=0,1,2$,
$\displaystyle{\frac{\partial}{\partial y_{i}}d_{N}(y,z)}$ exists, plus there
exists a non-negative function $d_{i,N}$, such that,

$\left\lvert\frac{\partial}{\partial y_{i}}d_{N}(y,z)\right\rvert\leq
d_{i,N},\leavevmode\nobreak\ \text{for all}\leavevmode\nobreak\
N,y,z,\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\sup_{N}\int\limits_{-\infty}^{\infty}\lvert z\rvert^{q}d_{i,N}(z)dz<\infty,$

for all $0\leq q\leq 4$, then the RMS error of the estimator (24) is
$O\left(M^{-1/2}+N^{-1}\right)$. Therefore, in order to achieve the RMS error
of $O(\epsilon)$ we need $M=O(\epsilon^{-2})$ and $N=O(\epsilon^{-1})$,
leading to the total computational complexity of $O(\epsilon^{-3})$. Authors
in [24] developed an adaptive sampling technique to deal with high
computational complexity previously discussed. Their approach was based on the
likelihood that an additional sample will result in a negative estimate of
$\widehat{Z}_{N+1}$ having estimated that $\widehat{Z}_{N}>0$ for given $Y$.
More specifically they showed that,

$\mathbf{P}\left[\widehat{Z}_{N+1}\leq
0\lvert\widehat{Z}_{N}\right]\leq\frac{\sigma^{2}}{\left(N\widehat{Z}_{N}(Y)+\mu\right)^{2}}\approx\frac{\sigma^{2}}{N^{2}\mu^{2}},$

where $\mu=\mathbf{E}[X\lvert Y]$ and $\sigma^{2}=\text{Var}[X\lvert Y]$.
Therefore, if
$\displaystyle{N\geq\frac{\epsilon^{-1/2}\sigma}{\lvert\mu\rvert}}$, then the
probability that
$H\left(\widehat{Z}_{N}(Y)\right)=H\left(\widehat{Z}_{N+1}(Y)\right)\approx
H\left(\mathbf{E}[X\lvert Y]\right)$ is equal to $1-\epsilon$. Based on these
observations, the authors in [24] introduced two algorithms, the first being
based on the minimization of the total number of samples for all inner Monte
Carlo samplers with respect to given tolerance $\epsilon$, and the second
being iterative, estimating $\lvert\mu\rvert$ and $\sigma$ after every
iteration, for given value of $Y$, using $N$ samples further adding more inner
samples till $\displaystyle{\frac{N\mu}{\sigma}}$ exceeds some error margin
threshold. Under these two algorithms it was observed that the overall
computational complexity is $O(\epsilon^{-5/2})$ [10]. The authors in [10]
introduced the above algorithms in the realm of multilevel simulation, wherein
they used multilevel estimator in order to achieve an approximation to the
outer expectation, while making use of the sample size in the inner
expectation as the discretization parameter. More specifically,

$\widetilde{\eta}\coloneqq\sum\limits_{l=0}^{L}\frac{1}{M_{l}}\sum\limits_{m=1}^{M_{l}}H\left(\widehat{Z}_{N_{l}}^{f,l,m}(y^{l,m})\right)-H\left(\widehat{Z}_{N_{l-1}}^{c,l,m}(y^{l,m})\right),$
(24)

where,

$\widehat{Z}_{N_{l}}^{f,l,m}(y)=\frac{1}{N_{l}}\sum\limits_{n=1}^{N_{l}}x^{f,l,m,n}(y),$
(25)

with $\\{x^{.,l,m,n}(y)\\}$ being the i.i.d samples of the random variable
$X$, given $Y=y$. Also, $H\left(\widehat{Z}_{-1}^{c,0,\dots}(y)\right)\equiv
0$. Now under the assumptions 3.1, it can be proved that [10],

$\bigg{\lvert}\mathbf{E}\left[H\left(\widehat{Z}_{N_{l}}(Y)\right)-H\left(\mathbf{E}[X\lvert
Y]\right)\right]\bigg{\rvert}=O\left(N_{l}^{-1}\right).$

Further, under the assumption that there exists constants $\delta_{0}$ and
$\rho_{0}$ such that, $\displaystyle{\rho(\delta)\leq\rho_{0}}$, for all
$\delta\in[0,\delta_{0}]$ where $\delta$ is the random variable with density
$\rho$, the authors in [10] proved that, if $X$ and $Y$ are the two random
variables, satisfying the stated assumption, then,

$\text{Var}\left[H\left(\widehat{Z}_{N}(Y)\right)-H\left(\mathbf{E}[X\lvert
Y]\right)\right]=O(N^{-1/2}).$ (26)

The above result determines the strong convergence property necessary to
analyze the full potential of the MLMC estimator, in this scenario. However,
if $\displaystyle{N_{l}=N_{0}2^{l}}$, then with standard MLMC complexity
analysis it is easy to determine that the computational complexity required to
achieve RMS error of $\epsilon$, we need $O\left(\epsilon^{-5/2}\right)$
computational complexity. To cater to this high computational demand, even in
the framework of MLMC, the authors undertook the adaptive approach developed
in [24] and extended it to the framework of MLMC.

Giles extended the studies carried out by authors in [24] to multilevel
paradigm with an aim to reduce the overall computational cost to
$O\left(\epsilon^{-2}\lvert\log\epsilon)\rvert^{2}\right)$. In addition to the
assumptions stated above, it is further assumed that,

$\sup_{y}\mathbf{E}[\sigma^{-q}\lvert X-\mathbf{E}[X\lvert Y]\rvert^{q}\lvert
Y=y]<\infty,\quad 2<q<\infty.$ ($C_{3}$)

Thus, under the above stated assumptions, it was proved in Lemma 2.5 (for the
perfect adaptive sampling) and Theorem 2.7 of [10], that if the maximum number
of sample paths is restricted to,

$N=\bigg{\lceil}\max\left(O\left(\epsilon^{-1}\right),C^{2}\frac{\sigma^{2}}{\lvert\mu\rvert^{2}}\right)\bigg{\rceil},$
(27)

then the further number of sample path of various level of resolutions are
given by,

$N_{l}=\bigg{\lceil}N_{0}4^{l}\max\left(2^{-l},\min\left(1,\left(C^{-1}N_{0}^{1/2}2^{l}\frac{\lvert\mu\rvert}{\sigma}\right)^{-r}\right)\right)\bigg{\rceil},$
(28)

with $C$ being some confidentiality constant and
$\displaystyle{1<r<2-\frac{2}{q}}$ for the perfect adaptive sampling and
$\displaystyle{1<r<2-\frac{\sqrt{4q+1}-1}{q}}$ when the values of
$\lvert\mu\rvert$ and $\sigma$ is approximated. Therefore,

$\text{Var}\left[H\left(\widehat{Z}_{N}(Y)\right)-H\left(\mathbf{E}[X\lvert
Y]\right)\right]=O\left(2^{-l}\right),$ (29)

thereby leading to the overall computational complexity of the desired order.
In a detailed discussion carried out in Section 4 of [10], it was proved
(pertaining to the calculation of VaR and CVaR) that in order to achieve the
overall computational cost of $O(\epsilon)$ RMS error, the required
computational complexity is
$O\left(\max(\epsilon^{-2}\lvert\log\epsilon)\rvert,K\epsilon^{-2})\right)$
for the estimation of VaR and CVaR, respectively. The numerical test on a
model problem undertaken shows the efficacy of the algorithm constructed.
Readers are directed to the referred paper for detailed discussion on the
proofs of the above stated results. It may be noted that the computational
complexity increases with an increase in the portfolio size, $K$. A random
sub-sampling approach, extending it to a multilevel framework, thereby
addressing the dependency on the portfolio size, to achieve the desired RMS
error was recently introduced in [21].

##  4 Conclusion

In this paper, we gave a brief overview of the recent trends in the paradigm
of the multilevel algorithm concerning the importance sampling, in the case of
option pricing and an adaptive sampling approach while determining the VaR and
CVaR for large portfolios. The algorithms discussed serves as the improvement
in the computational efficiency of the standard multilevel estimators, each
having its merits and shortcomings. As discussed in Section 2, the importance
sampling algorithm combined with multilevel estimators significantly decreases
variance at various resolution levels. However, the decrease in variance comes
at the cost of increased computational complexity in either case and an
increase in the sensitivity to approximate the optimal parameter. As for
developing the MLMC based algorithm for efficient risk estimation discussed in
Section 3, the adaptive sampling approach introduced in this paradigm leads to
a significant improvement in the overall computational complexity to achieve
the desired root mean squared error. However, the dependence of computational
complexity on the size of the portfolio is a subtle shortcoming of the
discussed algorithm. Overall, the presented ideas have substantially
contributed to the research and development of the multilevel algorithm for
various applications encountered in financial engineering problems. However,
the scope to enrich the standard algorithm with non-standard variance
reduction techniques is still an exciting path for future research.

## References

  * [1] M. B. Giles, Multilevel monte carlo path simulation, Operations research 56 (3) (2008) 607–617.
  * [2] M. Giles, Improved multilevel monte carlo convergence using the milstein scheme, in: Monte Carlo and Quasi-Monte Carlo Methods 2006, Springer, 2008, pp. 343–358.
  * [3] M. B. Giles, L. Szpruch, Antithetic multilevel monte carlo estimation for multidimensional sdes, in: Monte Carlo and Quasi-Monte Carlo Methods 2012, Springer, 2013, pp. 367–384.
  * [4] M. Giles, L. Szpruch, Multilevel monte carlo methods for applications in finance, Recent Developments in Computational Finance: Foundations, Algorithms and Applications (2013) 3–47.
  * [5] V. Lemaire, G. Pagès, Multilevel richardson–romberg extrapolation, Bernoulli 23 (4A) (2017) 2643–2692.
  * [6] D. Belomestny, J. Schoenmakers, F. Dickmann, Multilevel dual approach for pricing american style derivatives, Finance and Stochastics 17 (4) (2013) 717–742.
  * [7] S. Heinrich, Multilevel monte carlo methods, in: International Conference on Large-Scale Scientific Computing, Springer, 2001, pp. 58–67.
  * [8] A. Kebaier, J. Lelong, Coupling importance sampling and multilevel monte carlo using sample average approximation, Methodology and Computing in Applied Probability 20 (2) (2018) 611–641.
  * [9] M. B. Alaya, K. Hajji, A. Kebaier, Improved adaptive multilevel monte carlo and applications to finance, arXiv preprint arXiv:1603.02959 (2016).
  * [10] M. B. Giles, A.-L. Haji-Ali, Multilevel nested simulation for efficient risk estimation, SIAM/ASA Journal on Uncertainty Quantification 7 (2) (2019) 497–525.
  * [11] M. B. Giles, L. Szpruch, Antithetic multilevel monte carlo estimation for multi-dimensional sdes without lévy area simulation, The Annals of Applied Probability 24 (4) (2014) 1585–1620.
  * [12] B. Arouna, Adaptative monte carlo method, a variance reduction technique (2004).
  * [13] M. B. Alaya, K. Hajji, A. Kebaier, Importance sampling and statistical romberg method, Bernoulli 21 (4) (2015) 1947–1983.
  * [14] H.-F. Chen, L. Guo, A.-J. Gao, Convergence and robustness of the robbins-monro algorithm truncated at randomly varying bounds, Stochastic Processes and their Applications 27 (1987) 217–231.
  * [15] H.F.Chen, Y. Zhu, Stochastic approximation procedures with randomly varying truncations, Science in China, Ser. A (1986).
  * [16] C. Andrieu, É. Moulines, P. Priouret, Stability of stochastic approximation under verifiable conditions, SIAM Journal on control and optimization 44 (1) (2005) 283–312.
  * [17] J. Lelong, Almost sure convergence of randomly truncated stochastic algorithms under verifiable conditions, Statistics & Probability Letters 78 (16) (2008) 2632–2636.
  * [18] V. Lemaire, G. Pagès, Unconstrained recursive importance sampling, The Annals of Applied Probability 20 (3) (2010) 1029–1067.
  * [19] P. E. Kloeden, E. Platen, Stochastic differential equations, in: Numerical Solution of Stochastic Differential Equations, Springer, 1992, pp. 103–160.
  * [20] D. Sinha, S. P. Chakrabarty, Multilevel richardson-romberg and importance sampling in derivative pricing, arXiv preprint arXiv:2209.00821 (2022).
  * [21] M. B. Giles, A.-L. Haji-Ali, Sub-sampling and other considerations for efficient risk estimation in large portfolios, arXiv preprint arXiv:1912.05484 (2019).
  * [22] M. B. Gordy, S. Juneja, Nested simulation in portfolio risk measurement, Management Science 56 (10) (2010) 1833–1848.
  * [23] D. Giorgi, V. Lemaire, G. Pagès, Limit theorems for weighted and regular multilevel estimators, Monte Carlo Methods and Applications 23 (1) (2017) 43–70.
  * [24] M. Broadie, Y. Du, C. C. Moallemi, Efficient risk estimation via nested sequential simulation, Management Science 57 (6) (2011) 1172–1194.


[4]Xingran Xu

1]School of Information, Hunan University of Humanities, Science and
Technology, Loudi 417000, China

2]School of Physics and Electronic Science, Zunyi Normal University, Zunyi
563006, China 3]Department of Physics, Zhejiang Normal University, Jinhua
321004, China [4]School of Science, Jiangnan University, Wuxi 214122, China

# Dynamic protected states in the non-Hermitian system

Lei Chen<EMAIL_ADDRESS>Zhen-Xia Niu<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [  *

###### Abstract

The non-Hermitian skin effect and nonreciprocal behavior are sensitive to the
boundary conditions, which are unique features of non-Hermitian systems. The
eigenenergies will become complex and all eigenstates are localized at the
boundary, which is distinguished from the Hermitian topologies. In this work,
we theoretically study the dynamic behavior of the propagation of Gaussian
wavepackets inside a non-Hermitian lattice and analyze the self-acceleration
process of bulk state or Gaussian wavepackets toward the system’s boundary.
The initial wavepackets will not only propagate toward the side where the
eigenstates are localized, but also their momentum will approach to a specific
value where the imaginary parts of energy dispersion are the maximum. In
addition, if the wavepackets cover this specific momentum, they will
eventually exhibit exponentially increasing amplitudes with time evolution,
maintaining the dynamic protected condition for an extended period of time
until they approach the boundary. We also take two widely used toy models as
examples in one and two dimensions to verify the correspondence of the non-
Hermitian skin effect and the dynamic protected state.

###### keywords:

skin effect, photonic lattice, non-Hermitian Hamiltonian, dynamic behavior

##  1 Introduction

For an isolated quantum system, to ensure the conservation of energy and
probability, a Hermitian Hamiltonian is necessary. When a quantum system
couples to its environment, hence, particle loss and decoherence leads to the
breakdown of Hermiticity and requires non-Hermitian descriptions [1, 2, 3, 4,
5]. It has been shown that non-Hermitian considerations not only provide
suitable descriptions of open systems, but also bring novel physics and
unprecedented phenomena and applications, as found in a variety of physical
realms including parity-time symmetry and spectral singularity [6, 7, 8, 9],
non-reciprocal and chiral transport [10], as well as non-Hermitian topology
and unconventional band theory [11, 12, 13, 14]. Within this context,
particularly fascinating phenomenon is the recently discovered non-Hermitian
skin effect (NHSE) [14, 15, 16, 17].

The NHSE, an anomalous localization of extensive eigenstates at the open
boundaries of a non-Hermitian system, is drastically different from the
dynamics of the extended Bloch waves in Hermitian systems. This phenomenon was
first discovered by Hatano and Nelson [18, 19] in the late 1900s, where they
proposed a 1D disordered tight-binding lattice model with nearest-neighbour
nonreciprocal hopping and showed that non-Hermiticity induced by the
nonreciprocity can prevent Anderson localization, opening up a mobility region
characterized by unidirectional transportation. Recently, similar phenomena
have also been observed in many seminal works on the topological properties of
non-Hermitian systems. They either revisit the Hatano-Nelson (HN) model to
realize versatile directional transports [20, 21], interpret NHSE from
different perspectives [22, 23, 24, 25, 26, 27], even to develop new theories,
methods, and material designs to describe NHSE related topological properties
and the associated non-Hermitian bulk-boundary correspondence (BBC) [11, 15,
28, 29, 30, 31], or extend NHSE to higher dimensions [32, 33, 34, 35], to
explore new NHSE features enriched by heterogeneous degrees of freedom [36,
37, 38, 39, 40], and experimentally demonstrating the NHSE [41, 42, 43, 44].

One of the significant features of the NHSE is that the eigenstates are
sensitive to the boundary conditions. The bulk properties of the open boundary
condition (OBC) Hamiltonian can be well approximated by the Bloch Hamiltonian
with periodic boundary conditions (PBCs), which is called the BBC in band
theory [45]. However, the eigenstates with the NHSE cannot be characterized by
the BBC because all of them are localized at the edge with OBCs in the finite-
size sysytem [46, 47, 48]. In other words, the eigenstates and eigenenergies
with the NHSE under different boundary conditions are completely different
from each other. To solve this conflict, the generalized Brillouin zone (GBZ)
[28, 49, 50] and biorthogonal BBC [15, 51] have been proposed. GBZ theory
takes the complex momentum and modifies the Bloch factor to obtain the OBC
eigenenergies [15, 11, 52, 14, 53]. Therefore, the corresponding eigenstates
are exponentially localized at the edges of the system. Biorthogonal BBC
theory uses both the left and right eigenstates to obtain the OBC
eigenenergies from the PBC Hamiltonian. Moreover, the topological invariants
of the system with the NHSE are also quite different. The winding number can
be defined from the complex energy spectrum, where the complex energy spectrum
under PBCs can form a loop and encircle all of the complex energy spectrum
under OBCs [37, 54].

The largest influence of the NHSE on the dynamic properties of the system is
reflected in the appearance of a directional (chiral) bulk flow and persistent
currents [55, 25, 23]. The wavefunctions inside the bulk have a nonreciprocal
behavior where their propagation has a specific direction, which can be
detected by experiments [56, 57]. In addition, a novel bulk dynamic signature
of the NHSE arises in the early stage of the time evolution, i.e., self
acceleration of the wavefunction [57, 58]. The complex energy spectrum of the
eigenenergies will make the density of the wavepackets exhibit gain or loss
for different momentums. The long time dynamics of the wavefunction are
dominated by the interference of Bloch modes with the largest imaginary part
of the energy, leading to a directional flow of excitation along the lattice
at a constant drift velocity in systems, and displaying the NHSE. However, the
imaginary parts of the eigenenergies can have more than one maximum (beyond
zero), where the flow of condensates will have a competitive behavior between
the maximum points of the eigenenergies.

In this work, we study the dynamic evolution of wavepackets with different
velocities in non-Hermitian lattices. We find that the wavepackets with the
maximum imaginary eigenenergies will be protected by the bulk until they reach
the boundary. We discuss these dynamic protected states in several models. The
paper is organized as follows. In Sec. 2, we introduce the simplest lattice
with the NHSE and discuss the influence of the finite size on the energy
spectrum and localization. In Sec. 3, we investigate the bulk protected states
under OBCs and give the analytical result. In Sec. 4, we take two widely used
models as examples to prove the universality in the non-Hermitian system. In
Sec. 5, we conclude with a summary of our main results and final remarks.

##  2  Finite size and the non-Hermitian skin effect

Figure 1: a) Complex energy spectrum of the HN model with different system
sizes and b) eigenenergies with PBCs. c) Log function of the density of the
eigenstates, and d) the scaling parameter $\kappa$ as a function of the system
size $L$ in Eq. (3). We take $J_{1}=1$ and $J_{2}/J_{1}=1.5$ for a), b) and
c).

In this section, we review the influence of the size of the system with NHSE.
The HN model is the simplest NHSE model, with

$H=\sum_{n}\left(J_{1}C^{\dagger}_{n+1}C_{n}+J_{2}C^{\dagger}_{n}C_{n+1}\right),$
(1)

where $C_{n}(C_{n}^{\dagger})$ is the annihilation (creation) operator at the
$n$ site index, $J_{1}$ and $J_{2}$ are the nearest-neighbor hopping energies
along different directions.

The eigenenergies of the non-Hermitian skin Hamiltonian are always different
with different boundary conditions. In particular, the periodic boundary
eigenenergies can encircle all open boundary eigenenergies As shown in Fig. 1
a), the HN model fundamentally gives purely real eigenenergies when the system
has a small size $L$. However, when the system size increases, imaginary parts
of the eigenenergies appear and approach the PBC result (blue signs in Fig.
1a)). Based on this, the topological invariants can easily be defined [24, 59]
by

$W=\frac{1}{2\pi
i}\oint_{C}\frac{d}{dk}\log\left[Det\left(\hat{H(k)}-E_{R}\hat{I}\right)\right],$
(2)

where $H(k)=\sum_{k}\left[(J_{1}+J_{2})\cos k-i(J_{1}-J_{2})\sin
k\right]C_{k}^{\dagger}C_{k}$, can be obtained by the Fourier transformation
of the real-space Hamiltonian (see Fig. 1b) ) and $E_{R}$ is the reference
energy. The eigenenergies are complex when the system size is large enough or
the system has PBCs. The winding number $W$ is dependent on the sign of
$J_{1}-J_{2}$. If $J_{1}=J_{2}$, the winding number is zero, all eigenenergies
are real, and the NHSE vanishes. However, $J_{1}\neq J_{2}$ will let the
topological invariants arise and all eigenstates localize at the side.

In the HN model, all wavefunctions are localized at the edge. The amplitudes
of the absolute value of the density will exponentially increase or descrease
(depending on the magnitude of $J_{1}$ and $J_{2}$) with the site index $n$,
as illustrated in Fig. 1 (c). The relationship between the density of the
first site and the last site can be described by

$|\psi(n)|^{2}=e^{\kappa n}|\psi(1)|^{2},$ (3)

where $\kappa$ is the scaling parameter and $n$ is the site index. In a
Hermitian lattice, the Bloch factor $|e^{ik}|=1$, however, GBZ theory states
that $|e^{ik}|=\beta\neq 1$, where all wavefunctions can have exponential
localization along the index. So the translation symmetry is broken in the
non-Hermitian lattice [60, 61] . As is illustrated in Fig. 1 (d), the scaling
paramters $\kappa$ is proportional to the inverse of the system size
$\kappa\sim L^{-1}$. Therefore, the localization of the density will vanish
when the system size goes to infinity.

##  3 Dynamic protected state

In this section, we will analyze the time evolution of the Gaussian
wavepackets and obtain the motion of the center mass in the initial stage of
evolution.

Figure 2: Time evolution of the Gaussian wavepackets with different initial
momenta in real space (first column) and momentum space (second column). The
parameters are $J_{2}/J_{1}$=$1.5$, $k_{0}$=$0$, and $-\pi/2$ or $\pi/2$ for
different rows.

In a Hermitian topological system, only the edge states are topologically
protected, and their propagation is chiral. However, all eigenstates of the
non-Hermitian Hamiltonian are localized, including the bulk and edge states.
The wavefunction is forbidden to propagate to one side and get enhanced to
propagate to the opposite side in the 1D case[62, 63].

When we consider the time evolution of the wavepackets of the system, we need
to expand the initial state with the eigenstates. The eigenenergies are
complex for a non-Hermitian system with PBCs or a sufficiently large size. The
wavepacket with momentum $k$ evolutes with time as follows,

$\ket{\psi(k,t)}=e^{-i{Re}(E)t}e^{{Im}(E)t}\ket{\psi(k,0)},$ (4)

where the real eigenenergies (${Re}(E)$) govern the motion of the wavepacket
and the imaginary eigenenergies describe the decay or growth rate of the
corresponding eigenstates. When ${Im}(E)<0$ or ${Im}(E)>0$, the corresponding
eigenstate will exponentially decay or grow with time. Therefore, the time
evolution of the wavefunction will have a peculiar dynamic behavior where all
wavefunctions in momentum space will become the same.

The time evolution of the wavefunction in momentum space can be obtained by
the Fourier transformation, with

$\psi(k,0)=\sum_{n=0}^{L-1}\psi(n,0)e^{-ikn}.$ (5)

As shown in Fig. 2, we place the Gaussian state with different initial
velocities $k_{0}$ in the center of the HN lattice . The wavepacket has
nonreciprocal behavior, and its center of mass only moves to one side. In
addition, if the initial momentum is not equal to the momentum $k_{m}=\pi/2$
where the imaginary eigenenergies have the maximum, then the momentum of the
wavepacket will approach $k_{m}$. Specifically, when we use $k_{0}=k_{m}$ as
the momentum of the initial state, the Gaussian state is protected until it
reaches the boundary which is unique compared to the Hermitian topological
lattice. Due to the maximum of the imaginary energies in the non-Hermitian
lattice with PBCs, only the wavefunction with momentum $k_{m}$ will be the
most enhanced along with the time evolution. Therefore, the propagation of the
wavepacket with OBCs in the bulk can still obey the rule of PBCs and survive
for a sufficiently long time until it reaches the boundary.

Figure 3: The mean momentums of the wavepackets changes with time at different
system sizes $L$ a), and the duration of the dynamic protected state as a
function of the system size b). Parameters are $J_{2}/J_{1}=1.5$ and
$k_{0}=0.4\pi$.

The duration of the wavepacket with OBCs linearly increases with the system
size, as illustrated in Fig. 3. When we set the momentum of the initial state
near $k_{m}$ ($k_{0}=0.4\pi$ ), the speed of wavepacket $k(t)$ will keep a
duration at $k_{m}$, and then maintain a long time unless reaching the
boundary. The initial state in momentum space can be described as

$\psi(k,0)=e^{(k-k_{0})^{2}/2\sigma^{2}},$ (6)

where $k_{0}$ is the center of the wavepacket in $k$ space and $\sigma$ is the
width. When the Gaussian state is in the bulk, we can use the PBC Hamiltonian
to explain the wavepacket dynamic properties , with

$\psi(k,t)=e^{-it\sum_{k}(J_{1}+J_{2})\cos k-(J_{1}-J_{2})\sin
k}e^{(k-k_{0})^{2}/2\sigma^{2}}.$ (7)

The first term is the oscillation frequency and can be ignored compared to the
exponential decay or gain due to the imaginary eigenenergies.

We consider a non-Hermitian Hamiltonian, denoted as
$\hat{H}=\sum_{k}E_{k}\ket{k}\bra{k}$, where $\ket{k}$ represents the right
eigenvector with eigenenergy $E_{k}$. An initial state $\ket{\psi(t=0)}$ is
governed by the evolution operator $U(t)$ associated with this Hamiltonian.
The state at time $t$, denoted as $\ket{\psi(k,t)}$, is determined by the
evolution operator $U(t)$ associated with the Hamiltonian. Mathematically, it
can be expressed as

$\displaystyle\ket{\psi(k,t)}$ $\displaystyle=$
$\displaystyle\hat{U}(t)\ket{\psi(k,0)}=e^{-i\hat{H}t}\ket{\psi(k,0)}$ (8)
$\displaystyle=$
$\displaystyle\prod_{k}e^{-iE_{k}t\ket{k}\bra{k}}\ket{\psi(k,0)}.$

Let us examine the variation in the mean momentum over time, given by
$\bra{\psi(k,t)}\hat{k}\ket{\psi(k,t)}$.

$\displaystyle k(t)$ $\displaystyle=$
$\displaystyle\bra{\psi(k,t)}\hat{k}\ket{\psi(k,t)}$ (9) $\displaystyle=$
$\displaystyle\bra{\psi(k,0)}\prod_{k_{2}}e^{iE^{*}_{k_{2}}t\ket{k_{2}}\bra{k_{2}}}\hat{k}\prod_{k_{1}}e^{-iE_{k_{1}}t\ket{k_{1}}\bra{k_{1}}}\ket{\psi(k,0)}$
$\displaystyle=$ $\displaystyle\int
dk^{\prime}\braket{\psi(k,0)}{k^{\prime}}e^{i(E^{*}_{k^{\prime}}-E_{k^{\prime}})t}k^{\prime}\braket{k^{\prime}}{\psi(k,0)}$
$\displaystyle=$ $\displaystyle\int
dk^{\prime}e^{2{Im}E_{k^{\prime}}t}k^{\prime}|\braket{k^{\prime}}{\psi(k,0)}|^{2},$

where we insert a unit operator $\int
dk^{\prime}\ket{k^{\prime}}\bra{k^{\prime}}$. In Hermitian systems,
$E_{k^{\prime}}$ is purely real, observing that $k(t)=k(0)$ remains unchanged
with time evolution. In non-Hermitian systems, when ${Im}(E_{k^{\prime}})>0$,
the contribution of $k^{\prime}$ becomes more significant, while
${Im}(E_{k^{\prime}})<0$, the contribution of $k^{\prime}$ diminishes. A
predictable consequence is that $k(t)$ eventually approaches $k_{m}$, which
corresponds to the maximum value of ${Im}(E_{k_{m}})$. Without loss of
generality, we can shift ${Im}(E_{k})$ such that the maximum point corresponds
to zero. We can obtain

$k(t)=\int
dk^{\prime}e^{2{Im}(\tilde{E}_{k^{\prime}})t}e^{2at}k^{\prime}|\braket{k^{\prime}}{\psi(k,0)}|^{2},$
(10)

where ${Im}(\tilde{E}_{k^{\prime}})={Im}(E_{k^{\prime}}-a)$, where $a$ is a
positive constant. We consider a simplified case where $a=0$, indicating that
only dissipative effects are present in the system. The contribution of all
values of $k$ except for the maximum point decreases. Therefore, as long as
$|\braket{k_{m}}{\psi(0)}|^{2}$ is nonzero, $\lim_{t\to\infty}k(t)$ remains
equal to $k_{m}$.

Therefore, the dynamic protected states are the wavepacket that will not
change with time during the propagation period. The momentum of the
wavepackets will approach to the $k_{m}$, which is the maximum point of the
imaginary parts of the eigenenergies of the non-Hermitian Hamiltonian with
NHSE.

##  4 Dynamic protected states of different non-Hermitian systems

In this section, we take two widely used theoretical models (the SSH model and
2D polariton system) with different symmetries as examples to check the
universality of the dynamic protected state in non-Hermitian systems and
discuss the influence of the skin effect. The time and energy dimensions are
normalized with $J_{1}t$ and $E/J_{1}$.

###  4.1 SSH model with chiral symmetry

Figure 4:  The real a) and the imaginary b) parts of the eigenenergies of the
SSH model with different spins. Eigenfunctions of different spins c), and the
density evolution of spin-up with the initial velocity $k_{0}=-0.5\pi$ d). The
other parameters are $t_{2}/t_{1}=1.2$ and $\gamma/t_{1}=4/3$.

The SSH model with different nearest-neighbor hopping energies in a 1D lattice
is widely studied to realize the NHSE [64, 65]. The Hamiltonian in momentum
$k$ space can be described as

$H=d_{x}\sigma_{x}+(d_{y}+i\frac{\gamma}{2})\sigma_{y},$ (11)

where $\gamma/2$ changes the hopping term in the unit cell with different
hopping strengths. $d_{x}=t_{1}+t_{2}\cos k$, $d_{y}=t_{2}\sin k$ and
$\sigma_{x,y}$ are the Pauli matrices. $t_{1}$ and $t_{2}$ are the hopping
between the different spins in the same unit cell and the hopping between the
nearest unit cells.

The real parts and imaginary parts of eigenvalues of Eq. (11) with the Fourier
transformation in momentum space are shown in Figs. 4 a)-b). All eigenstates
with OBCs are localized at the left side (shown in Fig. 4 c) ) when
$t_{1}<t_{2}$. The eigenvalues appear in pairs, with

$E_{\pm}(k)=\pm\sqrt{t_{1}^{2}+t_{2}^{2}+2t_{1}t_{2}\cos
k-\gamma^{2}/4+it_{2}\gamma\sin k}.$ (12)

The OBC eigenenergies are zero modes in the finite-size system, and GBZ theory
gives the non-Bloch factor
$|e^{ik}|=\sqrt{\left|\frac{t_{1}-\gamma/2}{t_{1}+\gamma/2}\right|}\neq 1$,
where $k$ becomes a complex number and makes the eigenstates localized.

The imaginary parts of eigenvalues of different spins have different maximum
points $k_{m}$ and $-k_{m}$ (red solid line and blue dashed line in Fig. 4 b))
with PBC. When we place the initial Gaussian state in the middle of the
lattice with OBC, the mean momentum of the wavepacket will have two peaks
$k_{m}$ and $-k_{m}$. However, the density of $|\psi(k,t)|^{2}$ are most
localized at $k_{m}$ rather than $-k_{m}$ as shown in Fig. 4 d). The above
results imply that the boundary condition plays a significant role in the non-
Hermitian system with skin effect. The eigenstates of the Hamiltonian with OBC
are localized at the left side (shown in Fig. 4 c) ), so dynamic behaviors of
a condensation in the bulk are also affected by the boundary condition and
propagate to the left side with momentum $k_{m}$.

###  4.2  Two-dimensional non-Hermitian system

Figure 5: Complex energy spectrum with PBCs and OBCs a), and real b) and
imaginary c) energy bands of PBC eigenenergies. The black lines are the OBC
energies, and the blue and green areas are filled with the PBC energies of the
upper and lower bands. Meanwhile, the green dots are the complex energies
along with the selected integration direction for defining the winding number
in (a). Parameters: $J_{2}/J_{1}$=1, $k_{p}=\pi/4$, $\gamma_{+}/J_{1}$=$0.1$,
$\gamma_{-}/J_{1}$=$0.8$, and $A/J_{1}$=$4$.
Figure 6: Time evolution of the Gaussian states in momentum space for
different spins (different rows) obtained with Hamiltonian Eq. (13).
Parameters: $J_{2}/J_{1}$=1, $k_{p}=\pi/4$, $\gamma_{+}/J_{1}$=$0.1$,
$\gamma_{-}/J_{1}$=$0.8$, and $A/J_{1}$=$4$.

In this subsection we will discuss dynamic behaviors and preotected states of
the higher-order NHSE in 2D polariton system. The previous work [54] have
revealed that all wavefunctions are localized at the corners in 2D tight-
binding polariton model. The effective Hamiltonian can be described as

$\displaystyle H_{m,n}$ $\displaystyle=$
$\displaystyle\sum_{m,n}\left[J_{1}\left(\hat{a}_{m,n}^{\dagger}\hat{a}_{m,n+1}+\hat{b}_{m,n}^{\dagger}\hat{b}_{m,n+1}+h.c.\right)+J_{2}\left(\hat{a}_{m+1,n}^{\dagger}\hat{a}_{m,n}+\hat{b}_{m+1,n}^{\dagger}\hat{b}_{m,n}+h.c.\right)\right]$
(13) $\displaystyle+$
$\displaystyle\sum_{m,n}\left[-i\gamma_{+}\left(\hat{a}_{m,n}^{\dagger}\hat{a}_{m,n}\right)-i\gamma_{-}\left(\hat{b}_{m,n}^{\dagger}\hat{b}_{m,n}\right)\right]$
$\displaystyle+$
$\displaystyle\sum_{m,n}\left[A\hat{a}_{m,n}^{\dagger}\hat{b}_{m,n}e^{-ik_{p}\left(m+n\right)}+A\hat{b}_{m,n}^{\dagger}\hat{a}_{m,n}e^{ik_{p}\left(m+n\right)}\right],$

where $\hat{a}_{m,n}$ ($\hat{a}^{\dagger}_{m,n}$) and $\hat{b}_{m,n}$
($\hat{b}^{\dagger}_{m,n}$) denote the annihilation (creation )of particles
with spin $\pm$ at site $(m,n)$. $J_{1,2}$ are the hopping energies in the $x$
and $y$ directions, $\gamma_{\pm}$ are the different decay (or gain) rates of
different spins, and $A$ describes the strength of the polarization splitting.

As shown in Fig. 5 a), the eigenenergies with OBC (black dots) are surrounded
by the PBC energies (green and blue areas). All wavefunctions are localized at
the off-side corners and the winding number is dependent on the selected
paths. Furthermore, only along with the path $k_{x}=k_{y}$ (red circles), all
OBC eigenenergies are encircled by the PBC eigenenergies in this path [54].
The eigenenergies of different spins are plotted with different colors (green
and blue areas) and they are separated from each other. The winding number for
each spin is opposite and their corresponding eigenstates are localized at the
opposite corners [54]. The system can be transformed into a two-band model
with PBC and the real and imaginary parts of the corresponding eigenenergies
as shown in Figs. 5 b) and c). The real parts are separated from each other,
but the imaginary parts have crossovers. The higer order NHSE is different
from the 1D case since the winding number is determined by the chosen
integrated path.

The imaginary parts of the eigenenergies of Eq. (13) have multiple maximum
points as illustrated in Fig. 5 c). To check the density distrubution in the
momentum space, we place two Gaussian states with zero initial speeds at the
middle of the lattice in Fig. 6. The density distribution of two spins will
gain velocities because of the self-accelaration effect [55]. However, all
eigenstates of the OBC Hamiltonian are localized at the off-side corners.
Here, the effective decays for different spins $\gamma_{\pm}$ have the crucial
effect for letting the density decay along with different directions are
different in the time evolution.

At the begining, the density distribution for different spins in momentum
space are localized at the different $k_{x}$ and $k_{y}$ as shown in Figs. 6
a1)-(a2) and b1)-b2). After the sufficient long time, the two-dimensional
Gaussian states in momentum space will remain stable as shown in Figs. 6
a3)-a5) and b3)-b5). In addition, the density distributions for different
spins are localized at the same $k_{x}$ and $k_{y}$ although the imaginary
parts of the eigenenergies of the PBC Hamiltonian are different and have
different maximum points. The eigenstates of the non-Hermitian system are
sensitive to the boundary condition and so do the dynamic protected states.

##  5 Discussion

In this work, we study the time evolution of the wavepackets in a non-
Hermitian lattice. The initial states are prepared as Gaussian states and have
different center-mass positions and initial velocities. The mean momentum of
Gaussian wavepacket will approach to the certain momentum, at which the
eigenenergies have the largest imaginary parts. This dynamic protected state
is stable in the bulk with OBC and will exponentially accelerate from the
initial speed. Furthermore, this state will keep the certain momentum for a
long time during the time evolution until the wavepackets reach the boundary.
We also take the 1D SSH and 2D polariton system models as examples whose
eigenenergies have more than one maximum points. The momentum of Gaussian
wavepackets will only have one momentum depending on which sides are the
eigenstates localized and also sensitive to the boundary condition.

On the one hand, our results provide a better understanding for the non-
Hermitian skin effect and self acceleration of the wavefunction. On the other
hand, the sensitivity and robustness of the evolution of bulk states with
different momentum in the non-Hermitian lattices are of great significant to
design new quantum devices, such as making the electronic sensor with a high-
level sensitivity and a strong robustness [66, 67], advancing non-Hermitian
topological magnonics and photonics to deal with the engineering of
dissipation and/or gain for non-Hermitian topological phases [68, 69, 70] ,
among other [71, 72].

##  6 Acknowledgements

This work is supported by the Fundamental Research Funds for the Central
Universities ( Grant No. JUSRP123027 ), National Natural Science Foundation of
China ( Grant No. 12264061 ), the Zhejiang Provincial Nature Science
Foundation under ( Grant No. LQ22A040006 ), and the Science Foundation of
Guizhou Science and Technology Department ( Grant number QKHJZ[2021]033 ).

##  7 Competing Interests

The authors declare no competing interests.

##  8 Author Contribution Statement

Lei Chen and Zhen-Xia Niu carried out the analytical calculation. Xingran Xu
did the numerical simulation and supervised the project.

##  9 Data Availability Statement

The datasets used and analysed during the current study are available from the
corresponding author on reasonable request.

## References

  *   * Carmichael [1993] Carmichael, H.J.: Quantum trajectory theory for cascaded open systems. Phys. Rev. Lett. 70(15), 2273 (1993) https://doi.org/10.1103/physrevlett.70.2273
  * Bender and Boettcher [1998] Bender, C.M., Boettcher, S.: Real spectra in non-hermitian hamiltonians having pt symmetry. Phys. Rev. Lett. 80, 5243–5246 (1998) https://doi.org/10.1103/PhysRevLett.80.5243
  * Bender et al. [1999] Bender, C.M., Boettcher, S., Meisinger, P.N.: Pt-symmetric quantum mechanics. J. Math. Phys. 40(5), 2201–2229 (1999) https://doi.org/10.1063/1.532860
  * Bender et al. [2002] Bender, C.M., Brody, D.C., Jones, H.F.: Complex extension of quantum mechanics. Phys. Rev. Lett. 89(27), 270401 (2002) https://doi.org/%****␣sn-article.bbl␣Line␣100␣****10.1103/physrevlett.89.270401
  * Ashida et al. [2020] Ashida, Y., Gong, Z., Ueda, M.: Non-hermitian physics. Adv. Phys. 69(3), 249–435 (2020) https://doi.org/10.1080/00018732.2021.1876991
  * Bender [2007] Bender, C.M.: Making sense of non-hermitian hamiltonians. Rep. Prog. Phys. 70(6), 947 (2007) https://doi.org/10.1088/0034-4885/70/6/r03
  * El-Ganainy et al. [2018] El-Ganainy, R., Makris, K.G., Khajavikhan, M., Musslimani, Z.H., Rotter, S., Christodoulides, D.N.: Non-hermitian physics and pt symmetry. Nature Physics 14(1), 11–19 (2018) https://doi.org/10.1038/nphys4323
  * Özdemir et al. [2019] Özdemir, Ş.K., Rotter, S., Nori, F., Yang, L.: Parity–time symmetry and exceptional points in photonics. Nature materials 18(8), 783–798 (2019) https://doi.org/10.1038/s41563-019-0304-9
  * Miri and Alu [2019] Miri, M.-A., Alu, A.: Exceptional points in optics and photonics. Science 363(6422), 7709 (2019) https://doi.org/10.1126/science.aar7709
  * Fruchart et al. [2021] Fruchart, M., Hanai, R., Littlewood, P.B., Vitelli, V.: Non-reciprocal phase transitions. Nature 592(7854), 363–369 (2021) https://doi.org/10.1038/s41586-021-03375-9
  * Kawabata et al. [2019] Kawabata, K., Shiozaki, K., Ueda, M., Sato, M.: Symmetry and topology in non-hermitian physics. Phys. Rev. X 9, 041015 (2019) https://doi.org/10.1103/PhysRevX.9.041015
  * Zhou and Lee [2019] Zhou, H., Lee, J.Y.: Periodic table for topological bands with non-hermitian symmetries. Phys. Rev. B 99, 235112 (2019) https://doi.org/10.1103/PhysRevB.99.235112
  * Yokomizo and Murakami [2019] Yokomizo, K., Murakami, S.: Non-bloch band theory of non-hermitian systems. Phys. Rev. Lett. 123, 066404 (2019) https://doi.org/10.1103/PhysRevLett.123.066404
  * Yao and Wang [2018] Yao, S., Wang, Z.: Edge states and topological invariants of non-hermitian systems. Phys. Rev. Lett. 121, 086803 (2018) https://doi.org/%****␣sn-article.bbl␣Line␣250␣****10.1103/PhysRevLett.121.086803
  * Kunst et al. [2018] Kunst, F.K., Edvardsson, E., Budich, J.C., Bergholtz, E.J.: Biorthogonal bulk-boundary correspondence in non-hermitian systems. Phys. Rev. Lett. 121, 026808 (2018) https://doi.org/10.1103/PhysRevLett.121.026808
  * Lee and Thomale [2019] Lee, C.H., Thomale, R.: Anatomy of skin modes and topology in non-hermitian systems. Phys. Rev. B 99, 201103 (2019) https://doi.org/10.1103/PhysRevB.99.201103
  * McDonald et al. [2018] McDonald, A., Pereg-Barnea, T., Clerk, A.: Phase-dependent chiral transport and effective non-hermitian dynamics in a bosonic kitaev-majorana chain. Phys. Rev. X 8(4), 041031 (2018) https://doi.org/10.1103/physrevx.8.041031
  * Hatano and Nelson [1996] Hatano, N., Nelson, D.R.: Localization transitions in non-hermitian quantum mechanics. Phys. Rev. Lett. 77, 570–573 (1996) https://doi.org/10.1103/PhysRevLett.77.570
  * Hatano and Nelson [1997] Hatano, N., Nelson, D.R.: Vortex pinning and non-hermitian quantum mechanics. Phys. Rev. B 56, 8651–8673 (1997) https://doi.org/10.1103/PhysRevB.56.8651
  * Longhi et al. [2015a] Longhi, S., Gatti, D., Valle, G.D.: Robust light transport in non-hermitian photonic lattices. Scientific reports 5(1), 13376 (2015) https://doi.org/10.1038/srep13376
  * Longhi et al. [2015b] Longhi, S., Gatti, D., Della Valle, G.: Non-hermitian transparency and one-way transport in low-dimensional lattices by an imaginary gauge field. Phys. Rev. B 92(9), 094204 (2015) https://doi.org/10.1103/physrevb.92.094204
  * Martinez Alvarez et al. [2018] Martinez Alvarez, V., Barrios Vargas, J., Foa Torres, L.: Non-hermitian robust edge states in one dimension: Anomalous localization and eigenspace condensation at exceptional points. Phys. Rev. B 97(12), 121401 (2018) https://doi.org/10.1103/physrevb.97.121401
  * Zhang et al. [2020] Zhang, K., Yang, Z., Fang, C.: Correspondence between winding numbers and skin modes in non-hermitian systems. Phys. Rev. Lett. 125, 126402 (2020) https://doi.org/10.1103/PhysRevLett.125.126402
  * Okuma et al. [2020] Okuma, N., Kawabata, K., Shiozaki, K., Sato, M.: Topological origin of non-hermitian skin effects. Phys. Rev. Lett. 124, 086801 (2020) https://doi.org/10.1103/PhysRevLett.124.086801
  * Longhi [2019] Longhi, S.: Probing non-hermitian skin effect and non-bloch phase transitions. Phys. Rev. Research 1, 023013 (2019) https://doi.org/10.1103/PhysRevResearch.1.023013
  * Song et al. [2019] Song, F., Yao, S., Wang, Z.: Non-hermitian skin effect and chiral damping in open quantum systems. Phys. Rev. Lett. 123, 170401 (2019) https://doi.org/10.1103/PhysRevLett.123.170401
  * Yokomizo and Murakami [2021] Yokomizo, K., Murakami, S.: Scaling rule for the critical non-hermitian skin effect. Phys. Rev. B 104(16), 165117 (2021) https://doi.org/10.1103/physrevb.104.165117
  * Yang et al. [2020] Yang, Z., Zhang, K., Fang, C., Hu, J.: Non-hermitian bulk-boundary correspondence and auxiliary generalized brillouin zone theory. Phys. Rev. Lett. 125, 226402 (2020) https://doi.org/10.1103/PhysRevLett.125.226402
  * Xu et al. [2021] Xu, K., Zhang, X., Luo, K., Yu, R., Li, D., Zhang, H.: Coexistence of topological edge states and skin effects in the non-hermitian su-schrieffer-heeger model with long-range nonreciprocal hopping in topoelectric realizations. Phys. Rev. B 103(12), 125411 (2021) https://doi.org/10.1103/physrevb.103.125411
  * Claes and Hughes [2021] Claes, J., Hughes, T.L.: Skin effect and winding number in disordered non-hermitian systems. Phys. Rev. B 103(14), 140201 (2021) https://doi.org/10.1103/physrevb.103.l140201
  * Okuma and Sato [2021] Okuma, N., Sato, M.: Non-hermitian skin effects in hermitian correlated or disordered systems: quantities sensitive or insensitive to boundary effects and pseudo-quantum-number. Phys. Rev. Lett. 126(17), 176601 (2021) https://doi.org/10.1103/physrevlett.126.176601
  * Yao et al. [2018] Yao, S., Song, F., Wang, Z.: Non-hermitian chern bands. Phys. Rev. Lett. 121, 136802 (2018) https://doi.org/10.1103/PhysRevLett.121.136802
  * Fu et al. [2021] Fu, Y., Hu, J., Wan, S.: Non-hermitian second-order skin and topological modes. Phys. Rev. B 103, 045420 (2021) https://doi.org/10.1103/PhysRevB.103.045420
  * Kim and Park [2021] Kim, K.-M., Park, M.J.: Disorder-driven phase transition in the second-order non-hermitian skin effect. Physical Review B 104(12), 121101 (2021) https://doi.org/10.1103/physrevb.104.l121101
  * Ghorashi et al. [2021] Ghorashi, S.A.A., Li, T., Sato, M., Hughes, T.L.: Non-hermitian higher-order dirac semimetals. Phys. Rev. B 104(16), 161116 (2021) https://doi.org/10.1103/physrevb.104.l161116
  * Okuma and Sato [2019] Okuma, N., Sato, M.: Topological phase transition driven by infinitesimal instability: Majorana fermions in non-hermitian spintronics. Phys. Rev. Lett. 123, 097701 (2019) https://doi.org/10.1103/PhysRevLett.123.097701
  * Xu et al. [2021] Xu, X., Xu, H., Mandal, S., Banerjee, R., Ghosh, S., Liew, T.C.H.: Interaction-induced double-sided skin effect in an exciton-polariton system. Phys. Rev. B 103, 235306 (2021) https://doi.org/10.1103/PhysRevB.103.235306
  * Lin et al. [2021] Lin, Z., Ding, L., Ke, S., Li, X.: Steering non-hermitian skin modes by synthetic gauge fields in optical ring resonators. Opt. Lett. 46(15), 3512 (2021) https://doi.org/%****␣sn-article.bbl␣Line␣625␣****10.1364/OL.431904
  * Sun et al. [2021] Sun, X.-Q., Zhu, P., Hughes, T.L.: Geometric response and disclination-induced skin effects in non-hermitian systems. Phys. Rev. Lett. 127(6), 066401 (2021) https://doi.org/10.1103/physrevlett.127.066401
  * Yuce [2021] Yuce, C.: Nonlinear non-hermitian skin effect. Phys. Lett. A 408, 127484 (2021) https://doi.org/10.1016/j.physleta.2021.127484
  * Zhang et al. [2021a] Zhang, K., Yang, Z., Fang, C.: Universal non-Hermitian skin effect in two and higher dimensions (2021)
  * Zhang et al. [2021b] Zhang, L., Yang, Y., Ge, Y., Guan, Y.-j., Chen, Q., Yan, Q., Chen, F., Xi, R., Li, Y., Jia, D., Yuan, S.-q., Sun, H.-x., Chen, H., Zhang, B.: Acoustic non-Hermitian skin effect from twisted winding topology (2021)
  * Liang et al. [2022] Liang, Q., Xie, D., Dong, Z., Li, H., Li, H., Gadway, B., Yi, W., Yan, B.: Dynamic signatures of non-hermitian skin effect and topology in ultracold atoms. Phys. Rev. Lett. 129(7), 070401 (2022) https://doi.org/10.1103/physrevlett.129.070401
  * Shang et al. [2022] Shang, C., Liu, S., Shao, R., Han, P., Zang, X., Zhang, X., Salama, K.N., Gao, W., Lee, C.H., Thomale, R., Manchon, A., Zhang, S., Cui, T.J., Schwingenschlögl, U.: Experimental identification of the second-order non-hermitian skin effect with physics-graph-informed machine learning. Adv. Sci. 9(36), 2202922 (2022) https://doi.org/10.1002/advs.202202922
  * Alase et al. [2016] Alase, A., Cobanera, E., Ortiz, G., Viola, L.: Exact solution of quadratic fermionic hamiltonians for arbitrary boundary conditions. Phys. Rev. Lett. 117, 076804 (2016) https://doi.org/10.1103/PhysRevLett.117.076804
  * Helbig et al. [2020] Helbig, T., Hofmann, T., Imhof, S., Abdelghany, M., Kiessling, T., Molenkamp, L.W., Lee, C.H., Szameit, A., Greiter, M., Thomale, R.: Generalized bulk–boundary correspondence in non-hermitian topolectrical circuits. Nat. Phys. 16(7), 747–750 (2020) https://doi.org/%****␣sn-article.bbl␣Line␣775␣****10.1038/s41567-020-0922-9
  * Liu et al. [2020] Liu, C.-H., Zhang, K., Yang, Z., Chen, S.: Helical damping and dynamical critical skin effect in open quantum systems. Phys. Rev. Research 2, 043167 (2020) https://doi.org/10.1103/PhysRevResearch.2.043167
  * Lee et al. [2019] Lee, C.H., Li, L., Gong, J.: Hybrid higher-order skin-topological modes in nonreciprocal systems. Phys. Rev. Lett. 123, 016805 (2019) https://doi.org/10.1103/PhysRevLett.123.016805
  * Xiao et al. [2020] Xiao, L., Deng, T., Wang, K., Zhu, G., Wang, Z., Yi, W., Xue, P.: Non-hermitian bulk–boundary correspondence in quantum dynamics. Nat. Phys. 16(7), 761–766 (2020) https://doi.org/10.1038/s41567-020-0836-6
  * Imura and Takane [2019] Imura, K.-I., Takane, Y.: Generalized bulk-edge correspondence for non-hermitian topological systems. Phys. Rev. B 100, 165430 (2019) https://doi.org/10.1103/PhysRevB.100.165430
  * Xu et al. [2020] Xu, X., Liu, H., Zhang, Z., Liang, Z.: The non-hermitian geometrical property of 1d lieb lattice under majorana’s stellar representation. J. Phys-Condens. Mat. 32(42), 425402 (2020) https://doi.org/10.1088/1361-648x/ab9fd4
  * Gong et al. [2018] Gong, Z., Ashida, Y., Kawabata, K., Takasan, K., Higashikawa, S., Ueda, M.: Topological phases of non-hermitian systems. Phys. Rev. X 8, 031079 (2018) https://doi.org/10.1103/PhysRevX.8.031079
  * Lieu [2018] Lieu, S.: Topological phases in the non-hermitian su-schrieffer-heeger model. Phys. Rev. B 97, 045106 (2018) https://doi.org/10.1103/PhysRevB.97.045106
  * Xu et al. [2022] Xu, X., Bao, R., Liew, T.C.H.: Non-hermitian topological exciton-polariton corner modes. Phys. Rev. B 106, 201302 (2022) https://doi.org/10.1103/PhysRevB.106.L201302
  * Longhi [2022] Longhi, S.: Non-hermitian skin effect and self-acceleration. Phys. Rev. B 105, 245143 (2022) https://doi.org/10.1103/PhysRevB.105.245143
  * Weidemann et al. [2022] Weidemann, S., Kremer, M., Longhi, S., Szameit, A.: Topological triple phase transition in non-hermitian floquet quasicrystals. Nature 601(7893), 354–359 (2022) https://doi.org/10.1038/s41586-021-04253-0
  * Lin et al. [2022] Lin, Q., Li, T., Xiao, L., Wang, K., Yi, W., Xue, P.: Topological phase transitions and mobility edges in non-hermitian quasicrystals. Phys. Rev. Lett. 129, 113601 (2022) https://doi.org/10.1103/PhysRevLett.129.113601
  * Parto et al. [2023] Parto, M., Leefmans, C., Williams, J., Nori, F., Marandi, A.: Non-abelian effects in dissipative photonic topological lattices. Nature Communications 14(1) (2023) https://doi.org/10.1038/s41467-023-37065-z
  * Zhang et al. [2022] Zhang, X., Zhang, T., Lu, M.-H., Chen, Y.-F.: A review on non-hermitian skin effect. Advances in Physics: X 7(1) (2022) https://doi.org/10.1080/23746149.2022.2109431
  * Wang et al. [2023] Wang, Y.-C., Jen, H.H., You, J.-S.: Scaling laws for non-hermitian skin effect with long-range couplings. Phys. Rev. B 108, 085418 (2023) https://doi.org/%****␣sn-article.bbl␣Line␣1000␣****10.1103/PhysRevB.108.085418
  * Qin et al. [2023] Qin, F., Ma, Y., Shen, R., Lee, C.H.: Universal competitive spectral scaling from the critical non-hermitian skin effect. Phys. Rev. B 107, 155430 (2023) https://doi.org/10.1103/PhysRevB.107.155430
  * Hu et al. [2023] Hu, Y.-M.R., Ostrovskaya, E.A., Estrecho, E.: Wave-packet dynamics in a non-hermitian exciton-polariton system. Phys. Rev. B 108, 115404 (2023) https://doi.org/10.1103/PhysRevB.108.115404
  * Malzard et al. [2018] Malzard, S., Cancellieri, E., Schomerus, H.: Topological dynamics and excitations in lasers and condensates with saturable gain or loss. Opt. Express 26(17), 22506 (2018) https://doi.org/10.1364/OE.26.022506
  * Liu et al. [2020] Liu, J.S., Han, Y.Z., Liu, C.S.: A new way to construct topological invariants of non-hermitian systems with the non-hermitian skin effect. Chin. Phys. B 29(1), 010302 (2020) https://doi.org/10.1088/1674-1056/ab5937
  * Fu et al. [2020] Fu, Z., Fu, N., Zhang, H., Wang, Z., Zhao, D., Ke, S.: Extended SSH model in non-hermitian waveguides with alternating real and imaginary couplings. Appl. Sci. 10(10), 3425 (2020) https://doi.org/10.3390/app10103425
  * Yuan et al. [2023] Yuan, H., Zhang, W., Zhou, Z., Wang, W., Pan, N., Feng, Y., Sun, H., Zhang, X.: Non-hermitian topolectrical circuit sensor with high sensitivity. Advanced Science 10(19) (2023) https://doi.org/10.1002/advs.202301128
  * Zhang et al. [2024] Zhang, X., Wu, C., Yan, M., Liu, N., Wang, Z., Chen, G.: Observation of continuum landau modes in non-hermitian electric circuits. Nature Communications 15(1) (2024) https://doi.org/10.1038/s41467-024-46122-0
  * Yu et al. [2024] Yu, T., Zou, J., Zeng, B., Rao, J.W., Xia, K.: Non-hermitian topological magnonics. Physics Reports 1062, 1–86 (2024) https://doi.org/10.1016/j.physrep.2024.01.006
  * Yan et al. [2023] Yan, Q., Zhao, B., Zhou, R., Ma, R., Lyu, Q., Chu, S., Hu, X., Gong, Q.: Advances and applications on non-hermitian topological photonics. Nanophotonics 12(13), 2247–2271 (2023) https://doi.org/10.1515/nanoph-2022-0775
  * Li et al. [2023] Li, A., Wei, H., Cotrufo, M., Chen, W., Mann, S., Ni, X., Xu, B., Chen, J., Wang, J., Fan, S., Qiu, C.-W., Alù, A., Chen, L.: Exceptional points and non-hermitian photonics at the nanoscale. Nature Nanotechnology 18(7), 706–720 (2023) https://doi.org/10.1038/s41565-023-01408-0
  * Zhao et al. [2019] Zhao, H., Qiao, X., Wu, T., Midya, B., Longhi, S., Feng, L.: Non-hermitian topological light steering. Science 365(6458), 1163–1166 (2019) https://doi.org/10.1126/science.aay1064
  * Gu et al. [2022] Gu, Z., Gao, H., Xue, H., Li, J., Su, Z., Zhu, J.: Transient non-hermitian skin effect. Nature Communications 13(1) (2022) https://doi.org/10.1038/s41467-022-35448-2