# Sequence disorder-induced first order phase transition in confined
polyelectrolytes

V. Stepanyan  Yerevan State University, Yerevan, Armenia  A. Badasyan
University of Nova Gorica, Nova Gorica, Slovenia  V. Morozov  Yerevan State
University, Yerevan, Armenia  Y. Mamasakhlisov<EMAIL_ADDRESS>Yerevan
State University, Yerevan, Armenia  Institute of Applied Problems of Physics,
Yerevan, Armenia  R. Podgornik<EMAIL_ADDRESS>School of Physical
Sciences, University of Chinese Academy of Sciences, Beijing, 100049, China
CAS Key Laboratory of Soft Matter Physics, Institute of Physics, Chinese
Academy of Sciences, Beijing, 100190, China  Kavli Institute for Theoretical
Sciences, University of Chinese Academy of Sciences, Beijing, 100049, China
Wenzhou Institute of the University of Chinese Academy of Sciences, Wenzhou,
Zhejiang 325011, China  Department of Physics, Faculty of Mathematics and
Physics, University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia

###### Abstract

We consider a statistical mechanical model of a generic flexible
polyelectrolyte, comprised of identically charged monomers with long range
electrostatic interactions, and short-range interactions quantified by a
disorder field along the polymer contour sequence, which is randomly quenched.
The free energy and the monomer density profile of the system for no
electrolyte screening are calculated in the case of a system composed of two
infinite plates with an intervening oppositely charged polyelectrolyte chain.
We show that the effect of the contour sequence disorder mediated by short-
range interactions leads to an enhanced localization of the polyelectrolyte
chain and a first order phase transition at a critical value of the inter-
surface spacing. This phase transition results in an abrupt change of the
pressure from negative to positive values, effectively eliminating
polyelectrolyte mediated bridging attraction.

##  I Introduction

Interactions involving biologically relevant heteropolymers such as nucleic
acids, polysaccharides and polypeptides Leckband and Israelachvili (2001) are
grosso modo composed of two components Rubinstein and Colby (2003), of which
the relatively generic and longer ranged one originates in the electrostatic
charges, being at the origin of the polyelectrolyte phenomenology Dobrynin and
Rubinstein (2005), and is often invoked as the primary cause of their solution
behavior Muthukumar (2023). The more specific, shorter ranged interactions
French  _et al._ (2010), related to non-universal chemical identity of the
various monomer units Elias (2012), cannot be quantified by a single parameter
analogous to the electrostatic charge and are standardly invoked when the
observed polyelectrolyte phenomenology cannot be understood solely based on
the electrostatic interactions Smiatek (2020).

The origin of the specific short-range interactions can be understood within
different theoretical frameworks as stemming from the differences in
polarizability of the monomer subunits in the context of van der Waals
interactions between polymers Lu  _et al._ (2015), from the differences in the
dissociation properties of ionizable moieties within a generalized charge-
regulated electrostatics Borkovec  _et al._ (2001), or within a detailed
description of water structuring models accounting for the solvent mediated
interactions Blossey and Podgornik (2022), to list just a few. While each of
these various theoretical frameworks of accounting for specific interactions
within the context of polymer solution theory leads to differing predictions
on the details of these interactions, they all agree that these interactions
are important to understand the behavior of polymers in aqueous solutions.

The diversity of these specific interactions, originating from different types
of monomer units, such as 4 canonical nucleotides in nucleic acids Bloomfield
_et al._ (2000) or 20 canonical amino acids in polypeptides Finkelstein and
Ptitsyn (2016), is therefore much more pronounced in terms of short range non-
electrostatic interactions than long range electrostatic effects French  _et
al._ (2010)and is closely related to the sequence of monomers along the chain
leading to interesting effects for polymer-polymer interaction as well as
polymer-substrate interaction. Gauging the relative importance of the generic
long-range Coulomb interaction vs. the sequence specific short-range
interaction is particularly important when trying to assess their role in the
complicated multi-level problem of RNA folding Fallmann  _et al._ (2017), when
formulating realistic models of RNA virus co-assembly Twarock  _et al._
(2018), or manipulating electrostatic repulsion in the background of sequence-
specific Watson-Crick base pair interaction of ssRNAs Adhikari  _et al._
(2020). In all these cases an interplay of both types of interaction needs to
be properly accounted for.

In what follows we will investigate a model hetero-polyelectrolyte chain,
characterized by both generic long-range Coulomb interactions between its
equally charged monomers, with superimposed specific short-range interactions
depending on the type of the monomer along its sequence. We will assume that
the monomer sequence can be modeled as a quenched, random sequence represented
by a single random variable characterized by a Gaussian distribution with zero
mean. We will write the free energy of the system within the self-consistent
field theory of a single chain confined between two planar surfaces within a
replica symmetric Anzatz and investigate the interactions between two planar
surfaces, charged oppositely from the chain monomers, as a function of the
separation between them. We will show that the model confined polymer chain
exhibits an additional localization induced by sequence disorder akin to the
Anderson localization and will assess the relative importance of the disorder
effects superimposed on an electrostatic background.

##  II Model and Methods

###  II.1 Model

We consider a generic flexible hetero-polyelectrolyte, e.g. ssRNA, see Fig. 1,
confined between two planar surfaces, comprised of equally charged monomers
interacting via long-range Coulomb interactions, as well as short range
"chemical" interactions dependent on the monomer sequence. While the charge of
the monomers is assumed to be the same for all the monomers, the short range
interactions are sequence specific, with this specificity being encoded in the
sequence disorder. This implies that instead of considering a detailed
sequence of monomers characterized by different chemical identities, we
characterize the strength of the short-range interactions by a random,
Gaussian distributed, variable. Similar systems without the sequence disorder
affecting the short-range interactions Borukhov  _et al._ (1999); Shafir  _et
al._ (2003); Siber and Podgornik (2008); Podgornik (1991); Andelman (1995) and
the problem of a polymer chain with excluded volume interactions in quenched
random media Muthukumar (1989); Hribar-Lee  _et al._ (2011); Bratko and
Chakraborty (1995) have been studied before. Our main goal here is to address
the interplay between the disordered short-range interactions and long-range
Coulomb interactions. In order to simplify the calculations we consider a salt
free system with the surfaces oppositely charged from the polymer chain.

Figure 1:  The polyelectrolyte negatively charged chain is confined between
two positively charged infinite plates with surface charge density $\sigma$.
The charge per monomed is $pe$, where $e$ is the electron charge, and the
short-range interaction between two parts of the chain is
$v_{\tau\tau^{\prime}}({\bf x})=v_{0}\xi_{\tau}\xi_{\tau^{\prime}}\delta({\bf
x})$, where $\xi_{\tau}$ is the random sequence disorder variable and $\tau$
is the position along the contour of the chain.

The charge per monomer is assumed to be $pe$, where $e$ is the electron charge
and $0<p<1$.The position of a monomer is characterized by a continuous
Cartesian coordinate ${\bf r}(\tau)$, where $\tau\in[0,N]$ marks the position
along the contour of the chain, and $N$ is the dimensionless contour length of
the chain. The random variable characterizing the strength of the short-range
interaction is $\xi_{\tau}$. $\\{\xi\\}$s are considered to be independent,
quenched random variables with the probability distribution

${\cal P}\\{\xi\\}=\prod_{\tau}p(\xi_{\tau}),$ (1)

where the sequence disorder is supposed to be Gaussian in order to make the
formalism tractable, therefore

$p(\xi_{\tau})=\frac{\mathit{e}^{-\frac{(\xi_{\tau}-\bar{\xi})^{2}}{2\xi^{2}}}}{\sqrt{2\pi\xi^{2}}}$
(2)

For the sake of simplicity we restrict ourselves to the vanishing average
value, $\bar{\xi}=0$, case.

The Hamiltonian of the system in this case has a general Edwardsian form with

$\displaystyle\beta{\cal H}[{\bf
r}(\tau)]=\frac{3}{2\ell^{2}}\int_{0}^{N}\,d\tau(\partial_{\tau}{\bf
r}(\tau))^{2}+$ (3)
$\displaystyle+\frac{\beta}{2}\int_{0}^{N}\,d\tau\int_{0}^{N}\,d\tau^{\prime}v_{\tau\tau^{\prime}}({\bf
r}(\tau)-{\bf r}(\tau^{\prime}))+\beta V[\\{{\bf r(\tau)}\\}],$

where $\ell$ is the Kuhn length, $V[\\{{\bf r}(\tau)\\}]$ describes all the
long-range interactions in the system, i.e., the Coulomb interactions between
monomers as well as the Coulomb interactions between monomers and the surface,
while

$v_{\tau\tau^{\prime}}({\bf x})=v_{0}\xi_{\tau}\xi_{\tau^{\prime}}\delta({\bf
x}),$ (4)

with a positive interaction constant $v_{0}>0$. Thus, the different monomers
will attract, but the similar ones will repel.

###  II.2 Partition function and free energy

We formulate the partition function of this three-dimensional polymer chain by
consistently including the contribution of the monomer sequence disorder and
analyzing the combined effect of randomly quenched short-range and homogeneous
long-range interactions on the confinement free energy. The partition function
of the model polyelectrolyte chain for every fixed sequence $\\{\xi\\}$ is
given by

$Z\\{\xi\\}=\int\,{\cal D}[{\bf r}(\tau)]~{}e^{-\beta{\cal H}[{\bf
r}(\tau)]}.$ (5)

In the limit of a very long chain, $N\rightarrow\infty$, the disorder part of
the free energy obeys the principle of self-averaging yielding Binder and
Young (1986); Parisi (2023):

${\cal F}=-k_{B}T\langle{\ln Z\\{\xi\\}}\rangle_{{\cal P}},$ (6)

where $\langle{...}\rangle_{{\cal P}}$ stands for the average with the
distribution function (Eq. 1). Self-averaging physically means that the
distribution of the free energy has a vary narrow peak in the vicinity of its
maximum, corresponding to the mean value of the free energy (Eq. 6).

The quenched free energy (Eq. 6) can now be estimated using the replica trick
Binder and Young (1986); Parisi (2023)

$-\beta{\cal F}=\lim_{n\rightarrow
0}\frac{\langle{Z\\{\xi\\}^{n}}\rangle_{{\cal P}}-1}{n},$ (7)

where $\beta=(k_{B}T)^{-1}$. In order to calculate the $n$-replica partition
function

$\displaystyle\langle{Z\\{\xi\\}^{n}}\rangle_{{\cal P}}=$
$\displaystyle\qquad=\int\,{\cal D}{\bf
r}e^{-\frac{3}{2\ell^{2}}\sum_{a=1}^{n}\int_{0}^{N}\,d\tau(\partial_{\tau}{\bf
r}^{a}(\tau))^{2}-\beta\sum_{a=1}^{n}V\\{{\bf r}^{a}\\}}\times$
$\displaystyle\quad\int\,{\cal D}\xi{\cal P}\\{\xi\\}e^{-\frac{\beta
v_{0}}{2}\int_{0}^{N}\,d\tau\int_{0}^{N}\,d\tau^{\prime}\xi_{\tau}\xi_{\tau^{\prime}}\sum_{a=1}^{n}\delta({\bf
r}^{a}(\tau)-{\bf r}^{a}(\tau^{\prime}))},$ (8)

and thus to estimate the free energy, several order parameters need to be
introduced, namely, the inter-replica overlap $q_{ab}({\bf x},{\bf
x}^{\prime})$, with $a<b$, the density of $a$-th replica monomers
$\rho_{a}({\bf x})$, and the density of the random variable $\xi_{\tau}$ of
$a$-th replica monomers $m_{a}({\bf x})$ as

$\displaystyle q_{ab}({\bf x},{\bf x}^{\prime})$ $\displaystyle=$
$\displaystyle\int_{0}^{N}\,d\tau\delta({\bf x}-{\bf r}^{a}(\tau))\delta({\bf
x}^{\prime}-{\bf r}^{b}(\tau))$ $\displaystyle\rho_{a}({\bf x})$
$\displaystyle=$ $\displaystyle\int_{0}^{N}\,d\tau\delta({\bf x}-{\bf
r}^{a}(\tau))$ $\displaystyle m_{a}({\bf x})$ $\displaystyle=$
$\displaystyle\int_{0}^{N}\,d\tau\xi_{\tau}\delta({\bf x}-{\bf r}^{a}(\tau)).$
(9)

Introducing next the Lagrange multipliers ${\hat{q}}_{ab}({\bf x},{\bf
x}^{\prime})$ and ${\hat{\rho}}_{a}({\bf x})$, conjugated to the first two
order parameters in (Eq. II.2), allows us to rewrite the $n$-replica partition
function as Garel  _et al._ (1997)

$\displaystyle\langle{Z\\{\xi\\}^{n}}\rangle_{{\cal P}}$ $\displaystyle\propto
e^{-\frac{nV}{2\tilde{v}}\ln(2\pi\beta v_{0})}\times$ (10)
$\displaystyle\times\int\,{\cal D}{\rho}\,{\cal D}{{\hat{\rho}}}\,{\cal
D}{q}\,{\cal
D}{\hat{q}}~{}e^{g(\rho,{\hat{\rho}},q,{\hat{q}})+\ln\zeta({\hat{\rho}},{\hat{q}})},$

where $\tilde{v}$ is the monomer volume. The explicit expression for the
functional
$g(\rho,{\hat{\rho}},q,{\hat{q}})+\ln\zeta({\hat{\rho}},{\hat{q}})$, as given
by (Eq. A), allows us to obtain the final form of the free energy functional
as :

$\displaystyle-\beta{\mathcal{F}}=\lim_{n\rightarrow
0}\frac{1}{n}\bigg{[}-\frac{nV}{2\tilde{v}}\ln(\beta v_{0})-\frac{1}{2}{\rm
Tr}\ln\hat{A}-\frac{1}{2}{\rm Tr}\ln\bigl{\\{}I+$
$\displaystyle\xi^{2}\hat{A}^{-1}\hat{B}\bigr{\\}}-\beta\sum_{a}W_{el}^{a}+i\sum_{a}\langle\rho_{a}|\hat{\rho}_{a}\rangle+i\sum_{a<b}\langle
q_{ab}|\hat{q}_{ab}\rangle+\ln\zeta\bigg{]},$ (11)

where the operators $\hat{A}$ and $\hat{B}$ have been defined as $n\times n$
operator matrices $\langle{\bf x}|{\hat{A}}_{ab}|{\bf
x}^{\prime}\rangle=\delta_{ab}\delta({\bf x}-{\bf x}^{\prime})[\frac{1}{\beta
v_{0}}+\xi^{2}\rho_{a}({\bf x})]$ and $\langle{\bf x}|{\hat{B}}_{ab}|{\bf
x^{\prime}}\rangle=(1-\delta_{ab})q_{ab}({\bf x},{\bf x}^{\prime})$. For
relevant details and derivations see the Appendix A .

###  II.3 Replica - symmetric solution

In order to estimate the free energy of the system we maximize the functional
$g(\rho,{\hat{\rho}},q,{\hat{q}})+\ln\zeta({\hat{\rho}},{\hat{q}})$ from (Eq.
10). We will restrict ourselves to the replica - symmetric order parameters:
$q_{ab}({\bf x},{\bf x}^{\prime})=q({\bf x},{\bf x}^{\prime})$,
$\hat{q}_{ab}({\bf x},{\bf x}^{\prime})=\hat{q}({\bf x},{\bf x}^{\prime})$,
$\rho_{a}({\bf x})=\rho({\bf x})$ and $\hat{\rho}_{a}({\bf x})=\hat{\rho}({\bf
x})$. See Appendix B for details.

Solving the system of the saddle - point equations (B.12,B.16,B.17) we get two
separate equations for the electrostatic and configurational degrees of
freedom, corresponding to the polymer Poisson-Boltzmann equation Markovich
_et al._ (2021) for the saddle-point electrostatic potential $\varphi({\bf
r})$

$\varepsilon\varepsilon_{0}\nabla^{2}\varphi({\bf r})=2ec_{0}\sinh(\beta
e\varphi({\bf r}))+peN\psi({\bf r})^{2}-\rho_{s}({\bf r})$ (12)

as well as the Edwards equation for the saddle-point polymer one-replica
density field $\psi({\bf r})$

$\displaystyle\mathcal{E}_{0}\psi({\bf
r})=\bigg{[}-\frac{\ell^{2}}{6}\nabla^{2}-\beta pe\varphi({\bf
r})+\frac{1}{2}\frac{1}{\mu+N\psi({\bf r})^{2}}$
$\displaystyle\bigg{(}\frac{1}{\tilde{v}}+\frac{\kappa}{(1-\kappa)^{2}}\frac{2\mu\psi^{2}({\bf
r})+N\psi({\bf r})^{4}}{\mu+N\psi({\bf r})^{2}}\bigg{)}\bigg{]}\psi({\bf r}),$
(13)

corresponding to the monomer density $\rho({\bf r})=N\psi({\bf r})^{2}$.
Furthermore, we introduced

$\kappa=\int d^{3}x\frac{N\psi({\bf x})^{4}}{\mu+N\psi({\bf
x})^{2}}~{}~{}~{}~{}{\rm and}~{}~{}~{}~{}\mu=(\beta v_{0}\xi^{2})^{-1}.$

Upon substitution of the saddle - point equations (see Appendix, Eq. B.12 \-
B.17) into (Eq. II.2), we remain with

$\begin{split}\mathcal{F}&=\frac{Nl^{2}}{6\beta}\int(\nabla\psi)^{2}d^{3}{\bf
x}-\frac{V}{2\beta\tilde{v}}\ln{\mu}+\frac{1}{2\beta}\bigg{(}\frac{\kappa}{1-\kappa}+\ln{(1-\kappa)}\bigg{)}+\\\
&+\frac{1}{2\beta\tilde{v}}\int d^{3}{\bf x}\ln{(\mu+N\psi^{2}({\bf
x}))}+\frac{\varepsilon\varepsilon_{0}}{2}\int d^{3}{\bf x}(\nabla\varphi({\bf
x}))^{2},\end{split}$ (14)

which is the form of the free energy that we will use in what follows.

###  II.4 No-salt regime and planar confinement geometry

Let us now consider the model polyelectrolyte chain confined between two
infinite plates of area $S$ separated by a distance $2al$ in the salt-free
case, corresponding to zero added salt concentration $c_{0}=0$. In other
words, $2a$ is the inter – plate distance in the units of the Kuhn segment
length $l$. We will further assume that the surface charge density is
described by

$\rho_{s}({\bf r})=\sigma(\delta(zl-al)+\delta(zl+al)),$ (15)

and that

$\displaystyle\varphi(z)$ $\displaystyle=$ $\displaystyle\varphi({\bf r})$
(16) $\displaystyle\Phi(z)$ $\displaystyle=$ $\displaystyle\psi({\bf
r})\sqrt{Sl},$ (17)

where due to electroneutrality $\sigma=\frac{\tau N}{2S}$ and $\tau=pe$. We
also introduce the dimensionless coordinates ${\bf r}=(xl,yl,zl)$ and
stipulate the normalization $\int_{-a}^{a}dz~{}\Phi^{2}(z)=1$. With these
definnitions the Poisson-Boltzmann equation (Eq. (12)) becomes simply the
Poisson equation in the form

$\frac{d^{2}}{dz^{2}}\varphi(z)=\frac{\sigma
l}{\varepsilon\varepsilon_{0}}\bigg{(}2\Phi^{2}(z)-\delta(z-a)-\delta(z+a)\bigg{)},$
(18)

and thus has an obvious solution

$\varphi(z)=\frac{\sigma
l}{\varepsilon\varepsilon_{0}}\bigg{[}\int_{-a}^{a}\\!\\!\\!|z-z^{\prime}|\Phi^{2}(z^{\prime})dz^{\prime}-m(z,a)\bigg{]}.$
(19)

with $m(z,a)\equiv\max\big{(}z-a,-a-z,0\big{)}$. In the same limit the Edwards
equation (Eq. (II.3)) can be cast into the following form

$\begin{split}\mathcal{E}_{0}\Phi(z)&=\bigg{[}-\frac{1}{6}\frac{d^{2}}{dz^{2}}-\beta\tau\varphi(z)+\frac{1}{2\tilde{v}}\frac{1}{\mu+\frac{2\sigma}{\tau
l}\Phi^{2}(z)}+\\\
&+\frac{\Phi^{2}(z)}{Sl}\frac{\kappa}{2(1-\kappa)^{2}}\frac{2\mu+\frac{2\sigma}{\tau
l}\Phi^{2}(z)}{\big{(}\mu+\frac{2\sigma}{\tau
l}\Phi^{2}(z)\big{)}^{2}}\bigg{]}\Phi(z),\end{split}$ (20)

where

$\kappa=\int_{-a}^{a}dz\frac{2\sigma\Phi^{4}(z)}{\mu\tau
l+2\sigma\Phi^{2}(z)}.$ (21)

As we consider the bounding plates to be of infinite extension, $S\to\infty$,
the Edwards equation is simplified and assumes the form

$\frac{1}{6}\Phi^{\prime\prime}(z)=-\mathcal{E}_{0}\Phi(z)-\beta\tau\varphi(z)\Phi(z)+\frac{1}{2\tilde{v}}\frac{\Phi(z)}{\mu+\frac{2\sigma}{\tau
l}\Phi^{2}(z)}.$ (22)

Following the linearization procedure for (Eq. (19)), suggested in Podgornik
(1991), we get for $z\in[0,a]$ an even simpler version of the Edwards equation
that can be written as

$\frac{1}{6}\Phi^{\prime\prime}(z)=-\bigg{(}\mathcal{E}_{0}+\frac{\beta\sigma\tau
l}{\varepsilon\varepsilon_{0}}z\bigg{)}\Phi(z)+\frac{1}{2\tilde{v}}\frac{\Phi(z)}{\mu+\frac{2\sigma}{\tau
l}\Phi^{2}(z)},$ (23)

coupled with normalization

$\begin{split}2\int_{0}^{a}dz~{}\Phi^{2}(z)&=1\end{split}$ (24)

and the boundary conditions for the ground

$\begin{split}\Phi^{\prime}(0)=0\;\;\;\Phi(a)&=0.\end{split}$ (25)

and the first excited states

$\begin{split}\Phi(0)=0\;\;\;\Phi(a)&=0.\end{split}$ (26)

Furthermore, we introduce the following notations: $\lambda_{B}={\beta
e^{2}}/{4\pi\varepsilon_{0}\varepsilon},\nu={2\sigma}/{\tau l}={N}/{Sl}$ and
$\alpha=12\pi p^{2}l^{2}\lambda_{B}$, where $\lambda_{B}$ is the Bjerrum
length and $\nu$ is the charge concentration. The final form of the Edwards
equation is then given in the form of a non-linear Schrödinger equation

$\Phi^{\prime\prime}(z)=-\bigg{(}6\mathcal{E}_{0}+\alpha\nu~{}z\bigg{)}\Phi(z)+\frac{3}{\tilde{v}}\frac{\Phi(z)}{\mu+\nu\Phi^{2}(z)},$
(27)

which is the equation that we will solve numerically in what follows.

###  II.5 Reduced Free Energy

Using (Eq. (14)) and the results from the previous section, we then get the
surface density of the free energy for the salt-free case in the following
form

$\begin{split}A&=\frac{\mathcal{F}}{S}=\frac{2\sigma}{\tau\beta}\mathcal{E}_{0}-\frac{\sigma}{\tilde{v}\beta\tau}\int_{-a}^{a}dz\frac{\Phi^{2}(z)}{\mu+\frac{2\sigma}{\tau
l}\Phi^{2}(z)}+\\\
&\frac{1}{2S\beta}\bigg{(}\frac{\kappa}{1-\kappa}+\ln{(1-\kappa)}\bigg{)}+\\\
&+\frac{l}{2\tilde{v}\beta}\int_{-a}^{a}dz\ln{\bigg{(}1+\frac{2\sigma}{\tau
l\mu}\Phi^{2}(z)\bigg{)}}+\\\
&\sigma\bigg{(}\varphi(a)+\int_{-a}^{a}\Phi^{2}(z)\varphi(z)dz\bigg{)}.\end{split}$
(28)

After taking the $S\to\infty$ limit and using the notations introduced in (Eq.
27) together with $\tilde{A}=A/k_{B}T\nu l$ we obtain:

$\begin{split}\tilde{A}&=\mathcal{E}_{0}+\frac{\alpha\nu}{12}a+\frac{\alpha\nu}{12}\int_{-a}^{a}dz\int_{-a}^{a}dz^{\prime}|z-z^{\prime}|\Phi^{2}(z^{\prime})\Phi^{2}(z)-\\\
&-\int_{-a}^{a}dz\frac{\Phi^{2}(z)}{2\tilde{v}\mu+2\tilde{v}\nu\Phi^{2}(z)}+\frac{1}{2\tilde{v}\nu}\int_{-a}^{a}dz\ln{\big{(}1+\frac{\nu}{\mu}\Phi^{2}(z)\big{)}}.\end{split}$
(29)

No additional analytic approximations are feasible and/or appropriate at this
point. The integrals indicated above can only be preformed numerically with
the numerical solutions of the non-linear Schrödinger equation (Eq. 27).

##  III Results

Figure 2:  Free Energy(upper panel) and pressure (bottom panel) without
disorder, in the pure state ($\xi=0.0$) as a function of $a$, one-half of the
inter-surface separation $2a$. The solid lines are the free energies
calculated numerically based on the Eq.(29). The red circles on the left pane
describe the results from Ref. Podgornik (1991).

First, we compare our numerical results with the previously reported
analytical results in Ref. Podgornik (1991), where the reduced free energy can
be shown to correspond exactly to (Eq. (29)) with $\xi=0$, including some
additional approximations. We can see from Fig.2 that these additional
approximations are irrelevant and the results are effectively the same. The
pressure between the two bounding surfaces, i.e., the force per unit area is
calculated from

$P=-\frac{1}{l}\frac{\partial\tilde{A}}{\partial a}.$ (30)

Using this equation we calculate numerically the pressure of the system that
is presented in Fig.2, and compared to the analytical results of Ref.
Podgornik (1991).

Figure 3:  Dimensionless monomer density profile $\Phi^{2}(z)$ for $a=0.8$.
Blue(solid) line corresponds to the pure system, $\xi=0$, while the disordered
system follows orange(dashed dotted) line - $\xi=10$, and red(dashed) line -
$\xi=100$.

The configurational behavior of the confined polyelectrolyte chain with
short–range quenched disorder interactions is governed by the interplay
between the inter–plate separation $2a$ and the standard deviation of disorder
$\xi$. The distribution of the dimensionless polyelectrolyte monomers
$\Phi^{2}(z)$ is presented in Figs. 3, 4, and 5. For the given value of $\xi$
the behavior of the monomer density distribution on increase of the
inter–plate separation is qualitatively the same both for the pure case,
corresponding to $\xi=0$, as well as the disordered cases, $\xi\neq 0$, up to
the critical inter-plate separation $a_{c}$: upon further increase of the
separation there is a phase transition from the density profile having two
maximums separated by the region of lower non-zero density to the two-maximum
density profile with complete depletion of density at the midpoint between the
two bounding surfaces. The critical separation $a_{c}$ depends on the
magnitude of the disorder dispersion $\xi$ as well as on the charge parameters
of the system. The main effect of the quenched sequence disorder at the
inter–plate separations corresponding to $a<a_{c}$ is an overall depletion of
the monomer density close to the bounding surfaces and accumulation at the
relative maximum located at the midpoint between the two bounding surfaces (
Fig. 3) or at the points of maximum in case of two-maximum density, leading to
an increase of the monomer density in the region of the density maximum ( Fig.
4), i.e., there is an additional localization of the polyelectrolyte chain on
increase of the disorder parameter $\xi^{2}$. However, for separations
corresponding to $a>a_{c}$, there is a complete depletion of the monomer at
the midpoint between the two bounding surfaces with preferential accumulation
actually at the positions of the relative maximum of the monomer density at
the left and right halves of the slit ( Fig. 5).

Figure 4:  Dimensionless mmonomer density profile $\Phi^{2}(z)$ for $a=1.5$,
still below the critical separation $a_{c}=1.525$. Blue(solid) line - pure
system, orange(dashed dotted) line - disordered one, $\xi=10$, red(dashed)
line - disordered one, $\xi=100$.

In the presence of disorder the reduced free energy (Eq. 29) is defined by the
eigenvalue $\mathcal{E}_{0}$, corresponding to the "wave function" $\Phi(z)$
of the non-linear Schrödinger equation (Eq. 27). The free energy behavior
corresponding to the ground as well as the first excited states is given in
Fig. 6, where the singularity of the first derivative of the free energy at
$a=a_{c}$ is clearly shown. Thus, the confined polyelectrolyte with
short–range sequence disorder exhibits a phase transition, accompanied by a
loss of the original ground state at the critical distance $a_{c}$. In fact at
the critical distance $a_{c}$ the ground state free energy actually coincides
with that of the first excited state.

Figure 5: Dimensionless monomer density profile below (blue solid line) and
above (orange dashed dotted line) the critical inter-surface separation
$a_{c}$. The phase transition of the chain results in a complete depletion of
the monomer density at the midpoint between the bounding surfaces, effectively
suppressing the polyelectrolyte bridging induced attraction between the
surfaces.

The sequence disorder phase transition at a critical inter-surface separation
$a=a_{c}$ is of a first order type, and is accompanied by an abrupt change of
the inter–plate pressure (Fig. 7), as well as the vanishing of the monomer
density at the midpoint between the surfaces (Fig. 5). The pressure behavior
for the pure ($\xi=0$) and disordered ($\xi=10$) cases (Fig. 7) are
qualitatively similar below the critical distance $a<a_{c}$, but diverge
sharply for $a>a_{c}$, when the inter-surface pressure behavior changes
drastically from the negative (attractive) pressure in the pure case to the
positive (repulsive) pressure in the disordered case. This furthermore implies
that above the critical distance $a_{c}$ the repulsion between the two
positively charged plates is no longer compensated by the polyelectrolyte
bridging attraction mediated by the negatively charged polymer chain. The
reason for this qualitative difference can be traced to a complete charge
separation of the polyelectrolyte density between the left and right halves of
the inter-surface region, promoted by an increased localization of the monomer
density as clearly shown in Fig. 5.

Figure 6: The free energy dependence on the one-half inter–plate distance $a$
for the ground and first excited states of the equation (27) for the disorder
parameter $\xi=10$. Orange dashed line - for the first excited state, blue
line - for the ground state. Greater than the critical distance $a_{c}\approx
1.525$ the original excited state becomes the new ground state.

##  IV Discussion

We developed a self-consistent field (SCF) formalism describing behavior of
the polyelectrolyte chain confined between two charged impenetrable bounding
surfaces. There is no added salt present in the system and the charges on the
polyelectrolyte chain are compensated solely by the charges on the surfaces.
The identity of the monomers along the polyelectrolyte chain is described by a
random variable, $\xi_{\tau}$, characterized by a Gaussian distribution of
zero mean and variance $\xi^{2}$. The quenched distribution of monomers along
the contour only affects the short range interactions, while the long-range
interactions are Coulombic, corresponding to a fixed charge per monomer.

Figure 7: Inter-surface pressure dependence on $a$. At the critical separation
$a_{c}=1.525$ there is a first order-like phase transition of the
polyelectrolyte chain characterized by a discontinuous jump in the interaction
pressure that changes from attractive to repulsive on increase of the inter-
surface separation beyond the critical value $a_{c}$.

To calculate the free energy and the density profile of the monomers we have
used the replica approach and for the replica - symmetric solution obtained a
phase transition of the first order at the critical distance $a_{c}$. We have
considered the polyelectrolyte chain in no-added salt regime, confined between
two oppositely charged, impenetrable planar surfaces with the specific aim to
asses the contribution of the quenched sequence disorder and the short-range
interactions between the monomers. The obtained results can be interpreted in
terms of a “quantum-like" Hamiltonian
${\cal{\hat{H}}}_{n}=-\frac{\ell^{2}}{6}\sum_{a}\nabla_{a}^{2}+\imath\sum_{a}{\hat{\rho}}_{a}({\bf
x}_{a})+\imath\sum_{a<b}{\hat{q}}_{ab}({\bf x}_{a},{\bf x}_{b})$, describing
an effective particle in the random external field
$\imath\sum_{a}{\hat{\rho}}_{a}({\bf
x}_{a})+\imath\sum_{a<b}{\hat{q}}_{ab}({\bf x}_{a},{\bf x}_{b})$. The
confinement of the polyelectrolyte chain inside the slit with charged walls
would thus correspond to the quantum-mechanical particle localized in a
potential well. Behavior of the polymer density $\Phi(z)^{2}$ presented in
Figs.(3,4) exhibits the additional localization induced by sequence disorder
and consequently by the random effective field in the “quantum-like"
Hamiltonian ${\cal{\hat{H}}}_{n}$, akin to the Anderson localization of the
polymer chain in a random medium, as was e.g. demonstrated in Ref. Shiferaw
and Goldschmidt (2001).

Besides the additional, disorder driven localization, the sequence disorder
also results in a phase transition of a first order type, accompanied by an
abrupt change of the inter–surface pressure and monomer density (Fig. 5), on
increase of the inter-surface separation beyond the critical separation
$a_{c}$. The pressure behavior for the pure and disordered cases are actually
qualitatively similar below the critical separations, $a<a_{c}$, but at
separations $a>a_{c}$, the inter-surface pressure behavior changes
drastically, flipping from a negative pressure in the pure case to a positive
one in the disordered case. Above the critical separation $a_{c}$ the
repulsion between the two positively charged plates ceases to be mediated by
the polyelectrolyte chain bridging mechanism, due to the complete separation
of the polyelectrolyte density between the left and right halves of the slit
as a consequence of the disorder-driven localization of the monomer density.

The work presented above introduces an interesting and important new feature
in the phenomenology of polyelectrolyte-driven interactions in charged
systems, where the polymer chains additionally exhibit short-range "chemical"
sequence disorder.

###### Acknowledgements.

The SCS of Armenia supported this work, grant No. 22AA-1C023 (V.S., Y.M.),
grant No. 21AG-1C038 (V.S.), grant No. 21T-1F307 (Y.M.). A.B. and Y.M.
acknowledge the partial financial support from Erasmus+ Project No.
(2023-1-SI01-KA171-HED-000138423). R.P. acknowledges funding from the Key
Project No. 12034019 of the National Natural Science Foundation of China.

## author declarations

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contribution

V. Stepanyan: Data curation (equal); Formal analysis (equal); Investigation
(equal); Methodology (equal); Soft- ware (equal); Validation (equal);
Visualization (equal).

A. Badasyan: Conceptualization (equal); Formal analysis (equal); Methodology
(equal); Validation (equal); Writing – original draft (equal).

V. Morozov: Conceptualization (equal); Formal analysis (equal); Methodology
(equal); Supervision (equal).

Y. Mamasakhlisov: Conceptualization (equal); Formal analysis (equal);
Investigation (equal); Methodology (equal); Supervision (equal); Writing –
original draft (equal)

R. Podgornik: Conceptualization (equal); Formal analysis (equal); Methodology
(equal); Supervision (equal); Validation (equal); Writing – original draft
(equal).

## Data Availability Statement

The data that support the findings of this study are available within the
article.

## References

  * Leckband and Israelachvili (2001) D. Leckband and J. Israelachvili, Quarterly Reviews of Biophysics 34, 105–267 (2001).
  * Rubinstein and Colby (2003) M. Rubinstein and R. H. Colby,  _Polymer Physics_ (Oxford University Press, New York, 2003).
  * Dobrynin and Rubinstein (2005) A. V. Dobrynin and M. Rubinstein, Progress in Polymer Science 30, 1049 (2005).
  * Muthukumar (2023) M. Muthukumar,  _Physics of Charged Macromolecules: Synthetic and Biological Systems_ (Cambridge University Press, 2023).
  * French  _et al._ (2010) R. H. French, V. A. Parsegian, R. Podgornik, R. F. Rajter, A. Jagota, J. Luo, D. Asthagiri, M. K. Chaudhury, Y.-m. Chiang, S. Granick, S. Kalinin, M. Kardar, R. Kjellander, D. C. Langreth, J. Lewis, S. Lustig, D. Wesolowski, J. S. Wettlaufer, W.-Y. Ching, M. Finnis, F. Houlihan, O. A. von Lilienfeld, C. J. van Oss, and T. Zemb, Rev. Mod. Phys. 82, 1887 (2010).
  * Elias (2012) H.-G. Elias,  _Macromolecules· 1: Volume 1: Structure and Properties_ (Springer Science & Business Media, 2012).
  * Smiatek (2020) J. Smiatek, Molecules 25 (2020), 10.3390/molecules25071661.
  * Lu  _et al._ (2015) B.-S. Lu, A. Naji, and R. Podgornik, The Journal of Chemical Physics 142, 214904 (2015), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4921892/13254056/214904_1_online.pdf .
  * Borkovec  _et al._ (2001) M. Borkovec, B. Jönsson, and G. J. Koper,  _Ionization processes and proton binding in polyprotic systems: Small molecules, proteins, interfaces, and polyelectrolytes_ (Springer, 2001).
  * Blossey and Podgornik (2022) R. Blossey and R. Podgornik, Europhysics Letters 139, 27002 (2022).
  * Bloomfield  _et al._ (2000) V. A. Bloomfield, D. M. Crothers, and I. Tinoco,  _Nucleic acids: structures, properties, and functions_ (University Science Books, 2000).
  * Finkelstein and Ptitsyn (2016) A. V. Finkelstein and O. Ptitsyn,  _Protein physics: a course of lectures_ (Elsevier, 2016).
  * Fallmann  _et al._ (2017) J. Fallmann, S. Will, J. Engelhardt, B. Grüning, R. Backofen, and P. F. Stadler, Journal of biotechnology 261, 97 (2017).
  * Twarock  _et al._ (2018) R. Twarock, R. J. Bingham, E. C. Dykeman, and P. G. Stockley, Current opinion in virology 31, 74 (2018).
  * Adhikari  _et al._ (2020) P. Adhikari, N. Li, M. Shin, N. F. Steinmetz, R. Twarock, R. Podgornik, and W.-Y. Ching, Physical Chemistry Chemical Physics 22, 18272 (2020).
  * Borukhov  _et al._ (1999) I. Borukhov, D. Andelman, and H. Orland, J. Phys. Chem. B 103, 5042 (1999).
  * Shafir  _et al._ (2003) A. Shafir, D. Andelman, and R. R. Netz, J. Chem. Phys. 119, 2355 (2003).
  * Siber and Podgornik (2008) A. Siber and R. Podgornik, Phys. Rev. E 78, 051915 (2008).
  * Podgornik (1991) R. Podgornik, J. Phys. Chem. 95, 5249 (1991).
  * Andelman (1995) D. Andelman,  _Structure and Dynamics of Membranes Chap. 12_ , edited by E. S. R. Lipowsky, Vol. 1B (Elsevier, Amsterdam, Handbook of Biological Physics, 1995).
  * Muthukumar (1989) M. Muthukumar, The Journal of Chemical Physics 90, 4594 (1989).
  * Hribar-Lee  _et al._ (2011) B. Hribar-Lee, M. Lukšič, and V. Vlachy, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 107, 14 (2011).
  * Bratko and Chakraborty (1995) D. Bratko and A. Chakraborty, Physical Review E 51, 5805 (1995).
  * Binder and Young (1986) K. Binder and A. P. Young, Rev. Mod. Phys. 58, 801 (1986).
  * Parisi (2023) G. Parisi, Rev. Mod. Phys. 95, 030501 (2023).
  * Garel  _et al._ (1997) T. Garel, H. Orland, and E. Pitard,  _Spin glasses and random fields_ , edited by A. Young (World Scientific, 1997).
  * Markovich  _et al._ (2021) T. Markovich, D. Andelman, and R. Podgornik, in  _Handbook of lipid membranes_ (CRC Press, 2021) pp. 99–128.
  * Shiferaw and Goldschmidt (2001) Y. Shiferaw and Y. Y. Goldschmidt, Phys. Rev. E 63, 051803 (2001).

Appendices

##  Appendix A The free energy functional

The quenched free energy (Eq. 6) can be estimated using the replica trick
Binder and Young (1986)

$-\beta{\cal F}=\lim_{n\rightarrow
0}\frac{\langle{Z\\{\xi\\}^{n}}\rangle_{{\cal P}}-1}{n},$ (A.1)

where $\beta=(k_{B}T)^{-1}$. We therefore need to calculate the $n$-replica
partition function as

$\displaystyle\langle{Z\\{\xi\\}^{n}}\rangle_{{\cal P}}=$
$\displaystyle\int\,{\cal D}{\bf
r}e^{-\frac{3}{2\ell^{2}}\sum_{a=1}^{n}\int_{0}^{N}\,d\tau(\partial_{\tau}{\bf
r}^{a}(\tau))^{2}-\beta\sum_{a=1}^{n}V\\{{\bf r}^{a}\\}}\times$
$\displaystyle\int\,{\cal D}\xi{\cal P}\\{\xi\\}e^{-\frac{\beta
v_{0}}{2}\int_{0}^{N}\,d\tau\int_{0}^{N}\,d\tau^{\prime}\xi_{\tau}\xi_{\tau^{\prime}}\sum_{a=1}^{n}\delta({\bf
r}^{a}(\tau)-{\bf r}^{a}(\tau^{\prime}))}.~{}~{}~{}~{}~{}$ (A.2)

The average over the distribution function (Eq. 1) in the equation (Eq. A) is
transformed as

$\displaystyle\int\,{\cal D}\xi{\cal P}\\{\xi\\}e^{-\frac{\beta
v_{0}}{2}\int_{0}^{N}\,d\tau\int_{0}^{N}\,d\tau^{\prime}\xi_{\tau}\xi_{\tau^{\prime}}\sum_{a=1}^{n}\delta({\bf
r}^{a}(\tau)-{\bf r}^{a}(\tau^{\prime}))}=$ $\displaystyle\int\,{\cal D}\xi
e^{-\frac{1}{2\xi^{2}}\int_{0}^{N}\,d\tau\xi_{\tau}^{2}}e^{-\frac{\beta
v_{0}}{2}\int\,d^{3}{\bf x}\sum_{a=1}^{n}m_{a}({\bf x})^{2}},~{}~{}~{}~{}$
(A.3)

where the field $m_{a}({\bf x})=\int_{0}^{N}\,d\tau\xi_{\tau}\delta({\bf
x}-{\bf r}^{a}(\tau))$, and $v_{0}$ is the interaction strength parameter. The
r.h.s. of the equation (Eq. A) is linearized using the Hubbard-Stratonovich
transformation as

$\displaystyle e^{-\frac{\beta v_{0}}{2}\int\,d^{3}{\bf
x}\sum_{a=1}^{n}m_{a}({\bf x})^{2}}=e^{-\frac{nV}{2}\ln(2\pi\beta v_{0})}$
$\displaystyle\int\,{\cal D}{\Pi}e^{-\frac{1}{2\beta
v_{0}}\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf
x})^{2}+\imath\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf x})m_{a}({\bf x})}.$
(A.4)

Insertion of (Eqs. A,A) into the (Eq. A) upon averaging over the $\\{\xi\\}$
variables yields the $n$-replica partition function as

$\displaystyle\langle{Z\\{\xi\\}^{n}}\rangle_{{\cal P}}\propto
e^{-\frac{nV}{2\tilde{v}}\ln(2\pi\beta v_{0})}\int\,{\cal
D}{\Pi}e^{-\frac{1}{2\beta v_{0}}\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf
x})^{2}}\times$ $\displaystyle\int\,{\cal D}{\bf
r}e^{-\frac{3}{2\ell^{2}}\sum_{a=1}^{n}\int_{0}^{N}\,d\tau(\partial_{\tau}{\bf
r}^{a}(\tau))^{2}-\beta\sum_{a=1}^{n}V\\{{\bf r}^{a}\\}}\times$ $\displaystyle
e^{-\frac{\xi^{2}}{2}\sum_{a,b}\int_{0}^{N}\,d\tau\Pi_{a}({\bf
r}^{a}(\tau))\Pi_{b}({\bf r}^{b}(\tau))}.$ (A.5)

For further consideration it is useful to introduce two order parameters,
viz., the inter-replica overlap $q_{ab}({\bf x},{\bf x}^{\prime})$ , with
$a<b$, and the density of $a$-th replica monomers $\rho_{a}({\bf x})$ as

$\displaystyle q_{ab}({\bf x},{\bf x}^{\prime})=\int_{0}^{N}\,d\tau\delta({\bf
x}-{\bf r}^{a}(\tau))\delta({\bf x}^{\prime}-{\bf r}^{b}(\tau))$
$\displaystyle\rho_{a}({\bf x})=\int_{0}^{N}\,d\tau\delta({\bf x}-{\bf
r}^{a}(\tau)).$ (A.6)

Using these order parameters (Eq. A), the $n$-replica partition function (Eq.
A) is transformed as

$\displaystyle\langle{Z\\{\xi\\}^{n}}\rangle_{{\cal P}}\propto\int\\!\\!{\cal
D}{\varphi}\\!\\!\int\\!\\!{\cal D}{c^{\pm}}\\!\\!\int\\!\\!{\cal
D}{\Pi}~{}e^{-\frac{1}{2\beta v_{0}}\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf
x})^{2}}\times$ $\displaystyle\int\,{\cal D}{\bf
r}e^{-\frac{3}{2\ell^{2}}\sum_{a=1}^{n}\int_{0}^{N}\,d\tau(\partial_{\tau}{\bf
r}^{a}(\tau))^{2}}\times$ $\displaystyle
e^{-\beta\sum_{a=1}^{n}W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})-\frac{\xi^{2}}{2}\sum_{a}\int\,d^{3}{\bf
x}\Pi_{a}({\bf x})^{2}\rho_{a}({\bf x})}\times$ $\displaystyle
e^{-\xi^{2}\sum_{a<b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf x}^{\prime}\Pi_{a}({\bf
x})\Pi_{b}({\bf x}^{\prime})q_{ab}({\bf x},{\bf x}^{\prime})},~{}~{}~{}~{}~{}$
(A.7)

where the electrostatic part of the free energy
$W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})$ is written as a sum of
electrostatic energies of polyelectrolyte segments, fixed charges residing on
the surface, confining polyelectrolyte and charges of salt ions, as well as
the entropies of the salt ions. As shown in Borukhov  _et al._ (1999); Shafir
_et al._ (2003); Siber and Podgornik (2008),

$\displaystyle W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})=\int\,d^{3}{\bf
x}\bigg{\\{}-\frac{\mathcal{E}\mathcal{E}_{0}}{2}(\nabla\varphi_{a}({\bf
x}))^{2}+$ $\displaystyle\varphi_{a}({\bf x})\bigg{[}ec_{a}^{+}({\bf
x})-ec_{a}^{-}({\bf x})-pe\rho_{a}({\bf x})+\rho_{surf}({\bf x})\bigg{]}+$
$\displaystyle\sum_{i=\pm}\bigg{[}k_{B}T(c_{a}^{i}({\bf x})\ln c_{a}^{i}({\bf
x})-$ $\displaystyle c_{a}^{i}({\bf x})-(c_{0}^{i}\ln
c_{0}^{i}-c_{0}^{i}))-\mu^{i}(c_{a}^{i}({\bf
x})-c_{0}^{i})\bigg{]}\bigg{\\}}.$ (A.8)

Here $c^{\pm}$ are the concentrations of $\pm$ monovalent salt ions, with
$c_{0}^{\pm}$ being their bulk concentrations, and $\mu^{\pm}$ their chemical
potentials, $\varepsilon\varepsilon_{0}$ is the permittivity of water, and
$\rho_{surf}$ is the charge distribution over the bounding surfaces, confining
polyelectrolyte.

Introducing the Lagrange multipliers ${\hat{q}}_{ab}({\bf x},{\bf
x}^{\prime})$ and the ${\hat{\rho}}_{a}({\bf x})$ conjugated with the order
parameters (A), the $n$-replica partition function reads Garel  _et al._
(1997)

$\displaystyle\langle{Z\\{\xi\\}^{n}}\rangle_{{\cal P}}\propto
e^{-\frac{nV}{2\tilde{v}}\ln(2\pi\beta v_{0})}\int\,{\cal D}{\Pi}\,{\cal
D}{\rho}\,{\cal D}{{\hat{\rho}}}\,{\cal D}{q}\,{\cal D}{\hat{q}}$
$\displaystyle
e^{G(\Pi,\rho,{\hat{\rho}},q,{\hat{q}})+\ln\zeta({\hat{\rho}},{\hat{q}})},$
(A.9)

where

$\displaystyle
G(\Pi,\rho,{\hat{\rho}},q,{\hat{q}},\varphi,c^{\pm})=-\frac{1}{2\beta
v_{0}}\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf x})^{2}-$
$\displaystyle-\beta\sum_{a=1}^{n}W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})-\frac{\xi^{2}}{2}\sum_{a}\int\,d^{3}{\bf
x}\Pi_{a}({\bf x})^{2}\rho_{a}({\bf x})-$
$\displaystyle-\xi^{2}\sum_{a<b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf
x}^{\prime}\Pi_{a}({\bf x})\Pi_{b}({\bf x}^{\prime})q_{ab}({\bf x},{\bf
x}^{\prime})+$ $\displaystyle+\imath\sum_{a}\int\,d^{3}{\bf x}\rho_{a}({\bf
x}){\hat{\rho}}_{a}({\bf x})+\imath\sum_{a<b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf
x}^{\prime}q_{ab}({\bf x},{\bf x}^{\prime}){\hat{q}}_{ab}({\bf x},{\bf
x}^{\prime})$ (A.10)

and

$\displaystyle\ln\zeta({\hat{\rho}},{\hat{q}})=-N\min_{\psi}\bigg{\\{}\int\,d^{3n}{\bf
x}\psi({\bf x}_{1},...,{\bf x}_{n})\times$
$\displaystyle\bigg{(}-\frac{\ell^{2}}{6}\sum_{a}\nabla_{a}^{2}+\imath\sum_{a}{\hat{\rho}}_{a}({\bf
x}_{a})+\imath\sum_{a<b}{\hat{q}}_{ab}({\bf x}_{a},{\bf x}_{b})\bigg{)}\times$
$\displaystyle\psi({\bf x}_{1},...,{\bf
x}_{n})-\mathcal{E}_{0}\bigg{(}\int\,d^{3n}{\bf x}\psi({\bf x}_{1},...,{\bf
x}_{n})^{2}-1\bigg{)}\bigg{\\}}.$ (A.11)

Here $\mathcal{E}_{0}$ is the ground state energy and $\psi({\bf
x}_{1},...,{\bf x}_{n})$ is the corresponding eigenfunction of the “quantum-
like" Hamiltonian
${\cal{\hat{H}}}_{n}=-\frac{\ell^{2}}{6}\sum_{a}\nabla_{a}^{2}+\imath\sum_{a}{\hat{\rho}}_{a}({\bf
x}_{a})+\imath\sum_{a<b}{\hat{q}}_{ab}({\bf x}_{a},{\bf x}_{b})$. Since
$v_{0}>0$, we have no trend of segregation between different kinds of monomers
and we expect large fluctuations of the fields $\\{\Pi\\}$. Thus, we can
consider the saddle-point approximation in the equation (Eq. A) over the
fields $\\{\rho,{\hat{\rho}},q,{\hat{q}}\\}$ only and to integrate over the
fields $\\{\Pi\\}$. Let us introduce
$g(\rho,{\hat{\rho}},q,{\hat{q}},\varphi,c^{\pm})=\ln\int\,{\cal
D}{\Pi}~{}e^{G(\Pi,\rho,{\hat{\rho}},q,{\hat{q}})}$, where

$\displaystyle\int\,{\cal
D}{\Pi}~{}e^{G(\Pi,\rho,{\hat{\rho}},q,{\hat{q}})}=e^{-\beta\sum_{a=1}^{n}W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})}\times$
$\displaystyle e^{\imath\sum_{a}\int\,d^{3}{\bf x}\rho_{a}({\bf
x}){\hat{\rho}}_{a}({\bf x})+\imath\sum_{a<b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf
x}^{\prime}q_{ab}({\bf x},{\bf x}^{\prime}){\hat{q}}_{ab}({\bf x},{\bf
x}^{\prime})}\times$ $\displaystyle\int\,{\cal
D}{\Pi}e^{-\frac{1}{2}\sum_{a,b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf
x}^{\prime}P_{ab}({\bf x},{\bf x}^{\prime})\Pi_{a}({\bf x})\Pi_{b}({\bf
x}^{\prime})},$ (A.12)

and $P_{ab}({\bf x},{\bf x}^{\prime})=\delta_{ab}\delta({\bf x}-{\bf
x}^{\prime})[\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf
x})]+\xi^{2}(1-\delta_{ab})q_{ab}({\bf x},{\bf x}^{\prime})$. Thus,
$g(\rho,{\hat{\rho}},q,{\hat{q}})$ can be obtained as

$\displaystyle
g(\rho,{\hat{\rho}},q,{\hat{q}},\varphi,c^{\pm})=-\frac{1}{2}\ln\det\frac{{\hat{P}}}{2\pi}-$
$\displaystyle\beta\sum_{a=1}^{n}W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})+\imath\sum_{a}\int\,d^{3}{\bf
x}\rho_{a}({\bf x}){\hat{\rho}}_{a}({\bf x})+$
$\displaystyle\imath\sum_{a<b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf
x}^{\prime}q_{ab}({\bf x},{\bf x}^{\prime}){\hat{q}}_{ab}({\bf x},{\bf
x}^{\prime}).$ (A.13)

Let us consider the term $\ln\det$ in the equation (Eq. A). These calculations
may be carried out by employing properties of block-matrices and the operator
algebra defined over the Hilbert space $\\{|{\bf x}\rangle\\}$. We shall use
the compact notation by defining the operators $\langle{\bf
x}|{\hat{A}}_{ab}|{\bf x}^{\prime}\rangle=\delta_{ab}\delta({\bf x}-{\bf
x}^{\prime})[\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf x})]$ as the $n\times
n$ operator matrix and $\langle{\bf x}|{\hat{P}}_{ab}|{\bf
x^{\prime}}\rangle=P_{ab}({\bf x},{\bf x^{\prime}})$ as an element of the
$n\times n$ operator matrix

${\hat{P}}={\hat{A}}+\xi^{2}{\hat{B}},$ (A.14)

where $\langle{\bf x}|{\hat{B}}_{ab}|{\bf
x^{\prime}}\rangle=(1-\delta_{ab})q_{ab}({\bf x},{\bf x}^{\prime})$ is the
$n\times n$ operator matrix. Therefore

$\displaystyle\ln\det{\hat{P}}={\rm Tr}\ln{\hat{P}}={\rm
Tr}\ln\bigg{\\{}{\hat{A}}+\xi^{2}{\hat{B}}\bigg{\\}}={\rm Tr}\ln{\hat{A}}+$
$\displaystyle{\rm
Tr}\ln\bigg{\\{}{\hat{1}}_{n}\otimes{\hat{1}}+\xi^{2}{\hat{A}}^{-1}{\hat{B}}\bigg{\\}},$
(A.15)

where $({\hat{1}}_{n})_{ab}=\delta_{ab}$ and $\langle{\bf x}|{\hat{1}}|{\bf
x}^{\prime}\rangle=\delta({\bf x}-{\bf x}^{\prime})$. Calculation of the first
term on the r.h.s. of the equation (Eq. A) is straightforward

$\displaystyle{\rm
Tr}\ln{\hat{A}}=\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k}\sum_{a}\int\,d^{3}{\bf
x}\langle{\bf x}|({\hat{A}}-{\hat{1}})^{k}_{aa}|{\bf x}\rangle=$
$\displaystyle=\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k}\sum_{a}\int\,d^{3}{\bf
x}\delta({\bf x}-{\bf x})[\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf
x})-1]^{k}=$ $\displaystyle=\sum_{a}\int\,d^{3}{\bf x}\delta({\bf x}-{\bf
x})\ln[\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf x})]\approx$
$\displaystyle\approx\frac{1}{\tilde{v}}\sum_{a}\int\,d^{3}{\bf
x}\ln\bigg{(}\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf x})\bigg{)},$ (A.16)

using $\delta(0)\approx\frac{1}{\tilde{v}}$, where $\tilde{v}$ is volume of
the monomer.

##  Appendix B Replica symmetric solution

In the replica-symmetric case we will consider only the Hartree-like replica-
symmetric wave functions

$\psi({\bf x}_{1},...,{\bf x}_{n})=\prod_{a=1}^{n}\psi({\bf x}_{a}),$ (B.1)

where $\psi({\bf x}_{a})$ is the one-replica wave function. In the replica
symmetric regime the expression (Eq. II.2) may be written (using (Eq. A)) as

$\displaystyle-\beta{\mathcal{F}}=\lim_{n\rightarrow
0}\frac{1}{n}\bigg{[}-\frac{n}{2\tilde{v}}\int\,d^{3}{\bf
x}\ln\bigg{(}1+\frac{\rho_{a}({\bf x})}{\mu}\bigg{)}-\frac{1}{2}{\rm
Tr}\ln\bigl{\\{}I+\hat{C}\bigr{\\}}-$ $\displaystyle-\beta
nW_{el}+in\langle\rho|\hat{\rho}\rangle+i\frac{n(n-1)}{2}\langle
q|\hat{q}\rangle+\ln\zeta\bigg{]},$ (B.2)

where $\hat{C}=\xi^{2}\hat{A}^{-1}\hat{B}$ and $\mu=(\beta
v_{0}\xi^{2})^{-1}$. With the replica symmetry we have
$\hat{C}=\hat{T}\otimes\hat{C^{\prime}}$ where $\hat{T}_{ab}=1-\delta_{ab}$
and $\langle
x|\hat{C^{\prime}}|x^{\prime}\rangle=\frac{q(x,x^{\prime})}{\mu+\rho(x)}$.
Taking the limit we get the reduced free energy

$\displaystyle f=-\beta{\mathcal{F}}=-\frac{1}{2\tilde{v}}\int\,d^{3}{\bf
x}\ln\bigg{(}1+\frac{\rho_{a}({\bf x})}{\mu}\bigg{)}-\frac{\Xi}{2}-$
$\displaystyle-\beta W_{el}+i\langle\rho|\hat{\rho}\rangle-\frac{i}{2}\langle
q|\hat{q}\rangle+\ln\zeta,$ (B.3)

where we defined

$\Xi\equiv\lim_{n\to 0}\frac{1}{n}{\rm
Tr}\ln\bigl{\\{}I+\hat{C}\bigr{\\}}=\lim_{n\to
0}\frac{1}{n}\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k}{\rm Tr}(\hat{C}^{k}).$
(B.4)

Thus, using ${\rm Tr}(\hat{T}\otimes\hat{C^{\prime}})={\rm Tr}(\hat{T}){\rm
Tr}(\hat{C^{\prime}})$, the terms of the series are written

$\displaystyle{\rm Tr}(\hat{C})^{k}=T_{k}\int
d^{3k}x\frac{q(x^{(1)},x^{(2)})q(x^{(2)},x^{(3)})...q(x^{(k)},x^{(1)})}{(\mu+\rho(x^{(1)}))(\mu+\rho(x^{(2)}))..(\mu+\rho(x^{(k)}))}$
$\displaystyle
T_{k}=\sum_{a_{1},..,a_{k}}\bigg{[}\prod_{i=2}^{k}\biggl{\\{}(1-\delta_{a_{i},a_{i-1}})\biggr{\\}}(1-\delta_{a_{k},a_{1}})\bigg{]}$
(B.5)

Where $T_{k}={\rm Tr}(\hat{T}^{k})$. To calculate $T_{k}$ first notice the
following recursive equation

$T_{k}=T_{k-1}(n-2)+T_{k-2}(n-1)$ (B.6)

with the initial conditions $T_{2}=n(n-1)$ and $T_{3}=(n-2)T_{2}$. Let us
define a generating function $P(x)=\sum_{i=0}^{\infty}x^{i}T_{i+2}$. Using

$\displaystyle x(n-2)P(x)=\sum_{i=1}^{\infty}x^{i}T_{i+1}(n-2)$ (B.7)
$\displaystyle x^{2}(n-1)P(x)=\sum_{i=2}^{\infty}x^{i}T_{i}(n-1)$

we can furthermore write

$P(x)\bigg{(}1-(n-2)x-(n-1)x^{2}\bigg{)}=T_{2}=n(n-1)$ (B.8)

and

$P(x)=\sum_{i=0}^{\infty}x^{i}\bigg{(}(n-1)^{i+2}-(-1)^{i+1}(n-1)\bigg{)}=\sum_{i=0}^{\infty}x^{i}T_{i+2}.$
(B.9)

Thus, we can take the limit

$\displaystyle\lim_{n\to 0}\frac{T_{k}}{n}=\lim_{n\to
0}\frac{(n-1)^{k}-(-1)^{k-1}(n-1)}{n}=~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle=\lim_{n\to
0}\bigg{(}k(n-1)^{k-1}-(-1)^{k-1}\bigg{)}=(-1)^{k-1}(k-1)~{}~{}~{}~{}~{}~{}$
(B.10)

Substituting (Eq. B) into $\Xi$ defined by (Eq. B.4) we get

$\Xi=\sum_{k=2}^{\infty}\frac{k-1}{k}\int d^{3k}x\frac{q({\bf x}^{(1)},{\bf
x}^{(2)})q({\bf x}^{(2)},{\bf x}^{(3)})...q({\bf x}^{(k)},{\bf
x}^{(1)})}{(\mu+\rho({\bf x}^{(1)}))(\mu+\rho({\bf x}^{(2)}))..(\mu+\rho({\bf
x}^{(k)}))}$ (B.11)

The saddle - point equations can be written in the following form

$0=\frac{\delta f}{\delta\hat{\rho}({\bf r})}\Rightarrow\rho({\bf
r})=N\psi({\bf r})^{2}$ (B.12) $0=\frac{\delta f}{\delta\hat{q}({\bf r},{\bf
r}^{\prime})}\Rightarrow q({\bf r},{\bf r}^{\prime})=N\psi({\bf
r})^{2}\psi({\bf r}^{\prime})^{2}$ (B.13) $0=\frac{\delta f}{\delta\rho({\bf
r})}\Rightarrow\mathrm{i}\hat{\rho}({\bf r})=-\beta pe\varphi({\bf
r})+\frac{1}{2\tilde{v}}\frac{1}{\mu+\rho(\bf
r)}+\frac{1}{2}\frac{\delta\Xi}{\delta\rho}$ (B.14) $0=\frac{\delta f}{\delta
q({\bf r},{\bf r}^{\prime})}\Rightarrow\mathrm{i}\hat{q}({\bf r},{\bf
r}^{\prime})=-\frac{\delta\Xi}{\delta q}$ (B.15) $0=\frac{\delta
f}{\delta\varphi({\bf
r})}\Rightarrow\varepsilon\varepsilon_{0}\nabla^{2}\varphi({\bf
r})=-(ec^{+}-ec^{-}-pe\rho({\bf r})+\rho_{s}({\bf r}))$ (B.16) $0=\frac{\delta
f}{\delta c^{\pm}({\bf r})}\Rightarrow c^{\pm}({\bf r})=c_{0}\exp({\mp\beta
e\varphi({\bf r})})$ (B.17) $\displaystyle 0=\frac{\delta f}{\delta\psi({\bf
r})}=-N\frac{\delta}{\delta\psi({\bf
r})}\bigg{\\{}\frac{\ell^{2}}{6}\langle\nabla\psi|\nabla\psi\rangle+$
$\displaystyle\mathrm{i}\langle\psi^{2}|\hat{\rho}\rangle--\frac{\mathrm{i}}{2}\langle\psi^{2}|\hat{q}|\psi^{2}\rangle-\mathcal{E}_{0}(\langle\psi|\psi\rangle-1)\bigg{\\}}=$
$\displaystyle-N\bigg{\\{}-\frac{\ell^{2}}{3}\nabla^{2}\psi({\bf
r})+2\mathrm{i}\psi({\bf r})\hat{\rho}({\bf r})-\mathrm{i}\psi({\bf
r})\int\,d^{3}{\bf x}\hat{q}({\bf r},{\bf x})\psi({\bf x})^{2}-$
$\displaystyle\mathrm{i}\psi({\bf r})\int\,d^{3}{\bf x}\hat{q}({\bf x},{\bf
r})\psi({\bf x})^{2}-2\mathcal{E}_{0}\psi({\bf r})\bigg{\\}}$ (B.18)

where

$\begin{split}&\frac{\delta\Xi}{\delta\rho({\bf
x})}=-\sum_{k=2}^{\infty}\frac{k-1}{(\mu+\rho({\bf x}))}\\\ &\int d^{3k-3}{\bf
x}\frac{q({\bf x},{\bf x}^{(1)})q({\bf x}^{(1)},{\bf x}^{(2)})...q({\bf
x}^{(k-1)},{\bf x})}{(\mu+\rho({\bf x}))(\mu+\rho({\bf
x}^{(1)}))..(\mu+\rho({\bf x}^{(k-1)}))}\end{split}$ (B.19)
$\begin{split}&\frac{\delta\Xi}{\delta q({\bf x},{\bf
x}^{\prime})}=\sum_{k=2}^{\infty}\frac{k-1}{(\mu+\rho({\bf x}))}\\\ &\int
d^{3k-6}{\bf x}\frac{q({\bf x}^{\prime},{\bf x}^{(1)})q({\bf x}^{(1)},{\bf
x}^{(2)})...q({\bf x}^{(k-2)},{\bf x})}{(\mu+\rho({\bf
x}^{\prime}))(\mu+\rho({\bf x}^{(1)}))..(\mu+\rho({\bf
x}^{(k-2)}))}\end{split}$ (B.20)

By further defining

$\kappa\equiv\int d^{3}x\frac{N\psi({\bf x})^{4}}{\mu+N\psi({\bf x})^{2}}$

and using (Eq. B.12, B.13) we get

$\frac{\delta\Xi}{\delta\rho({\bf
r})}=-\frac{\kappa}{(1-\kappa)^{2}}\frac{N\psi({\bf r})^{4}}{(\mu+N\psi({\bf
r})^{2})^{2}}$ (B.21) $\frac{\delta\Xi}{\delta q({\bf r},{\bf
r}^{\prime})}=\frac{1}{(1-\kappa)^{2}}\frac{N\psi({\bf r})^{2}\psi({\bf
r}^{\prime})^{2}}{(\mu+N\psi({\bf r})^{2})(\mu+N\psi({\bf r}^{\prime})^{2})}$
(B.22)

Substituting equations (Eqs. B.14,B.15) into the equation (Eq. B) we finally
obtain

$\displaystyle\frac{\ell^{2}}{6}\nabla^{2}\psi({\bf r})=\psi({\bf
r})\bigg{[}-\beta pe\varphi({\bf r})+\frac{1}{2\tilde{v}}\frac{1}{\mu+\rho(\bf
r)}+\frac{1}{2}\frac{\delta\Xi}{\delta\rho}\bigg{]}+$ $\displaystyle+\psi({\bf
r})\int\,d^{3}{\bf x}\frac{1}{2}\frac{\delta\Xi}{\delta q({\bf r},{\bf
x})}\psi({\bf x})^{2}+$ $\displaystyle\psi({\bf r})\int\,d^{3}{\bf
x}\frac{1}{2}\frac{\delta\Xi}{\delta q({\bf x},{\bf r})}\psi({\bf
x})^{2}--\mathcal{E}_{0}\psi({\bf r})$ (B.23)

Together with equations (Eqs. B.21,B.22) the latter equation is transformed
into

$\displaystyle\frac{\ell^{2}}{6}\nabla^{2}\psi({\bf r})=\psi({\bf
r})\bigg{[}-\beta pe\varphi({\bf
r})-\mathcal{E}_{0}+\frac{1}{2\tilde{v}}\frac{1}{\mu+N\psi({\bf r})^{2}}-$
$\displaystyle-\frac{\kappa}{(1-\kappa)^{2}}\frac{N\psi({\bf
r})^{4}}{2(\mu+N\psi({\bf r})^{2})^{2}}+$ $\displaystyle\int\,d^{3}{\bf
x}\frac{1}{(1-\kappa)^{2}}\frac{N\psi({\bf r})^{2}\psi({\bf
x})^{2}}{(\mu+N\psi({\bf r})^{2})(\mu+N\psi({\bf x})^{2}}\psi({\bf
x})^{2})\bigg{]}.$ (B.24)


# Quasi-stationary solutions in non-Hermitian systems

Cem Yuce<EMAIL_ADDRESS>Department of Physics, Faculty of Science,
Eskisehir Technical University, Turkey

###### Abstract

Eigenstates exhibit localization at an open edge in a non-Hermitian lattice
due to non-Hermitian skin effect. We here explore another interesting feature
of non-Hermitian skin effect and predict quasi-stationary solutions, which are
approximately time-independent. We show that the transition from such states
to eigenstates is dramatically non-perturbative. We consider a non-Hermitian
variant of the Su-Schrieffer-Heeger (SSH) model and predict non-topological
but robust quasi-stationary zero energy modes.

Non-Hermitan systems are ubiquitous in nature and has attracted great
attention in the last two decades. There are many ongoing experiments to check
theoretical predictions for non-Hermitian systems. Photonic lattices have
become the most popular area in the last decade in the community of non-
Hermitian physics. Introducing gain and loss to a photonics lattice is not the
only way to obtain a non-Hermitian lattice. Anisotropy of hopping amplitudes
between adjacent sites leads also to non-Hermiticity. This way of obtaining
non-Hermiticity is of particular interest as the spectrum strongly depends on
the boundary conditions. In the case of open boundary conditions (OBC), all
eigenstates are localized around the boundary and no extended modes in the
bulk of the lattice exist. This unique behavior is dubbed non-Hermitian skin
effect nhs1 , originally introduced in the study of topological phase of non-
Hermitian systems reva1 ; reva2 ; reva3 ; reva4 ; reva420 . This extreme
localization of all eigenstates implies that the bulk-boundary correspondence
fails nhs3 ; nhs2 ; nhs4 ; nhs5 ; nhs6 ; nhs7 ; nhs8 ; nhs9 ; cemAnnals ;
nhs10 ; nhs11 ; nhs12 ; nhs13 ; nhs14 ; nhs15 ; nhs16 ; nhs17 ; refek2 ;
refek3 ; refek4 ; refek5 ; refek6 ; refek7 ; refek8 ; refek9 ; refek10 ;
refek11 ; refek12 ; refek13 ; refek14  and some interesting consequences such
as light funneling funnel , in which any light excitation in the lattice moves
towards the open boundary and non-Bloch band theory of non-Hermitian systems
nhs16 . The non-Hermitian skin effect is not only theoretical curiosity as it
has recently been realized in some experiments nhs18 ; nhs19 ; nhs20 ; nhs21 .
Recently, we predict eigenstate clustering around exceptional points and show
that eigenstate localization due to non-Hermitian skin effect is a special
type of eigenstate clustering cluster .
In this paper, we predict another interesting consequence of non-Hermitian
skin effect. We find that very long and infinitely long lattices can have
dramatically different spectra as opposed to Hermitian systems. We explore
quasi-stationary states with continuum spectrum in a finite lattice with open
boundary conditions. They are not truly stationary but stay stationary in the
time scale of a typical experiment in photonics systems ourpaperc1 ;
ourpaperc2 . We discuss that the transition from the eigenstates to quasi-
stationary states is dramatically non-perturbative. These unique states can be
used in a variety of applications in non-Hermitian systems. We apply our idea
to a one dimensional non-Hermitian topological system and predict non-
topological but robust quasi-stationary zero energy modes.

##  I Quasi-stationary solutions

Consider a non-Hermitian tight-binding lattice with asymmetric hopping. The
complex field amplitude $\displaystyle{\psi_{j}}$ at site $j$ is governed

$\psi_{j+1}+{\gamma}~{}\psi_{j-1}=E~{}\psi_{j}$ (1)

where $\displaystyle{j=1,2,...,N}$ and $\displaystyle{N+1}$ is the total
number of lattice sites, $\displaystyle{0\leq\gamma<1}$ is the non-Hermitian
degree and $\displaystyle{E}$ are energy eigenvalues. This is a second order
linear homogeneous recurrence relation. The most general solution of Eq. (1)
for $\displaystyle{\gamma\neq 0}$ is then given by

$\psi_{j}=\Delta_{+}^{j}~{}c_{1}+\Delta_{-}^{j}~{}c_{2}$ (2)

where $\displaystyle{\Delta_{\mp}=\frac{E\mp\sqrt{E^{2}-4\gamma}}{2}}$ and
$\displaystyle{c_{1,2}}$ are arbitrary constants. On the other hand, Eq. (1)
becomes a first order linear homogeneous recurrence relation for
$\displaystyle{\gamma=0}$ and its solution is given by
$\displaystyle{\psi_{j}=E^{j-1}~{}\psi_{1}}$. The energy eigenvalues and the
arbitrary constants $\displaystyle{c_{1,2}}$ are determined from the boundary
conditions and normalization condition. Applying OBC,
$\displaystyle{\psi_{0}=\psi_{N+1}=0}$, yields $\displaystyle{c_{2}=-c_{1}}$
and

$\Delta_{+}^{N+1}-\Delta_{-}^{N+1}=0$ (3)

This equation can readily be solved:
$\displaystyle{\Delta_{+}=e^{i\frac{2n\pi}{N+1}}~{}\Delta_{-}}$, where
$\displaystyle{n=1,..,N}$. We then obtain the corresponding energy eigenvalues
for OBC

$E_{n}=2~{}\sqrt{\gamma}\cos\frac{n\pi}{N+1}$ (4)

Note that the energy eigenvalues for PBC are given by
$\displaystyle{E_{k}=e^{ik}+{\gamma}~{}e^{-ik}}$ corresponding to the extended
eigenstates $\displaystyle{\psi_{j}=\psi_{0}~{}e^{ikj}}$. Let us now compare
the spectra for OBC and PBC. Firstly, the eigenvalues are real for OBC while
they are complex for PBC. Secondly, an exceptional point occur for OBC at
$\displaystyle{\gamma=0}$ since all eigenstates coalesce to a localized
exceptional eigenstate with zero energy,
$\displaystyle{\psi_{j}=\psi_{1}~{}\delta_{j,1}}$. However, no exceptional
points occurs for PBC. It is interesting to consider a long lattice with a
small value of $\displaystyle{\gamma}$. In this case, the OBC and PBC lattices
differ only at the edges but such a slight difference has nonlocal and non-
perturbative consequences on both the eigenstates and eigenvalues This is
known as non-Hermitian skin effect.
We have derived Eq. (4) by assuming that $\displaystyle{N}$ is a finite
number. Surprisingly, there is another solution of (3) which appears only when
$\displaystyle{N\rightarrow\infty}$ (semi-infinite lattice for which energy
eigenvalues form a continuum band). The solution reads
$\displaystyle{\Delta_{\mp}<1}$ for which the OBC is also satisfied
$\displaystyle{\psi_{\infty}=0}$. If we solve the equation
$\displaystyle{\Delta_{\mp}<1}$, we obtain a continuum band inside an ellipse
in the complex plane

$\frac{E_{R}^{2}}{(1+\gamma)^{2}}+\frac{E_{I}^{2}}{(1-\gamma)^{2}}<1$ (5)

where $\displaystyle{E_{R},E_{I}}$ are real and imaginary parts of energy
eigenvalues, respectively. It is interesting to see that the PBC spectrum
formula $\displaystyle{E_{k}=e^{ik}+{\gamma}~{}e^{-ik}}$ is equivalent to
$\displaystyle{\frac{E_{R}^{2}}{(1+\gamma)^{2}}+\frac{E_{I}^{2}}{(1-\gamma)^{2}}=1}$.
This means that the PBC energy spectrum makes a closed loop enclosing the OBC
spectrum. We plot them in Fig.1 for $\displaystyle{\gamma=0}$ and
$\displaystyle{\gamma=0.6}$. The light blue region is for OBC while the red
curve is for PBC. The elliptic band is compressed by increasing
$\displaystyle{\gamma}$ and flattened to a line on the $E_{R}$ axis in the
Hermitian limit $\displaystyle{\gamma=0}$. This means that PBC and OBC spectra
for the semi-infinite lattice overlap only for the Hermitian case.
Let us also compare (4) and (5). The former one becomes
$\displaystyle{-2\sqrt{\gamma}<E<2\sqrt{\gamma}}$ in the limit
$\displaystyle{N\rightarrow\infty}$, which predicts a real valued continuum
band in between $\displaystyle{\mp 2\gamma}$. This band is quite different
from what Eq. (5) predicts. This shows that our solutions, which appear only
for infinitely long lattice are novel. Furthermore, this novel solution is
unique to non-Hermitian systems. In the Hermitian case
$\displaystyle{\gamma=1}$, the expressions (4) and (5) yield the same results:
$\displaystyle{-2<E<2}$, which is exactly what PBC predicts
$\displaystyle{E_{k}=2\cos{k}}$. As a result, we say that mathematically
extending the boundary of long non-Hermitian lattices to infinity may not be
trivial. In Hermitian system, such an extension is always trivial as it
slightly changes the spectrum. However, the spectrum of a very long lattice
can be dramatically different from the spectrum of the infinitely long lattice
in non-Hermitian systems. We stress that one can not obtain these novel
solutions directly by numerically solving the eigenvalue equation since the
corresponding matrix is infinite dimensional. To this end, we mention another
important difference between the finite and semi-infinite open lattices. The
exceptional point at $\displaystyle{\gamma=0}$, which exist for a finite
lattice disappears for the semi-infinite lattice. In other words, there exist
only one solution with zero energy
$\displaystyle{\psi_{j}=\psi_{1}\delta_{j,1}}$ for a finite lattice at
$\gamma=0$ while infinitely many solutions
$\displaystyle{\psi_{j}=E^{j-1}~{}\psi_{1}}$ with complex valued energy
$\displaystyle{|E|<1}$ for the semi-finite lattice at
$\displaystyle{\gamma=0}$.
We conclude that we use (4) for the spectrum for a finite lattice no matter
how long the lattice is, while we use (5) for the spectrum for the semi-
infinite lattice. The latter formula can’t be obtained from the former one in
the limit $\displaystyle{N\rightarrow\infty}$. A real physical system has a
finite number of lattice sites, so we should use the former solution.
Fortunately, the latter one can be utilized for a sufficiently long lattice.
Assume that $N$ is large. Then the right open boundary condition is not
exactly but approximately satisfied. For example, suppose
$\displaystyle{\gamma=0}$ and $\displaystyle{N=100}$ at which the right open
boundary condition is exactly satisfied only when $E=0$. However,
$\displaystyle{\psi_{101}\approx 3.9\times 10^{-31}}$ at
$\displaystyle{E=0.5}$ and $\displaystyle{\psi_{101}\approx 1.6\times
10^{-10}}$ at $\displaystyle{E=0.8}$. These are very small numbers and these
solutions can be considered as quasi-stationary solutions. They eventually
change their forms since they are not true eigenstates. Fortunately, these
quasi-stationary solutions stay almost stationary in the time scale of a
typical non-Hermitian experiment.
Let us now qualitatively discuss quasi-stationary solutions. They appear in a
sufficiently long and highly nonreciprocal latices where non-Hermitian skin
effect occurs. In such a lattice, all eigenstates exhibit localization at
either edge. If the localization occurs at the left edge, then we say that all
eigenstates have extremely low densities at the right edge of the sufficiently
long lattice. The key idea here is that the right boundary condition is
automatically satisfied to a very good approximation. This has an interesting
implication. Lifting the constraint at the right boundary leads to continuum
of energy values since the energy values don’t need to take some certain
values to satisfy the right boundary condition. They are not true eigenstates
but quasi-stationary solutions as OBC is approximately satisfied thanks to the
non-Hermitian skin effect.
Quasi-stationary solutions can also occur for non-Hermitian systems with gain
and loss provided that non-Hermitian skin effect occurs. Let us illustrate our
idea on an another model for further understanding of quasi-stationary
solutions. Consider a tight-binding non-Hermitian lattice with alternating
gain and loss. The field amplitude satisfies

$\psi_{j+1}+\gamma~{}\psi_{j-1}+i~{}(-1)^{j}~{}V_{0}~{}\psi_{j}=E~{}\psi_{j}$
(6)

where $\displaystyle{i}$ is the complex number, $j=1,2,..,N$, and the real
number $\displaystyle{V_{0}}$ is the gain/loss strength. Since the lattice is
open, $\displaystyle{\psi_{0}=\psi_{N+1}=0}$.
We start with the case $\displaystyle{\gamma=0}$ at which analytical solutions
are available. In this case, there are two exceptional points of order $N/2$
and the eigenstates coalesce to two eigenstates
$\displaystyle{\psi_{j,+}=\\{1,0,0,...,0\\}^{T}}$ and
$\displaystyle{\psi_{j,-}=\frac{1}{\sqrt{5}}\\{1,-2,0,...,0\\}^{T}}$
corresponding to the eigenvalues $E_{\mp}={\mp}~{}i~{}V_{0}$. Let us find the
unnormalized novel solutions for the semi-infinite lattice recursively. Assume
$\displaystyle{\psi_{1}=1}$. Then $\displaystyle{\psi_{2}=E+iV_{0}}$,
$\displaystyle{\psi_{3}=(E+iV_{0})(E-iV_{0})}$ and so on. The field amplitude
at infinity becomes

$\psi_{\infty}=\left((E+iV_{0})(E-iV_{0})\right)^{\infty}$ (7)

Note that this recursive solution is reduced to the eigenstates if we replace
$E=E_{\mp}={\mp}iV_{0}$. To obtain the energy eigenvalues for the semi-
infinite lattice, we require $\displaystyle{|E^{2}+V_{0}^{2}|<1}$. This
condition is satisfied if $E$ is inside a continuum band in the complex plane

$4E_{R}^{2}E_{I}^{2}+(E_{R}^{2}-E_{I}^{2}+V_{0}^{2})^{2}<1$ (8)

The continuum band is inside a circle of radius $1$ in complex plane at
$V_{0}=0$. As $V_{0}$ is increased, the continuum band is elongated on $E_{I}$
axis until $V_{0}=1$ at which the band looks like the $8$ shape. For
$V_{0}>1$, the band is divided into two daughter bands, well separated in
$E_{I}$ axis. To study whether PBC spectrum encloses the OBC bands, we perform
numerical calculation for various values of $V_{0}$ for very large number of
$N$. We numerically find that PBC spectrum satisfy
$\displaystyle{4E_{R}^{2}E_{I}^{2}+(E_{R}^{2}-E_{I}^{2}+V_{0}^{2})^{2}=1}$,
which proves our statement. Note that the two exceptional point for a finite
lattice disappear for the semi-infinitely long lattice.
We have so far studied the fully asymmetric hopping at which analytical
solutions are available. Consider now the case with $\displaystyle{\gamma\neq
0}$. In this case, obtaining the spectrum analytically is challenging. Here,
we follow another approach to find the continuum band. We start to obtain the
PBC spectrum numerically for a large lattice (suppose $N=10^{5}$). The PBC
spectrum is a closed loop in complex domain and it is in fact supposed to
determine the boundary for the continuum band. To check this statement, we
arbitrarily choose many $E$ values inside or outside the closed loop. We make
iterations for Eq. (6) with the initial value $\displaystyle{\psi_{1}=1}$ at
these arbitrary values of $E$. We numerically see that
$\displaystyle{\psi_{j}}$ goes to zero at the edge if $E$ is inside the closed
loop, while it diverges $E$ is outside the closed loop, This proves our above
statement.

Figure 1:  The plots show the elliptic continuum bands in the complex plane
for the semi-infinite lattice at $\displaystyle{\gamma=0}$ (a) and
$\displaystyle{\gamma=0.6}$ (b). The red loop is for PBC and the blue shaded
are is for OBC. The PBC spectrum encloses the OBC spectrum. The band becomes a
line on $E_{R}$ axis in the Hermitian limit $\displaystyle{\gamma=0}$. In this
case, the PBC and OBC spectra overlap.

##  II Non-topological robust zero-energy edge states

We have discussed that quasi-stationary states for a finite lattice can be
quite different from the eigenstates. Therefore quasi-stationary states bring
us new degrees of freedom in non-Hermitian systems. Here we present an example
where zero-energy quasi-stationary states appear in
$\it{topologically~{}trivial}$ region of a 1D non-Hermitian system. We show
that they are surprisingly immune to the same kind of disorder as topological
zero-energy eigenstates.
Let us consider a non-reciprocal extension of the well-known SSH model, which
is a one dimensional dimerized tight-binding lattice with non-reciprocal
alternating hopping amplitudes. The Hamiltonian reads

$\hat{H}=\sum_{n=1}^{N-1}T_{n}|n><n+1|+\gamma~{}T_{n-1}|n+1><n|$ (9)

where $\displaystyle{T_{n}=1+(-1)^{n}\delta}$ are the alternating hopping
amplitudes with $\displaystyle{-1<\delta<1}$ and
$\displaystyle{0\leq\gamma<1}$ is the non-Hermitian degree and the total
number of lattice sites $\displaystyle{N}$ is assumed to be an even number.
Note that non-Hermitian skin effect occurs and all eigenstate become localized
at the left edge when $\displaystyle{\gamma}$ is a small number.
In the Hermitian system, $\displaystyle{\gamma=1}$, it is well known that the
point $\displaystyle{\delta=0}$ is the topological phase transition point and
topological zero energy modes appear when $\displaystyle{\delta>0}$. In the
non-Hermitian case, two topological zero-energy modes are still available.
However, both of them as well as bulk eigenstates are localized around the
left edge due to the non-Hermitian skin effect. In our recent paper cemAnnals
, we showed that a zero-energy mode at the right edge can still be available,
but it tends to move rapidly to the left edge.
One can obtain zero-energy quasi-stationary states by solving the
corresponding recursion equation for the above Hamiltonian. They are given by

$\psi_{E=0}=A\sum_{n=1}^{N/2}\left(\frac{-\gamma~{}(1-\delta)}{~{}~{}~{}~{}(1+\delta)}\right)^{n-1}|2n-1>$
(10)

where $\displaystyle{A}$ is the normalization constant and the density is non-
zero only for odd numbers of lattice sites. This zero energy solution
satisfies the OBC for the semi-infinite lattice if
$\displaystyle{\left|\frac{\gamma~{}(1-\delta)}{~{}~{}(1+\delta)}\right|<1}$.
Then we obtain the domain for the existence of the zero-energy modes

$-\frac{1-\gamma}{1+\gamma}<\delta<1$ (11)

This solution is exact for the semi-infinite lattice while it is quasi-
stationary for a sufficiently long lattice. We recover the topological zero-
energy modes when $\displaystyle{0<\delta<1}$ and we further predict non-
topological zero-energy quasi-stationary states when
$\displaystyle{-\frac{1-\gamma}{1+\gamma}<\delta\leq 0}$. In Fig. 2, we plot
the energy eigenvalues as a function of $\displaystyle{\delta}$ at
$\displaystyle{\gamma=0.2}$. As can be seen, no zero-energy modes (in blue)
appear in the topologically trivial region $\displaystyle{\delta<0}$. However,
zero-energy quasi-stationary states (in red) exist in between
$\displaystyle{-0.67<\delta<0}$. They are almost stationary for a long time.
For example, we numerically find that the quasi-stationary zero energy state
at $\displaystyle{\delta=-0.4}$ is almost constant up to
$\displaystyle{t=23}$.
A question arises. Are these zero-energy quasi-stationary modes in the
topologically trivial region robust? Let us study robustness of these non-
topological zero-energy quasi-stationary modes against weak hopping amplitude
disorder, which can be introduced in our system by having randomized weak
forward hopping amplitudes $\displaystyle{t_{n}\rightarrow t_{n}+\delta
t_{n}}$ and $\displaystyle{\gamma\rightarrow\gamma+\delta\gamma_{n}}$, where
$\displaystyle{|\delta t_{n}|<<1}$ and
$\displaystyle{|\delta\gamma_{n}|<<\gamma}$ are real-valued random set of
constants. We numerically see that zero-energy quasi-stationary modes survive
in the presence of the disorder, which shows that they are robust. The
topological modes and quasi-stationary modes are robust against the same kind
of disorder. We are then tempted to say that zero energy quasi-stationary
modes might be topological, too. This notion needs for a further study as
topological invariants introduced in the literature so far does not predict
the region (11). This is an open problem.

Figure 2:  The plot shows the energy eigenvalues (in blue) as a function of
$\displaystyle{\delta}$ at $\displaystyle{\gamma=0.2}$ and
$\displaystyle{N=40}$. There exists no topological zero modes in topologically
trivial region $\displaystyle{\delta<0}$. However, zero energy quasi-
stationary states (in red) exist in topologically trivial region. These zero-
energy quasi-stationary states are robust to weak hopping amplitude disorder
as they are for topological zero modes for $\displaystyle{\delta>0}$.

##  III Conclusion

In this paper, we predict novel quasi-stationary solutions that are non-
perturbatively different from eigenstates but preserve their forms up to a
large time. Eigenvalue equations are not compatible with finding such
solutions since they normally appear for infinitely long lattices. They have
no analogue in Hermitian systems as they can only be seen in non-Hermitian
systems. We show that very long and infinitely long lattices can have
dramatically different spectra in non-Hermitian systems as opposed to
Hermitian system. We discuss that the non-Hermitian skin effect is responsible
for the existence of quasi-stationary solutions. We show that PBC spectrum
determines the boundary of OBC spectrum in complex plane.
We apply our idea to a topological system and find non-topological but robust
zero-energy quasi-stationary states. They are robust against the same kind of
disorder as topological zero-energy eigenstates.
Our formalism is constructed in 1D and higher dimensional generalization of
quasi-stationary solutions is straightforward. We think that they can find
some applications in non-Hermitian systems.
This study is supported by Eskisehir Technical University Scientific Research
Projects Commission under the grant no. 20ADP168.

## References

  * (1) Shunyu Yao and Zhong Wang, Phys. Rev. Lett. 121, 086803 (2018).
  * (2) Ananya Ghatak, and Tanmoy Das, J. Phys. Condens. Matter 31, 263001 (2019).
  * (3) V. M. Martinez Alvarez, J. E. Barrios Vargas, M. Berdakin, L. E. F. Foa Torres, Eur. Phys. J. Special Topics 227, 1295 (2018).
  * (4) Emil J. Bergholtz, Jan Carl Budich, Flore K. Kunst, arXiv:1912.10048 (2019).
  * (5) Luis E F Foa Torres, J. Phys. Mater. 3 014002 (2020).
  * (6) C. Yuce, Phys. Lett. A 379, 1213 (2015).
  * (7) C. Yuce, Phys. Lett. A 384, 126094 (2020).
  * (8) Shunyu Yao, Fei Song, and Zhong Wang, Phys. Rev. Lett. 121, 136802 (2018).
  * (9) Fei Song, Shunyu Yao, and Zhong Wang, Phys. Rev. Lett. 123, 170401 (2019).
  * (10) P. Wang, L. Jin, and Z. Song, Phys. Rev. A 99, 062112 (2019).
  * (11) Kazuki Yokomizo and Shuichi Murakami, Phys. Rev. Lett. 123, 066404 (2019).
  * (12) Motohiko Ezawa, Phys. Rev. B 99, 121411(R) (2019).
  * (13) Ken-Ichiro Imura and Yositake Takane, Phys. Rev. B 100, 165430 (2019).
  * (14) Nianzu Fu, Ziwei Fu, Huaiyuan Zhang, Qing Liao, Dong Zhao and Shaolin Ke , Optical and Quantum Electronics 52, 61 (2020).
  * (15) C Yuce, Ann. Phys. 415, 168098 (2020).
  * (16) Linhu Li, Ching Hua Lee, Sen Mu, Jiangbin Gong, Nat. Commun. 11, 5491 (2020).
  * (17) Ching Hua Lee and Ronny Thomale Phys. Rev. B 99, 201103(R) (2019).
  * (18) Ching Hua Lee, Linhu Li, and Jiangbin Gong, Phys. Rev. Lett. 123, 016805 (2019).
  * (19) Stefano Longhi, Phys. Rev. Research 1, 023013 (2019).
  * (20) Hui Jiang, Li-Jun Lang, Chao Yang, Shi-Liang Zhu, and Shu Chen, Phys. Rev. B 100, 054301 (2019).
  * (21) Jong Yeon Lee, Junyeong Ahn, Hengyun Zhou, and Ashvin Vishwanath, Phys. Rev. Lett. 123, 206404 (2019).
  * (22) Flore K. Kunst, Guido van Miert, and Emil J. Bergholtz, Phys. Rev. B 99, 085427 (2019).
  * (23) Yifei Yi, Zhesen Yang, arXiv e-prints, arXiv:2003.02219 (2020).
  * (24) Zongping Gong, Yuto Ashida, Kohei Kawabata, Kazuaki Takasan, Sho Higashikawa, and Masahito Ueda, Phys. Rev. X 8, 031079 (2018).
  * (25) Nobuyuki Okuma, Kohei Kawabata, Ken Shiozaki, and Masatoshi Sato, Phys. Rev. Lett. 124, 086801 (2020).
  * (26) C. Yuce, H. Ramezani, Phys. Rev. A 100, 032102 (2019).
  * (27) Ya-Nan Zhang, Shuang Xu, Hao-Di Liu and Xue-Xi Yi 60, 355 (2021).
  * (28) Ziwei Fu, Nianzu Fu, Huaiyuan Zhang, Zhe Wang, Dong Zhao and Shaolin Ke, Appl. Sci. 10, 3425 (2020)
  * (29) Shuai Li, Min Liu, Fuli Li and Bo Liu, Phys. Scr. 96, 015402 (2021).
  * (30) Zhi-Xu Zhang, Rong Huang, Lu Qi, Yan Xing, Zhan-Jun Zhang, and Hong-Fu Wang, Ann. Phys. (Berlin) 533, 2000272(2021.
  * (31) Yanxia Liu and Shu Chen, Phys. Rev. B 102, 075404 (2020).
  * (32) Boran Zhou, Rui Wang, and Baigeng Wang, Phys. Rev. B 102, 205116 (2020).
  * (33) Ling-Zhi Tang, Guo-Qing Zhang, Ling-Feng Zhang, Dan-Wei Zhang, arXiv:2101.05505 (2021).
  * (34) Thomas Benjamin Smith and Alessandro Principi, Physica E 126, 114423 (2021).
  * (35) Rong Huang, Yu Yan, Zhi-Xu Zhang, Lu Qi, Hong-Fu Wang and Shou Zhang, Commun. Theor. Phys. 72, 105101 (2020).
  * (36) Qi-Bo Zeng, Yan-Bin Yang, and Rong Lu, Phys. Rev. B 101, 125418 (2020).
  * (37) S. Weidemann, et. al., Science 368, 311 (2020).
  * (38) Xueyi Zhu, Huaiqiang Wang, Samit Kumar Gupta, Haijun Zhang, Biye Xie, Minghui Lu, and Yanfeng Chen, Phys. Rev. Research 2, 013280 (2020).
  * (39) Tobias Helbig, Tobias Hofmann, et. al., Nat. Phys. 16, 747 (2020).
  * (40) Lei Xiao, Tianshu Deng, Kunkun Wang, Gaoyan Zhu, Zhong Wang, Wei Yi, and Peng Xue, Nat. Phys. 16, 761 (2020).
  * (41) Ananya Ghatak, Martin Brandenbourger, Jasper Van Wezel, and Corentin Coulais, Proc. Natl. Ac. Sc. 117, 29561 (2020).
  * (42) C Yuce, Phys. Rev. A 102, 032203 (2020).
  * (43) C. Yuce, H. Ramezani, Opt. Lett. 46, 765 (2021).
  * (44) C. Yuce, H. Ramezani, arXiv:2009.12880 (2020).