# Extreme statistics of the excitations in the random transverse Ising chain

István A. Kovács<EMAIL_ADDRESS>Department of Physics and
Astronomy, Northwestern University, Evanston, IL 60208-3112, USA  Tamás Pető
<EMAIL_ADDRESS>Institute of Theoretical Physics, Szeged University,
H-6720 Szeged, Hungary  Ferenc Iglói<EMAIL_ADDRESS>Wigner Research
Centre for Physics, Institute for Solid State Physics and Optics, H-1525
Budapest, P.O. Box 49, Hungary  Institute of Theoretical Physics, Szeged
University, H-6720 Szeged, Hungary

###### Abstract

In random quantum magnets, like the random transverse Ising chain, the low
energy excitations are localized in rare regions and there are only weak
correlations between them. It is a fascinating question whether these
correlations are completely irrelevant in the sense of the renormalization
group. To answer this question, we calculate the distribution of the
excitation energy of the random transverse Ising chain in the disordered
Griffiths phase with high numerical precision by the strong disorder
renormalization group method and - for shorter chains - by free-fermion
techniques. Asymptotically, the two methods give identical results, which are
well fitted by the Fréchet limit law of the extremes of independent and
identically distributed random numbers. Given the finite size corrections, the
two numerical methods give very similar results, but they differ from the
correction term for uncorrelated random variables. This fact shows that the
weak correlations between low-energy excitations in random quantum magnets are
not entirely irrelevant.

##  I Introduction

Many-body systems in the presence of quenched disorder have unusual dynamical
properties due to rare-region effects. In these systems - due to extreme
fluctuations of strong couplings - domains are formed, which can remain
locally ordered even in the paramagnetic phase. The relaxation time, $\tau$,
associated with turning the spins in such domains can be extremely large and
it has no upper limit in the thermodynamic limit. This type of Griffiths
singularities are responsible for non-analytical behaviour of several physical
quantities (susceptibility, specific heat, auto-correlation function) in the
so called Griffiths phase, which is an extended part of the paramagnetic
phasegriffiths .

In random quantum systems Griffiths-effects are stronger than in classical
ones, which is manifested in power-law decay of the auto-correlation function,
as well as power-law singularities of the susceptibility and the specific heat
at low temperaturesmccoy ; vojta . In random quantum magnets with discrete
symmetry, such as in the random transverse Ising model (RTIM) the low-energy
excitations are localised and their properties can be successfully studied by
the so called strong disorder renormalization group (SDRG) methodim . As
initiated be Ma, Dasgupta and Humdh  and further developed by Fisherfisher
the SDRG is a local renormalization method, in which quantum and disorder
fluctuations are treated at the same time and degrees of freedom with a large
excitation energy are successively eliminated. In these random quantum magnets
the SDRG method is expected to provide asymptotically exact results not only
at the critical point, the properties of which are governed by an infinite
disorder fixed pointdaniel_review , but in the Griffiths phase as well, at
least for the dynamical singularities.

In phenomenological descriptions it is often assumed that the localized
excitations in random quantum magnets are independentthill_huse . For example
distribution of low-energy excitations in these systems are well approximated
by the Fréchet distributionjli , which represents the limit law of the
extremes of independent and identically distributed (iid) random numbers.
Recently extreme value statistics (EVS) has been applied to several problems,
we can mention earthquakes, tsunamis, extreme flooding, big wildfire, extremes
of climate, stock market risks in finance, sport records, etcfisher_tippett ;
gumbel ; gnedenko ; weibull ; galambos ; ev . In a mathematical point of view
complete understanding of extreme value statistics is known for uncorrelated
random numbers, in which case also the convergence to the limit laws is
derived by mathematicalhaan  and renormalization groupprl ; pre ; bertin ;
zarfaty  (RG) methods. The iid limit distributions also apply to weakly
correlated random numbers, while new types of limit distributions appear for
the strongly correlated casesmps .

It is an intriguing problem to what extent the localised excitations in random
quantum magnets are independent? Whether the weak correlations between the
rare regions are completely irrelevant or these are manifested in some
effects, such as in the form of finite-size corrections? In this paper we make
an effort to answer this question and consider the paradigmatic model, the
RTIM in one dimension (1D) and calculate numerically the distribution of the
first energy gap by the asymptotically exact SDRG method with high accuracy.
For moderate $L$ system sizes, we use also free-fermion techniquesJW ; lieb61
; pfeuty79  to calculate the gaps of the random samples exactly. The
distribution of the gaps for finite $L$ are compared with an appropriate
Fréchet distribution and their difference is analysed through finite-size
scaling. To check the potential form of numerical errors the same type of
numerical test is repeated for iid random numbers, too.

The rest of the paper is organised in the following way. In Sec.II the RTIM
model is introduced and in Sec.III the methods to calculate its excitation
energy are described. In Sec.IV the distributions of gaps of the RTIM are
calculated and the finite-size corrections are compared with the analytical
results of uncorrelated variables. To test the numerical accuracy we repeat
this analysis for uncorrelated Kesten variableskesten73 . Our paper concludes
with a discussion in the final section. Related results about the distribution
of extremes of uncorrelated variables are recapitulated in the Appendix.

##  II Model and known results

Here, we consider the RTIM in 1D defined by the Hamiltonian

$\hat{H}=-\sum_{i=1}^{L-1}J_{i}\sigma_{i}^{z}\sigma_{i+1}^{z}-\sum_{i=1}^{L}h_{i}\sigma^{x}_{i}\;,$
(1)

in terms of the $\sigma_{i}^{x,z}$ Pauli matrices at site $i$ and the nearest
neighbour couplings, $J_{i}>0$ and the transverse fields, $h_{i}>0$, are taken
from the distributions, $\pi_{1}(J)$ and $\pi_{2}(h)$, respectively. Generally
we use open boundary conditions (OBC-s) and work at zero temperature, $T=0$.

In the thermodynamic limit, $L\to\infty$, the control-parameter is defined
asfisher :

$\delta=\frac{[\ln h]_{\rm av}-[\ln J]_{\rm av}}{\rm{var}(h)+\rm{var}(J)}\;,$
(2)

where $[\cdots]_{\rm av}$ denotes averaging over quenched disorder and
$\rm{var}(x)$ stands for the variance of $x$. For $\delta<0$ ($\delta>0$) the
system is in the ordered ferromagnetic (disordered paramagnetic) phase and at
$\delta=0$ there is a random quantum critical point. According to SDRG
calculationsfisher ; ijl01 ; Igloi  and numerical resultsyoung_rieger ;
bigpaper  the critical behaviour of the system is controlled by an infinite
disorder fixed pointdaniel_review . For example, the energy scale, defined by
the lowest gap, $\varepsilon$, and the length-scale are related as

$\ln\varepsilon\sim L^{1/2},\quad\delta=0\;.$ (3)

In the paramagnetic phase this relation is in a power-law form

$\varepsilon\sim L^{-z},\quad\delta>0\;,$ (4)

which is due to Griffiths singularities. Here, the dynamical exponent,
$z=z(\delta)$, depends on the distance from the critical point and is given by
the positive root of the equation Z ; ijl01 ; Igloi :

$\left[\left(\frac{J}{h}\right)^{1/z}\right]_{\rm av}=1\;,$ (5)

which in the vicinity of the critical point diverges as: $z\approx
1/(2\delta)$.

The distribution of the first gap has been calculated analytically through the
SDRG methodfisher_young :

$P_{L}(\varepsilon;z)=\frac{1}{z}u^{1/z-1}\exp\left(-u^{1/z}\right)\;,$ (6)

in terms of $u=u_{0}\varepsilon(L/\xi)^{z}$, $\xi$ being the correlation
length and $u_{0}$ is a constant, defined by the standardisation. This
relation is valid for $L\gg\xi$ and for $\delta\ll 1$, in which limit
$\xi\sim\delta^{-2}$. We observe that Eq.(6) is just the Fréchet distribution,
see in Eq.(26).

An approximate form of the distribution of the low-energy excitations can be
obtained from the assumption that these excitations are localised and are due
to extreme fluctuations of say $n\ll L$ consecutive strong bondsthill_huse .
The probability to find such a strongly connected cluster in the system is
given by: $P_{L}(n)\sim L\exp(-an)$. At the same time the excitation energy
due to such a cluster is exponentially small: $\varepsilon\sim\exp(-bn)$.
Combining these expressions we obtain $P_{L}(\ln\varepsilon)\sim
L\varepsilon^{1/z}$, with $z=b/a$ being the dynamical exponent in agreement
with Eq.(4). Then, the cumulated distribution of the relaxation times,
$\tau\sim 1/\varepsilon$ is obtained in this approach:

$\mu(\tau)\approx 1-A\tau^{-1/z},\quad\tau\gg 1\;.$ (7)

If we assume that the excitations are uncorrelated following the reasoning in
the Appendix we arrive at the distribution in Eq.(6), which is calculated by
the SDRG approach in the given limits.

In this paper, we aim to study this problem in more detail consider the
following points. i) To calculate the distribution function through numerical
iteration of the SDRG approach and to check if the result in Eq.(6) is valid
in the entire Griffiths phase. ii) To confirm the results of the SDRG
calculation with the exact gaps obtained through free-fermion techniques. iii)
To check the form of the finite-size corrections of the two methods and to
compare with the analytical result for the extremes of uncorrelated random
variables.

##  III Methods to calculate the gap in the RTIM

The energy scale of the model is given by the lowest excitation energy of the
Hamiltonian in Eq.(1). This is calculated for finite chains by the
asymptotically exact SDRG method through iteration and for shorter chains by
free-fermion techniques.

###  III.1 The SDRG method

In the SDRG procedureim  we perform consecutive decimation steps, each time
considering the local excitations, say at position $i$. These excitations
correspond to either couplings or sites, having the value of the associated
gaps: $2J_{i}$ and $2h_{i}$, respectively. These gaps are sorted in descending
order and the largest one, denoted by $\Omega$, which sets the energy-scale in
the problem, is eliminated. Then, between remaining degrees of freedom, new
terms in the Hamiltonian are generated through perturbation calculation. This
procedure is successively iterated, during which $\Omega$ monotonously
decreases. At the fixed point, with $\Omega^{*}=0$, one makes an analysis of
the distribution of the different parameters and calculates the scaling
properties. In the following, we describe the elementary decimation steps.

####  III.1.1 Strong-coupling decimation

In this case, the largest local term in the Hamiltonian is a coupling, say
$\Omega=J_{i}$, connecting sites $i$ and $i+1$ and the two-site Hamiltonian is
given by

$\hat{H}_{cp}=-J_{i}\sigma^{z}_{i}\sigma^{z}_{i+1}-h_{i}\sigma^{x}_{i}-h_{i+1}\sigma^{x}_{i+1}\;.$
(8)

The spectrum of $\hat{H}_{cp}$ contains four levels, the lower two being
separated from higher two by a gap of $2J_{i}$. We omit the higher two levels,
corresponding to merging the two strongly coupled sites into a spin cluster in
the presence of a (renormalized) transverse field $\tilde{h}$, the value of
which is given by second-order perturbation calculation

$\tilde{h}=\frac{h_{i}h_{i+1}}{J_{i}}\;.$ (9)

####  III.1.2 Strong-transverse-field decimation

In this case, the largest local term is a transverse field, say $h_{i}$, and
due to its large value this site does not contribute to the longitudinal
magnetisation and therefore it is eliminated. The renormalised coupling
between the remaining sites is given by

$\tilde{J}=\frac{J_{i-1}J_{i}}{h_{i}}\;,$ (10)

which is calculated by second order perturbation method.

To calculate the smallest gap of a given sample, $\varepsilon$, we perform
$(L-1)$ decimation steps up to the last spin cluster having an effective
transverse field $\tilde{h}=\varepsilon/2$.

###  III.2 Free-fermion technique

In this method, $\hat{H}$ is expressed in terms of spinless free fermions JW ;
lieb61 ; pfeuty79 . In the first step the spin operators $\sigma_{i}^{x,y,z}$
are mapped to fermion creation (annihilation) operators $c_{i}^{\dagger}$
($c_{i}$) by using the Jordan-Wigner transformation JW :
$c^{\dagger}_{i}=a_{i}^{+}\exp\left[\pi\imath\sum_{j=1}^{i-1}a_{j}^{+}a_{j}^{-}\right]$
and
$c_{i}=\exp\left[\pi\imath\sum_{j=1}^{i-1}a_{j}^{+}a_{j}^{-}\right]a_{i}^{-}$,
where $a_{j}^{\pm}=(\sigma_{j}^{x}\pm\imath\sigma_{j}^{y})/2$, and the Ising
Hamiltonian in Eq.(1) is written in a quadratic form

$\displaystyle{\hat{H}}=-\sum_{i=1}^{L}h_{i}\left(2c^{\dagger}_{i}c_{i}-1\right)-\sum_{i=1}^{L-1}J_{i}(c^{\dagger}_{i}-c_{i})(c^{\dagger}_{i+1}+c_{i+1})\;.$
(11)

In the second step, the Hamiltonian in Eq.(11) is diagonalized through a
canonical transformation lieb61 , in terms of the new fermion creation
(annihilation) operators $\eta_{k}^{\dagger}$ ($\eta_{k}$)

${\hat{H}}=\sum_{k=1}^{L}\epsilon_{k}\left(\eta_{k}^{{\dagger}}\eta_{k}-\frac{1}{2}\right)\;.$
(12)

The energies of free fermionic modes, $\epsilon_{k}$, are given by the
eigenvalues of a $2L\times 2L$ tridiagonal matrix

$T=\begin{pmatrix}0&h_{1}&&&&&\cr h_{1}&0&J_{1}&&&&\cr
0&J_{1}&0&h_{2}&&&\cr&&h_{2}&0&\ddots&&\cr&&&\ddots&\ddots&J_{L-1}&\cr&&&&J_{L-1}&0&h_{L}\cr&&&&&h_{L}&0\cr\end{pmatrix}$
(13)

and we consider only the $\epsilon_{k}\geq 0$ part of the
spectrumigloiturban96 . The smallest gap of ${\hat{H}}$ in Eq.(1) is given by
$\varepsilon={\rm min}|\epsilon_{k}|$.

##  IV Numerical results

In the numerical calculation, the parameters of the Hamiltonian, $J_{j}$ and
$h_{j}$ are taken from box-like distributions

$\displaystyle\begin{split}\pi_{1}(J)&=\begin{cases}1&\hskip
18.49411pt\text{for }0<J\leq J_{0}\,,\\\ 0&\hskip
18.49411pt\text{otherwise.}\end{cases}\\\
\pi_{2}(h)&=\begin{cases}1/h_{0}&\text{for }0<h\leq h_{0}\,,\\\
0&\text{otherwise,}\end{cases}\end{split}$ (14)

and set the energy-scale with $J_{0}=1$. In this case the critical point of
the RTIM is located at $h_{0}=1$ and in the disordered phase, $h_{0}>1$, the
dynamical exponent satisfies the equation:

$(1-z^{-2})h_{0}^{1/z}=1\;,$ (15)

see in Eq.(5). In the vicinity of the critical point $z$ diverges as
$z\approx\tfrac{1}{\ln h_{0}}$, $h_{0}\to 1^{+}$, whereas for large $h_{0}$ it
approaches $1$ as $z\approx 1+\tfrac{1}{2h_{0}}$, $h_{0}\to\infty$.

In the actual calculations we considered three points of the paramagnetic
phase, $h_{0}=2$ ($z=1.747655$), $h_{0}=3$ ($z=1.33542$) and $h_{0}=4$
($z=1.2112289$). In the free-fermion calculation of the gaps the lengths of
the systems were $L=16,24,32,48$ and $64$, whereas with the SDRG method larger
systems up to $L=512$ are treated. In all cases $10^{10}$ independent samples
are investigated.

###  IV.1 SDRG gaps and the Fréchet distribution

We start by analysing the data collected by the SDRG method and presenting the
distributions of the log-gaps in Fig.1 for different sizes. To test their
relation with the Fréchet distribution, we calculated in each case the best
fit of the analytical function in Eq.(32), having the dynamical exponent,
$z_{L}\equiv\gamma_{L}$ and the position of the maximum, $x_{0}=x_{0}(L)$ as
fit parameters. These are shown in Fig.1, too.

Figure 1: (Color online) Distribution of the energy gaps calculated by the
SDRG algorithm in the Griffiths phase with $h_{0}=4$ (upper panel), $h_{0}=3$
(middle panel) $h_{0}=2$ (lower panel) for finite chains of lengths
$L=32,64,128,256$ and $512$, from left to right. The best fit of the Fréchet
distribution is indicated by full lines.

As seen in Fig.1, the Fréchet distribution describes the numerical data well
and the difference between the numerical points and the fitting curve is
decreasing for increasing values of $L$. At a given value of $h_{0}$ the
fitted value of the dynamical exponent, $z_{L}$, has a small variation with
the size, and approaches the exact asymptotic value, as given in Eq. (15). At
the value $h_{0}=4$ this is illustrated in Table 1, where the second column
shows the fitted, finite size values of $z_{L}$, while the third column shows
the difference from the exact value. As a first step, we used a power-law form
$z-z_{L}\sim L^{-a}$ to fit the finite-size corrections, and estimates for the
exponent $a$ are calculated through two-point fit. These are listed in the
fourth column of the table and have a slow convergence, which generally
indicates a strong correction to scaling term. Having this possibility in mind
we have used another functional form:

$z-z_{L}\sim\frac{\ln^{\omega}L}{L}\;.$ (16)

In this case the effective exponents for $\omega$ are also calculated through
two-point fit and are presented in the fifth column of the table. This type of
fitting turned out to be more stable, the effective exponents seem to converge
to $\omega\approx 1.33$.

$L$ | $z_{L}$ | $z-z_{L}$ | $a$ | $\omega$
---|---|---|---|---
32 | 1.102788 | 0.108441 |  |
64 | 1.142667 | 0.068562 | 0.6614 | 1.4081
128 | 1.169523 | 0.041706 | 0.7172 | 1.3724
256 | 1.186572 | 0.024657 | 0.7583 | 1.3405
512 | 1.196923 | 0.014306 | 0.7854 | 1.3388
Table 1: Finite-size estimates of the dynamical exponent, $z_{L}$, at
$h_{0}=4$ calculated by the SDRG algorithm from the best fit of the Fréchet
form in Eq.(32) and its difference from the asymptotically exact value
$z=1.2112288988$. Exponent of a power-law fit, $a$ and that of the logarithmic
correction, $\omega$. (See text.)

We have repeated the same analysis at the other two points of the Griffiths
phase. In both cases the log-correction form is found to provide the better
fit, with the correction exponents $\omega\approx 1.83$ at $h_{0}=2$ and
$\omega\approx 1.50$ at $h_{0}=3$.

Most importantly, at a given value of $h_{0}$ in Fig.1 the distributions are
shifted with increasing size, and according to Eq.(4) the position of the
maximum is expected to follow the rule: $x_{0}(L)\approx const.+z_{L}\ln L$.
Comparing the position of the maximum of the distribution at two sizes one can
obtain estimates for the dynamical exponent as:
$z(L,2L)=(x_{0}(2L)-x_{0}(L))/\ln 2$. We have checked that generally
$z_{L}<z(L,2L)<z_{2L}$ and these estimates approach the asymptotic exact value
with the same type of corrections as noticed for the case of $z_{L}$ in the
previous paragraph.

We note that analysis of the data obtained by free-fermion calculation of the
gap gives similar results, but due to smaller values of $L$ the asymptotic
region of the effective exponents, $z_{L}$ is more remote.

###  IV.2 Finite-size corrections to the gap distributions

|  |
---|---|---
|  |
|  |
Figure 2: (Color online) Finite-size corrections to the Fréchet distribution
for the inverse gaps, $\tau=1/\varepsilon$ in the RTIM at different sizes $L$
calculated through SDRG iterations (left panels), and by the free-fermion
method (middle panels) with the first standardisation in Eq.(23) and compared
with the analytical results with $\gamma=z(h_{0})$ and $\gamma^{\prime}=-1$
(dashed lines). In the right panels the maximum of uncorrelated Kesten random
variables are shown, having the same scaling exponents, see text. The scaling
exponents $z(h_{0})$ correspond to $h_{0}=4$ (first row), $h_{0}=3$ (second
row) and $h_{0}=2$ (third row), see in Eq.(15).

We considered the gaps calculated by SDRG iteration for $L=24,32,48,64,96$ and
$128$ and for comparison we calculated those by free-fermion techniques as
well, for shorter chains with $L=16,24,32,48$ and $64$, except for $h_{0}=4$,
where the corrections are the smallest and we went up to $L=48$. We analysed
finite-size scaling of distributions in two different ways. First, we used the
_first standardisation convention_ in Eq.(23) and calculated the difference
with the Fréchet extreme distribution, with an effective $z_{L}$, calculated
from the relations below Eq.(32). This difference is then rescaled by a factor
$\Delta z=z-z_{L}$, and the results for $h_{0}=2,3$ and $4$ are drawn in
Fig.2. In these figures the analytical results calculated for iid random
numbers in Eq.(30) with $\gamma=z$ and $\gamma^{\prime}=-1$ are also
presented. Here, the correction to scaling exponent, $\gamma^{\prime}=-1$
corresponds to the expected scaling form in Eq.(16), which is obtained through
numerical analysis of the data in Table 1. According to Fig.2, we can draw the
following conclusions. i) The scaled finite-size difference of the
distribution function seems to approach a limiting curve for large $L$, which
depends on $z(h_{0})$. ii) At a given $z(h_{0})$ the limiting curves are
similar (if not identical) for the SDRG and the free-fermion data. iii) The
convergence to this limit curve is slow, much slower than that of the iid
random numbers, see later in Sec.IV.3 and the figures in the third column of
Fig.2. This slow convergence is probably related to the logarithmic correction
to the scaling of the dynamical exponent. iv) Finally, the expected limit
curve of the numerical distribution differs from that of the analytical result
calculated for iid random numbers (having identical parameter, $z$). Even
though the overall shapes of the curves are similar, there are noticeable
differences. In particular, the low-energy part of the numerical curves are
stronger represented in the numerical curves, which can be interpreted as the
reduction of the value of the gap due to small, but relevant correlations
between the rare regions.

We have repeated the analysis with the distribution of the log-gaps and using
the second standardisation condition in Eq.(24). The results about the finite-
size corrections of the distributions are shown in Fig.3 together with the
analytical curves for iid random numbers. In Fig.3 the distributions of the
log-gaps have a faster finite-size convergence, than those of the gaps in
Fig.2. The curves in Fig.3 are almost indistinguishable for larger sizes. In
addition, we notice a difference between the shape of the numerical curves and
the analytical results, the former being somewhat shifted to the right around
zero. This indicates a reduction of the value of the gap due to small, but
relevant correlations between the rare regions.

|  |
---|---|---
|  |
|  |
Figure 3: (Color online) The same as in Fig.2 for the log-variable and using
the second standardisation condition in Eq.(24).

###  IV.3 Numerical test for uncorrelated Kesten variables

In order to test the possible convergence of the finite-size corrections for
uncorrelated variables we have repeated the analysis in the previous
subsection for a parent distribution generated by Kesten random numbers. Here,
we remark that Kesten-type random variables are defined askesten73 :

$u_{m}=1+\sum_{i=1}^{m}\prod_{j=1}^{i}s_{j}+\dots\;,$ (17)

where the $s_{j}$-s are iid random numbers. It is known, that in the limit of
$m\to\infty$ there is a limit distribution, $\rho_{\infty}(u)$, provided $[\ln
s]_{\rm av}<0$. This limit distribution has a power-law tail

$\rho_{\infty}(u)\sim u^{-(1+\alpha)};\quad u\gg 1\;,$ (18)

where the exponent is the positive root $\alpha>0$ of the equation

$\left[s^{\alpha}\right]_{\rm av}=1\;.$ (19)

Note that the relation for $\alpha$ is analogous to the equation for the
dynamical exponent, $z$ of the RTIM, see in Eq.(5), and so we have the
relation $\alpha=1/z$. To have a direct relation with the RTIM calculations we
set $s_{j}=J_{j}/h_{j}$ where $J_{j}$ and $h_{j}$ are taken from the
distributions in Eq.(14). In the numerical calculation, we have $m=64$, and
considered the maximum of a set of $L=16,24,32,48$ and $64$ Kesten numbers,
denoted by $\tau$, so that we have $m\geq L$ in each case. We have checked
that at such value of $m$ the truncation of the series in Eq.(17) has
negligible error.

We have analysed the distribution of the maximum values at the three points
$h_{0}=2$, $h_{0}=3$ and $h_{0}=4$ as done previously for the gaps of the
RTIM. Using the first standardisation condition in Eq.(23) the finite-size
corrections are shown in the third column of Fig.2. Here, the analytical
results are obtained with the exponents: $\gamma=z(h_{0})$ and
$\gamma^{\prime}=-1$, the latter follows from the analytical results in
Ref.[clnp, ]. It is seen in Fig.2 that there is an overall good agreement
between the numerical and analytical results. There is some size-dependence of
the corrections, which changes sign between $h_{0}=2$ and $h_{0}=3$. In
comparison to the distribution of the gaps in the RTIM, the agreement with the
analytical results is much better and the finite-size corrections are smaller.

We have repeated the analyses of the data by considering log-variables, as
given in Eqs.(31) and (32) and using the second standardisation condition in
Eq.(24). Results of the numerical analysis in this case are given in the third
column of Fig.3. As seen in this figure, there are some finite-size
dependences of the scaled numerical data, but the expected asymptotic curves
agree well with the analytical results. We note that practically no finite-
size dependence of the correction term is observed for a pure power parent
distribution in Fig.4 in the Appendix.

##  V Discussion

In this paper we have considered a paradigmatic model of random quantum
magnets, the random transverse Ising model in 1D and studied the distribution
of low-energy excitations in the paramagnetic Griffiths phase, with extensive
numerical methods. We have considered a large set of random samples
($10^{10}$) and the calculation is performed by the approximate, but
asymptotically correct SDRG method (up to $L=512$) and for comparison we also
used the free-fermion method for shorter chains (up to a size $L=64$).
Analysing the distribution of the gaps we have demonstrated with high
precision that – in agreement with previous expectations – the limit
distribution in the thermodynamic limit is in the Fréchet form. The Fréchet
distribution depends on the value of the dynamical exponent, $z$ and we have
shown that a powerful method of calculation is through fitting a Fréchet curve
to the numerical gap distributions. In this way, effective, finite-size
estimates are obtained for the dynamical exponent $z_{L}$, which then are
extrapolated to $L\to\infty$. According to the numerical data this convergence
is in the form: $z-z_{L}\sim\ln^{\omega}L/L$, where the exponent of the
logarithm, $\omega$, depends on the distance from the critical point.

More interestingly, we have systematically studied the finite-size corrections
to the limit law and showed that the difference between the numerical
distribution and the asymptotic Fréchet form scales with $z-z_{L}$, and this
scaled difference, $P_{1}(\varepsilon;z)$ is a unique function of the value of
the gap, $\varepsilon$. We have performed this type of analysis for the gap,
using the first standardisation condition in Eq.(23), as well as for the log-
gap, when the second standardisation condition in Eq.(24) was used. In both
cases the asymptotic form of the scaled finite-size correction function are
found similar (if not identical) for the data with the free-fermion
calculation and that from the SDRG iteration. We can thus conclude that the
SDRG method provides not only the correct asymptotic form of the gap
distribution function, but the correct finite-size correction function as
well.

The measured finite-size correction functions are also compared with the
analytical results of iid random numbers, having the same decay exponent
($\gamma=z$ for the gaps and $\gamma=0$ for the log-gaps) and correction to
scaling exponent, $\gamma^{\prime}=-1$. The two curves are found to have
similar shape, but there are also differences in the asymptotic forms. We have
checked that the observed differences are larger than the statistical error of
the calculation. For this purpose, we have analysed with identical methods the
same set of uncorrelated Kesten random variables having the same
characteristic exponent, $z$. For these iid random numbers the asymptotic
finite-size corrections are found to be well described by the analytical
results.

The observed difference in the finite-size corrections between the numerical
curves and the analytical iid results indicates that the weak correlations
between low-energy excitations in the RTIM are not completely irrelevant. This
is probably related to the fact that the linear extension of rare regions,
$n$, with extreme fluctuations of the strong couplings, scales as $n\sim\ln
L$, see the reasoning above Eq.(7). Since $n$ is not limited, the finite-size
dynamical exponents contain a logarithmic multiplicative factor as given in
Eq.(16).

Our investigations have concluded on the RTIM in 1D, but the observed results
are possibly valid for other random quantum systems as well with localised
excitations. For example, we can mention the RTIM in higher dimensions2d ;
2dki ; ddRG , the random quantum Pottssenthil  and Ashkin-Teller modelscarlon-
clock-at  and generally random quantum magnets with short range interactions
and with discrete symmetry. Similar conclusions apply to some non-equilibrium
processes, such as to the random asymmetric exclusion processigloipartial  and
the random contact processigloicontact  as well.

###### Acknowledgements.

This work was supported by the National Research Fund under Grants No.
K128989, No. K115959 and No. KKP-126749. F.I. is indebted to Katalin Ozogány
for sending her PhD thesis Ref.[ozogany, ] and for valuable correspondence.

## Appendix - Extreme statistics of uncorrelated variables

Let us have a process of iid random numbers $u_{1},u_{2},\dots,u_{N}$ and
consider the maximum value $u={\rm max}(u_{1},u_{2},\dots,u_{N})$. Each random
variable is distributed by the same, so called _parent_ distribution function
$\rho(u)$ and the cumulative (or integrated) parent distribution is given by
$\mu(u)=\int_{-\infty}^{u}\rho(t){\rm d}t$. The cumulative distribution of the
maximum value, $u$ is easy to compute

$M_{N}^{\textrm{max}}(u)={\rm Prob}[u_{1}\leq u,u_{2}\leq u,\dots,u_{N}\leq
u]=\mu^{N}(u)\;.$ (20)

According to extreme value statistics (EVS) $M_{N}^{\textrm{max}}(u)$ has a
limit distribution for large $N$ and large $x$ in terms of the scaling
combination $v=(u-b_{N})/a_{N}$ as

$M_{N}(v)=M_{N}^{\textrm{max}}(a_{N}v+b_{N})\to M(v)\;,$ (21)

and similarly for the extreme density $P_{N}(v)={\rm d}M_{N}/{\rm d}v$

$P_{N}(v)=a_{N}P_{N}^{\textrm{max}}(a_{N}v+b_{N})\to P(v)\;.$ (22)

Here, the $a_{N}$ and $b_{N}$ are free up to an additive constant, the value
of which is fixed by different standardisation conditions. For the analytical
calculation the condition:

$M(0)=P(0)=1/e\;,$ (23)

is convenient to use, whereas for analysing numerical data it is often better
to require

$\int_{-\infty}^{\infty}vP(v)=0,\quad\int_{-\infty}^{\infty}v^{2}P(v)=1\;,$
(24)

provided the second moment of the distribution exists. In the paper, we refer
to Eq.(23) and Eq.(24) as _first_ and _second standardisation condition_ ,
respectively.

According to EVS theory, the limit distributions can be of three forms,
depending on the large $u$ tail of the parent distribution. In a
renormalization group treatmentprl ; pre  with the first standardisation in
Eq.(23), the limit distributions are the fixed-point solutions and are given
in the form

$\displaystyle M(v;\gamma)$ $\displaystyle=$
$\displaystyle\exp\left[-(1+\gamma v)^{-1/\gamma}\right]\;,$ (25)
$\displaystyle P(v;\gamma)$ $\displaystyle=$ $\displaystyle(1+\gamma
v)^{-1/\gamma-1}\exp\left[-(1+\gamma v)^{-1/\gamma}\right]\;,$ (26)

for $1+\gamma v\geq 0$ and the different universality classes depend on the
value of $\gamma$. For $\gamma>0$ it corresponds to a parent distribution

$\mu(u)\approx 1-Au^{-1/\gamma},\quad u\gg 1\;,$ (27)

which represents the Fréchet universality class. For $\gamma=0$ the asymptotic
approach of the cumulative parent distribution is faster than a power and
represents the Gumbel universality class. Finally, for $\gamma<0$ the
asymptotic value is approached at a finite upper border with a power
$-1/\gamma$ and is given as the Weibull distribution.

Finite-size corrections in the RG treatment are given in the form

$\displaystyle M_{N}(v)$ $\displaystyle\approx$ $\displaystyle
M(v;\gamma)+\varepsilon_{N}M_{1}(v;\gamma,\gamma^{\prime})$ (28)
$\displaystyle=$ $\displaystyle
M(v;\gamma)+\varepsilon_{N}P(v;\gamma)\psi(v;\gamma,\gamma^{\prime})\;,$

where

$\varepsilon_{N}=\gamma-\gamma_{N}\propto N^{\gamma^{\prime}}\;,$ (29)

and

$\psi(v;\gamma,\gamma^{\prime})=\frac{(1+\gamma v)+\gamma^{\prime}v-(1+\gamma
v)^{\gamma^{\prime}/\gamma+1}}{\gamma^{\prime}(\gamma^{\prime}+\gamma)}\;.$
(30)

Consequently, the correction term to the cumulative distribution
$M_{1}(v;\gamma,\gamma^{\prime})$ and that of the density,
$P_{1}(v;\gamma,\gamma^{\prime})={\rm d}M_{1}(v;\gamma,\gamma^{\prime})/{\rm
d}v$ depend on the decay exponent of the parent distribution $\gamma$, as well
as the correction exponent $\gamma^{\prime}$.

In practical analysis of the data with a parent distribution in Eq.(27), it is
convenient to use logarithmic variables, $x=\ln u$, so that
$e^{x-x_{0}}=1+\gamma v$. Since the parent distribution of $x$ is exponential,
the extreme distribution is given by the Gumbel form

$\tilde{M}(x,\gamma)=\exp\left(-\exp\left(\frac{-x+x_{0}}{\gamma}\right)\right)\;,$
(31)

and similarly for the probability density

$\tilde{P}(x,\gamma)=\frac{1}{\gamma}\exp\left(\frac{-x+x_{0}}{\gamma}\right)\exp\left(-\exp\left(\frac{-x+x_{0}}{\gamma}\right)\right)\;,$
(32)

so that $\tilde{M}(x_{0},\gamma)=1/e$ and
$\tilde{P}(x_{0},\gamma)=1/(e\gamma)$. In this case, the finite-size
correction term is given by $\tilde{P}_{1}(x;0,\gamma^{\prime})$. In the
numerical examples, we have $\gamma=0$ and $\gamma^{\prime}=-1$, and the
correction to the density in the second standardisation condition is given
byozogany

$\displaystyle\tilde{P}_{1}(x;0,-1)$
$\displaystyle=ae^{-\tilde{x}-e^{-\tilde{x}}}\left[(e^{-\tilde{x}}-1)(1-x/a-e^{-\tilde{x}})\right.$
(33) $\displaystyle\left.-a^{-2}+e^{-\tilde{x}}\right]\;,$

with $\tilde{x}=ax+b$, $a=\pi/\sqrt{6}$ and $b=\gamma_{E}=0.5772156649$ being
the Euler-Mascheroni constant. As an illustration, we consider a pure power
parent distribution $\mu(u)=1-u^{-1/\gamma},~{}u\geq 1$ and in Fig.4 show the
numerically calculated finite-size corrections to the distributions of the
log-extremes, denoted by $\ln(\tau)$, for $L=16,24,32,48$ and $64$ with
$\gamma=z(h_{0}=3)$ in the second standardisation condition. These are to be
compared with the analytical result in Eq.(33), since in this case
$\gamma^{\prime}=-1$. The agreement is almost perfect.

Figure 4: (Color online) Finite-size corrections to the distributions of the
log-extremes for a pure power parent distribution with the second
standardisation condition, see text.

## References

  * (1) R.B. Griffiths, Phys. Rev. Lett. 23, 17 (1969).
  * (2) B. McCoy, Phys. Rev. Lett. 23, 383 (1969).
  * (3) T. Vojta, J. Phys. A: Math. Gen. 39, R143 (2006).
  * (4) For reviews, see: F. Iglói and C. Monthus, Physics Reports 412, 277, (2005); Eur. Phys. J. B 91, 290 (2018).
  * (5) S.-K. Ma, C. Dasgupta, and C. K. Hu, Phys. Rev. Lett. 43, 1434 (1979); C. Dasgupta and S.-K. Ma, Phys. Rev. B 22, 1305 (1980).
  * (6) D.S. Fisher, Phys. Rev. Lett. 69, 534 (1992); Phys. Rev. B 51, 6411 (1995).
  * (7) D.S. Fisher, Physica A 263, 222 (1999).
  * (8) M.J. Thill and D.A. Huse, Physica A 15, 321 (1995).
  * (9) R. Juhász, Y.-C. Lin and F. Iglói, Phys. Rev. B 73, 224206 (2006).
  * (10) R. A. Fisher and L. H. C. Tippett, Proc. Cambridge Phil. Soc. 24, 180 (1928).
  * (11) E. J. Gumbel, Statistics of Extremes (Dover, New York, 1958).
  * (12) B. V. Gnedenko, Ann. Math. 44, 423 (1943).
  * (13) W. Weibull, J. Appl. Mech. 18, 293 (1951).
  * (14) J. Galambos, The Asymptotic Theory of Extreme Order Statistics (Wiley, New York, 1978).
  * (15) M. R. Leadbetter, G. Lindgren, H. Rootzén, Extremes and related properties of random sequences and processes (Springer, New York, 2012),
  * (16) L. de Haan and S. Resnick, Ann. Probab. 24, 97 (1996); L. de Haan and U. Stadtmüller, J. Austral. Math. Soc. 61, 381 (1996).
  * (17) G. Györgyi, N. R. Moloney, K. Ozogány, and Z. Rácz, Phys. Rev. Lett. 100, 210601 (2008).
  * (18) G. Györgyi, N. R. Moloney, K. Ozogány, Z. Rácz and M. Droz, Phys. Rev. E. 81, 041135 (2010).
  * (19) E. Bertin and G. Györgyi, J. Stat. Mech., P08022 (2010).
  * (20) L. Zarfaty, E. Barkai, and D. A. Kessler, arXiv:2006.13677.
  * (21) S. N. Majumdar, A. Pal and G. Schehr, Physics Reports 840, 1 (2020).
  * (22) P. Jordan and E. Wigner, Z. Phys. 47 631 (1928).
  * (23) E. Lieb, T. Schultz and D. Mattis, Ann. Phys. (N.Y.) 16, 407 (1961); S. Katsura, Phys. Rev. 127, 1508 (1962); P. Pfeuty, Ann. Phys. (Paris) 57, 79 (1970).
  * (24) P. Pfeuty, Phys. Lett. 72A 245 (1979).
  * (25) H. Kesten, Acta Math. 131, 208 (1973).
  * (26) F. Iglói, R. Juhász and P. Lajkó, Phys. Rev. Lett. 86, 1343 (2001).
  * (27) F. Iglói, Phys. Rev. B 65, 064416 (2002).
  * (28) A. P. Young and H. Rieger, Phys. Rev. B 53, 8486 (1996).
  * (29) F. Iglói and H. Rieger, Phys. Rev. B 57, 11404 (1998).
  * (30) F. Iglói and H. Rieger, Phys. Rev. E 58, 4238 (1998).
  * (31) D. S. Fisher and A. P. Young, Phys. Rev. B 58, 9131 (1998).
  * (32) F. Iglói and L. Turban, Phys. Rev. Lett. 77, 1206 (1996).
  * (33) C. de Calan, J.M. Luck, Th. M. Nieuwenhuizen and D. Petritis, J. Phys. A 18, 501 (1985).
  * (34) O. Motrunich, S.-C. Mau, D.A. Huse and D.S. Fisher, Phys. Rev. B 61, 1160 (2000); Y.-C. Lin, N. Kawashima, F. Iglói and H. Rieger, Progress in Theor. Phys. 138, (Suppl.) 479 (2000); D. Karevski, Y-C. Lin, H. Rieger, N. Kawashima and F. Iglói, Eur. Phys. J. B 20 267 (2001); Y-C.Lin, F. Iglói and H. Rieger, Phys. Rev. Lett. 99, 147202 (2007); R. Yu, H. Saleur and S. Haas, Phys. Rev. B 77, 140402 (2008);
  * (35) I. A. Kovács and F. Iglói, Phys. Rev. B 80, 214416 (2009), Phys. Rev. B 82, 054437 (2010).
  * (36) I. A. Kovács and F. Iglói, Phys. Rev. B 83, 174207 (2011), J. Phys. Condens. Matter 23, 404204 (2011).
  * (37) T. Senthil and S. N. Majumdar Phys. Rev. Lett. 76, 3001 (1996).
  * (38) E. Carlon, P. Lajkó, and F. Iglói Phys. Rev. Lett. 87, 277201 (2001).
  * (39) R. Juhász, L. Santen and F. Iglói, Phys. Rev. Lett. 94, 010601 (2005).
  * (40) J. Hooyberghs, F. Iglói, and C. Vanderzande Phys. Rev. Lett. 90, 100601 (2003); Phys. Rev. E 69, 066140 (2004).
  * (41) K. Ozogány, PhD thesis - Loránd Eötvös University (2012).


# A machine learning approach to drawing phase diagrams of topological lasing
modes

Stephan Wong  School of Physics and Astronomy, Cardiff University, Cardiff
CF24 3AA, UK  Jan Olthaus  Institut für Festkörpertheorie, Universität
Münster, 48149 Münster, Germany  Thomas K. Bracht  Institut für
Festkörpertheorie, Universität Münster, 48149 Münster, Germany  Doris E.
Reiter  Institut für Festkörpertheorie, Universität Münster, 48149 Münster,
Germany  Condensed Matter Theory, Department of Physics, TU Dortmund, 44221
Dortmund, Germany  Sang Soon Oh<EMAIL_ADDRESS>School of Physics and
Astronomy, Cardiff University, Cardiff CF24 3AA, UK

###### Abstract

Identifying phases and analyzing the stability of dynamic states are
ubiquitous and important problems which appear in various physical systems.
Nonetheless, drawing a phase diagram in high-dimensional and large parameter
spaces has remained challenging. Here, we propose a data-driven method to
derive the phase diagram of lasing modes in topological insulator lasers. The
classification is based on the temporal behaviour of the topological modes
obtained via numerical integration of the rate equation. A semi-supervised
learning method is used and an adaptive library is constructed in order to
distinguish the different topological modes present in the generated parameter
space. The proposed method successfully distinguishes the different
topological phases in the Su-Schrieffer-Heeger (SSH) lattice with saturable
gain. This demonstrates the possibility of classifying the topological phases
without needing for expert knowledge of the system and may give valuable
insight into the fundamental physics of topological insulator lasers via
reverse engineering.

##  I Introduction

Over the last few years, significant research efforts have been made on
photonic topological insulator (PTI) lasers. While the efforts have been
concentrated on the spatial stability of the topologically protected edge
modes, namely on the existence of such topological edge modes in non-Hermitian
PTIs [1, 2, 3, 4, 5, 6, 7, 8, 9], the temporal stability has not been the
focus of interest so far [10, 11, 12, 12, 13, 14, 15]. Due to the non-linear
nature of PTI lasers, the temporal stability is an important characteristic to
take into account for experimental demonstrations and real-life applications.
Although the spatial stability of the topological modes, _i.e._ , its
robustness, may be guaranteed in active non-Hermitian PTIs, unstable behaviour
may be present in the time domain. In this regard, the temporal dynamics of
the topologically protected modes have been studied [10, 12, 13], mainly using
linear stability analysis. It is, however, a challenging task to apply the
same approach to more complex structures because of the high-dimensional phase
space and parameter space as well as the lack of analytical solutions [11].

Machine learning (ML) can be advantageous for the theoretical study of the
stability of PTIs, which requires repetitive numerical simulations for several
varying parameters. ML is a data-based method which can be implemented with
different strategies, and the most appropriate ML strategy depends on the
dataset under study. For instance, a supervised learning strategy relies on
labelled data, a dataset of input-output pairs. This has been utilized in
topological photonics to draw topological phase diagrams [16], calculate
topological invariants [5], or explore topological band structures [17]. On
the other hand, an unsupervised learning strategy consists of extracting
information from the dataset for which we do not have labels. This is used for
dimensional reductions by keeping only the main features of the high-
dimensional structure of the dataset or for clustering problems from which the
data is divided into different types [18]. For instance, this has been
successful in obtaining the phase transition in the Ising model [19], and
clustering Hamiltonians that belong to the same symmetry classes [20].

In the unsupervised learning strategy, modal decomposition is a common and
successful method which reduces the analysis of very high-dimensional data to
a set of relatively few modes. Among the modal decomposition methods,
principal component analysis (PCA) is a method which derives the eigenmodes or
the main features based on their variance in the data [21]. These eigenmodes
can then be utilized as a basis to represent the dataset [22]. This reduced-
order model method has been extended to identify distinct non-linear regimes
[23, 24, 25, 26, 27, 28, 29] by constructing a library composed of
representatives of these regimes: this is known as representation
classification. Nevertheless, the preliminary identification of the regimes
composing the library and the construction of the library is a manual process
and requires expert knowledge of the complex system.

In this paper, we propose a representation classification method to study the
spatio-temporal dynamics of non-linear topological systems. The results will
be based on the phase diagram of the Su-Schrieffer-Heeger (SSH) lattice [30]
with a domain wall and with saturable gain [12, 13]. To remove the necessity
of the required expertise on the complex system, we present an algorithm which
constructs an appropriate library of the different phases automatically. For
this goal, we propose two approaches: a top-down approach in which the library
has numerous phases that are merged into the equivalent phases, and a bottom-
up approach in which the library is completed on the fly to get the most
accurate classification. Via reverse engineering, our proposed method can be
used as a tool to find novel topological lasing modes in more complicated
settings. For given rate equations of a lasing system, one would only need to
integrate the differential equations in the desired parameter space region and
then apply the adaptive representation classification to obtain the phase
diagram.

##  II Results

###  II.1 Phase diagram of the SSH Model

Figure 1:  Phase diagram of the domain-wall-type SSH lattice with saturable
gain. (a) Schematic of the domain-wall-type SSH lattice considered. The
vertical red line is a guide to the eye for the domain wall between the two
SSH lattices, SSH 1 and SSH 2. (b) Phase diagram of the domain-wall-type SSH
lattice with saturable gain and linear loss on the A and B sublattices. The
white and grey areas correspond to the non-oscillating and oscillating
topological phases, respectively. (c) Representative total intensity $I_{A}$
(and $I_{B}$) of the A (and B) sublattice in blue (and orange) for the non-
oscillating and oscillating topological lasing mode on the top and bottom
panel, respectively. Spatio-temporal dynamics of the (d) non-oscillating and
(e) oscillating topological lasing modes. The non-oscillating and oscillating
topological modes displayed are chosen at
($\gamma_{AB},g_{A}-\gamma_{AB})=(0.48,0.06)$ and $(0.16,0.44)$, respectively.

As a toy model, we will consider the domain-wall-type SSH lattice with
saturable gain [Fig. 1(a)]. The system has a domain wall, at the A site $n=0$,
which separates two SSH lattices, namely lattices composed of two sites per
unit cell, A and B, and characterized by intra- and inter-unit cell couplings,
$t_{\text{intra}}$ and $t_{\text{inter}}$, respectively.
$t_{\text{intra}}=t_{1}$ and $t_{\text{inter}}=t_{2}$
($t_{\text{intra}}=t_{2}$ and $t_{\text{inter}}=t_{1}$) for the lattice on the
left (right) side of the domain wall, _i.e._ , the sites with $n<0$ ($n\geq
0$). The dynamics of the system
$\psi(t)=\left(\psi_{-N}(t),\ldots,\psi_{-1}(t),\psi_{0}(t),\psi_{1}(t),\ldots,\psi_{N}(t)\right)\equiv
x(t)$, with $N_{s}=2N+1$ sites, reads, for $n=-N,\ldots,N$:

$\displaystyle\begin{split}i\frac{d\psi_{n}}{dt}&=i\left(\frac{g_{n}}{1+|\psi_{n}|^{2}}-\gamma_{n}\right)\psi_{n}\\\\[4.30554pt]
&\qquad+t_{\text{intra,n+1}}\psi_{n+1}+t_{\text{inter,n-1}}\psi_{n-1}\end{split}$
(1)

where $\psi_{n}$ is the amplitude of the $n$-th site. $g_{n}$ and $\gamma_{n}$
are the linear gain and linear loss at the $n$-th site, respectively. Using
explicitly the amplitudes $a_{p}$ and $b_{p}$ of the $A$ and $B$ sites on the
$p$-th unit cell, respectively, in
$\psi(t)=\left(\ldots,a_{p}(t),b_{p}(t),\ldots\right)$, Eq. (1) can be re-
written as:

$\displaystyle i\frac{da_{p}}{dt}$
$\displaystyle=i\left(\frac{g_{A}}{1+|a_{p}|^{2}}-\gamma_{A}\right)a_{p}+t_{\text{intra}}b_{p}+t_{\text{inter}}b_{p-1},$
(2) $\displaystyle i\frac{db_{p}}{dt}$
$\displaystyle=i\left(\frac{g_{B}}{1+|b_{p}|^{2}}-\gamma_{B}\right)b_{p}+t_{\text{intra}}a_{p}+t_{\text{inter}}a_{p+1},$
(3)

where $g_{\sigma}$ and $\gamma_{\sigma}$ are the linear gain and linear loss
at the site $\sigma=A,B$.

In the passive setting ($g_{A}=g_{B}=0$, $\gamma_{A}=\gamma_{B}=0$), this
configuration, with $t_{1}>t_{2}$, is known to give a single topologically
protected zero-energy (non-oscillating) mode localized at the domain wall and
with non-vanishing amplitudes only on the A sites [31]. This is due to the
bulk-boundary correspondence at the domain wall between trivial and non-
trivial topological SSH lattices. Indeed, an SSH lattice is topologically
trivial (non-trivial) if the intra-unit cell coupling is lower (greater) than
the inter-unit cell coupling.

In the active setting, it has been shown that the topological phase of the
system depends on the gain and coupling parameters [12, 13]. As parameters, we
use the values from Ref. [13] with $t_{1}=1$, $t_{2}=0.7$, $g_{B}=0$ and
$\gamma_{A}=\gamma_{B}\equiv\gamma_{AB}$. Figure 1(b) shows the phase diagram,
by varying $g_{A}$ and $\gamma_{AB}$ in the parameter space, for a lattice
composed of $N_{s}=21$ sites ($N=10$). In this configuration, the system has
two distinct topological phases: a non-oscillating phase [white area in Fig.
1(b)] and an oscillating phase [grey area in Fig. 1(b)]. The dynamics of the
two topological phases can be visualized by plotting the total intensity
$I_{A}=\sum_{p}|a_{p}|^{2}$ ($I_{B}=\sum_{p}|b_{p}|^{2}$) of the A (B) sites
in Fig. 1(c) as in Ref. [12, 13]. Alternatively, more details can be
understood by plotting the space-time dynamics of the topological modes as
shown in Figs. 1(d) and 1(e). The non-oscillating phase is similar to the
zero-energy mode in the passive SSH lattice. We can see in Fig. 1(d) that the
mode is localized at the interface and has the majority of its amplitudes on
the A sublattice. On the other hand, the system with saturable gain exhibits a
new topological phase with no counterpart in the passive setting. The new
topological mode is characterized by an edge mode at the domain wall with an
oscillating behaviour of the amplitudes on the A and B sites, as shown in Fig.
1(e).

The classification of the new topological phases in non-linear systems
requires, so far, an expert knowledge of the given non-linear systems, for
example, the known results derived in Ref. [12]. In fact, the phase diagram in
Fig. 1(b) has been obtained solely by the fast Fourier transform of the time
series in the parameter space. Thus, the main aim of this paper is to develop
a tool to explore the topological phases of PTI lasers in more complicated
settings for which we have little knowledge.

The phase diagram shown in Fig. 1(b) will serve as a reference for our
proposed method. The dataset we will utilize throughout this paper is composed
of about $1000$ samples which are randomly generated from the same coupling
and gain parameters’ range as in Fig. 1(b). The coupled-mode equations [Eqs.
(2) and (3)] are integrated using the fourth-order Runge-Kutta method and with
$a_{p}(t=0)=b_{p}(t=0)=0.01$, $\forall p$, as an initial condition. Although
the integration has been performed using a fixed time step $dt=0.01$ until a
final time at $t=1400$, only $2000$ time snapshots are uniformly retrieved in
order to keep the time series at a reasonable size. For the parameters given
above, this sample rate leaves about $10$ time steps per period for the
oscillating regime case [Fig. 1(c)]. The phase diagram is then obtained solely
from the time series of the states within the given parameter space.

###  II.2 Representation classification method

To classify topological lasing modes based on their distinct non-linear
regimes, we use a representation classification method [23, 27, 28]. The
general idea of representation classification relies on the assumption, and
common situations, that the dynamics of a high-dimensional system evolves on a
low-dimensional attractor such as fixed points or periodic orbits [32]. The
low-dimensional structure of the attractor allows for a reduced-order model
that accurately approximates the underlying behaviour of the system: the
dynamics of the complex system can thus be written using a basis that spans
the low-dimensional space. Representation classification consists of
constructing a library of appropriate basis, representative of the dynamical
regimes of interest, and only then employ a filtering strategy to identify the
regime corresponding to a given unknown time series. In the following, we will
use the term “regime” to denote the different dynamical behaviours or the
different topological phases in the non-linear SSH lattice with saturable
gain. Besides, for convenience, we will plot only the total intensity on the A
($I_{A}=\sum_{p}|a_{p}|^{2}$) and B ($I_{B}=\sum_{p}|b_{p}|^{2}$) sublattices
to represent the given regimes. Nevertheless, the time series of the complex
amplitudes at each site will be considered for the construction of the
library.

As is common in complex dynamical systems, the dynamics of a system close to
an attractor lie in a low-dimensional space. This means that a given spatio-
temporal dynamics, denoted by the vector $x(t)$, can be approximately written
in terms of a basis $\Phi=\\{\phi_{i}\\}_{i=1,\ldots,D}$ spanning the low
$D$-dimensional space, namely:

$x(t)\approx\sum_{i=1}^{D}\phi_{i}\beta_{i}(t)=\Phi\beta(t)$ (4)

where $\beta_{i}$ are the weighted coefficients in the above linear
combination of basis states $\phi_{i}$. Using the terminology used in the
literature [23, 27, 28], $x(t)$ will, in the following, be referred to the
state measured at time $t$.

However, one of the main characteristics of non-linear systems is the
drastically different dynamical behaviours with respect to the system’s
parameters. Therefore, the reduced-order modelling strategy using a single
representative basis, _i.e._ , corresponding to a single regime, is bound to
fail. Instead of finding a global basis, we here construct a set of local
bases, _i.e._ , construct a library composed of the bases of each non-linear
regime of interest:

$\mathcal{L}=\\{\Phi_{1},\cdots,\Phi_{J}\\}=\\{\phi_{j,i}\\}_{j=1,\ldots,J\,,\,i=1,\ldots,D},$
(5)

where $J$ is the number of regimes, $\Phi_{j}$’s are the bases of each of the
dynamical regime $j$, and $\phi_{j,i}$’s are the corresponding basis states.
This is the supervised learning part of the method, from which the data-driven
method attempts to capture the dynamics of the system in the reduced-order
model. Therefore, the library $\mathcal{L}$ contains the representative basis
of each regime of interest, and corresponds to an overcomplete basis that
approximates the dynamics of the system across the given parameter space. A
better approximation of $x(t)$, instead of using Eq. (4) for a single basis
regime, then reads:

$x(t)\approx\sum_{j=1}^{J}\sum_{i=1}^{D}\phi_{j,i}\beta_{j,i}(t)=\sum_{j=1}^{J}\Phi_{j}\beta_{j}(t)$
(6)

where $\beta_{j,i}$ are the weighted coefficients in the above linear
combination in the overcomplete basis library $\mathcal{L}$. It is worth
noting that the library modes $\phi_{j,i}$ are not orthogonal to each other,
but instead orthogonal in groups of modes for each different regime $j$.

Throughout this paper, the bases used for constructing the library
$\mathcal{L}$ will be generated by using a time-augmented dynamical mode
decomposition (aDMD) method [28] to consider both the spatial and temporal
behaviours (see Supplementary section SI for additional information).

Here, we use a classification scheme based on a simple hierarchical strategy
[28]. The regime classification approach is fundamentally a subspace
identification problem, where each regime is represented by a different
subspace. Given the state $x(t_{i}:t_{i+N_{w}})$ measured within the time
window $[t_{i},t_{i+N_{w}}]$, with $N_{w}$ the time step window size, the
correct regime $j^{*}$ is identified as the corresponding subspace in the
library $\mathcal{L}$ closest to the measurement in the $L2$-norm sense [28].
In other words, the classification strategy is to find the subspace that
maximises the projection of the measurement onto the regime subspace:

$j^{*}=\arg\max_{j=1,\ldots,J}\lVert P_{j}x(t_{i}:t_{i+N_{w}})\rVert_{2},$ (7)

where $P_{j}$ is a projection operator given by:

$P_{j}=\Phi_{j}\Phi_{j}^{+}$ (8)

with $\Phi_{j}^{+}$ being the pseudo-inverse of $\Phi_{j}$,
$\lVert\cdot\rVert_{2}$ the $L2$-norm of a vector, $\lVert
v\rVert_{2}:=\sqrt{\sum_{i}|v_{i}|^{2}}$, and $\arg\max$ the function that
returns the index of the maximum value.

For the 1D system exemplified in this paper, the data collection is not too
expensive because $N_{s}$ is reasonable. However, if $N_{s}$ is very large,
sparse measurement might be desirable and a slight change in the methodology
is then needed as explained in Supplementary section SII.

###  II.3 Phase Diagrams

In the following we will draw phase diagrams using different approaches to the
representation classification. In the phase diagrams we will mark the
different identified regimes by the color of the dots, where the (dark or
light) blue dots always mark the oscillating regime, the green dots the non-
oscillating regime, the red dots the transient regime and the orange dots the
transition regime. The white and grey areas are overlays of the referenced
phase diagram obtained in Fig. 1. The aDMD bases have been generated with
$N_{w}=25$ from the time series starting at the $1800$-th time step.

####  II.3.1 Fixed library

Figure 2:  Representation classification based on a fixed library. Phase
diagram obtained from the library composed of (a) two regimes (one non-
oscillating and one oscillating) and (b) four regimes (three non-oscillating
and one oscillating). The library is constructed by the magenta dots located
at ($\gamma_{AB},g_{A}-\gamma_{AB})=(0.48,0.06)$ and $(0.16,0.44)$, and for
(b) additionally at $(0.31,0.11)$ and $(0.17,0.24)$.

Figure 2(a) displays the phase diagram derived from a library basis made of
the two topological regimes known from Fig. 1. These two topological modes,
used for the construction of the library, are represented by the two magenta
dots. They are randomly chosen from the known regimes’ region in Fig. 1(b) and
the details of the resulted phase diagrams is dependent from that random
choice. The remaining coloured dots in the plot represent the identified
regime $j^{*}$ [Eq. (7)] of each sample depicted by the dots in the parameter
space. However, we can see that the phase diagram fails to correctly predict
all the dynamical behaviours. Indeed, we observe that many time series are not
correctly identified. Using a different choice of topological regimes in the
parameter space to construct the library could be a solution to find a better
phase diagram, but our attempts only showed marginal improvement of the
agreement. Testing every possible choices in the dataset with no guarantee of
finding the accurate phase diagram is therefore not a practical solution.
Instead, using four bases in the library instead of two bases, or equivalently
considering four regimes from the given parameter space, the phase diagram in
Fig. 2(b) shows better results: The identified oscillating and non-oscillating
regimes are more separated and have a better fitting with the referenced phase
diagram, even though they belong to a distinct regime index $j^{*}$. Merging
the three oscillating regimes present in the library would then give a more
accurate derived phase diagram. Therefore, Figure 2 suggests that increasing
the number of bases in the library $\mathcal{L}$ and merging some of them
might help to get closer to the desired phase diagram, as we will see in the
later sections.

####  II.3.2 Top-down adaptive library

In the previous section, the construction of the library $\mathcal{L}$ was a
manual process from which we already knew the different dynamical regimes.
This, however, requires prior knowledge of the complex system considered. The
strategy, here, is to adaptively construct the library based on the given data
samples. Here, we employ a top-down approach in which we start with too many
samples for the construction of the library, and then reduce the library size
by merging some of them. Based on some measures in the decision process, this
automated construction of the library thus removes the manual construction of
the regimes.

The underlying assumption of the classification scheme is based on the
dissimilarity between the subspace of different regimes. We thus propose to
consider regimes that are similar as equivalent regimes. This would, for
instance, help us to merge the three phases in the non-oscillating region in
Fig. 2(b), and consider them as a single regime. In other words, the regimes
$i$ and $j$ are said to be equivalent, denoted by $i\sim j$, if the fraction
of information retained after the projection onto each other,
$\gamma_{ij}\in[0,1]$, is high enough:

$\gamma_{ij}>\gamma_{th},$ (9)

where $\gamma_{th}\in[0,1]$ is the hyper-parameter which decides the threshold
value for merging different regimes and $\gamma_{ij}$ is the subspace
alignment given by:

$\gamma_{ij}:=\frac{\lVert P_{i}P_{j}\rVert_{F}^{2}}{\lVert
P_{i}\rVert_{F}\lVert P_{j}\rVert_{F}},$ (10)

with $\lVert\cdot\rVert_{F}$ the Frobenius norm of a matrix, $\lVert
M\rVert_{F}:=\sqrt{\sum_{ij}|M_{ij}|^{2}}$. Importantly, the relation Eq. (9)
is numerically computed in such a way that the transitivity property of the
equivalence relation is satisfied, namely that if $i\sim j$ and $j\sim k$ then
$i\sim k$. The relation Eq. (9) is then indeed an equivalence relation because
the reflexive ($i\sim i$) and symmetric ($i\sim j\Rightarrow j\sim i$)
property of the relation is automatically satisfied from the definition of
$\gamma_{ij}$ [Eq. (10)]. Supplementary section SIII gives more details on the
top-down library generation principle.

Figure 3:  Representation classification based on a top-down adaptive library.
(a) Phase diagram obtained using the top-down classification strategy with an
initial library composed of $J=60$ regimes randomly chosen and
$\gamma_{th}=0.75$. (b) Representative total intensity $I_{A}$ (and $I_{B}$)
of the A (and B) sublattice in blue (and orange) for the different regimes.
The black vertical dotted line indicates the starting time from which the
bases are generated.

The top-down representation classification strategy is to classify the time
series according to a large library of bases, and only then merge the
equivalent identified regimes via the calculated alignment subspace
$\gamma_{ij}$ and the equivalence relation Eq. (9). Figure 3 shows the phase
diagram obtained from the top-down algorithm that started with a library
composed of $J=60$ regimes randomly chosen, along with the representative of
each phase. We observe in Fig. 3(a) that the derived phase diagram is able to
distinguish the non-oscillating [top panel of Fig. 3(b)] and oscillating
regimes [middle panel of Fig. 3(b)]. In addition, a third regime corresponding
to a transient regime [bottom panel of Fig. 3(b)] is found close to the
$\gamma_{AB}=0$ or $g_{A}-\gamma_{AB}=0$ axis. This transient regime indicates
that a longer simulation time might be needed to be considered either in the
oscillating or non-oscillating regimes.

However, we can see that the derived phase diagram is still failing in the low
$\gamma_{AB}$ and low $g_{A}-\gamma_{AB}$ region (bottom-left region of the
present phase diagram), where some time series are interpreted as non-
oscillating instead of oscillating regime. This shows the limitation of this
method where the initially constructed library may lack some of the paths that
may connect similar bases. For example, regimes $i$ and $k$ might not be
similar enough to be considered as equivalent directly [Eq. (9)] but are both
equivalent to the regime $k$, _i.e._ , $i\sim j$ and $j\sim k$, which is
missing in the library. The natural workaround would be to increase the
initial library size and ensure that the regimes in the library have no
missing paths, as we will see in the next section.

Figure 4:  $\gamma_{th}$-hyper-parameter dependency. (a) Library size against
$\gamma_{th}$. Phase diagrams derived using the top-down classification
strategy with (a) $\gamma_{th}=0.55$ and (b) $\gamma_{th}=0.95$. The initial
library is composed of $J=60$ regimes randomly chosen.

The hyper-parameter $\gamma_{th}$ is an important quantity in the algorithm
since it dictates which regimes are equivalent or not. A low threshold
$\gamma_{th}$ will easily merge regimes while a high $\gamma_{th}$ will barely
reduce the size of the library as depicted in Fig. 4(a). The threshold is here
arbitrarily chosen based on Fig. 4, and based on the refinement of the desired
library. For example, we can see in Fig. 4(b) that the derived phase diagram
with $\gamma_{th}=0.55$ has two different phases. For this coarse threshold,
the transient regime is identified but the distinct non-oscillating and
oscillating phases are merged together into a single phase. On the other hand,
with the same library as in Fig. 4(b), Figure 4(c) displays the obtained phase
diagram for a finer threshold value $\gamma_{th}=0.95$. The plot shows that
the algorithm separates the parameter space into several regimes which can be
grouped into four main regimes. In addition to the non-oscillating, the
oscillating and the transient regimes, there is a regime corresponding to the
transition between the two topological phases. Besides, this finer description
allows us to see distinct sets of modes in the oscillating parameter space
region [dark blue and light blue dots in Fig. 4(c)]. Nevertheless, we observe
again that the initial library misclassifies some of the non-oscillating time
series most likely because of some missing paths, as said previously.

####  II.3.3 Bottom-up adaptive library

We propose an alternative and dual approach which considers fewer samples in
the library. The core idea of this bottom-up approach is then to add samples
on the fly during the classification of the given sample if the library is not
good enough.

Here, the library is considered to be good enough if the maximal projection of
the measurement onto the regimes’ subspace is high enough. In other words, the
library is said to be good enough if the worst relative reconstruction error,
$\epsilon$, is low enough:

$\epsilon<\epsilon_{th},$ (11)

where $\epsilon_{th}$ is the hyper-parameter which decides the threshold
quality of the library and

$\epsilon:=\max_{j=1,\ldots,J}\frac{\lVert P_{j}y(t)-y(t)\rVert_{2}}{\lVert
y(t)\rVert_{2}},$ (12)

with $\lVert\cdot\rVert_{2}$ the $L_{2}$-norm of a vector. Supplementary
section SIV gives more details on the bottom-up library generation principle.

The advantage of this bottom-up approach is the full exploration of the
parameter space region and the automatic construction of a library based on
its quality. This method does not suffer from the randomly chosen samples used
to construct the library, and the library composition is not restricted to a
narrow parameter space region. Using a good enough library quality, the
algorithm should therefore be able to sort the missing paths issue in the top-
down method.

Figure 5:  Representation classification based on a bottom-up adaptive
library. (a) Phase diagram obtained using the bottom-up classification
strategy with a starting library composed of a single regime randomly chosen,
$\epsilon_{th}=0.005$ and $\gamma_{th}=0.75$. (b) Representative total
intensity $I_{A}$ (and $I_{B}$) of the A (and B) sublattice in blue (and
orange) for the different regimes.

The bottom-up representation classification scheme consists of classifying the
time series according to a given library or adding this sample into the
library if the library is not good enough, and only then merging the different
phases obtained into groups of equivalent regimes using the top-down method.
Figure 5 depicts the phase diagram derived from the bottom-up classification
algorithm with a starting library composed of a single regime. Similarly to
the top-down approach in Fig. 3(a), we observe three distinct regimes
corresponding to the non-oscillating, oscillating and transient regimes [Fig.
3(b)]. Nevertheless, the obtained phase diagram now better predicts the
regimes. The misclassifications of the non-oscillating and oscillating
regimes, which were due to missing paths in the library, are now reduced, and
only very few dots are not correctly identified due to being close to the
topological transition boundary. Likewise, the transient points are
indications of longer simulations needed because of the long transient time.

Figure 6:  $\epsilon_{th}$-hyper-parameter dependency. (a) Library size
against $\epsilon_{th}$. The inset is a zoom-in of the plot. Phase diagram
derived using the bottom-up classification strategy with $\gamma_{th}=0.95$
and (b) $\epsilon_{th}=0.005$ and (c) $\epsilon_{th}=0.05$. The initial
library is composed of a single regime randomly chosen.

Along with the $\gamma_{th}$ hyper-parameter, the threshold hyper-parameter
$\epsilon_{th}$ is an important parameter since it tells us whether we want to
add or not a given sample into the library. We observe in Fig. 6 that a low
threshold $\epsilon_{th}$ will add many samples to the library, whereas a high
$\epsilon_{th}$ will not add samples to the library at all. The threshold
value $\epsilon_{th}$ is, again, arbitrarily chosen but with a preference for
a high-quality library, _i.e._ , low $\epsilon_{th}$, in order to avoid
missing paths. For example, we can see in Fig. 6(c) the phase diagram derived
using $\epsilon_{th}=0.05$ (and $\gamma_{th}=0.95$), namely with a library
that gives less than $5\%$ of the reconstruction error of the measurement.
This set of hyper-parameter gives four main regimes that seem to correspond to
the non-oscillating, oscillating, transition and transient regimes. Yet, there
is some misidentification of the two topological phases most likely because of
missing paths of the obtained library. On the other hand, with a better
library quality, here $\epsilon_{th}=0.005$ (and $\gamma_{th}=0.95$), the
missing paths are retrieved and the derived phase diagram correctly predicts
the topological phases [Fig. 6(b)]. Figure 6(b) shows that the different
regimes, obtained previously with a lower library quality, are now better
defined. The non-oscillating and oscillating regimes are well located in their
respective parameter space region, and the transition points follow the
transition boundary between the two topological phases. In addition, the
bottom-up representation classification is predicting distinct oscillating
modes [dark blue and light blue dots in Fig. 6(b)].

The presence of distinct oscillating modes is an example of new insights given
by the data-driven classification method. Indeed, the complex values of the
amplitudes of $x(t)$ are, here, taken into account instead of solely the total
intensity of each sublattice $A$ and $B$ as in Ref. [12, 13]. This allows for
a finer description of the dynamic pattern based on the whole lattice with the
relative phase difference of the sites or the absolute value of amplitudes.

##  III Conclusion

We have proposed a data-driven approach to identify and classify topological
phases of dynamical systems. By utilizing the representation classification
strategy based on the aDMD, we have successfully drawn the phase diagrams of
the domain-wall-type SSH lattice with saturable gain. To avoid manual
labelling in the classification, we have proposed two automatic library
construction schemes: top-down and bottom-up approaches that merge similar
phases in a library or adaptively construct a library according to its
quality, respectively. While the bottom-up method is preferred due to the
missing paths issue, both methods allow exploration of parameter space without
any expert knowledge of the complex non-linear system. Our approach is
advantageous in doing the reverse engineering to find new phases because only
a small desired parameter space region is required for solving the laser rate
equation, applying the bottom-up representation classification strategy and
building a starting library. Our method opens the door for finding novel
topological lasing modes by providing new insights into the dynamics of
coupled lasers in more complicated settings.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Code availability

The code that support the findings of this study are available from the
corresponding author upon reasonable request.

## acknowledgments

We acknowledge the support of the Sêr Cymru II Rising Star Fellowship
(80762-CU145 (East)) which is part-funded by the European Regional Development
Fund (ERDF) via the Welsh Government.

## Author contributions

S.W., D.R. and S.S.O. conceived the project. S.W developed the theoretical
model and performed all the numerical simulations. All authors contributed to
the discussions and the preparation of the manuscript.

## Competing interests

The authors declare that they have no competing interests.

## Additional information

Correspondence and requests for materials should be addressed to Sang Soon Oh.

##  SI Basis generation methods

This supplementary section review few methods for generating a basis from a
time series. There are many ways to generate the low-order model of a given
dynamical behaviour. An important quantity is the so-called data matrix $X$
built from the data at hand. The data matrix is a $(N_{s}\times N_{t})$-matrix
that collects the $N_{t}$ data snapshots $x(t_{i})$ into columns:

$X=\left[x(t_{0}),x(t_{1}),\ldots,x(t_{N_{t}})\right].$ (S1)

Here, the vector $x(t_{i})$ is chosen to be the complex-valued amplitudes of
the modes at the A and B sites. Other “observables”, such as the absolute
values or the total intensity per sublattices, can be used. This may give
different (better or worse) results and is left for a future study. The bases
can then be constructed using dimensional reduction techniques on the data
matrix. Here, we will cover different methods in order to highlight their
importance and limitations for the classification scheme used.

###  SI.1 Proper orthogonal decomposition

Figure S1:  POD method. Singular values of the (a) non-oscillating and (b)
oscillating regimes. The red dots correspond to the singular values
accumulating $99\%$ of the total variance of the data. (c) Real and (d)
imaginary parts of the field profile of the first POD mode for the non-
oscillating regime. (e) Real and (f) imaginary parts of the field profile of
the (top) first and (bottom) second POD mode for the oscillating regime. The
POD bases have been generated from the time series starting at the $1800$-th
time step.

Proper orthogonal decomposition (POD) [24] is a commonly used tool for
dimensional reduction of physical systems. This decomposition relies on the
singular value decomposition (SVD) of the data matrix, given by:

$X=U\Sigma V^{\dagger}$ (S2)

where $U$ and $V^{\dagger}$ are $(N_{s}\times N_{s})$ and $(N_{s}\times
N_{t})$ unitary matrices, respectively. $\Sigma$ is a diagonal $(N_{0}\times
N_{0})$-matrix $\text{diag}=(\sigma_{1},\ldots,\sigma_{N_{0}})$, with
$N_{0}=\min(N_{s},N_{t})$. The diagonal entries of $\Sigma$ are the so-called
singular values and are ordered in ascendant order
$\sigma_{1}>\sigma_{2}>\ldots>\sigma_{N_{0}}\geq 0$. The SVD gives us two
orthonormal bases $U$ and $V^{\dagger}$ since the matrices $U$ and
$V^{\dagger}$ are unitary matrices. The columns of $U$ are ordered according
to the variance $\sigma_{i}$ they capture in the data matrix and are called
the singular vectors: these are the POD modes that are used in the basis
$\Phi$. Moreover, the POD modes are complex because of the complex data $X$.

For a low-dimensional attractor, the POD basis can be safely truncated at a
cut-off value $r$ while retaining the main information of the data matrix.
Explicitly, the SVD reads:

$X_{im}=\sum_{n=0}^{N_{0}}U_{in}\sigma_{n}V^{\dagger}_{nm},$ (S3)

and keeping only the $r$ highest terms in the decomposition [Eq. S3], we have
the approximation:

$X_{im}\simeq\sum_{n=0}^{r}U_{in}\sigma_{n}V^{\dagger}_{nm}.$ (S4)

This is re-written, in a matrix form, as:

$X\simeq U_{r}\Sigma_{r}V^{\dagger}_{r}$ (S5)

where $U_{r}$, $\Sigma_{r}$ and $V^{\dagger}_{r}$ are the truncated matrix of
$U$, $\Sigma$ and $V^{\dagger}$, respectively. Although the cut-off value $r$
can be chosen based on different criteria [33], $r$ is typically chosen so
that the POD modes retain a certain amount of the variance (or energy)
$\sigma_{X}$ in the data, namely:

$\sum_{i=0}^{r}\sigma_{i}>\sigma_{X}.$ (S6)

Figure S1 displays the POD method of the non-oscillating and oscillating
regimes in the domain-wall SSH lattice with saturable gain [Fig. 1]. The
truncation has been chosen such that $\sigma_{X}=99\%$ of the total variance
is retained. We observe in Figs. S1(a) and S1(b) the normalised singular
values, and that a single POD mode is retained for the zero-mode-like, whereas
three POD modes are needed for the oscillating regime, as marked by the red
dots. The real and imaginary parts of the field profile of the corresponding
first few POD modes are plotted in Figs. S1(c) and S1(d) and Figs. S1(e) and
S1(f) for the zero-mode-like non-oscillating and oscillating regimes,
respectively. One can see that the main spatial feature of the zero-mode is
captured in the single POD mode obtained after truncation, where the majority
of its amplitudes are on the A sublattice. On the other hand, the POD modes of
the oscillating dynamical regime also capture part of the information with
some finite amplitudes on the A and B sublattices.

Importantly, in this decomposition [Eq. S3], SVD is implicitly doing a space-
time separation of the data matrix, where the POD modes $U$ contain the
spatial information while $V$ have the temporal information at each spatial
grid point. Therefore, the POD modes give a static basis and do not explicitly
model the temporal dynamics of the time series. This method will therefore
most likely fail to identify the correct dynamical regime in the
classification step.

###  SI.2 Dynamical mode decomposition

Figure S2:  DMD method. (a) Singular values of the oscillating regime in the
DMD. The red dots correspond to the singular values accumulating $99\%$ of the
total variance of the data. (b) Plot of the logarithm of the DMD values
$\ln(\lambda)$ in the complex plane. The size of the open circle is
proportional to their corresponding singular values. (c) Real and (d)
imaginary parts of the field profile of the (top) first and (bottom) second
DMD mode for the oscillating regime. The DMD basis has been generated from the
time series starting at the $1800$-th time step.

Dynamical mode decomposition (DMD) [34, 35, 24] is an alternative to the POD
method for learning the dynamics of non-linear systems. DMD is an explicitly
temporal decomposition that takes the sequences of snapshots into account, and
is able to derive the spatio-temporal patterns of the data matrix $X$. The
dynamics of the system are taken into account by considering a linear matrix
$A$ that maps a data matrix $X_{1}$, starting at some time steps $(t_{1})$, to
the data matrix $X_{2}$, starting at the next time step ($t_{2}$). The matrix
$A$ is thus defined as:

$X_{2}=AX_{1},$ (S7)

and the corresponding data matrices are given by:

$X_{1}=\left[x(t_{1}),x(t_{2}),\ldots,x(t_{N_{t}-1})\right]$ (S8)

and

$X_{2}=\left[x(t_{2}),x(t_{3}),\ldots,x(t_{N_{t}})\right].$ (S9)

Interestingly, Equation S7 is similar to a linear stability analysis
formulation for discrete maps if we think of the stability matrix as the
linear matrix $A$. The DMD method thus consists of solving the following
eigenvalues problem:

$A\Phi=\Phi\Lambda$ (S10)

where the columns of $\Phi$ are the DMD modes $\phi_{i}$ and the corresponding
DMD eigenvalues $\lambda_{i}$ are the diagonal entry of $\Lambda$. The DMD
modes $\phi_{i}$ give us the spatial eigenmodes while their corresponding
eigenvalues $\lambda_{i}$ have their temporal information. Using a change of
units from the data snapshots, observed at every $\Delta t$, to units in time,
the eigenvalues are complex-valued scalars:

$\frac{\ln\left(\lambda_{i}\right)}{\Delta t}=\mu_{i}+i\omega_{i}$ (S11)

where $\mu_{i}$ gives the growth (decay) rate if $\mu_{i}>0$ ($\mu_{i}<0$) and
$\omega_{i}$ the oscillation frequency of the DMD modes $\phi_{i}$.

However, the size of the data matrix usually makes the eigendecomposition not
feasible. The goal, here, is therefore to approximate the eigenvalues and
eigenvectors of A, using only the data matrices $X_{1}$ and $X_{2}$. The idea
is to start by the truncated SVD of $X_{1}=U_{r}\Sigma_{r}V_{r}^{\dagger}$ in
which Eq. S7 becomes:

$X_{2}=AU_{r}\Sigma_{r}V_{r}^{\dagger}.$ (S12)

Then the linear matrix $A$ is reduced by considering its projection onto the
truncated POD subspace:

$A_{r}:=U_{r}^{\dagger}AU_{r}=U_{r}^{\dagger}X_{2}V_{r}\Sigma^{-1}.$ (S13)

The eigenvalue problem for $A_{r}$ is solved with:

$A_{r}W=W\Lambda,$ (S14)

from which we have the relation:

$\Phi=X_{2}V\Sigma^{-1}W.$ (S15)

The key feature of the DMD method is that it decomposes the data into a set of
coupled spatio-temporal modes. The DMD resembles a mixture of the POD in the
spatial domain and the discrete Fourier transform (DFT) in the time domain.
Figure S2 shows the DMD results for the oscillating regime. We can, indeed,
see in Fig. S2(a) that the singular values are similar to that of the POD.
Besides, we observe in Figs. S2(c) and S2(d) that the field profile of the DMD
modes closely resembles the POD modes in Fig. S1. The largest DMD modes not
only look similar to the POD modes, but they also contain the oscillation
frequencies from $\omega_{i}$, as in DFT. The DMD even goes beyond DFT by
giving an estimate of the growth (decay) rate in time via $\mu_{i}>0$
($\mu_{i}<0$). This can be seen by plotting the DMD modes, scaled by their
contribution in the decomposition $\sigma_{i}$, in the frequency plane of
$\lambda_{i}$. We can see in Fig. S2(b) that the dynamical regime has a single
DMD mode with $\omega_{i}=0$ akin to the offset of the oscillation amplitudes,
and two DMD modes with opposite $\omega_{i}\neq 0$ corresponding to the
oscillating behaviours. All the above three DMD modes have vanishing growth or
decay rate $\mu_{i}=0$.

###  SI.3 Time-augmented dynamical mode decomposition

Figure S3:  aDMD method. (a) Singular values of the oscillating regime in the
aDMD. The red dots correspond to the singular values accumulating $99\%$ of
the total variance of the data. (b) Plot of the logarithm of the aDMD values
$\ln(\lambda)$ in the complex plane. The size of the open circle is
proportional to their corresponding singular values. (c) Real and (d)
imaginary parts of the “field profile” of the (top) first and (bottom) second
aDMD mode for the oscillating regime. Here, by “Site n”, we mean the $n$-th
entry of the eigenvector. The aDMD basis has been generated with $N_{w}=25$
from the time series starting at the $1800$-th time step.

Although the DMD gives the temporal behaviours of the non-linear system, the
temporal information is not fully incorporated into the DMD basis $\Phi$ since
only the DMD modes are used. Exploiting the time evolution in the dynamical
regime requires the use of DMD modes along with their eigenvalues. The idea is
therefore to incorporate the dynamic information by augmenting the basis
$\Phi$ [28]. This time-augmented DMD will be denoted by aDMD in the remaining
of this chapter. Using the defining relation of the eigenvalues $\lambda_{i}$
as similar to a time evolution operator, _i.e._ multiplying by $\lambda$ is
the same as shifting by one time step, we have for a given DMD mode
$\phi_{i}$, its evolution given by $\lambda_{i}^{N_{w}}\phi_{i}$ at $N_{w}$
time step ahead in time. Therefore, the time-augmented basis vector reads:

$\left[\begin{array}[]{c}\phi_{i}\\\ \lambda_{i}\phi_{i}\\\ \vdots\\\
\lambda_{i}^{N_{w}}\phi_{i}\end{array}\right].$ (S16)

By considering a time window $N_{w}$, the time-augmentation of the DMD basis
provides us with the dynamical information of the non-linear regime. Figure S3
shows the aDMD results for the same oscillating regime as in Fig. S2. We can
see that the singular values and aDMD eigenvalues plots [Figs. S3(a) and
S3(b)] are the same as in the DMD algorithm [Figs. S2(a) and S2(b)] whereas
the “field profile” of the aDMD modes [Figs. S3(c) and S3(d)] carry some
temporal evolution information. In particular, we observe in the top panel of
Figs. S3(c) and S3(d) the static behaviours of the aDMD mode with
$\omega_{i}=0$. On the other hand, we can see, in the bottom panel of Figs.
S3(c) and S3(d), one of the first aDMD modes with $\omega_{i}\neq 0$ featuring
some oscillating behaviour in time. The size of the basis mode is larger than
the plain DMD, and can exhibit its time evolution. Nevertheless, the graphs do
not exactly plot the temporal evolution of the DMD modes since the first
$N_{s}$ entry of the basis state is for the $N_{s}$ sites; the next $N_{s}$
for again the $N_{s}$ sites but at the next time step, etc.

###  SI.4 Classification results from different decomposition methods

Figure S4:  Phase diagram derived using different decomposition methods. Phase
diagrams obtained from a library composed of two regimes (one non-oscillating
and one oscillating) from which the bases have been generated using the (a)
POD, (b) DMD, (c) aDMD with $N_{w}=5$ and (d) aDMD with $N_{w}=25$. The green
and blue dots correspond respectively to the identified non-oscillating and
oscillating regimes. The magenta dots represent the regimes used for the
construction of the library. These are located at
($\gamma_{AB},g_{A}-\gamma_{AB})=(0.48,0.06)$ and $(0.16,0.44)$. The white and
grey areas are overlays of the referenced phase diagram obtained in Fig. 1.
The bases have been generated from the time series starting at the $1800$-th
time step.

The classification results [Eq. (7)] from the decomposition methods reviewed
is shown in Fig. S4. We observe in Figs. S4(a) and S4(b) that the phase
diagrams fail to correctly predict the distinct dynamical regimes. This is
expected since the POD or DMD modes do not contain enough information about
the temporal behaviours. Besides, the classification for these diagrams is
based on a single snapshot ($N_{w}=0$). Thus, it is expected the
classification fails to capture the correct dynamics since a single snapshot
only relies on the spatial pattern of the regime. On the other hand, we can
see in Figs. S4(c) and S4(d) that the derived phase diagrams have better
accuracy when the bases are time-augmented, or equivalently when the
classification scheme uses several snapshots. By increasing the time window in
the classification, the derived phase diagram is even better as illustrated in
Figs. S4(c) and S4(d) for $N_{w}=5$ and $N_{w}=25$, respectively. The phase
diagram will get improved until the time window is large enough to capture the
dynamic behaviour.

##  SII Sparse measurement

This section detailed the slight change in the methodology in case of sparse
sensing. Sparse sensing is often desirable since the measurement and the data
collection can be expensive for a complex system if the space grid is too
fine, _i.e._ if $N_{s}$ is very large. The compressed measurement $y(t)$ is
derived from the full-state measurement $x(t)$ and the measurement matrix $C$:

$y(t)=Cx(t),$ (S17)

where $C$ is a matrix of size ($N_{p}\times N_{s}$) with $N_{p}$ the number of
measurements. Although the measurement matrix $C$ can be represented by some
advance and complex mapping [36], here we focus on point-wise measurements,
namely the $C_{ij}$ entry in the matrix measurement corresponds to the $i$-th
measurement at the $j$-th spatial grid point. Therefore the compressed basis
is given by:

$\Theta=C\Phi$ (S18)

where $\Phi$ is the basis obtained from the full-state data collection. The
library of bases for the $J$ distinct dynamical regimes is similarly re-
written as:

$\mathcal{L}=\\{\Theta_{1},\ldots,\Theta_{J}\\}.$ (S19)

Nevertheless, the size of the current SSH lattice is, here, reasonable and
allows us to choose the matrix measurement as the identity matrix
$C=1_{N_{s}}$. We will thus use the full-state instead of sparse measurements,
but retain the $\Theta$ and $y(t)$ notations to keep the general formalism.

##  SIII Top-down library generation principle

Figure S5:  Principle of the adaptively top-down library generation. (a)
Subspace alignment matrix for an initial library composed of $J=60$ regimes.
(b) Subspace alignment matrix from the same library as in (a) but grouped into
equivalent regimes. (c) Subspace alignment matrix from the reduced library
after the equivalent regimes are merged with $\gamma_{th}=0.75$ in (b). (d)
Representative total intensity $I_{A}$ (and $I_{B}$) of the A (and B)
sublattice in blue (and orange) for the different regimes. The black vertical
dotted line indicates the starting time from which the bases are generated.
The aDMD bases have been generated with $N_{w}=25$ from the time series
starting at the $1800$-th time step.

This section supplement the top-down library generation principle. The general
idea of the top-down construction of the library is to start with a library
made of a high number of bases, generated from the time series randomly chosen
in the given parameter space region, and then merge them into groups of
equivalent regimes. Figure S5 illustrates the top-down algorithm. Starting
with a library composed of $J$ bases (here $J=60$), the subspace alignment
$\gamma_{ij}$ is computed [Fig. S5(a)] and then grouped into equivalent
regimes according to Eq. (9) in the main text [Fig. S5(b)]. Each
representative of the regimes is then randomly selected within each group
[Fig. S5(c)]. We plot in Fig. S5(d) the time series of the representative of
each regime. These regimes are the non-oscillating and oscillating regimes, as
well as a third regime which may correspond to a transient regime. The
vertical dashed line in the time series represents the initial time used for
constructing the bases in the library.

##  SIV Bottom-up library generation principle

This section supplement the bottom-up library generation principle. Figure S6
illustrates the bottom-up methodology proposed. We start with a single sample
in the library, randomly chosen in the parameter space region [Fig. S6(a)].
The library is then adaptively constructed according to the relative
reconstruction error $\epsilon$ [Fig. S6(b)]. Finally, with the large library
at hand, the top-down approach is used to reduce the library size by merging
equivalent regimes [Fig. S6(c)]. The representative of the regimes is plotted
in Fig. S6(d) and corresponds to the non-oscillating, oscillating and
transient regimes, respectively.

Figure S6:  Principle of the adaptively bottom-up library generation. (a)
Total intensity $I_{A}$ (and $I_{B}$) on the A (and B) sublattice in blue (and
orange) for the single and randomly chosen regime composing the library. (b)
Subspace alignment matrix from the adaptively constructed library with the
single initial regime in (a) and with $\epsilon_{th}=0.01$. (c) Subspace
alignment matrix from the reduced library after the equivalent regimes are
merged with $\gamma_{th}=0.75$ in (b). (d) Representative total intensity
$I_{A}$ (and $I_{B}$) of the A (and B) sublattice in blue (and orange) for the
different regimes. The black vertical dotted line in (a) and (d) indicates the
starting time from which the bases are generated. The aDMD bases have been
generated with $N_{w}=25$ from the time series starting at the $1800$-th time
step.

## References

  * Schomerus [2013] H. Schomerus, Topologically protected midgap states in complex photonic lattices, Opt. Lett. 38, 1912 (2013).
  * Poli  _et al._ [2015] C. Poli, M. Bellec, U. Kuhl, F. Mortessagne, and H. Schomerus, Selective enhancement of topologically induced interface states in a dielectric resonator chain, Nat. Commun. 6, 6710 (2015).
  * Bahari  _et al._ [2017] B. Bahari, A. Ndao, F. Vallini, A. El Amili, Y. Fainman, and B. Kanté, Nonreciprocal lasing in topological cavities of arbitrary geometries, Science 358, 636 (2017).
  * Leykam  _et al._ [2017] D. Leykam, K. Y. Bliokh, C. Huang, Y. D. Chong, and F. Nori, Edge modes, degeneracies, and topological numbers in non-Hermitian systems, Phys. Rev. Lett. 118, 28 (2017).
  * Zhang  _et al._ [2018] P. Zhang, H. Shen, and H. Zhai, Machine learning topological invariants with neural networks, Phys. Rev. Lett. 120, 066401 (2018).
  * Parto  _et al._ [2018] M. Parto, S. Wittek, H. Hodaei, G. Harari, M. A. Bandres, J. Ren, M. C. Rechtsman, M. Segev, D. N. Christodoulides, and M. Khajavikhan, Edge-mode lasing in 1D topological active arrays, Phys. Rev. Lett. 120, 113901 (2018).
  * Takata and Notomi [2018] K. Takata and M. Notomi, Photonic topological insulating phase induced solely by gain and loss, Phys. Rev. Lett. 121, 213902 (2018).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and topology in non-Hermitian physics, Phys. Rev. X 9, 041015 (2019).
  * Shao  _et al._ [2020] Z.-K. Shao, H.-Z. Chen, S. Wang, X.-R. Mao, Z.-Q. Yang, S.-L. Wang, X.-X. Wang, X. Hu, and R.-M. Ma, A high-performance topological bulk laser based on band-inversion-induced reflection, Nat. Nanotechnol. 15, 67 (2020).
  * Longhi [2018] S. Longhi, Non-Hermitian gauged topological laser arrays, Ann. Phys. 530, 1800023 (2018).
  * Longhi  _et al._ [2018] S. Longhi, Y. Kominis, and V. Kovanis, Presence of temporal dynamical instabilities in topological insulator lasers, EPL 122, 14004 (2018).
  * Malzard  _et al._ [2018] S. Malzard, E. Cancellieri, and H. Schomerus, Topological dynamics and excitations in lasers and condensates with saturable gain or loss, Opt. Express 26, 22506 (2018).
  * Malzard and Schomerus [2018] S. Malzard and H. Schomerus, Nonlinear mode competition and symmetry-protected power oscillations in topological lasers, New J. Phys. 20, 063044 (2018).
  * Wong and Oh [2021] S. Wong and S. S. Oh, Topological bulk lasing modes using an imaginary gauge field, Phys. Rev. Res. 3, 033042 (2021).
  * Gong  _et al._ [2020] Y. Gong, S. Wong, A. J. Bennett, D. L. Huffaker, and S. S. Oh, Topological insulator laser using valley-Hall photonic crystals, ACS Photonics 7, 2089 (2020).
  * Araki  _et al._ [2019] H. Araki, T. Mizoguchi, and Y. Hatsugai, Phase diagram of a disordered higher-order topological insulator: A machine learning study, Phys. Rev. B 99, 085406 (2019).
  * Peano  _et al._ [2021] V. Peano, F. Sapper, and F. Marquardt, Rapid exploration of topological band structures using deep learning, Phys. Rev. X 11, 021052 (2021).
  * Pedregosa  _et al._ [2011] F. Pedregosa, R. Weiss, M. Brucher, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and É. Duchesnay, Scikit-learn: Machine learning in Python, J. Mach. Learn. Res. 12, 2825 (2011).
  * Wang [2016] L. Wang, Discovering phase transitions with unsupervised learning, Phys. Rev. B 94, 195105 (2016).
  * Scheurer and Slager [2020] M. S. Scheurer and R.-J. Slager, Unsupervised machine learning and band topology, Phys. Rev. Lett. 124, 226401 (2020).
  * Jolliffe [2005] I. Jolliffe, Principal component analysis, in  _Encycl. Stat. Behav. Sci._ (John Wiley & Sons, Ltd, Chichester, UK, 2005).
  * Wright  _et al._ [2009] J. Wright, A. Yang, A. Ganesh, S. Sastry, and Yi Ma, Robust face recognition via sparse representation, IEEE Trans. Pattern Anal. Mach. Intell. 31, 210 (2009).
  * Brunton  _et al._ [2014] S. L. Brunton, J. H. Tu, I. Bright, and J. N. Kutz, Compressive sensing and low-rank libraries for classification of bifurcation regimes in nonlinear dynamical systems, SIAM J. Appl. Dyn. Syst. 13, 1716 (2014).
  * Proctor  _et al._ [2014] J. L. Proctor, S. L. Brunton, B. W. Brunton, and J. N. Kutz, Exploiting sparsity and equation-free architectures in complex systems, Eur. Phys. J. Spec. Top. 223, 2665 (2014).
  * Fu  _et al._ [2014] X. Fu, S. L. Brunton, and J. Nathan Kutz, Classification of birefringence in mode-locked fiber lasers using machine learning and sparse representation, Opt. Express 22, 8585 (2014).
  * Brunton  _et al._ [2016] B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz, Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition, J. Neurosci. Methods 258, 1 (2016).
  * Bright  _et al._ [2016] I. Bright, G. Lin, and J. N. Kutz, Classification of spatio-temporal data via asynchronous sparse sampling: Application to flow around a cylinder, Multiscale Model. Simul. 14, 823 (2016).
  * Kramer  _et al._ [2017] B. Kramer, P. Grover, P. Boufounos, S. Nabi, and M. Benosman, Sparse sensing and DMD-based identification of flow regimes and bifurcations in complex flows, SIAM J. Appl. Dyn. Syst. 16, 1164 (2017).
  * Ozan  _et al._ [2021] D. E. Ozan, A. Iannelli, M. Yin, and R. S. Smith, Regularized classification and simulation of bifurcation regimes in nonlinear systems, IFAC-PapersOnLine 54, 352 (2021).
  * Su  _et al._ [1979] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Solitons in polyacetylene, Phys. Rev. Lett. 42, 1698 (1979).
  * Asbóth  _et al._ [2016] J. K. Asbóth, L. Oroszlány, and A. Pályi,  _Lect. Notes Phys._, Lecture Notes in Physics, Vol. 919 (Springer International Publishing, Cham, 2016) p. 166.
  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65, 851 (1993).
  * Gavish and Donoho [2014] M. Gavish and D. L. Donoho, The optimal hard threshold for singular values is 4/$\sqrt{3}$, IEEE Trans. Inf. Theory 60, 5040 (2014).
  * Schmid [2010] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, J. Fluid Mech. 656, 5 (2010).
  * Tu  _et al._ [2014] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, On dynamic mode decomposition: Theory and applications, J. Comput. Dyn. 1, 391 (2014).
  * Candès [2006] E. J. Candès,  _Proc. Int. Congr. Math. Madrid, August 22–30, 2006_ (European Mathematical Society Publishing House, 2006) pp. 1433–1452.