1 INTRODUCTION

Modulational Instability and Dynamical Growth Blockade in the

Nonlinear Hatano-Nelson Model
Stefano Longhi
Prof. Stefano Longhi
Dipartimento di Fisica, Politecnico di Milano,Piazza L. da Vinci 32, I-20133 Milano, Italy
& IFISC (UIB-CSIC), Instituto de Fisica Interdisciplinar y Sistemas Complejos, E-07122 Palma de Mal-
lorca, Spain
Email Address: [email protected]
Keywords: non-Hermitian skin effect, modulation instability, discrete nonlinear Schrödinger equation
The Hatano-Nelson model is a cornerstone of non-Hermitian physics, describing asymmetric hopping dynamics on a one-dimensional
lattice, which gives rise to fascinating phenomena such as directional transport, non-Hermitian topology, and the non-Hermitian skin
arXiv:2501.01226v1 [nlin.PS] 2 Jan 2025

effect. It has been widely studied in both classical and quantum systems, with applications in condensed matter physics, photonics,
and cold atomic gases. Recently, nonlinear extensions of the Hatano-Nelson model have opened a new avenue for exploring the inter-
play between nonlinearity and non-Hermitian effects. Particularly, in lattices with open boundary conditions, nonlinear skin modes
and solitons, localized at the edge or within the bulk of the lattice, have been predicted. In this work, we examine the nonlinear ex-
tension of the Hatano-Nelson model with periodic boundary conditions and reveal a novel dynamical phenomenon arising from the
modulational instability of nonlinear plane waves: growth blockade. This phenomenon is characterized by the abrupt halt of norm
growth, as observed in the linear Hatano-Nelson model, and can be interpreted as a stopping of convective motion arising from self-
induced disorder in the lattice.

1 Introduction
The Hatano-Nelson (HN) model is a non-Hermitian quantum mechanical model introduced by Naomichi
Hatano and David Nelson in the 1990s [1, 2, 3]. Originally studied in the context of Anderson localiza-
tion in random media under the influence of an imaginary gauge field [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15], the HN model has become fundamental to understanding non-Hermitian topology [16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]–a rapidly expanding field in non-Hermitian physics that in-
vestigates topological phases of matter in systems where energy is not conserved [30]. The HN model
exhibits a nontrivial point-gap topology and a strong sensitivity of the energy spectrum on the bound-
ary conditions, resulting in the non-Hermitian skin effect [18, 21, 26, 27, 28, 30, 31, 32]– a unique phe-
nomenon where the bulk states of a system become localized at the boundaries, defying conventional
bulk-boundary correspondence seen in Hermitian systems. Introducing many-body effects and nonlin-
earity into the HN model leads to a variety of novel effects, which are attracting a great attention re-
cently [33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56]. In par-
ticular, nonlinear extension of the original HN model incorporating a Kerr-type (focusing or defocusing)
nonlinearity, has been introduced in recent works [51, 52, 53, 54, 55, 56], leading to the predictions and
very recent observation of nonlinear skin effect and skin solitons in systems with open boundary con-
ditions. The Hermitian limit of the nonlinear HN model is the discrete nonlinear Schrödinger Equation
(DNLSE), which is a fundamental model in nonlinear dynamics and mathematical physics describing the
evolution of wave-like excitations on discrete lattices with nonlinearity in a variety of physical systems
ranging from optics to Bose-Einstein condensates and lattice systems (see e.g. [57, 58, 59, 60, 61, 62, 63]
and references therein). The DNLSE can be regarded as a mean field model of the many-body Bose-
Hubbard model [64]. A fundamental phenomenon in nonlinear wave dynamics is the modulational in-
stability (MI) [58, 59, 63, 65]. It refers to the exponential growth of small perturbations in a plane wave
(spatially-homogeneous) solution, leading rather generally to the formation of spatial localized structures
such as solitons, breathers or localized patterns, or complex spatio-temporal dynamics.
In this work we investigate the phenomenon of modulational instability in the nonlinear Hatano-Nelson
model under periodic boundary conditions and unveil a novel dynamical phenomenon arising from the
modulational instability of nonlinear plane waves: growth blockade. This phenomenon is characterized
by the abrupt halt of intensity growth, as observed in the linear Hatano-Nelson model, and can be ex-
plained as a stopping of convective motion arising from self-induced disorder in the lattice. Our results

1
 2 MODULATIONAL INSTABILITY IN THE NONLINEAR HATANO-NELSON MODEL

indicate that combination of nonlinearity and non-Hermitian characteristics opens up new avenues for
studying wave propagation, localization, and instabilities in both classical and quantum non-Hermitian
systems.

2 Modulational instability in the nonlinear Hatano-Nelson model

Before considering the nonlinear HN model, let us first briefly review the concept of modulational insta-
bility in the DNLSE [65], to which the nonlinear HN reduces in the Hermitian limit. The DNLSE on a
one-dimensional lattice reads [58, 59, 60]
dψn
i = κψn+1 + κψn−1 + χ|ψ|2 ψn (1)
dt
for the excitation amplitude ψn at the n-th lattice site, where κ is the hopping amplitude and χ is the
strength of the cubic (Kerr) nonlinearity. We consider a finite lattice comprising N sites with periodic
boundary conditions ψn+N = ψn . Under such periodic boundaries, the DNLSE admits of a family of
nonlinear plane waves, which are the nonlinear extensions of the usual Bloch eigenmodes in the linear
limit χ = 0. Such plane waves are given by
ψn (t) = A exp[iqn − iω(q)t − iθ(t)] ≡ ϕn (t) (2)
where q is the wave number, which is quantized as q = ql = 2l(π/N ) (l = 0, 1, 2, .., N − 1),
ω(q) = 2κ cos(q) (3)
is the dispersion relation of the linear lattice, and
θ(t) = χ|A|2 t
is the additional nonlinear-induced phase shift. The stability of the plane wave solution against small
perturbations can be performed by a standard linear stability analysis [60, 65], looking for a solution to
Eq.(1) of the form
ψn (t) = [A + un (t)] exp[iqn − iω(q)t − iθ(t)] (4)
where un (t) is a small perturbation. Assuming |un (t)| ≪ |A|, substitution of Eq.(4) into Eq.(1) yields the
following linearized equation for the perturbation amplitudes un (t)
dun
i = κ exp(iq)un+1 + κ exp(−iq)un−1 + (−ω(q) + χ|A|2 )un + χA2 u∗n . (5)
dt
The most general solution to Eq.(5) is an arbitrary superposition of waves
un (t) = X exp[iQn − iλ(Q)t] + Y ∗ exp[−iQn + iλ∗ (Q)t] (6)
where Q is the spatial wave number of the perturbation and (X, Y ), λ(Q) are the eigenvectors and corre-
sponding eigenvalues (Lypaunov exponents) of the 2 × 2 linearized matrix
 
2κ cos(q + Q) − ω(q) + χ|A|2 χA2
M (Q) = . (7)
−χA∗2 −2κ cos(q − Q) + ω(q) − χ|A|2
This yields q
2
λ(Q) = −2κ sin q sin Q ± 4κ cos q sin2 (Q/2) − χ|A|2 − χ2 |A|4 . (8)
An instability arises whenever Im(λ(Q)) > 0, with a MI gain given by g(Q) = Im(λ(Q)). Clearly, for
a focusing nonlinearity χ > 0 the plane waves with wave number q such that cos q > 0 undergo MI,
while those with wave number q such that cos q < 0 are linearly stable. The reverse occurs in the defo-
cusing regime χ < 0. The main effects of MI is to destabilize the homogeneous intensity pattern of the

2
 2 MODULATIONAL INSTABILITY IN THE NONLINEAR HATANO-NELSON MODEL

Figure 1: Onset of modulational instability in the DNLSE. (a) Numerically-computed temporal evolution of the amplitudes
|ψn (t)| on a pseudocolor map in a lattice comprising N = 12 sites with periodic boundary conditions for hopping rate
κ = 1 and nonlinear parameter χ = 1. The initial condition is the plane wave ψn (t) = A exp(iqn) with amplitude A = 1
and spatial wave number q = 0, perturbed by adding a small random noise. Note that the plane wave is modulationally
unstable. (b) Same as (a), but for χ = −1. In this case the modulational instability is absent and the plane wave is stable.
(c) Same as (a) but for q = 2π/3. The plane wave is stable. (d) Same as (c) but for χ = −1. The plane wave is unstable.

plane wave leading rather generally to irregular (chaotic-like) dynamics [65] or the creation of more reg-
ular localized solutions such as train of solitons and breathing modes [60]. A few examples of the onset
of MI in the DNLSE, leading to irregular spatio-temporal dynamics, are shown in Fig.1. The results are
obtained by numerical integration of Eq.(1) using an accurate variable-step fourth-order Runge-Kutta
method for a few different initial conditions ψn (0) in a lattice comprising N = 12 sites with periodic
boundary conditions. A focusing (χ = 1) or defocusing (χ = −1) nonlinearity has been considered. As
an initial condition, we assumed a nonlinear plane wave with wave number q and amplitude A, slightly
perturbed by a small random noise, namely

ψn (0) = A(1 + rn ) exp(iqn)

where rn is a complex random variable whose real and imaginary parts are uniformly distributed in the
interval (−δ, δ) (typically we set δ = 10−3 ). The results clearly show the different MI regions depending
on the sign of nonlinearity χ, and the existence of stable plane wave solutions.
Let us now consider the nonlinear extension of the Hatano-Nelson model incorporating Kerr-type non-
linearity, as introduced in recent studies [51, 52, 54, 55], which is basically obtained from the DNLSE (1)
by making the left/right hopping amplitudes asymmetric via an imaginary gauge field h [1]. The nonlin-
ear HN model reads
dψn
i = κ exp(h)ψn+1 + κ exp(−h)ψn−1 + χ|ψ|2 ψn (9)
dt
where h ≥ 0 is the imaginary gauge field. Clearly, in the Hermitian limit h = 0 the above equation re-
duces to the DNLSE (1), whereas in the linear limit χ = 0 it reproduces the disorder-free HN model
[1, 2, 16]. The linear HN model displays a nontrivial point-gap topology, characterized by a topologi-
cal winding number [16], and the non-Hermitian skin effect (NHSE) under open or semi-infinite bound-
ary conditions, i.e. the localization of bulk eigenstates at the lattice edge [18, 27]. The central result in
the band theory of NH systems is that the energy spectra σ(HP BC ), σ(HOBC ) and σ(HSIBC ) of the sys-
tem Hamiltonian H, corresponding to periodic boundary conditions (PBC), open boundary condition

3
 2 MODULATIONAL INSTABILITY IN THE NONLINEAR HATANO-NELSON MODEL

(OBC), and semi-infinite boundary conditions (SIBC), are distinct, which implies the emergence of the
NHSE, topological NH edge states and the need for a non-Bloch band theory [27, 28, 29]. Specifically,
for the linear HN model: (i) σ(HP BC ) is a closed loop (an ellipse) in complex energy plane described by
the Bloch Hamiltonian
H(q) = κ exp(iq + h) + κ exp(−iq − h)
(−π ≤ q < π is the Bloch wave number); (ii) σ(HOBC ) is the segment (−2κ, 2κ) on theSreal energy axis,
which is always topological trivial in terms of a point gap; (iii) σ(HSIBC ) = σ(HP BC ) B, where B is
the interior of the PBC energy spectrum loop such that for E ∈ B the winding number W (E), defined
by [27] Z π
1 d
W (E) = dq log det {H(q) − E}
2πi −π dq
is non vanishing. If W (E) < 0, then E is an eigenvalue of HSIBC of multiplicity |W (E)|, and the corre-
sponding (right) eigenvectors are exponentially localized at the left edge. Such a result provides a bulk-
boundary correspondence for the linear HN model, relating the appearance of skin edge states in a semi-
infinite lattice to the topology of the PBC energy spectrum [16, 27]. In the nonlinear extension of the
HN model, the system dynamics are no longer governed by linear superposition principles but instead
exhibit nonlinear effects like feedback and self-interactions. In particular, it is not possible anymore to
classify edge (skin) states in terms of a topological winding number since under open boundary condi-
tions the inclusion of Kerr-type nonlinearity alters the NHSE-induced localization leading to new bound-
ary and bulk phenomena investigated in some recent works [51, 52, 54, 55]. While in the linear model
the NHSE results in the localization of states at one edge as established by the topological winding num-
ber W (E) of the PBC energy spectrum, nonlinearity can modify or redistribute the eigenstate localiza-
tion [54]. Depending on the strength and nature of the nonlinearity, it can either enhance or suppress
the NHSE due to the interplay between nonlinear self-trapping and linear convective drift, potentially
leading to the emergence of new localized states in the bulk, such as skin solitons [55], which clearly demon-
strate a breakdown of the bulk-edge correspondence observed in the linear HN model. Interestingly, un-
der open boundary conditions the nonlinear HN model can be mapped onto the Hermitian DNLSE with
spatially-inhomogeneous nonlinearity via the non-unitary gauge transformation ψn (t) → ψn (t) exp(−hn)
[55].
Unlike recent studies [51, 52, 53, 54, 55], here we focus our attention to a fully different regime, corre-
sponding to a lattice with periodic boundary conditions where such a mapping is prevented. As we will
show below, in this case the nonlinearity drives the system into an Anderson-type disordered phase trig-
gered by the modulational instability of nonlinear plane waves. In fact, under periodic boundary condi-
tions, i.e. ψn+N (t) = ψn (t), Eq.(9) admits of the following family of exact nonlinear plane wave solutions
ψn (t) = A exp[iqn − iω(q)t − iθ(t)] ≡ ϕn (t) (10)

where q is the wave number, which is quantized as q = ql = 2l(π/N ) (l = 0, 1, 2, .., N − 1),

ω(q) = 2κ cosh(h + iq) (11)
is the dispersion relation of plane waves of the linear HN model, and
Z t
2 χ|A|2
θ(t) = χ|A| dτ exp[2Im(ω(q))τ ] = [exp(2Im(ω(q))t) − 1] (12)
0 2Im(ω(q))
is the nonlinear time-varying phase shift accumulated during the evolution. Clearly, in the Hermitian
(conservative) limit h = 0 such plane waves reduce to Eq.(2) of the DNLSE. For h > 0, the frequency
ω(q) is complex, corresponding to wave damping for q < 0 and wave amplification for q > 0, with the
largest growth rate attained at q = π/2. Unlike the Hermitian limit of the DNSE, for h > 0 the plane
solutions are exponentially growing or damped along the propagation, with an exponential growth of the
phase term θ(t) according to Eq.(12). Only for q = 0, π the growth rate Im(ω) vanishes and θ(t) varies

4
 3 DYNAMICAL GROWTH BLOCKADE

linearly with time t like in the Hermitian limit. Hence, the usual concept of ”linear stability” and ”mod-
ulational instability” of nonlinear plane wave solutions for the DNSLE, established by the behavior of
Lyapunov exponents (eigenvalues) of linearized equations [Eq.(8)], should be properly extended when h
is nonvanishing. To test the stability of the plane wave solution (10), let us consider the perturbed solu-
tion
ψn (t) = ϕn (t) [1 + un (t)] (13)
where un (t) is a small perturbation, i.e. |un (t)| ≪ 1. Substitution of Ansatz (13) into Eq.(9) and after
linearization yields the following evolution equation for the perturbations un (t)
dun
i = −ω(q)un + κ exp(h + iq)un+1 + κ exp(−h − iq)un−1 + θ̇(un + u∗n ) (14)
dt
where θ̇ = (dθ/dt) = χ|A|2 exp[2Im(ω(q))t]. The plane wave solution (10) is said to be linearly stable
whenever, for any initial perturbation un (0) at time t = 0 with |un (0)| ≪ 1, as t → ∞ |un (t)| remains
much smaller than one. The solve the linearized equation (14), let us introduce the Ansatz
un (t) = F (t) exp(iQn) + G∗ (t) exp(−iQn), (15)
where Q is the spatial wave number of the perturbation. Substitution Eq.(15) into Eq.(14) yields the fol-
lowing coupled equations for the complex amplitudes F (t) and G(t)
dF
i = [ω(q + Q) − ω(q)]F + θ̇(F + G) (16)
dt
dG
−i = [ω ∗ (q − Q) − ω ∗ (q)]G + θ̇(G + F ) (17)
dt
To establish the asymptotic behavior of the amplitudes F and G as t → ∞, we should distinguish three
cases, depending on the value of the wave number q of the nonlinear wave mode.
(i) For a wave number −π < q < 0, corresponding to Im(ω(q)) < 0, the nonlinear wave is damped,
i.e. ψn (t) → 0 as t → ∞, and in the large t limit the evolution equations for the amplitudes F (t) and
G(t) [Eqs.(16) and (17)] are uncoupled and their growth rates are given by the imaginary parts of ω(q ±
Q). Clearly, there are intervals of wave number Q of spatial perturbations such that Im(ω(q ± Q)) >
Im(ω(q)). This means that, as t → ∞, the perturbation amplitude |un (t)| exponentially grows from a
small initial value, indicating that the nonlinear plane wave with wave number q is unstable.
(ii) For q = 0, π, the frequency ω(q = 0) = ±2κ cosh(h) of the nonlinear plane wave is real, θ̇ is constant
in time and a standard linear stability analysis based on eigenvalue analysis, as in the DNLSE case, can
be performed in Eqs.(16) and (17). The analysis indicates that the nonlinear wave is modulationally un-
stable is some intervals of the perturbation wave numbers Q.
(iii) For a wave number 0 < q < π, corresponding to Im(ω(q)) > 0, the nonlinear wave is amplified
and it turns out that at any wave number Q of spatial perturbations the amplitudes G(t), F (t) are un-
bounded. This result is nontrivial and involves a detailed asymptotic analysis of the non-autonomous
system Eqs.(16) and (17), which is presented in Appendix A.
The above results indicate that, unlike the NLDSE, in the nonlinear HN model any nonlinear plane wave
undergoes MI, regardless the value of the wave number q. Remarkably, such a result holds even for plane
waves displaying an exponential growth, i.e. with Im(ω(q)) > 0, which are unstable to any spatial wave
number perturbation.

3 Dynamical growth blockade

The linear stability analysis of plane wave solutions in the nonlinear HN model reveals that these modes
are always modulationally unstable, unlike the conservative limit of the DNSLE (h = 0). However, this
analysis does not provide insight into the system’s behavior as the instability evolves. To address this,
numerical simulations of Eq.(9) are required. Assuming as an initial condition a plane wave with spa-
tial wave number q, slightly perturbed by a small random fluctuations of the amplitudes, an instability

5
 3 DYNAMICAL GROWTH BLOCKADE

Figure 2: Illustration of the growth blockade effect. (a) Numerically-computed temporal evolution of the amplitudes
|ψn (t)| on a pseudocolor map in a lattice comprising N = 12 sites with periodic boundary conditions for parameter val-
ues κ = 1 (hopping rate), h = 0.1 (imaginary gauge field) and for a nonlienar coefficient χ = 1 (upper panels) and χ = −1
(lower panels). The initial condition is the plane wave ψn (t) = A exp(iq0 n) with amplitude A = 1 and spatial wave num-
ber q0 = π/2, P perturbed by adding a small random noise. (b) Corresponding temporal evolution of the total excitation
2
intensity P (t) = n |ψn (t)| (solid curve). The dashed curve is the behavior of P (t) of the nonlinear plane wave solution,
given by the exponential law P (t) = N |A|2 exp {2Im(ω(q0 ))t}. Note that the exponential growth of the total excitation
intensity is arrested after an initial transient, and P (t) settles down to a nearly constant value with small fluctuations for
times t >∼ 10. The arrest of the secular intensity growth is associated to the emergence of an irregular spatial distribution
of ψn (t).

clearly sets in leading to complex spatio-temporal patterns. Examples of numerical results are shown in
Figs.2,3 and 4 for both positive and negative sign of the nonlinear coefficient χ. The figures illustrate
the typical evolution of the wave amplitude |ψn (t)| on a pseudocolor map, and corresponding behavior
of total intensity P (t), defined by
XN
P (t) = |ψn (t)|2 , (18)
n=1

as obtained by numerical integration of Eq.(9) using an accurate variable-step fourth-order Runge-Kutta

method for a few different initial conditions ψn (0). A lattice comprising N = 12 sites with periodic
boundary conditions has been assumed in the numerical simulations, and both the focusing (χ = 1)
and defocusing (χ = −1) nonlinearities have been considered. As an initial condition, we assumed a
nonlinear plane wave with wave number q and amplitude A, slightly perturbed by a small random noise,
namely
ψn (0) = A(1 + rn ) exp(iqn) (19)
where rn is a complex random variable whose real and imaginary parts are uniformly distributed in the
interval (−δ, δ) (typically we set δ = 10−3 ). A general characteristic of the clean HN model under pe-
riodic boundary conditions –corresponding to the linear limit, χ = 0, in Eq. (9) – is that, for a wide
range of initial conditions (including random ones), the total excitation intensity in the system P (t) dis-
plays a secular growth in time, with a largest growth rate given by 2κ sinh h, which is the largest value
of Im(ω(q)) attained at q = π/2. One might naively expect a similar behavior in the nonlinear regime

6
 3 DYNAMICAL GROWTH BLOCKADE

Figure 3: Same as Fig.2, but for q0 = 0.

Figure 4: Same as Fig.2, but for q0 = −π.

7
 3 DYNAMICAL GROWTH BLOCKADE

(χ ̸= 0), since the nonlinear term is conservative and does not directly affect the power balance in the
system. However, this is not the case: after an initial transient, a blockade of the secular intensity growth
is observed, as shown in Figs.2,3 and 4. Interestingly, when the initial condition is the plane wave corre-
sponding to the wave number q = π/2 with the largest growth rate 2κ sinh h, the temporal evolution
of P (t) shows an initial exponential growth of the nonlinear plane wave mode, followed by an abrupt
halt of growth, as shown in Fig.2(b). A natural question arises: what is the physical origin of the growth
blockade in the nonlinear HN model observed in the numerical simulations? Broadly speaking, this growth
arrest can be understood as a stopping of convective motion due to self-induced disorder in the system
originating from the modulation instability. In fact, in the clean (disorder-free) linear HN model the growth
of P (t) is the result of a convective instability in the bulk of the lattice [26, 66, 67, 68], associated to a
biased drift of the excitation: preventing the biased drift, for example by an edge boundary or by strong
impurities or strong enough on-site potential disorder in the lattice, halt the secular growth of the total
excitation intensity P (t). For example, when strong on-site potential disorder Vn is added to the linear
HN model, i.e. when we consider the model
dψn
i = κψn+1 + κψn−1 + Vn ψn , (20)
dt
owing to Anderson localization the drift along the lattice is prevented and the corresponding growth of
excitation is arrested. Such an effect is well understood in the linear HN model with disorder and related
to the dubbed non-Hermitian Anderson localization transition predicted in the early studies of the HN
model [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. In the nonlinear HN model, there is not any on-site potential dis-
order in the system nor edges, however a kind of on-site potential disorder is self-induced by the modu-
lational instability of the nonlinear plane wave solution. In fact, the MI makes the local intensity term
Vn = |ψn |2 in Eq.(9) basically equivalent to a self-induced spatially-inhomogeneous and irregular on-
site potential, which thus prevents drift and amplification of small-amplitude excitation along the non-
linear lattice. To check the significance of such a physical interpretation, let ψn (t) = ϕn (t) be a spatially-
inhomogeneous (and irregular) solution to Eq.(9), as observed in Fig.2 after the initial transient growth,
and let us consider the evolution of a small perturbation un (t) added to ϕn (t). After linearization, the
perturbation un (t) satisfies the linearized system
dun
i = κ exp(h)un+1 + κ exp(−h)un−1 + 2Vn (t)un + Wn (t)u∗n (21)
dt
where we have set
Vn (t) ≡ χ|ϕn (t)|2 , Wn (t) ≡ χϕ2n (t) (22)
which act as effective potentials for the the perturbation wave amplitude un (t). For an irregular solution
ϕn (t), as observed in the numerical simulations of Fig.1 after the initial transient growth, the two effec-
tive potentials are spatially irregular and time varying. To provide qualitative insights on the growth
blockade effect, we can neglect rapid time oscillations of the effective potentials and replace them in Eqs.(21)
and (22) by the their local time average
Z t+∆t Z t+∆t Z t+∆t Z t+∆t
χ χ 2 χ χ
Vn = dτ Vn (τ ) = dτ |ϕn (τ )| , W n = dτ Wn (τ ) = dτ ϕ2n (τ ) (23)
∆t t ∆t t ∆t t ∆t t
where ∆t is the time interval of averaging. For large enough ∆t, the two averaged potentials are slow-
ing varying with time t and can be assumed to be locally constant, as shown as an example in Fig.5(a).
Under such approximation, the growth of perturbations un (t) in time is established by the position in
complex plane of the 2N eigenvalues λ of the following 2N × 2N matrix
 
A B
S=
−B ∗ −A

8
 4 CONCLUSION

R t+∆t
Figure 5: (a) Numerically-computed behavior of the local averaged potential V n (t) = (χ/∆t) t |ϕn (τ )|2 dτ versus time
t, corresponding to the solution ϕn (t) shown in Fig.2 with χ = 1, on a pseudo color map. Average time interval is ∆t = 1.
(b) Detailed behavior of V n for time t = 19. (c) Eigenvalue spectrum of the matrix S with potential V n plotted in panel
(b). Note that the imaginary part of all eigenvalues vanishes.

where the N × N matrix A is given by

 
V1 κ exp(h) 0 ... 0 0 κ exp(−h)
κ exp(−h) V 2 κ exp(h) ... 0 0 0 
 
A=  0 0 0 ... 0 0 ,

 0 0 0 ... κ exp(−h) V N −1 κ exp(h) 
κ exp(h) 0 0 ... 0 κ exp(−h) VN

B is the N × N diagonal matrix defined by Bn,m = W n δn,m , and B ∗ is the element-wise complex con-
jugate of B. A typical eigenvalue spectrum λ in complex plane is shown in Fig.5(c), indicating that all
states are marginally stable, i.e. the imaginary part of λ is vanishing. Hence small-amplitude pertur-
bations cannot secularly grow once the self-induced potentials Vn and Wn have been created. Such an
analysis is of course not rigorous, since it disregard time fluctuations of the self-induced disordered po-
tentials, however it is capable of catching the main physics underlying the growth blockade effect as ob-
served in numerical simulations.

4 Conclusion
In this work, we have investigated modulational instability in the nonlinear extension of the Hatano-
Nelson model under periodic boundary conditions, unveiling a novel dynamical phenomenon: growth
blockade. In contrast to the unbounded growth of the total excitation intensity seen in the linear HN
model – where an unsaturated convective instability dominates – in the nonlinear regime, we observe a
suppression of this expected growth after an initial transient phase. This growth arrest can be under-
stood as a self-induced suppression of transport in the lattice, driven by dynamically induced disorder,
which is a hallmark of the HN model with strong disorder or edges. This mechanism is fundamentally
different from the growth clamping of modulation instability of planes waves seen in other conservative
or dissipative systems, where the instability is stopped by some gain saturation mechanisms in the sys-
tem leading to pattern formation or localized structures [69].
Our results highlight the critical role played by the interplay between nonlinearity and non-reciprocal
transport in shaping the system’s behavior. The nonlinear HN model introduces entirely new dynami-
cal regimes that are absent in both the linear HN model and more conventional nonlinear systems like
the discrete nonlinear Schrödinger equation. Specifically, the combination of Kerr-type nonlinearity and
non-reciprocal transport inherent in the HN model produces a range of complex behaviors, including the
nonlinear skin effect and the formation of skin solitons in systems with open boundary conditions, which

9
 A APPENDIX A: STABILITY ANALYSIS OF NONLINEAR PLANE WAVES

have been recently predicted and observed in related works. The observed growth blockade in periodic
settings adds another layer to our understanding of these systems.
This work serves as a foundation for further exploration of nonlinear non-Hermitian systems, with po-
tential applications in fields such as optics, photonics, and quantum non-Hermitian many-body physics
in the mean-field regime.
Conflict of Interest
The author declares no conflict of interest.
Data Availability Statement The data that support the findings of this study are available from the
corresponding author upon reasonable request.
Acknowledgements
The author acknowledges the Spanish State Research Agency, through the Severo Ochoa and Maria de
Maeztu Program for Centers and Units of Excellence in R&D (Grant No. MDM-2017- 0711).

A Appendix A: Stability analysis of nonlinear plane waves

In this Appendix it is shown that the nonlinear wave given by Eq.(4), for a spatial wave number 0 < q <
π corresponding to an amplified wave (Im(ω(q)) > 0), is modulationally unstable for any spatial wave
number Q of the perturbation. To this aim, let us rewrite the coupled equations (16) and (17), describ-
ing the evolution of the amplitudes F (t) and G(t) of the perturbation, after introduction of the stretched
time variable T , defined via the relation
Z t
exp[2Im(ω(q))t] − 1
T = dτ exp[2Im(ω(q))τ ] = . (24)
0 2Im(ω(q))
Taking into account that d/dt = exp[2Im(ω(q))t](d/dT ), Eqs.(16) and (17) take the form
dF
i = ϵ1 (T )F + σ(F + G) (25)
dT
dG
i = ϵ2 (T )F − σ(F + G) (26)
dT
where we have set σ ≡ χ|A|2 and
ω(q + Q) − ω(q)
ϵ1 (T ) = (27)
1 + 2Im(ω(q))T
ω ∗ (q) − ω ∗ (q − Q)
ϵ2 (T ) = . (28)
1 + 2Im(ω(q))T
Let us note that ϵ1,2 (T ) ∼ 1/T → 0 as T → ∞ for any value of the spatial perturbation wave number Q.
Let us first consider the case Q = 0, i.e. the long-wavelengh limit (Q → 0) of the perturbation. In this
case ϵ1,2 (T ) = 0 and one obtains the autonomous linear system
dF
i = χ|A|2 (F + G) (29)
dT
dG
i = −χ|A|2 (G + F ). (30)
dT
which can be solved by standard methods. Since the matrix of such a linear system displays two vanish-
ing eigenvalues corresponding to an exceptional point, for rather arbitrary initial conditions it turns out
that the amplitudes F (T ) and G(T ) are unbounded and grow as ∼ T . In physical time t, according to
Eq.(24) this means that the amplitudes F and G of the perturbation grow as ∼ exp[2Im(ω(q))t], indicat-
ing that the nonlinear plane wave solution ϕn (t) is unstable.

10
 A APPENDIX A: STABILITY ANALYSIS OF NONLINEAR PLANE WAVES

Let us now consider the case of a spatial perturbation with Q ̸= 0. In this case ϵ1,2 (T ) are distinct and
complex, vanishing as T → ∞, i.e. the exceptional point of the linear system is attained asymptotically
in time. To study the asymptotic behavior of the solutions in this case, we note that Eqs.(26) and (29)
can be regarded as the equations of a two-level system described by the time-dependent non-Hermitian
Hamiltonian  
ϵ1 (T ) + σ σ
H(T ) = (31)
−σ ϵ2 (T ) − σ
for which a non-adiabatic asymptotic analysis can be performed, similar to the one developed to study
exceptional point encircling (see e.g. [70, 71, 72, 73]). To this aim, let us indicate by |u1,2 (T )⟩ and |u†1,2 (T )⟩
the instantaneous right and left eigenvectors, respectively, of H(T ), and by λ1,2 (T ) the corresponding in-
stantaneous eigenvalues, i.e. H|u1,2 (T )⟩ = λ1,2 (T )|u1,2 (T )⟩ and H† |u†1,2 (T )⟩ = λ∗1,2 (T )|u†1,2 (T )⟩, with the
normalization conditions ⟨u†i |uk ⟩ = δi,k (i, k = 1, 2). It can be readily shown that, in the large T limit, at
leading order one has
√ p Θ
λ1,2 (T ) = ± σ ϵ1 (T ) − ϵ2 (T ) ∼ ± √ (32)
T
for the instantaneous eigenvalues, where
s
√ ω(q + Q) − ω(q) − ω ∗ (q − Q) + ω ∗ (q)
Θ≡ σ ,
2Im(ω(q))

and
λ∗1,2 + σ
   
σ 1
|u1,2 ⟩ = , |u†1,2 ⟩ = (33)
λ1,2 − σ 2σλ∗1,2 σ
for the corresponding right and left eigenvectors. Note that, according to Eq.(32), the two eigenvalues λ1
and λ2 = −λ1 are complex with opposite imaginary parts. Without loss of generality, we will assume λ1
as the dominant eigenvalue, i.e. Im(λ1 ) > 0 and Im(λ2 ) < 0, corresponding to Im(Θ) > 0. Let us write
the solution to Eqs.(25) and (26) in the adiabatic basis using the Ansatz
   Z T   Z T 
F (T )
= α1 (T ) exp −i dtλ1 (t) |u1 (T )⟩ + α2 (T ) exp −i dtλ2 (t) |u2 (T )⟩. (34)
G(T )

Substitution of Eq.(34) into Eqs.(25) and (26) yields the following coupled equations for the adiabatic
amplitudes α1 (T ) and α2 (T )
 Z T 
dα1 † du1 † du2
= −α1 ⟨u1 | ⟩ − α2 ⟨u1 | ⟩ exp i dt(λ1 − λ2 ) (35)
dT dT dT
 Z T 
dα2 † du2 † du1
= −α2 ⟨u2 | ⟩ − α1 ⟨u2 | ⟩ exp −i dt(λ1 − λ2 ) (36)
dT dT dT
which are exact at this stage. Using Eq.(33), Eqs.(35) and (36) ca be cast in the form

λ˙1 λ˙2
 Z T 
dα1
= − α1 − α2 exp i dt(λ1 − λ2 ) (37)
dT 2λ1 2λ1
λ˙2 λ˙1
 Z T 
dα2
= − α2 − α1 exp −i dt(λ1 − λ2 ) (38)
dT 2λ2 2λ2

where the dot indicates the derivative with respect to T , i.e. λ̇1,2= (dλ1,2 /dT ). Note that, using Eq.(32)
one has asymptotically
√ 
 Z T  
exp −i dt(λ1 − λ2 ) ∼ exp −4iΘ T (39)

11
REFERENCES REFERENCES

√
which grows in time as ∼ exp[4Im(Θ) T ]. To determine the asymptotic behavior of amplitudes α1,2 (T )
in the large T limit, one can make the assumption, to be justified a posteriori, that the second term√on
the right hand side of Eq.(37) is negligible, i.e. that α2 (T ) can grow but slower than ∼ exp[4Im(Θ) T ].
−1/2
Under this assumption Eq.(37) can be solved yielding α1 (T ) = λ1 , which secularly grows in time as ∼
T 1/4 . From Eq.(38), one can then calculate the asymptotic long-time behavior of α2 (T ) by neglecting the√
first term on the right hand side, and the corresponding
√ solution grows in time as ∼ (1/T 1/4 ) exp[4Im(Θ) T ].
Such a growth is smaller than ∼ exp[4Im(Θ) T ], consistently with the initial trial assumption. To sum
up, the asymptotic analysis indicates that in the large T limit the behavior of F (T ), G(T ) is dominated
by the first term√on the right hand side of Eq.(34), and they grow in time asymptotically as |F |, |G| ∼
T 1/4 exp[2Im(Θ) T ]. This proves that the nonlinear plane wave solution is modulationally unstable for
Q ̸= 0 as well.

References
[1] N. Hatano, D. R. Nelson, Phys. Rev. Lett. 1996, 77, 570.
[2] N. Hatano, D.R. Nelson Phys. Rev. B 1997, 56, 8651.
[3] N. Hatano, D. R. Nelson, Phys. Rev. B 1998, 58, 8384.
[4] P. W. Brouwer, P. G. Silvestrov, C. W. J. Beenakker, Phys. Rev. B 1997, 56, R4333(R).
[5] N. M. Shnerb, D. R. Nelson, Phys. Rev. Lett. 1998, 80, 5172.
[6] P. G. Silvestrov, Phys. Rev. B 1998, 58, R10111(R).
[7] I. Ya. Goldsheid, B. A. Khoruzhenko, Phys. Rev. Lett. 1998, 80, 2897.
[8] I. V. Yurkevich, I. V. Lerner, Phys. Rev. Lett. 1999, 82, 5080.
[9] A. V. Kolesnikov, K. B. Efetov, Phys. Rev. Lett. 2000, 84, 5600.
[10] J. Heinrichs, Phys. Rev. B 2001, 63, 165108.
[11] L. G. Molinari, J. Phys. A 2009, 42, 265204.
[12] F. Hebert, M. Schram, R. T. Scalettar, W. B. Chen, Z. Bai, Eur. Phys. J. B 2011, 79, 465.
[13] S. Longhi, D. Gatti, G. Della Valle, Sci. Rep. 2015, 5, 13376.
[14] B.T. Torosov, G. Della Valle, and S. Longhi, Phys. Rev. A 2014, 89, 063412.
[15] S. Longhi, D. Gatti, G. Della Valle, Phys. Rev. B 2015, 92, 094204 .
[16] Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, M. Ueda, Phys. Rev. X 2018, 8,
031079.
[17] H. Shen, B. Zhen, L. Fu, Phys. Rev. Lett. 2018, 120, 146402.
[18] S. Yao, Z. Wang, Phys. Rev. Lett. 2018, 121, 086803.
[19] S. Yao, F. Song, Z. Wang, Phys. Rev. Lett. 2018, 121, 136802.
[20] F.K. Kunst, E. Edvardsson, J.C. Budich, E.J. Bergholtz, Phys. Rev. Lett. 2018, 121, 026808.
[21] C.H. Lee, R. Thomale, Phys. Rev. B 2019, 99, 201103(R).
[22] H. Zhou , J. Y. Lee, Phys. Rev. B 2019, 99, 235112.
[23] F. Song, S. Yao, Z. Wang, Phys. Rev. Lett. 2019, 123, 246801..

12
REFERENCES REFERENCES

[24] K. Yokomizo, S. Murakami, Phys. Rev. Lett. 2019, 123, 066404.

[25] K. Kawabata, K. Shiozaki, M. Ueda, M. Sato, Phys. Rev. X 2019, 9, 041015.
[26] S. Longhi, Phys. Rev. Research 2019, 1, 023013.
[27] N. Okuma, K. Kawabata, K. Shiozaki, M. Sato, Phys. Rev. Lett. 2020, 124, 086801.
[28] E.J. Bergholtz, J.C. Budich, F.K. Kunst, Rev. Mod. Phys. 2021, 93, 015005.
[29] S. Longhi, Phys. Rev. Lett. 2022, 128, 157601,
[30] N. Okuma, M. Sato, Ann.Rev.Condensed Matter Phys. 2023, 14 83.
[31] X. Zhang, T. Zhang, M.H. Lu, Y.F. Chen, Adv. Phys.: X 2022, 7 .2109431.
[32] R. Lin, T. Tai, L. Li, C.H. Lee, Front. Phys 2023,.18, 53605.
[33] R. Hamazaki, K. Kawabata, M. Ueda, Phys. Rev. Lett. 2019, 123, 090603.
[34] T. Liu, J.J. He, T. Yoshida, Z.-L. Xiang, F. Nori, Phys. Rev. B 2020, 102, 235151.
[35] L.-J. Zhai, S. Yin, G.-Y. Huang, Phys. Rev. B 2020, 102, 064206.
[36] S.-B. Zhang, M. M. Denner, T. Bzdusek, M. A. Sentef, T. Neupert, Phys. Rev. B 2022, 106,
L121102.
[37] K. Kawabata, K. Shiozaki, S. Ryu, Phys. Rev. B 2022, 105, 165137.
[38] K. Suthar, Y.-C. Wang, Y.-P. Huang, H. H. Jen, J.-S. You, Phys. Rev. B 2022, 106, 064208.
[39] B. Dóra, C.P. Moca, Phys. Rev. B 2022, 106, 235125.
[40] T. Yoshida, Y. Hatsugai, Phys. Rev. B 2022, 106, 205147.
[41] R. Shen, C.H. Lee, Commun. Phys. 2022, 5, 238.
[42] W.N. Faugno, T. Ozawa, Phys. Rev. Lett. 2022, 129, 180401.
[43] S. Longhi, Phys. Rev. B 2023, 108, 075121.
[44] S. Longhi, Ann. Phys. 2023, 535, 2300291.
[45] T. Orito, K.-I. Imura, Phys. Rev. B 2023, 108, 214308.
[46] M. Zheng, Y. Qiao, Y. Wang, J. Cao, S. Chen, Phys. Rev. Lett. 2024, 132, 086502.
[47] Y. Qin, L. Li, Phys. Rev. Lett. 2024, 132, 096501.
[48] T. Yoshida, S.-B. Zhang, T. Neupert, N. Kawakami, Phys. Rev. Lett. 2024, 133, 076502.
[49] B.H. Kim, J.-H. Han, M.J. Park, Commun. Phys. 2024, 7, 73.
[50] M. Yang, L. Yuan, C.H. Lee, 2024 arXiv:2410.01258v1.
[51] C. Yuce, Phys. Lett. A 2021, 408, 127484.
[52] M. Ezawa, Phys. Rev. B 2022, 105, 125421.
[53] Z.-X. Zhang, J. Cao, J.-Q. Li, Y. Zhang, S. Liu, S. Zhang, H.-F. Wang, Phys. Rev. B 2023, 108,
125402.
[54] B.M. Manda, R.Carretero-Gonzalez, P.G., Kevrekidis, V. Achilleos, Phys. Rev. B 2024, 109,
094308.

13
REFERENCES REFERENCES

[55] I. Komis., Z.H. Musslimani, and K.G. Makris, Opt. Lett. 2023, 48, 6525..
[56] S. Wang, B. Wang, C. Liu, C. Qin, L. Zhao, W. Liu, S. Longhi, P. Lu, 2024 arXiv:2409.19693.
[57] S. Flach, C.R. Willis, Phys. Rep. 1998, 295, 181.
[58] P.G. Kevrekidis, K.O. Rasmussen, A.R. Bishop, Int. J. Mod. Phys. B 2001, 15, 2833.
[59] P.G. Kevrekidis, The Discrete Nonlinear Schrödinger Equation (Springer Tracts in Modern Physics
232, Springer, 2009).
[60] F. Lederer, G.I. Stegeman, D.N. Christodoulides, G. Assanto, M. Segev, Y. Silberberg, Phys. Rep.
2008, 463, 1
[61] H. S. Eisenberg, Y. Silberberg, R. Morandotti, A. R. Boyd, J. S. Aitchison, Phys. Rev. Lett. 1998,
81, 3383.
[62] A. Trombettoni, A. Smerzi, Phys. Rev. Lett. 2001, 86, 2353.
[63] D.K. Campbell, S. Flach, Yu.S. Kivshar, Phys. Today 2004, 57, 43.
[64] E. Picari, A. Ponno, L. Zanelli, Ann. Henri Poincaré 2022, 23 1525.
[65] Y.S. Kivshar, M. Peyrard, Phys. Rev. A 1992, 46, 3198 (1992).
[66] M. Santagiustina, P. Colet, M. San Miguel, D. Walgraef, Phys. Rev. Lett. 1997, 79, 3633.
[67] P. Colet, D., Walgraef, M. San Miguel, Eur. Phys. J. B 1999, 11, 517.
[68] S. Longhi, Phys. Rev. A 2013, 88, 052102.
[69] M. C. Cross, P. C. Hohenberg, Rev. Mod. Phys. 1993, 65, 851.
[70] R. Uzdin, A. Mailybaev, N. Moiseyev, J. Phys. A: Math. Theor. 2011, 44, 435302.
[71] M. V. Berry, R. Uzdin, J. Phys. A: Math. Theor. 2011, 44, 435303.
[72] H.Wang, L. J. Lang, Y. D. Chong, Phys. Rev. A 2018, 98, 012119.
[73] S. Longhi, L. Feng, Phys. Rev. B 2023, 107, 085122.

14