# Quasidisorder Induced Topology

M. F. Madeira  Departamento de Física and CeFEMA, Instituto Superior Técnico,
Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal  P. D.
Sacramento  Departamento de Física and CeFEMA, Instituto Superior Técnico,
Universidade de Lisboa, Av. Rovisco Pais, 1049-001 Lisboa, Portugal

(August 28, 2024)

###### Abstract

We study the effects of quasidisorder and Anderson disorder on a two
dimensional topological superconductor with an applied external magnetic
field. The cases of a $p$-wave superconductor and a noncentrosymmetric
superconductor with mixed $p$ and $s$-wave pairings and Rashba spin-orbit
coupling are studied. We show that, for a perpendicular magnetic field, the
introduction of quasidisorder leads to the appearance of topological phases in
new regions, characterised by an integer value of the Chern number. For a
parallel magnetic field, we identify regimes with the appearance of new
Majorana flat bands and also new unidirectional Majorana edge states, as
quasidisorder is introduced. We show that the Majorana flat bands have a
quantized Berry phase of $\pi$ and identify it as a topological invariant. Two
topological transitions are identified and the values of the critical
exponents $z$ and $\nu$ are obtained. The fractal nature of the eigenstates is
discussed both for Anderson disorder and Aubry-André disorder.

††preprint: APS/123-QED

##  I Introduction

The search and study of topological properties of matter has proved fruitful
in recent years in research in materials science and condensed matter physics.
Superconductors have long been a focus of interest due to their promising
applications. Superconductors with intrinsic topological properties, in
particular, have recently attracted theoretical and experimental interest due
to phenomena associated with surface or edge Majorana modes, which appear from
an interplay between topology and bulk-boundary correspondence [1, 2, 3].
These Majorana zero modes emerge with non-Abelian exchange statistics and are
sought after due to their promising expected applications in quantum
computing, being candidates for the building blocks of a quantum qubit [4, 5].

It has been theoretically predicted that Majorana states appear as flat
dispersion bands in gapless superconducting phases, such as in the
$d_{xy}$+$p$-wave pairing noncentrosymmetric superconductor in two dimensions
with preserved time-reversal symmetry [6, 7], or for a $p$-wave topological
superconductor in two dimensions, with broken time reversal symmetry by an
applied magnetic field parallel to the two dimensional plane of the system
[8]. Flat bands also emerge on the surface of three dimensional
noncentrosymmetric superconductors, with spin-orbit coupling and which
preserve time-reversal symmetry [9, 10]. It is predicted that flat bands can
increase the critical temperature for superconductivity, and even give rise to
room-temperature superconductivity [11, 12]. Similar behavior has been found
when one has finite-size systems (with increased fluctuations of the density
of states) [13, 14], non-homogeneous order parameters [15, 16, 17, 18, 19,
20], or fractal (critical) states [21, 22, 23, 24, 25, 26, 27, 28, 29] with
corresponding spatial fluctuations of the amplitude of the wave functions. The
difference between an isolated flat band and a flat band with band touchings
has also been recently discussed [30]. It was shown that isolated flat bands
are not needed to achieve a higher superconducting temperature, and that band
touchings can actually increase it. Flat electronic bands can also be found in
some Kagome-type superconductors [31]. A growing interest has been seen in
these types of materials, AV3Sb5 (with A$=$K, Rb, Cs), which can host exotic
quantum properties, displaying topological phases, an unconventional charge
density wave, and evidence of time-reversal symmetry breaking [32, 33, 34,
35].

The study of perturbations in condensed matter systems, namely through the
introduction of disorder, is a central issue. On one hand, introducing
disorder can destroy some phases and their properties, preventing their
experimental observation. In this sense, the study of their robustness becomes
crucial. On the other hand, disorder can by itself lead to new phenomena or
stabilize previously existing phases. One type of disorder that has been
attracting interest in the research field is quasiperiodic disorder. These
systems are somewhat in between periodic and truly random systems, and exhibit
interesting phenomena, in transport [36, 37], topological properties [38, 39,
40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51], and critical behaviour [36,
52, 53, 54]. It is possible to realize these types of systems in experimental
setups of ultracold atoms [55, 56], in optical lattices [57, 58] or in
photonics systems [59]. In addition to systems subject to quasiperiodic
potentials, as in the Aubry-André model [60], there has been growing interest
in Moiré systems in which two incommensurate lattices are connected, or in
which layers of lattices are put in contact and rotated, such as the 2d
twisted bilayer graphene [61, 62, 63, 64, 65]. In such systems, a superlattice
potential is created from proximity coupling between the two lattices, which,
depending on the angle of rotation between the two, may exhibit
quasiperiodicity.

An example of the study of coexisting quasidisorder and superconductivity,
which is significant in the context of this work, is the one dimensional
Kitaev chain with Aubry-André modulation [66, 67, 68, 69]. Without
superconductivity the model has a topological nature revealed by its mapping
to a $2d$ quantum Hall system [70], maintaining a topological nature as we add
superconducting pairing. In general, the mappings involve a corresponding
model in a higher dimension in the form of some parent Hofstadter generalized
Hamiltonian. Topology in quasicrystals may be understood considering mappings
to higher dimensions, typically of the types $1d$ to $2d$ and $2d$ to $4d$. In
$1d$ with no superconductivity the model is self-dual (position and momentum
space) and there is a single transition from an extended state phase to a
phase where all the states are localized. At the transition point the system
has critical states. Generalized models show the existence of mobility edges,
such that there is a separation as a function of energy between extended and
localized states [71, 72, 73, 74, 75, 76, 77, 78] and the existence of hidden
dualities leads to a rich class of systems, where such edges appear [79]. The
introduction of $p$-wave pairing in the Aubry-André model leads to the
appearance of a finite extent region of critical (fractal) states, between the
regions of extended and localized states. Remarkably, the transitions between
localised and critical regimes have been studied and were found to deviate
from the known Aubry-André universality class [68, 69].

In this work we study a model of a two-dimensional superconductor with spin
triplet $p$-wave pairing, or mixed $p$ and $s$-wave pairings with Rashba spin-
orbit coupling, in the presence of a time reversal symmetry breaking magnetic
field. Some materials which are candidates for realizing triplet pairing
superconductivity include Sr2RuO4 [80], UPt3 [81] and CuxBi2Se3 [82]. In the
presence of $s$-wave pairing and Rashba spin orbit coupling, the model
describes a noncentrosymmetric superconductor, of which are examples
$\text{CePt}_{3}\text{Si}$ [83], $\text{CeIrSi}_{3}$ [84] and
$\text{CeRhSi}_{3}$ [85]. In the noncentrosymmetric regime the breaking of
inversion symmetry allows for the mixture of spin-triplet and spin-singlet
pairings. This mixing is expected to lead to novel phenomena such as higher
than usual values of the upper critical field [86, 87].

The clean model has been studied, in both the centrosymmetric and the
noncentrosymmetric regimes, and is known to possess diverse topological
properties. If time-reversal symmetry is preserved, the model displays gapless
Majorana edge states and is characterised by a $\mathbf{Z}_{2}$ invariant. The
observed properties when time-reversal symmetry is broken by an external
magnetic field are found to be very dependent on its direction in relation to
the two-dimensional superconducting plane. If the magnetic field is such that
it is perpendicular to the plane of the superconductor, the model has a rich
phase diagram indexed by the Chern number [88]. When the magnetic field is
parallel to the plane of the system, interesting phenomena, such as Majorana
flat bands or Majorana unidirectional states, appear on phases with a gapless
bulk [8, 89].

The effect of disorder may be considered in different ways. One possibility is
to consider a non-homogeneous magnetic field, achieved by inserting magnetic
impurities in the clean superconductor [90] which may give rise to or change
topological properties in the system. Examples include the addition of chains
of magnetic adatoms [91, 92], islands of magnetic impurities [93] or fully
random distributions of impurities [94]. Another possibility is to consider
potential scattering impurities on the superconductor in the presence of a
constant magnetic field, either perpendicular or parallel to the system. We
are interested in studying the effects of quasidisorder in these regimes.
Besides Aubry-André disorder, we will also consider Anderson disorder as a
comparison to the effects of quasi-periodicity.

Anderson localization does not require full randomness. If differences in
potential between sites are large enough compared to hoppings, one may expect
a transition to localized states. In addition to full randomness, a
quasidisordered potential leads to localization if the disorder amplitude is
large enough [95]. One expects that the Aubry-André quasiperiodic potential
should affect the long-range nature of states, and in particular topological
states that are by themselves of long-range nature. Aubry-André is expected to
be naturally of a multifractal nature. Anderson and Aubry-André are different
and, in particular, critical states due to Anderson appear at the transition
to localization while in the Aubry-André added to the Kitaev $1d$ model one
finds phases with this behavior (or in $2d$ a mixture of critical states in
the crossover to localization). Multifractality probes long distances and
therefore one expects that it may enhance superconductivity due to Chalker
scaling [96, 97], as expected and observed with other inhomogeneities.
Multifractal wave functions have larger spatial overlap and stronger state to
state correlations for states with similar energies.

As stated previously, quasiperiodicity may also lead to topological properties
[98, 99]. A $2d$ topological insulator plus quasiperiodic potential shows a
transition from a trivial insulator to a topological insulator. Flat
topological bands and eigenstate criticality have also been shown as a result
of a quasiperiodic perturbation in the context of the Bernevig-Hughes-Zhang
model plus $2d$ quasiperiodic potential [100].

The presence of gapless states in a system may also be associated with long-
distance behavior. For instance, nodal points of Weyl semimetals may lead to
interesting behavior in the presence of disorder. It has been shown that they
survive the presence of moderate disorder [101]. On the other hand, in the
case of gapless states of the form of nodal loops, any amount of disorder
mixes states. Disorder-driven multifractality has been shown in Weyl nodal
loops [102]. In the case of magic angle semimetals quasiperiodicity
generically leads to flat bands in nodal, semi-metallic structrures. A
transition from a Weyl semimetal to metal driven by quasiperiodic potential
has been found in $3d$ [52, 103].

While the influence of disorder, either Anderson or quasidisorder has been
extensively considered in the case of one-dimensional systems, including in
the presence of superconductivity, it is interesting to consider their effects
on a two-dimensional $p$-wave superconductor, and in particular in the
presence of a magnetic field. In the clean system the topology is influenced
by the orientation of the magnetic field and, in particular, the gapped or
gapless nature of the states may lead to different responses to disorder. As
stated before, the difference of symmetry classes plus disorder gives rise to
new universality classes. Also, topology may be induced by quasiperiodicity,
which leads to the expectation of new universality classes (beyond the usual
classification), as found in the one-dimensional case. In particular, one may
expect interesting effects with the interplay of quasiperiodicity due to the
presence of critical bulk states, and the existence of Majorana flat bands.
The long-range nature of the quasiperiodic potential and the intrinsic long-
range nature of the gapless states may lead to an interesting competition. A
distinction between Anderson disorder (with moderate intensity) and
quasidisorder is therefore interesting to consider, as shown in non-
superconducting systems, where for instance nodal points and nodal loops are
affected differently by Anderson disorder, or on a semimetal where imposing a
quasiperiodic potential leads to flat bands.

The rest of the paper is organized as follows. Section II introduces the model
of the Hamiltonian and the topological properties of the clean system are
discussed, first under a perpendicular and second under a parallel magnetic
field, respectively in subsections II.1 and II.2. In subsection II.2 we derive
the regions where the model is topological, and show that the topological
regions are characterized by a Berry phase of $\pi$. In section III we present
the results for the disordered model under a perpendicular magnetic field. We
show that the introduction of Aubry-André disorder leads to the appearance of
topological phases in new regions. In IV we present the results for the
disordered model under a parallel magnetic field. First we discuss the
localization properties of the system in real space under different types of
disorder, using the inverse participation ratio (IPR). We then turn to a mixed
space description and discuss the evolution of the system as Anderson or
Aubry-André disorder are introduced. We show that the introduction of Aubry-
André disorder leads to the appearance of new regimes: for the $p$-wave
superconductor, new gapless regimes with Majorana flat bands appear, and for
the noncentrosymmetric superconductor, new regimes with unidirectional edge
states appear. We then obtain the Berry phase using twisted boundary
conditions and show it is quantized to a value of $\pi$ for the quasidisorder
induced flat bands. Identifying it as a topological invariant, we study two
topological transitions and obtain the critical exponents $z$ and $\nu$, which
we find to deviate from the known universality classes. Finally, using the IPR
we study the nature of the eigenfunctions distinguishing between localized,
single-fractal and multifractal regimes in the thermodynamic limit for both
Anderson and Aubry-André disorder. We conclude in section V. Three appendices
discuss some further results on the disorder driven transitions under a
perpendicular magnetic field in Appendix A, the influence of the
dimensionality of the quasidisorder potential in Appendix B and the energy
spectra and density of states for the disordered noncentrosymmetric
superconductor in Appendix C.

##  II Model Hamiltonian

In momentum space, the Bogoliubov-de Gennes (BdG) Hamiltonian matrix of the
two dimensional model is written as

$\mathcal{H}(\mathbf{k})=\left(\begin{array}[]{cc}\xi(\mathbf{k})+\mathbf{B}\cdot\bm{\sigma}&\Delta(\mathbf{k})\\\
\Delta^{\dagger}(\mathbf{k})&-\xi^{T}(-\mathbf{k})-\mathbf{B}\cdot\bm{\sigma}^{*}\end{array}\right)$
(1)

in a basis
$(\bm{c}_{\mathbf{k}}^{\dagger},\bm{c}_{-\mathbf{k}})=(c_{\mathbf{k}\uparrow}^{\dagger},c_{\mathbf{k}\downarrow}^{\dagger},c_{-\mathbf{k}\uparrow},c_{-\mathbf{k}\downarrow})$
with $c_{\mathbf{k}\sigma}^{\dagger}$ ($c_{\mathbf{k}\sigma}$) the creation
(annihilation) operator for an electron with momentum $\bm{k}=(k_{x},k_{y})$
and spin projection $\sigma$. In the BdG Hamiltonian,
$\xi(\bm{k})=\epsilon_{\bm{k}}\sigma_{0}+\mathbf{s}\cdot\bm{\sigma}$, where
$\epsilon_{\mathbf{k}}=\left[-2t\left(\cos k_{x}+\cos
k_{y}\right)-\mu\right]\sigma_{0}$ is the kinetic term, with $t$ the nearest-
neighbour hopping integral and $\mu$ the chemical potential,
$\mathbf{s}\cdot\bm{\sigma}=-\alpha(-\sin k_{y},\sin
k_{x},0)\cdot\bm{\sigma}=-\alpha\left[-\sin k_{y}\sigma_{x}+\sin
k_{x}\sigma_{y}\right]$ is the Rashba spin-orbit term with $\mathbf{s}$ the
spin-orbit vector. The term $\mathbf{B}\cdot\bm{\sigma}$ describes the Zeeman
coupling of the electrons with an external magnetic field $\mathbf{B}$ and
$\hat{\Delta}(\mathbf{k})=\left[\Delta_{s}+\mathbf{d}(\mathbf{k})\cdot\bm{\sigma}\right]\left(i\sigma_{y}\right)$
is the superconducting gap function. The pairing vector is chosen as
$\mathbf{d}=d(-\sin k_{y},\sin k_{x},0)$, so that $d$ is the $p$-wave pairing
amplitude and $\Delta_{s}$ is the $s$-wave pairing amplitude. The simultaneous
existence of $s$ and $p$-wave terms is possible with a nonzero spin-orbit
term, which breaks the parity symmetry.

The case of study is that of a system with periodic boundary conditions along
the $x$ direction and open boundary conditions in the $y$ direction, such as
in a cylinder geometry. Thus we can also write the Hamiltonian in a mixed
space, $(k_{x},y)$, where a Fourier transform to the reciprocal space is only
done in the $x$ direction. In this case, for each value of $k_{x}$ the
Hamiltonian matrix has a dimension $(4\times N_{y})\times(4\times N_{y})$,
where $N_{y}$ is the number of sites in $y$. It is also of interest to write
the Hamiltonian in real space. In this case the Hamiltonian matrix has
dimension $(4\times N)\times(4\times N)$ with $N=N_{x}\times N_{y}$ the total
number of sites and $N_{x}$, $N_{y}$ the number of sites in the $x$ and $y$
directions, respectively.

When $\mathbf{B}=0$, the system respects the time-reversal symmetry (TRS)
$\mathcal{T}=(\sigma_{0}\otimes\mathrm{i}\sigma_{y})$ and the particle-hole
symmetry (PHS) $\mathcal{P}=(\sigma_{x}\otimes\sigma_{o})$ such that

$\begin{split}&\mathcal{P}\mathcal{H}(\mathbf{k})\mathcal{P^{\dagger}}=-\mathcal{H}^{*}(-\mathbf{k}),\\\
&\mathcal{T}\mathcal{H}(\mathbf{k})\mathcal{T}^{\dagger}=\mathcal{H}^{*}(-\mathbf{k}),\end{split}$
(2)

and $\mathcal{T}^{2}=-1$, $\mathcal{P}^{2}=1$. Therefore the Hamiltonian
belongs to the DIII symmetry class, and if $|d|>|\Delta_{s}|$ the system has a
nontrivial $\mathbf{Z}_{2}$ number, displaying gapless counterpropagating
Majorana edge states [88, 8].

For $\mathbf{B}\neq 0$ the time-reversal symmetry is broken. The system
exhibits different topological properties whether the applied magnetic field
is perpendicular or parallel to the plane of the system, as will be now
discussed.

###  II.1 Perpendicular Magnetic Field

Let us first consider the case in which the external magnetic field is
perpendicular to the plane of the system, $\mathbf{B}=(0,0,B_{z})$. We have a
gap closing point if one of the equations is satisfied [88]:

$\begin{split}(-4t-\mu)^{2}+\Delta^{2}_{s}=B^{2}_{z},\\\
\mu^{2}+\Delta^{2}_{s}=B^{2}_{z},\\\
(4t-\mu)^{2}+\Delta^{2}_{s}=B^{2}_{z}.\end{split}$ (3)

Eqs. 3 define the boundaries between regions in which the system has different
topological properties. At the gap closing points the $D$ class system with
broken time reversal symmetry undergoes topological transitions between gapped
phases with different Chern numbers. The phase diagram of the system (indexed
by the Chern number) is presented in Fig. 1(a) for $t=1$, $\Delta_{s}=0$, and
$d>0$.

Figure 1: Phase diagram for a) Chern number and b) winding number
$I(k_{y}=0,\pi)$ as a function of $\mu$ and $B_{z}$, for $t=1$,
$\Delta_{s}=0$, $d>0$.

The regimes with a Chern number of zero and $B_{z}<2$, $0<|\mu|<4t$ exhibit
edge states, besides having $C=0$. This can be explained by one additional
topological invariant. It can be defined noting that the Hamiltonian obeys a
particle hole symmetry $\mathcal{P}=(\sigma_{x}\otimes\sigma_{0})$ with

$\mathcal{P}\mathcal{H}(\mathbf{k})\mathcal{P}^{\dagger}=-\mathcal{H}^{*}(-\mathbf{k}).$
(4)

For the values $k_{y}=0$ and $k_{y}=\pi$, the Hamiltonian obeys
$\mathcal{H}^{*}(-\mathbf{k})=\mathcal{H}(\mathbf{k})$ and thus anticommutes
with $\mathcal{P}$, $\\{\mathcal{H}(\mathbf{k}),\mathcal{P}\\}=0$. Therefore
the basis which diagonalizes $\mathcal{P}$ anti-diagonalizes the Hamiltonian.
A winding number $I(k_{y})$ can then be defined as [88]

$\begin{split}I\left(k_{y}\right)=\frac{1}{4\pi
i}\int_{-\pi}^{\pi}dk_{x}\operatorname{tr}\big{[}q^{-1}(k_{x})\partial_{k_{x}}q(k_{x})-\\\
q^{\dagger-1}(k_{x})\partial_{k_{x}}q^{\dagger}(k_{x})\big{]},\quad
k_{y}=0,\pi,\end{split}$ (5)

with

$q(k_{x})=\left(\begin{array}[]{cc}-\epsilon_{\mathbf{k}}-B_{z}+id\sin
k_{x}&\Delta_{s}-i\alpha\sin k_{x}\\\ -\Delta_{s}+i\alpha\sin
k_{x}&-\epsilon_{\mathbf{k}}+B_{z}+id\sin k_{x}\end{array}\right)$ (6)

the anti-diagonal block of the Hamiltonian matrix. The values of $I(0)$ and
$I(\pi)$ inside each phase are represented in Fig. 1(b).

The invariant $I(k_{y})$ loses its meaning if a finite magnetic field in the
$y$ direction, $B_{y}$, is applied. However, we found that this is not true
for the Chern number. Fig. 2 shows phase diagrams indexed by the Chern number
as a function of $B_{z}$ and $B_{y}$ for three different values of $\mu$. In
this case the Chern number depends only on the value of
$\sqrt{B^{2}_{y}+B^{2}_{z}}$. Also note that the diagrams only concern values
of $B_{z}>0$, excluding the points where $B_{z}=0$ and $B_{y}\neq 0$. In Fig.
3 we present the phase diagram of the system as a function of $\mu$ and
$B_{z}$ for constant values of $B_{y}$.

Figure 2: Phase diagrams for ($B_{z}>0,B_{y}$), indexed by the Chern number,
obtained numerically for a) $\mu=0$, b) $\mu=1$ and c) $\mu=-3.5$ for
$\Delta_{s}=0$. Figure 3:  Phase diagram indexed by the Chern number as a
function of $\mu$ and $B_{z}$ (with $B_{z}>0$), for $t=1$, $\Delta_{s}=0$,
$d>0$, and a) $B_{y}=2$ and b) $B_{y}=4.5$.

###  II.2 Parallel Magnetic Field

Now let us consider the case in which the applied magnetic field is parallel
to the system, $\mathbf{B}=(B_{x},B_{y},0)$. This could be realized, for
instance, by threading a wire through the center of the superconductor in a
cylindrical geometry.

Taking first the $s$-wave term $\Delta_{s}$ and the spin-orbit term $\alpha$
to be zero, the eigenvalues of the Hamiltonian are given by

$E(\mathbf{k})=\pm\sqrt{z_{1}\pm 2\sqrt{z_{2}}},$ (7)

with

$\begin{split}z_{1}=\mathbf{d}\cdot\mathbf{d}+\epsilon^{2}_{\mathbf{k}}+\mathbf{B}\cdot\mathbf{B},\\\
z_{2}=\epsilon^{2}_{\mathbf{k}}(\mathbf{B}\cdot\mathbf{B})+(\mathbf{B}\cdot\mathbf{d})^{2}.\end{split}$
(8)

The gap closing points are solutions of the equation $z_{1}=2\sqrt{z_{2}}$,
which is equivalent to the two equations being simultaneously satisfied:

$\begin{split}&\mathbf{d}\cdot\mathbf{d}+\epsilon^{2}_{\mathbf{k}}=\mathbf{B}\cdot\mathbf{B},\\\
&(\mathbf{B}\cdot\mathbf{B})(\mathbf{d}\cdot\mathbf{d})=(\mathbf{B}\cdot\mathbf{d})^{2}.\end{split}$
(9)

Eqs. 9 simplify if we consider the magnetic field aligned with one of the
axes. Let us then take the magnetic field aligned with the $y$ direction,
$\mathbf{B}=(0,B_{y},0)$. In this case, the second equation simplifies to
$\sin{k_{y}}=0$ which implies the bulk gap will close at
$k_{y,0}=n\pi,n\in\mathbb{Z}$, provided there are values of $k_{x}$ that
satisfy the equations

$d^{2}\sin^{2}{k_{x}}+(-2t(\cos{k_{x}}\pm 1)-\mu)^{2}=B^{2}_{y}.$ (10)

When the $p$-wave superconductor is in a gapless phase, and for a certain
range of magnetic field, Majorana flat bands (MFBs) will appear in the system.
This will be discussed next.

When finite spin-orbit $\alpha$ and $s$-wave pairing $\Delta_{s}$ terms are
also considered, the flat bands will (for certain values of the magnetic
field) acquire a slope, giving origin to unidirectional Majorana edge states
(MESs). The appearance of such states is only possible with a gapless bulk,
where a counter-propagating bulk current is created to cancel the edge current
[8].

####  II.2.1 Flat bands: winding number and Berry phase quantization

When the system is subject to an applied magnetic field, it no longer respects
time-reversal symmetry. If the applied field has a generic form
$\mathbf{B}=(B_{x},B_{y},0)$ we can, however, take $k_{x}$ as a fixed
parameter of the Hamiltonian and find a set of symmetries that are only
satisfied in the $y$ direction. It is found that the Hamiltonian respects the
symmetries:

$\begin{split}\mathcal{T}_{k_{y}}^{-1}\mathcal{H}(k_{x},k_{y})\mathcal{T}_{k_{y}}=\mathcal{H}(k_{x},-k_{y}),\\\
\mathcal{P}_{k_{y}}^{-1}\mathcal{H}(k_{x},k_{y})\mathcal{P}_{k_{y}}=-\mathcal{H}(k_{x},-k_{y}),\\\
\end{split}$ (11)

where $\mathcal{T}_{k_{y}}=(\sigma_{z}\otimes\sigma_{z})K$ and
$\mathcal{P}_{k_{y}}=(\sigma_{y}\otimes\sigma_{y})K$ are, respectively,
defined as a ”time-reversal-like” symmetry and a ”particle-hole-like” symmetry
[8] with $\mathcal{T}_{k_{y}}^{2}=\mathcal{P}_{k_{y}}^{2}=1$ ($K$ is the
complex conjugate operator). From these we can define a third chiral-like
symmetry $\mathcal{S}_{k_{y}}=\mathcal{T}_{k_{y}}\mathcal{P}_{k_{y}}$:

$\begin{split}\mathcal{S}_{k_{y}}^{-1}\mathcal{H}(k_{x},k_{y})\mathcal{S}_{k_{y}}=-\mathcal{H}(k_{x},k_{y}).\\\
\end{split}$ (12)

Since we have that $\mathcal{T}_{k_{y}}^{2}=\mathcal{P}_{k_{y}}^{2}=1$, the
Hamiltonian belongs to the BDI symmetry class and, since the problem is
effectively reduced to one dimension, the system can be characterized by an
integer topological invariant. We can then write the Hamiltonian in the basis
where $\mathcal{S}_{k_{y}}$ is diagonal, in which the Hamiltonian takes an
anti-diagonal form. From here it is possible to obtain a winding number
$\mathcal{W}$ at each value of $k_{x}$. It can be shown [8] that the winding
number is calculated as

$\mathcal{W}(k_{x})=\frac{i}{\pi}\left[\log{\left(\frac{sgn(\mathcal{M}(k_{y}=0))}{sgn(\mathcal{M}(k_{y}=\pi))}\right)}\right]$
(13)

with

$\begin{split}&\mathcal{M}\left(k_{x},k_{y}\right)=\\\ &\left[\mu+2t\left(\cos
k_{x}+\cos
k_{y}\right)\right]^{2}+d^{2}\sin^{2}k_{x}-B_{y}^{2}+B_{x}^{2}.\end{split}$
(14)

In the regimes with $|\mathcal{W}|=1$ the system has a topological nature and
Majorana flat bands appear, as is shown in Fig. 4. These are protected by the
chiral symmetry $\mathcal{S}_{k_{y}}$ as defined in Eq. 12.

Figure 4: Energy spectrum, absolute value of the winding number $\mathcal{W}$
and Berry phase $\gamma$ normalized by $2\pi$, as a function of $k_{x}/\pi$.
The values of the parameters are $t=1$, $d=1/6$, $\mu=-3.5$ and a) $B_{y}=d$,
b) $B_{y}=3.5d$.

The existence of topological flat bands may also be identified by a non-
trivial Berry phase. In general, the Berry phase can take any real value. In
the presence of certain symmetry constraints, the Berry phase can become
quantized to 0 or $\pi$ and carry topological information (at the value of
$\pi$). This quantization can happen in the presence of inversion or chiral
symmetries, also leading to the quantization of polarization [104]. As the
problem is reduced to one dimension, we can obtain a Berry phase $\gamma_{B}$
at each value of $k_{x}$, given by:

$\gamma_{B}(k_{x})=i\int_{0}^{2\pi}dk_{y}\langle\Psi(k_{x},k_{y})\mid\frac{\partial}{\partial
k_{y}}\Psi(k_{x},k_{y})\rangle$ (15)

with $\Psi$ the ground-state wavefunction. The calculation is done numerically
by discretizing the Brillouin zone [104, 105, 106, 107] in the $y$ direction.
As is shown in Fig. 4, we have found that in the regimes with
$|\mathcal{W}|=1$, the Berry phase is also quantized to a value of $\pi$.

####  II.2.2 Domain of flat band existence: topological and gapless regions

From Eq. 13 it is found that $|\mathcal{W}|=1$ in the regimes where
$\mathcal{M}\left(k_{x},k_{y}=0\right)$ and
$\mathcal{M}\left(k_{x},k_{y}=\pi\right)$ have opposite signs. This is only
possible if $|B_{y}|>|B_{x}|$, thus this is a necessary condition for the
appearance of MFBs. The flat band regions can be summarized in (with
$\tilde{B}^{2}=B_{y}^{2}-B_{x}^{2}$):

  * •

$(1)$ $\mu\geq 2t$

$\mathcal{D}_{+}>\tilde{B}^{2}>\mathcal{D}_{-}$ (16)

  * •

$(2)$ $\mu\leq-2t$

$\mathcal{D}_{-}>\tilde{B}^{2}>\mathcal{D}_{+}$ (17)

  * •

$(3)$ $-2t<\mu<2t$

$(\mathcal{D}_{+}>\tilde{B}^{2}>\mathcal{D}_{-})\vee(\mathcal{D}_{-}>\tilde{B}^{2}>\mathcal{D}_{+})$
(18)

where

$\mathcal{D}_{\pm}=\left[\mu+2t\left(\cos k_{x}\pm
1\right)\right]^{2}+d^{2}\sin^{2}k_{x}.$ (19)

Eqs. 16, 17 and 18 define the regions where the superconductor is in a
nontrivial regime with $|\mathcal{W}|=1$, for a certain value of $k_{x}$.
Furthermore, since MFBs can only appear in a gapless phase, the equations also
define the regions where the bulk is gapless, as a function of the in-plane
magnetic field.

Figure 5: a) Domain of existence of Majorana flat bands (shaded region) for
$B_{y}$ vs. $k_{x}$ for the parameters $t=1$, $d=1/6$, $\mu=-3.5$. b) Closeup
of a) in the region $B_{y}\in$ [$-1,1$] and $k_{x}\in$ [$-1.5,1.5$].

Note that the chiral-like symmetry that protects the flat bands is broken by
either a non-zero $s$-wave pairing term $\Delta_{s}$ or a non-zero spin-orbit
term $\alpha$. A finite perpendicular magnetic field $B_{z}$ is also found to
break the chiral-like symmetry, leading to the absence of flat bands. If the
flat band includes the point $k_{x}=0$, the addition of a finite $B_{z}$ will
lead to the appearance of bands with a finite slope that cross at zero energy
at $k_{x}=0$. Otherwise, the bands will be lifted to finite energy.

Figure 6: Phase diagrams indexed by the Chern number $C$ for a system with
20x20 sites, for several values of disorder strength $\lambda$ and
perpendicular magnetic field $B_{z}$, obtained for an average over 10 disorder
configurations. For Aubry-André disorder, each random disorder configuration
is obtained by selecting a random value of $\phi$. The first row with panels
a)-c) concerns the case of Anderson disorder (2d), the second row with panels
d)-f) concerns the case of Anderson disorder (1d along $y$, uniform along
$x$), the third row with panels g)-i) concerns Aubry-André disorder (1d along
$y$, uniform along $x$), and the fourth row with panels j)-l) concerns Aubry-
André disorder (2d). The values of the parameters are $t=1$, $d=0.6$ and
$\mu=0$ (left), $d=0.6$ and $\mu=1$ (middle), $d=1/6$ and $\mu=3d-4t=-3.5$
(right).

##  III Disordered model under a perpendicular magnetic field

We first want to investigate the effects of quasidisorder and disorder on the
system subject to an applied magnetic field in the perpendicular direction,
$\mathbf{B}=(0,0,B_{z})$. Here we limit ourselves to the study of the system
in real space and, to classify the topological nature of the system, the Chern
number is obtained numerically [108]. We consider four different types of
disorder potentials:

  1. 1.

Anderson disorder (2d), where the disorder term is random at each site and
varies with uniform probability within an interval:

$\Lambda(x,y)\in[-\lambda,\lambda].$ (20)

  2. 2.

Anderson disorder (1d along $y$, uniform along $x$), where the potential is of
the same type as described above but varies only along the $y$ direction,
being uniform along the $x$ direction:

$\Lambda(x,y)=\Lambda(y)\in[-\lambda,\lambda].$ (21)

  3. 3.

Aubry-André disorder (1d along $y$, uniform along $x$), where the disorder
term is a quasiperiodic potential of the form:

$\Lambda(x,y)=\Lambda(y)=\lambda\cos(2\pi\beta f(x,y)+\phi)$ (22)

with $f(x,y)$ a function of the lattice sites, $\beta=\frac{\sqrt{5}-1}{2}$
the inverse golden ratio, and $\phi$ a phase between $0$ and $2\pi$. Here we
take $f(x,y)=y$, so that the considered quasiperiodic potential is uniform in
the $x$ direction.

  4. 4.

Aubry-André disorder (2d), where the disorder term is a sum of two
quasiperiodic potentials of the form:

$\Lambda(x,y)=\lambda\cos(2\pi\beta x+\phi)+\lambda\cos(2\pi\beta y+\phi)$
(23)

so that disorder potentials are introduced in both the $x$ and $y$ directions.

In Fig. 6 we show the phase diagrams indexed by the Chern number, for three
different values of $\mu$ and $d$ (with $t=1$ in all cases) and for a system
with size $20\times 20$.

When Anderson disorder is introduced in the system (first row), the
topological regimes are destroyed as the disorder strength is increased. There
is, however, some difference in robustness as a function of the magnetic
field. This is noticeable in Figs. 6(a) and (b), where we see that the
robustness of the topological phases increases with the increase of $B_{z}$.
In Fig. 6(a) and for a small region of magnetic field (for $B_{z}>4$) we
observe reentrant topology as disorder is increased, as in Fig. 6(b), for
lower values of magnetic field ($B_{z}<1$).

The second row of the figure is obtained when disorder is considered with
uniformity in the $x$ direction. Unexpectedly, the topological regions are to
be less robust if compared with the previous case where Anderson disorder was
considered with no modulation. Small traces of induced topology are observed
for $B_{z}<1$ in Fig. 6(e) and $B_{z}<0.5$ in Fig. 6(f).

For Aubry-André disorder uniform in the $x$ direction (third row) we obtain
phase diagrams with well defined boundaries, and with induced topological
regions. Here, the topological phases show an interesting and unexpected
response to the increase of quasidisorder. There is a clear difference in
robustness for different values of $B_{z}$, which originates the seemingly
effect of “peaks” and “valleys” in the phase diagram, respectively at more
robust and more vulnerable values of $B_{z}$. Induced topology is visible in
panels g)-i), with topological transitions to finite values of $C$ happening
at low and high values of the magnetic field with the increase of disorder.

The last row of Fig. 6 concerns the case of two-dimensional Aubry-André
disorder. The introduction of disorder leads to the appearance of new
topological regions, where several are characterized by values of $C$ that are
not seen in the clean system, in the range of $[-4,4]$. Also, some regions
appear where the Chern number oscillates within an interval between two
integer values, without tending clearly to one of them.

We may argue that by adding disorder, local fluctuations of $\mu$ may lead to
changes of the Chern number. This is particularly seen in the presence of
quasidisorder. This suggests that the long-range quasiperiodicity resonates
more with the calculation of the Chern number, that reflects the global
structure of the states. Further details on the effect of disorder are shown
in Appendix A.

Figure 7: Average IPR of the whole system as a function of disorder strength,
for a) Anderson disorder (2d), b) Anderson disorder (1d along $y$, uniform
along $x$), c) Aubry-André disorder (1d along $y$, uniform along $x$), d)
Aubry-André disorder (2d). The IPR is averaged over all eigenstates of the
system in a given disorder configuration, and averaged over 10 disorder
configurations. In a) fits are done to functions of the form
$IPR=C_{1}\exp{C_{2}\lambda}$ in the range $\lambda\in[1.5,3]$, giving the
values $(C_{1},C_{2})=(5.6\times 10^{-4},1.13)$ for $B_{y}=0.5d$ and
$B_{y}=3.5d$, and $(C_{1},C_{2})=(3.6\times 10^{-4},1.22)$ for $B_{y}=4d$,
$\alpha=0.2d$, $\Delta_{s}=0.5d$. In b) fits are done to functions of the form
$IPR=C_{1}\lambda+C_{2}$ in the range $\lambda\in[0.5,3]$, giving the values
$(C_{1},C_{2})=(4.4\times 10^{-3},5\times 10^{-5})$ for $B_{y}=0.5d$ and
$B_{y}=3.5d$, and $(C_{1},C_{2})=(2.8\times 10^{-3},1.4\times 10^{-4})$ for
$B_{y}=4d$, $\alpha=0.2d$, $\Delta_{s}=0.5d$.

##  IV Disordered model under a parallel magnetic field

We now introduce disorder on the system with an applied parallel magnetic
field in the $y$ direction, $\mathbf{B}=(0,B_{y},0)$. In section A we consider
the system in real space, where the disorder term takes the same form as in
the previous section (Eqs. 20, 21, 22,23).

From section B onwards we study the system in a mixed $(k_{x},y)$ space. In
this case the disorder term is either of the form of Eq. 21 or Eq. 22, with
the potential varying in the $y$ direction and being the same for all $k_{x}$,
as:

  1. 1.

Anderson disorder:

$\Lambda(y)\in[-\lambda,\lambda],$ (24)

  2. 2.

Aubry-André disorder:

$\Lambda(y)=\lambda\cos(2\pi\beta y+\phi),$ (25)

where, as before, $\beta=\frac{\sqrt{5}-1}{2}$ is the inverse golden ratio and
$\phi$ is a phase between $0$ and $2\pi$.

###  IV.1 Localization properties in real space

Here we briefly consider the system in real space and study its localization
properties in three different regimes. We fix the parameter values as $t=1$,
$d=1/6$, $\mu=-3.5$ and consider three different cases: the case of a $p$-wave
superconductor for which a magnetic field $B_{y}=0.5d$ is added, such that the
system is in a phase with a gapped bulk but gapless edge states; a $p$-wave
superconductor with an added magnetic field of $B_{y}=3.5d$, where the system
has a gapless bulk and is in the MFB regime; and a case of the
noncentrosymmetric superconductor, with $B_{y}=4d$ and added $s$-wave pairing
and spin-orbit terms, $\Delta_{s}=0.3d$ and $\alpha=0.2d$, where the system
has a gapless bulk and unidirectional MESs.

To quantify the effects of disorder on the system’s localization we use the
inverse participation ratio, IPR. For a given eigenstate labeled by $m$, the
IPR is defined as:

$\text{IPR}_{m}=\sum_{i}\left|\psi_{i}^{m}\right|^{4},$ (26)

with $\psi_{i}^{m}$ the wavefunction of the eigenstate $m$ at a site $i$. For
perfectly localized states we have that $\text{IPR}_{m}\sim 1$ and for
delocalized states $\text{IPR}_{m}\sim 1/N$. In Fig. 7 we present results for
the average IPR as a function of disorder for a system of size $N=N_{x}\times
N_{y}=41\times 41$ [109] and for the same types of disorder as before: Fig.
7(a) Anderson disorder (2d), Fig. 7(b) Anderson disorder (1d along $y$,
uniform along $x$), Fig. 7(c) Aubry-André disorder (1d along $y$, uniform
along $x$), and Fig. 7(d) Aubry-André disorder (2d).

From observation of Figs. 7(a)-(d) we find four qualitatively different
behaviours. In Fig. 7(a) (Anderson disorder) we see that the IPR shows an
exponential-like behaviour for $\lambda>1.5$. A fit of the form
$\text{IPR}=C_{1}\exp{C_{2}\lambda}$ is done in the range $\lambda\in[1.5,3]$,
and is presented in Fig. 7(a) in dashed lines. We find that for $\lambda>1.5$
the IPR follows an exponential behaviour closely, while for $\lambda<1.5$
there is a deviation from it. As disorder is increased, the low energy states
become increasingly localized inside the bulk. From inspection of the
wavefunctions we observed that the edge states quickly lose their structure
for low values of disorder, although they do not become as quickly localized
as the remaining bulk states. Accordingly, the IPR of these low energy states
shows a slower increase than what is observed in Fig. 7(a).

Figure 8: Energy spectra evolution with a) Anderson disorder and b) Aubry-
André disorder, for $B_{y}=0.5d$.

For Anderson disorder along $y$ and with $x$ uniformity, we find a different
localization behaviour. The IPR grows linearly with the increase of disorder,
although with some fluctuations and a deviation for $\lambda<0.5$. A fit of
the form $\text{IPR}=C_{1}\lambda+C_{2}$ is done to the range
$\lambda\in[0.5,3]$ and presented in red dashed lines. The change in behaviour
in relation to 7(a) is a result of imposing periodicity in the $x$ direction
on the disorder term, which unexpectedly causes the behaviour of the IPR to
become linear. The IPR of the low energy states follows a similar behaviour to
what is seen for the average IPR. By increasing disorder the low energy edge
states are removed from the edges and localize inside the bulk, while
remaining periodic in the $x$ direction.

In Fig. 7(c) (Aubry-André disorder with $x$ uniformity) we see a threshold
behaviour where a transition happens around $\lambda=2$. For $\lambda<2$
(approximately) there is a slow increase of the IPR, while for $\lambda>2$ the
IPR greatly increases. This resembles some known results: in the one-
dimensional Aubry-André model, where the system undergoes an extended-
localized transition at $\lambda=2t$, after which the average IPR shows a
marked increase; for a one dimensional $p$-wave superconductor with an Aubry-
André potential this transition point changes to $\lambda$ = 2(t+d) with $d$
the $p$-wave pairing amplitude (when the chemical potential is taken as zero)
[41, 42]. For values of $\lambda$ before the transition, we observed that some
bulk states acquire a critical like behaviour in the $y$ direction, while
remaining periodic in the $x$ direction. The low energy states are more robust
to disorder if compared with the Anderson disorder cases, and are only removed
from the system at the transition: after the threshold value of $\lambda$
there are no edge states in the system.

Fig. 7(d) concerns the case of two-dimensional Aubry-André disorder. The IPR
shows again a threshold behaviour, and as in c) a transition is seen slightly
below $\lambda=2$. However, the transition between two different regimes is
abrupt in the IPR, and more closely resembles that of the one dimensional
Aubry-André chain. Also, unlike cases a)-c), the IPR follows the same
behaviour for the three regimes considered. By a closer inspection of the IPR
we see that this is only true for values of disorder over $\lambda=0.7$, as
for $\lambda<0.7$ the noncentrosymmetric regime shows a consistently lower
IPR, as before. As disorder is increased for $\lambda>2$, states localize
along both the $x$ and $y$ directions.

The comparison between $1d$ Aubry-André disorder and $2d$ Aubry-André
disorder, as well as a comparison between a perpendicular and a parallel
magnetic field, is detailed in Appendix B, with particular emphasis on the
existence of critical states and the apparent absence of a transition between
extended and critical states, in contrast to what is found in the one-
dimensional case.

Although results are not explicitly shown, the effect of edge disorder was
also briefly studied, extending previous results obtained for a time-reversal
invariant system [110]. We considered both Anderson and Aubry-André disorder
potentials which were introduced locally at the edges at $y=0$ and $y=N_{y}$,
varying along the $x$ direction (along the edge) also for a system of size
$N=41\times 41$. We found that the bulk states and the system as a whole are
almost unaffected by edge disorder, and the average IPR of the system remains
nearly constant. However, the edge states are affected, and their behaviour
depends on the type of disorder introduced. For Anderson disorder, the states
localize continuously along the edges, while for Aubry-André disorder there is
also a threshold behaviour, similarly to what is presented in Fig. 7(c).

Figure 9:  Energy spectra evolution with a) Anderson disorder and b) Aubry-
André disorder, for $B_{y}=d$.

###  IV.2 Energy spectra evolution and density of states

We now consider the system in a mixed $(k_{x},y)$ space, with finite width
along $y$ and OBC. We fix the parameter values as $t=1$, $d=1/6$, $\mu=-3.5$,
and $B_{y}=0.5d$ or $B_{y}=d$ (such that the system describes a $p$-wave
superconductor) and obtain the evolution of the energy spectra for several
values of disorder strength $\lambda$ for both Anderson and Aubry-André
disorder. Since the values of $t$, $d$ and $\mu$ will be kept constant we will
now omit them.

In Fig. 8 we show the energy spectra for $B_{y}=0.5d$ with a) Anderson
disorder and b) Aubry-André disorder. The clean system has gapless edge states
and the bulk gap is not closed by $B_{y}$. As Anderson disorder is increased,
the edge states lose their structure and the bulk gap is closed. Accordingly,
there is an increase in the density of states at $E=0$ and around zero energy
as it can be seen in Fig. 10(a).

Introducing quasidisorder, as seen in panel b), leads to a closing of the bulk
gap with the appearance of new Majorana flat bands. As disorder is increased,
the flat band then splits in two and disappears as a gap opens in the system
for around $\lambda=1.8$. The appearance of MFBs leads to an increase of the
density of states at zero energy, as can be seen in Fig. 10(b) for the value
of $\lambda=1.4$. At higher values of disorder, the system is gapped and the
DOS at $E=0$ goes to zero. The reopening of the gap contrasts with what was
found for Anderson disorder, where the bulk remains gapless as disorder is
increased. We observed that the edge states inside the quasidisorder induced
flat bands appear localized at both edges simultaneously. While the edge
states of the clean system are localized symmetrically on both edges, the flat
band states lose this symmetry and localize more near one of the edges if
quasidisorder is present. Near the edge on which a given state appears less
localized, there is also a deviation from the edge, and the state mostly
localizes on the subsequent sites in $y$.

Figure 10:  Density of states evolution with a) Anderson disorder and b)
Aubry-André (AA) disorder, for $B_{y}=0.5d$ and with c) Anderson disorder and
d) Aubry-André disorder for $B_{y}=d$.

In Fig. 9 the clean system with $B_{y}=d$ is in a gapless phase with both edge
states and a range of $k_{x}$ supporting Majorana flat bands. As Anderson
disorder is increased, the bulk remains gapless and there is a sharp increase
in the density of states at zero energy, as the bulk states come from finite
energies to lower energies. The sharp peak in the DOS observed at $E=0$ is
reminiscent of the characteristic behaviour of a two-dimensional disordered
superconductor with broken time-reversal invariance in the thermal metal
regime [111] in which the density of states displays a logarithmic divergence
at zero energy.

In Fig. 9(b) when Aubry-André disorder is introduced, the edge states appear
to be robust up until around $\lambda\approx 0.8$. However, the MFBs which are
present at $\lambda=0$ are more robust if compared with the edge states, with
the band staying at zero energy but the initial range of $k_{x}$ hosting flat
bands decreasing as $\lambda$ increases. Simultaneously, flat bands appear for
new values of $k_{x}$, as is can be seen in the figure for $\lambda=1.4$, and
accordingly, the density of states at zero energy increases. At higher values
of disorder there is a collapse of states to lower energies and the density of
states exhibits a peak at $E=0$ which is reminiscent of the behaviour found
for Anderson disorder for the same parameter values. Contrary to what is
observed in for a lower magnetic field, there is no opening of the bulk gap
for larger values of $\lambda$. When quasiperiodic disorder is introduced, a
gap will only open for larger values of $\lambda$ if the bulk was gapped prior
to introducing disorder, as in Fig. 8, otherwise the bulk will remain gapless.

Let us now consider the addition of finite values of $\alpha$ and
$\Delta_{s}$. The addition of finite values of spin-orbit coupling and
$s$-wave pairing potential breaks the chiral-like symmetry
$\mathcal{S}_{k_{y}}$ (defined in Eq. 12) that protects the flat bands. As a
result, the latter are lifted to a finite energy and the spectrum acquires a
tilt. For certain regimes of $B_{y}$, the noncentrosymmetric superconductor in
the clean system shows unidirectional edge states. In such regimes, the
addition of Aubry-André disorder leads to the appearance of ”flipped”
unidirectional states in the system. This can be seen in Fig. 11 for
$\lambda=1.4$.

Figure 11: Energy spectra evolution with Aubry-André disorder for $B_{y}=4d$,
$\alpha=0.2d$, $\Delta_{s}=0.3d$ and a) $\lambda=0$, b) $\lambda=1.4$.

###  IV.3 Topological nature of quasidisorder induced flat bands

We want to investigate if the Majorana flat bands that arise in the presence
of a quasiperiodic potential have a topological nature, such as is the case of
the flat bands in the ordered system. Since the Berry phase was found to be
quantized to a value of $\pi$ in the clean system in the region of flat bands,
we calculate it here for the disordered case. The Berry phase $\gamma_{B}$ is
obtained in real space using twisted boundary conditions. Considering a
twisted boundary phase $\theta_{y}$ we have:

$\gamma_{B}(k_{x})=i\int_{0}^{2\pi}d\theta_{y}\langle\Psi(k_{x},\theta_{y})\mid\frac{\partial}{\partial\theta_{y}}\Psi(k_{x},\theta_{y})\rangle$
(27)

where $\Psi$ denotes the ground-state many body wavefunction, which is given
by the Slater determinant of the single particle wavefunctions. We can
represent the ground state wavefunction by an $M\times N$ matrix
$\mathbf{\Psi}^{\theta_{y}}$ where $N$ is the number of sites in $y$ and $M$
is the number of occupied states (negative energy states). Numerically, the
twist variable is discretized into $L$ points between $0$ and $2\pi$, such
that $\theta_{y}$ is constrained to take the values
$\theta_{y,n}=\frac{2\pi}{L}n$, with $n$ an integer that goes from $0$ to
$L-1$. A link variable can then be defined as
$U(\theta_{y,n})=\text{det}\left[\mathbf{\Psi^{\dagger}}_{\theta_{y,n}}\mathbf{\Psi}_{\theta_{y,n+1}}\right]$,
and the Berry phase is obtained as

$\gamma_{B}=-\mathrm{i}\sum_{n=1}^{L}\log{U(\theta_{y,n})}.$ (28)

Figure 12: Energy spectrum and Berry phase $\gamma$ normalized by $2\pi$, as a
function of $k_{x}$. The values of the parameters are $t=1$, $d=t/6$,
$\mu=3d-4t$ and a) $B_{y}=0.5d$, $\lambda=1.4$ b) $B_{y}=d$, $\lambda=1.6$,
with $\lambda$ the strength of the quasiperiodic Aubry-André potential.

We find that at the values of $k_{x}$ where Majorana flat bands appear the
Berry phase is quantized to $\pi$, as shown in Fig. 12. We have also found
that, for the considered system sizes, the values of $k_{x}$ at which the
Berry phase is quantized to $\pi$ are independent of the phase $\phi$ in the
Aubry-André potential.

To quantify the induced bands at zero energy and study the transition to a
$\pi$-quantized Berry phase, we use the concept of Majorana pair density,
defined as [112]:

$\rho_{\gamma}=\frac{N_{\gamma}}{N_{k}},$ (29)

where $N_{k}$ is the number of discrete points of $k_{x}$ taken inside the
interval $\left[-\pi,\pi\right]$, and $N_{\gamma}$ is the number of such
points which support MFBs at the edges. Numerically it is more convenient to
consider the number of $k_{x}$ points for which the Berry phase is quantized
to $\pi$, $N_{\pi}$, since it was found that $N_{\pi}=N_{\gamma}$. A
transition from $\rho_{\gamma}=0$ to $\rho_{\gamma}\neq 0$ then signals a
transition from a trivial to a topological regime ($\pi$-quantized Berry
phase). Fig. 13(a) shows the evolution of $\rho_{\gamma}$ as a function of
quas-idisorder strength for the case $t=1$, $d=t/6$, $\mu=-3.5$ and
$B_{y}=0.5d$, and for the range $\lambda\in[1,2]$. A transition
$\rho_{\gamma}=0\rightarrow\rho_{\gamma}\neq 0$ occurs between $\lambda=1.22$
and $\lambda=1.23$ at a certain critical value $\lambda_{C,1}$. The value of
$\rho_{\gamma}$ grows until $1.49\pm 0.01$ when the flat band splits in two
and the behaviour of $\rho_{\gamma}$ changes, with an abrupt change in the
sign of the second derivative. A second transition occurs between $1.79$ and
$1.8$, at a critical value $\lambda_{C,2}$, where $\rho_{\gamma}$ becomes
zero.

Figure 13: a) Values of $\rho_{\gamma}$ for the case $t=1$, $d=t/6$,
$\mu=-3.5$ and $B_{y}=0.5d$ vs. quasidisorder strength $\lambda$. Obtained for
a system with $76$ sites in $y$. b) Value of the DOS at $E=0$ for the same
parameter values as in a), vs. quasidisorder strength $\lambda$, and the
contribution for $\rho(E=0)$ which comes from the Majorana flat bands in the
corresponding regime.

In Fig. 13(b) the density of states at zero energy $\rho(E=0)$ (normalized by
the system size) is shown, for the same parameters as in Fig. 13(a) and for
$N_{y}=76$, along with the corresponding contribution for the zero energy
density of states which comes from the MFB, $\rho(E=0)_{\gamma}$. Inside the
topological phase, which is highlighted, we can see that the finite value of
$\rho(E=0)$ observed in the system with OBC comes almost entirely from the
presence of flat bands.

$N_{y}$ | $\lambda_{C,1}$ | $\lambda_{C,2}$
---|---|---
76 | $1.225\pm 0.005$ | $1.800\pm 0.005$
100 | $1.215\pm 0.005$ | $1.775\pm 0.005$
175 | $1.230\pm 0.005$ | $1.805\pm 0.005$
200 | $1.220\pm 0.005$ | $1.800\pm 0.005$
400 | $1.230\pm 0.005$ | $1.805\pm 0.005$
800 | $1.225\pm 0.005$ | $1.805\pm 0.005$
Table 1: Values of the critical points $\lambda_{C,1}$ and $\lambda_{C,2}$,
for the parameter values $t=1$, $d=t/6$, $\mu=-3.5$ and $B_{y}=0.5d$ and for
the system sizes $\\{76,100,175,200,400,800\\}$.

Table 1 shows the values of $\lambda_{C,1}$ and $\lambda_{C,2}$ for several
system sizes, obtained for random values of the phase $\phi$ in the Aubry-
André potential, where the uncertainty is taken as the minimum interval
considered between values of $\lambda$. Is is found that the values of the
critical points show little variation with the system size, and we also found
that the critical points are independent of $\phi$ for the system sizes
considered.

###  IV.4 Scaling of the density of states: critical exponents

####  IV.4.1 A detour to the clean system

Let us first briefly consider the clean system, without disorder. For the
clean case, it is possible to obtain the values of the dynamical exponent $z$
and of the critical exponent $\nu$ analytically, for the transition that
occurs as $B_{y}$ is increased, corresponding to a transition from a winding
number of 0 to 1 or a Berry phase of 0 to $\pi$. Here we consider the case of
$\mu<-2t$ (such that the topological phase is within the region described by
Eq. 17). At the topological transition to a gapless phase, the gap closing
points in $k_{x}$, $k_{x,0}$, are given by

$k_{x,0}=\pm\arccos{\left[-\frac{2(t\mu+2t^{2})}{-d^{2}+4t^{2}}\right]}+2n\pi,n\in\mathbb{Z}.$
(30)

The values of $k_{y}$ for which the gap closes are given by $k_{y,0}=n\pi$,
$n\in\mathbb{Z}$ (general solution). In this case the transition happens at
$k_{y,0}=2n\pi$, $n\in\mathbb{Z}$. The gap closes at a critical value of the
magnetic field, $B_{y_{C}}$, which, fixing $k_{y}=k_{y,0}$, is defined from
the value of $k_{x,0}$ as

$B^{2}_{y_{C}}=\left[\mu+2t\left(\cos
k_{x,0}+1\right)\right]^{2}+d^{2}\sin^{2}k_{x,0}.$ (31)

We can now first expand the expressions for the bulk energy around $k_{x,0}$
to find the dependence of the energy on $k_{x}$. We only need to consider the
first positive energy band, $E_{+}(k_{x})$. Taking $k_{y}=k_{y,0}$ and
expanding around $k_{x}=k_{x,0}$ we find

$\begin{split}E_{+}(k_{x})\propto(k_{x}-k_{x,0}),\end{split}$ (32)

implying a value of the dynamical exponent $z=1$ for the transition. We can
now take $k_{x}=k_{x,0}$ and see how the gap closes as a function of $B_{y}$.
We find

$\begin{split}E_{+}(k_{x}=k_{x,0})=\left||B_{y_{C}}|-|B_{y}|\right|.\end{split}$
(33)

Near a quantum phase transition as a critical point $\lambda_{C}$ is
approached, the gap behaves as $\Delta\sim|\lambda-\lambda_{C}|^{z\nu}$,
therefore at $k_{x}=k_{x,0}$ the gap vanishes linearly, with an exponent
$z\nu=1$. Since $z=1$, this implies $\nu=1$, and

$z=1,\quad\nu=1.$ (34)

####  IV.4.2 Quasidisorder: numerical calculation of the critical exponents

Around a critical point, the density of states $\rho(E)$ follows [113]

$\rho(E)=\delta^{(D-z)\nu}f(|E|\delta^{-z\nu}),$ (35)

with $D$ the dimension of the system (here $D=2$),
$\delta=\frac{|\lambda-\lambda_{C}|}{\lambda_{C}}$ the normalized distance to
the critical point $\lambda_{C}$, and $f$ a scaling function. Right at the
critical point, when $\delta=0$, the DOS behaves as

$\rho(E)\sim|E|^{\frac{D}{z}-1}.$ (36)

From the behaviour of the density of states near the phase transition and
using Eqs. 35 and 36 it is possible to obtain the values of the critical
exponents numerically.

Here we study a system with $N_{y}=800$ sites in $y$ and consider the obtained
critical values $\lambda_{C,1}=1.225$ and $\lambda_{C,2}=1.805$ (as shown in
table 1 for this system size). A fit of the form of Eq. 36 for the density of
states at the critical points, done in the interval E $\in[0.005,0.025]$,
gives the values of the critical exponents $z=1.27\pm 0.04$ for the first
transition and $z=1.23\pm 0.03$ for the second transition. To determine the
value of $\nu$ we take values of $\lambda$ inside the topological (gapless)
phase, $\lambda>1.225$ and $\lambda<1.805$, and obtain the density of states
close to zero energy. For small values of $\delta$ and close to zero energy a
collapse of the scaled values of the density of states according to Eq. 35 is
expected.

Figure 14: Density of states for $E\in[0.005,0.025]$ and several values of
$\lambda$ close to the critical values, for a) $\lambda_{C,1}=1.225$ and b)
$\lambda_{C,2}=1.805$, scaled according to Eq. 35 for a) $z=1.27$ and
$\nu=0.95$ and b) $z=1.23$ and $\nu=1.00$.

In Fig. 14 we show the results for the scaled density of states for: a) values
close to the first transition at $\lambda_{C,1}=1.225$, and b) values close to
the second transition at $\lambda_{C,2}=1.805$. The density of states shows a
collapse for a) $z=1.27$ and $\nu=0.95$ and b) $z=1.23$ and $\nu=1.00$.

The quantum phase transitions in the disordered regime are therefore in a
different universality class than that of the clean case, which was found to
behave with $z=\nu=1$. The obtained values also differ significantly from the
known results for the Anderson or the Aubry-André transitions in one
dimension, the first belonging to an universality class with with $\nu=2$ and
$z=2/3$, and the second case with critical exponents $\nu=1$ and $z=2.375$
[114]. Recent results show that for a one dimensional system with $p$-wave
superconductivity subject to an Aubry-André potential the quasidisorder driven
transitions also deviate from the normal Aubry-André class. For the localized-
critical transition line and when the $p$-wave pairing term is finite, the
correlation length exponent has been obtained as $\nu=0.997$ and the dynamical
exponent as $z=1.373$ in [68], and as $\nu=1.000$, $z=1.388$ in [69]. Note,
however, that the referred results are for $D=1$ while we are studying a two
dimensional system, and concern systems with no applied magnetic field.
Nevertheless, one could say that the aforementioned results make it so that
deviations from the known universality classes are also expected for
transitions in the system at study. Up to numerical errors, the values of
$\nu$ obtained for the disordered driven transitions coincide with that of the
Aubry-André transition; nevertheless the value of $z$ deviates from that of
the known classes, which suggests these transitions belong to new universality
classes. The identified transitions, where Marojana flat bands appear as a
result of a quasidisorder induced gap closing, and the subsequent opening of
the bulk gap, are found to happen for other values of the imposed parameters.
Considering the values of the parameters $\mu$, $t$ and $d$ are such that the
topological regions of the superconductor are described by Eq. 17, then as
long as $B_{y}<B_{y,C}$ (when the bulk is gapless) with $B_{y,C}$ defined as
in Eq. 31) the same type of transitions will take place with the increase of
$\lambda$.

$\lambda_{C}$ | $z$ | $\nu$
---|---|---
1.225 | $1.27\pm 0.04$ | $0.95\pm 0.05$
1.805 | $1.23\pm 0.03$ | $1.00\pm 0.05$
Table 2: Values of $z$ and $\nu$ obtained numerically for the topological
transitions for $N_{y}=800$.

###  IV.5 Fractal Analysis

Figure 15: Results of $\tau$ vs. $q$, for several values of disorder strength,
$\lambda$, for $k_{x}=0.02\pi$ and $k_{x}=0.2\pi$, for a) Anderson disorder
and b) Aubry-André disorder. In all cases, the IPR is averaged for the states
within the energy range $E\in[0.05,1]$.

One of the effects of Anderson transitions is the emergence of
multifractality, which is characterized by fluctuations of eigenstates. These
fluctuactions are manifested in the generalized inverse participation ratio.
For a given eigenstate labeled by $m$, the generalized IPR is defined as:

$\text{IPR}(q)_{m}=\sum_{i}\left|\psi_{i}^{m}\right|^{2q},$ (37)

where, as before, $\psi_{i}^{m}$ is the wavefunction of the eigenstate $m$ at
a site $i$. At criticality, the generalized IPR behaves as [115]

$\text{IPR}(q)\sim L^{\tau(q)}$ (38)

where $L$ is the system size and the exponent $\tau(q)$ is defined in terms of
a generalized dimension $D(q)$ as $\tau(q)=D(q)(q-1)$. In a metallic phase,
$D(q)=d$ and for an insulating phase $D(q)=0$. Wavefunction multifractality is
characterized by a $q$ dependent value of $D(q)$, whereas the cases of a
constant $D(q)$ are single fractals [116].

Here we want to make a simple fractal analysis of the system both for disorder
and quasidisorder. We take $k_{x}$ at fixed values, such that system is
reduced to an effective one dimension. The IPR as a function of $q$ is
calculated and averaged within the energy range $E\in[0.05,1]$. We fix the
parameters $t=1$, $d=t/6$, $\mu=3d-4t$ and $B_{y}=0.5d$ and consider both the
cases of Aubry-André disorder and Anderson disorder. The following
subintervals of $L$ are considered, to which a fit of an equation of the form
of Eq. 38 is done:

$\begin{split}&L_{1}=\\\
&\\{75,100,150,175,200,255,275,400,475,600,675,800\\},\\\
&L_{2}=\\{150,175,200,255,275,400,475,600,675,800\\},\\\
&L_{3}=\\{200,255,275,400,475,600,675,800\\},\\\
&L_{4}=\\{275,400,475,600,675,800\\}.\end{split}$ (39)

The obtained results are presented in Figs. 15 and 16.

####  IV.5.1 Anderson disorder

Fig. 15(a) shows the values of $\tau(q)$ for $k_{x}=0.02\pi$ and
$k_{x}=0.2\pi$, for several values of $\lambda$ and considering the system
size interval $L_{1}$. One thing that can be immediately noticed is that for
the clean system, $\lambda=0$, the values of $\tau(q)$ closely follow the line
$\tau(q)=(q-1)$, indicating that $D(q)$ is $q-$independent and equal to $1$.
This is the expected behaviour of the clean system (taking a fixed $k_{x}$
where the system is reduced to one dimension) and reveals that the bulk states
are extended in the $y$ direction. For higher values of disorder, $\tau(q)$
approaches the line $\tau(q)=0$, where $D(q)=0$, suggesting the states are
localized. For other values of disorder strength, starting at $\lambda=0.1$,
$\tau(q)$ does not follow a behaviour characteristic either of $D(q)=1$ or
$D(q)=0$. In order to take a conclusion, it is necessary to evaluate $\tau(q)$
as the system size tends to infinity. To do this, the subintervals of $L$ in
Eq. 39 are considered, to which a fit of equation of the form of Eq. 38 is
done. The results are presented in Fig. 16(a). We find that for the values
$\lambda=0.1$ and above, as larger values of $L$ are considered, the curves
$\tau(q)$ approach $\tau(q)=0$. This confirms a localization of the bulk
states in the thermodynamic limit for small values of disorder.

Figure 16: a) Values of $\tau$ at different values of $q$ and Anderson
disorder strength $\lambda$, for $k_{x}=0.02\pi$ and $k_{x}=0.2\pi$. b) Values
of $\tau$ at different values of $q$ and quasidisorder strength $\lambda$, for
$k_{x}=0.02\pi$ and $k_{x}=0.2\pi$. In all cases, the IPR is averaged for the
states within the energy range $E\in[0.05,1]$.

####  IV.5.2 Aubry-André disorder

Fig. 16(b) shows the values of $\tau(q)$ for $k_{x}=0.02\pi$ and
$k_{x}=0.2\pi$, for several values of quasidisorder strength $\lambda$ for the
size interval $L_{1}$. Unlike the previous case with Anderson disorder, we see
that the results differ for each $k_{x}$, and that for some values of disorder
strength $\tau(q)$ follows the line $q-1$ closely until some value of $q$
where the behaviour suddenly changes. In Fig. 16(b) we show, as before, values
of $\tau(q)$ fitted for the considered size intervals $L_{1}$, $L_{2}$,
$L_{3}$ and $L_{4}$. or lower values of $q$, $\tau$ remains at the values
defined by the Eq. $\tau(q)=D(q)(q-1)$ with $D(q)=1$. However, at higher
values of $q$, this behaviour changes. Contrary to the case with Anderson
disorder, there is no clear tendency of $\tau(q)$ at increased system sizes,
and the behaviour also depends on the value of $q$. This deviation from the
$D(q)=1$ line is verified as soon as disorder is introduced, and suggests the
system is in a multifractal regime. Accordingly, we see the appearance of
critical bulk states in the system. From inspection of Fig. 15(b) and of the
corresponding values of $\tau(q)$ at larger system sizes, we identify a
transition to a localized phase around $\lambda\in[2.0,2,1]$.

##  V Conclusions

In this work we studied a two-dimensional topological superconductor in the
presence of a magnetic field. We introduced disorder and quasidisorder in the
system with the aim of studying the effects on topological and localization
properties. Considering previous results on other systems such as insulators,
semimetals and one-dimensional superconductors, and the results we found on
the effect of quasidisorder in two-dimensional superconductors, we may expect
that the results may be general considering gapless systems or topological
systems (gapped or gapless), in which regions displaying similar topological
and localization properties may be found.

The system was first studied under a perpendicular magnetic field $B_{z}$.
Four types of disorder were considered: Anderson disorder (two-dimensional),
Anderson disorder (one-dimensional along $y$, uniform along $x$), Aubry-André
disorder (one-dimensional along $y$, uniform along $x$) and Aubry-André
disorder (two-dimensional). We observed that the response of the topological
phases of the system differs depending on the type of disorder, and that
quasidisorder induces topological phases in new regions of $B_{z}$,
characterized by integer values of the Chern number $C$. The critical points
at these phase boundaries were shown to become sharper as the system size
increases, allowing us to conclude that the obtained phase diagrams apply to
bigger system sizes.

The real space system was also briefly studied when a parallel magnetic field
is applied in the $y$ direction. We studied the cases of bulk disorder, with
the same four different spatial configurations. For two-dimensional Anderson
disorder, we found that the average IPR of the system increases with an
exponential behaviour as a function of $\lambda$ for $\lambda>1.5$. When
uniformity in the $x$ direction is imposed in the Anderson disorder potential,
we found that the IPR shows a linear increase as a function of $\lambda$, for
$\lambda>0.5$. For Aubry-André disorder, the behaviour of the average IPR of
the system reveals the existence of two different regimes. In the first, the
average IPR shows a slow increase with $\lambda$, and in the second the IPR
greatly increases. The transition between the two regimes is located around
$\lambda=2$.

We then studied the system in a mixed ($k_{x},y$) space with an applied
parallel magnetic field. The clean superconducting system is known to possess
flat bands in the gapless regime. At the corresponding values of $k_{x}$ these
have a winding number $\mathcal{W}$ of 1, which is defined from reducing the
two dimensional system to an effective one dimension. We showed that these are
also characterized by a $\pi-$quantized Berry phase at the same values of
$k_{x}$.

We showed that the introduction of quasidisorder induces new gapless phases in
parameter regimes where they were absent in the clean case. For the $p$-wave
system subject to a parallel magnetic field this leads to new regimes with
Majorana flat bands. This is not only true for phases with a gapless bulk but
also for gapped phases, where quasidisorder closes the bulk gap and Majorana
flat bands appear. We then obtained the Berry phase with twisted boundary
conditions and concluded the quasidisorder induced flat bands also have a
quantized Berry phase of $\pi$. For the noncentrosymmetric superconductor with
added $s$-wave superconducting pairing and Rashba spin orbit coupling, we
showed that new regimes with unidirectional Majorana edge states appear. In
particular, we showed that for a phase where right-moving unidirectional edge
states were present in the system, the introduction of quasidisorder leads to
the appearance of edge modes in the opposite moving direction, and for a
certain quasidisorder range these modes coexist in the system.

Figure 17: Phase diagrams of a system with 20x20 sites indexed by the Chern
number $C$, for several values of Aubry-André disorder strength $\lambda$ and
perpendicular magnetic field $B_{z}$, obtained for the average over 20
disorder configurations. The values of the parameters are $t=1$, $d=1/6$ and
$\mu=4.5$ . Figure 18: Values of the Chern number $C$ vs. one-dimensional
Aubry-André disorder strength $\lambda$ for the system sizes $20\times 20$,
$30\times 30$ and $41\times 41$ and for a) $t=1$, $\mu=0$, $d=0.6$,
$B_{z}=0.4$, b) $\mu=1$, $d=0.6$, $B_{z}=1.1$, c) $\mu=-3.5$, $d=1/6$,
$B_{z}=0.3$. The results were averaged over 20 disorder configurations. Figure
19: Values of the Chern number $C$ vs. two-dimensional Aubry-André disorder
strength $\lambda$ for the system sizes $20\times 20$, $30\times 30$ and
$41\times 41$ and for a) $t=1$, $\mu=0$, $d=0.6$, $B_{z}=0.4$, b) $\mu=1$,
$d=0.6$, $B_{z}=1.1$, c) $\mu=-3.5$, $d=1/6$, $B_{z}=0.3$. The results were
averaged over 10 disorder configurations.

The identification of the quasidisorder induced topologically non-trivial flat
bands with a quantized Berry phase of $\pi$ allowed us to study in detail two
topological transitions, for the $p$-wave superconductor with a parallel
applied magnetic field $B_{y}$. The two critical points were identified and
studied by obtaining the density of induced Majorana bound states in relation
to $k_{x}$ points. We found that the values of the critical points show almost
no variation with the system size for systems bigger than $76$ sites along
$y$. The values of the dynamical critical exponents and correlation length
critical exponents were obtained as $z=1.27\pm 0.04$ and $\nu=0.95\pm 0.05$
for the first critical point and $z=1.23\pm 0.03$, $\nu=1.00\pm 0.05$ for the
second critical point, which puts these transitions in novel universality
classes. We then investigated the multifractal nature of the wavefunctions by
calculating the values of $\tau(q)$ from the IPR values at several values of
disorder, at the same parameter values as those in which the topological
transition was studied. From the behaviour as the thermodynamic limit is
approached, we concluded that the introduction of quasidisorder induces
multifractality in the system. A transition to a localized regime was
identified for $\lambda_{AA}\in[2.0,2.1]$. The same analysis was made for the
system with Anderson disorder. The behaviour of $\tau(q)$ as the system size
tends to infinity suggests that the introduction of Anderson disorder will
drive the system to a localized phase (in the thermodynamic limit).

We have also shown that the average inverse participation ratio is not very
sensitive to the magnitude of $B_{z}$ or $B_{y}$. Although $B_{y}$ leads to
gapless behavior and $B_{z}$ in general leads to gapped behavior, and although
each magnetic field direction leads to different topological properties and
symmetries in the clean system, the localization properties are similar and
the existence of critical states is also similar. It seems that a magnetic
field in the $y$ direction, $B_{y}$, leads to a more localized behavior in the
presence of a quasidisordered potential. For both magnetic field directions we
found critical states and no mobility edges were found. We found a crossover
as a function of $\lambda$ with a mixture of extended and critical states that
grow in number as quasidisorder increases. In this context, Aubry-André along
$1d$ in the two dimensional system or along $2d$ does not lead to
qualitatively different results (in the sense that the crossover in
localization is seen for both cases), besides the differences in induced
topology.

###### Acknowledgements.

We acknowledge partial support from FCT through the Grant UID/CTM/04540/2019.
M.F.M. acknowledges partial support through the grant (1801P.01102.1.01)
QMSP2021 - CEFEMA - IST-ID.

##  Appendix A Additional details on the disordered model under a
perpendicular magnetic field

In Figs. 6(h) and 6(i), we observe that for low magnetic fields the increase
of quasidisorder induces topological phases in regions for which the Chern
number was zero. In the clean system, however, these regions correspond to a
phase that is topological and characterized by a finite value of $I(k_{y})$.
In Fig. 17 we show a phase diagram for the parameters $d=1/6$ and $\mu=4.5$.
In this case, the topological invariant $I(k_{y}=0,\pi)$ is trivial for low
values of magnetic field ($B_{z}<0.5$) when $C=0$. Contrarily to what is
observed in Fig. 6(i), there is no induced topological region with $C=-1$.
This thus suggests that the reentrant regions observed in Fig. 6(h) and Fig.
6(i) can possibly be related with the topological nature of the phases
characterized by $I(k_{y})$ and with $C=0$.

To see how the different critical values for quasiperiodic disorder scale with
the system size, three transitions for the phase diagrams obtained with Aubry-
André disorder along $y$ with uniformity along $x$ (third row in Fig. 6) and
Aubry-André in two dimensions (fourth row in Fig. 6) were considered, at fixed
values of $B_{z}$ and $\mu$. The results are presented in Figs. 18 and 19 for
the system sizes $20\times 20$, $30\times 30$ and $41\times 41$. We found that
within the considered system size range, the transitions become sharper as the
size increases, thus suggesting that the obtained phase diagrams in Fig. 6
should apply to larger systems.

Figure 20:  Average participation ratio, APR, for a $p$-wave superconductor in
the presence of Aubry-André quasidisorder for a) one-dimensional system and b)
two-dimensional system, as a function of the superconducting pairing term $d$
and disorder strength $\lambda$.  Figure 21:  Average participation ratio,
APR, for a $2d$ system with Anderson or Aubry-André (AA) quasidisorder, for
$d=t/6,\mu=-3.5$, in a a) perpendicular or b) parallel magnetic field.

##  Appendix B Participation ratio of $2d$ Aubry-André quasidisorder

In Fig. 7 we considered the average inverse participation ratio regarding both
disorder or quasidisorder along one spatial direction (the $y$ axis) and
disorder or quasidisorder in the plane, comparing a set of values of parallel
magnetic field. In this Appendix we carry out a more extensive analysis. We
want to focus our attention on the regime of increasing disorder strength,
from small values to larger values as localization takes place, in particular
on the possible separation between extended, critical and localized regimes.
The inverse participation ratio is particularly useful to study the transition
to the localized regime, but is not as revealing in the extended-critical
regimes. In this Appendix we will consider instead the participation ratio,
which is given by the inverse of Eq. 26. It is of the order of one for
extended states, and becomes of the order of the inverse of the system size in
the localized regime.

Figure 22:  Average participation ratio as a function of $\lambda$ for
perpendicular magnetic field, $B_{z}$, and a) $1d$ and b) $2d$ quasidisorder
and parallel magnetic field, $B_{y}$, and c) $1d$ and $2d$ quasidisorder.
Figure 23:  Participation ratios of the eigenstates for
$\lambda=1,d=t/6,\mu=-3.5$, as a function of energy for a) $B_{z}=2$ and b)
$B_{y}=2$.

In the case of a one-dimensional $p$-wave superconductor (Kitaev model) in the
presence of an Aubry-André potential and with no magnetic field, it is known
that the three regimes of extended, critical and localized states are present,
as one changes the amplitude of the pairing, $d$, and the quasidisorder
strength, $\lambda$ [42]. The average participation ratio (APR), considering a
single disorder configuration, is shown in Fig. 20(a), where clear transitions
are shown separating the various regimes (here we only consider a disorder
configuration, but the results are characteristic of a larger set of disorder
configurations). The average participation ratio has different plateaus as the
parameters change. Considering a two-dimensional $p$-wave superconductor with
two-dimensional quasidisorder, there is no clear transition between the
regimes and one finds crossovers as $\lambda$ increases. The system sizes
considered in Fig. 20(b) are small ($21\times 21$), but a smooth transition
showing the decreasing of the average participation ratio seems to indicate no
clear separation of a critical regime before the localized phase takes place.

In Fig. 21 we consider the average participation ratio for the two-dimensional
case in the presence of a perpendicular or parallel magnetic field, for the
cases of Anderson disorder and quasidisorder. These results highlight the
extended/critical regimes at lower values of $\lambda$. Anderson disorder
behaves similarly for the two magnetic field directions, and quasidisorder
leads to higher values of the average participation ratio as disorder
increases, with a sharper transition to the localized regime, as shown in Fig.
7. In the case of the parallel magnetic field quasidisorder has a stronger
localization effect. A difference with respect to the perpendicular magnetic
field is the existence of gapless states, more sensitive to the long-range
disorder associated with the Aubry-André potential. In the case of
perpendicular magnetic field, the system remains gapped (except at the
transitions between the various topological regimes) and therefore is expected
to be less sensitive to the quasidisorder potential.

The crossover behavior is clearly seen in Fig. 22, independently of the
magnetic field direction ($B_{z}$ or $B_{y}$). Also, the consideration of
quasidisorder along the $y$ direction and periodic along $x$ or quasidisorder
that is fully two-dimensional leads to similar results. Some differences are
visible for small magnetic fields or small values of $\lambda$. Except for
these regions, the average participation ratio is quite independent of the
amplitude of the magnetic field, but the effect of a parallel magnetic field
is more significant, as discussed.

Figure 24:  Wave function amplitudes for the parameters of Fig. 23, as a
function of space location in a $21\times 21$ system, for states with a)
$B_{z}=2$ and participation ratio $PR=0.781$, b) $B_{z}=2$ and $PR=0.186$, c)
$B_{y}=2$ and $PR=0.643$, d) $B_{y}=2$ and $PR=0.169$.

In order to have a better understanding of the possible existence of critical
states in the regime prior to the transition to localization, we show a
typical case in Fig. 23 of the participation ratios of the various
eigenstates, as a function of their energies, for perpendicular and parallel
magnetic fields. The parameters are chosen so that we are in an intermediate
regime, where the average participation ratio is in the crossover between
fully extended and localized states. There are significant fluctuations
between states with high participation ratios (characteristic of extended
states) and low participation ratios (characteristic of intermediate,
critical, states) but still larger than values that correspond to the
localized regime. The results do not show a mobility edge, and the states mix
throughout the energy range. Also, as $\lambda$ increases, we have found that
the percentage of critical-like states increases, explaining the crossover
behavior. The extended versus critical character nature of the states is shown
in Fig. 24, where a few wave functions are shown (for the system with periodic
boundary conditions in both directions), characteristic of extended and
critical states coexisting in the same energy spectrum.

##  Appendix C Energy spectra evolution for the noncentrosymmetric
superconductor with Aubry-André and Anderson disorder in $(k_{x},y)$ space

Figure 25: Energy spectra evolution with a) Anderson disorder and b) Aubry-
André disorder, for $B_{y}=d$, $\alpha=0.2d$ and $\Delta_{s}=0.5d$.

Here we present results for the evolution of the energy spectra for the
noncentrosymmetric superconductor, with spin orbit coupling $\alpha$ and mixed
$p$ and $s$-wave pairings, subject to Anderson and Aubry-André disorder in the
$(k_{x},y)$ space. We take the same parameter values for $t$, $d$ and $\mu$ as
in section IV.2, and consider two different cases: $B_{y}=d$, $\alpha=0.2d$,
$\Delta_{s}=0.5d$ and $B_{y}=4d$, $\alpha=0.2d$ and $\Delta_{s}=0.3d$.

In Fig. 25 we present the evolution of the energy spectrum for a) Anderson
disorder and b) Aubry-André disorder, for $B_{y}=d$, $\alpha=0.2d$,
$\Delta_{s}=0.5d$. As Anderson disorder is introduced in the system, the edge
states are destroyed and the considered energy range gets filled with bulk
states, but the tilt of the spectrum is preserved. As a result, the flat bands
which were previously lifted due to the introduction of finite values of
$\alpha$ and $\Delta_{s}$ do not collapse to zero energy. For high values of
$\lambda$ the density of states exhibits two peaks which result from the
inclination of the bulk energy spectrum (Fig. 27 a) ).

Figure 26: Energy spectra evolution with a) Anderson disorder and b) Aubry-
André disorder, for $B_{y}=4d$, $\alpha=0.2d$, and $\Delta_{s}=0.3d$.

In Fig. 25(b), as Aubry-André disorder is introduced, we see an evolution that
is similar to Fig. 9(b), but instead of new flat band regimes, new
unidirectional edge states appear. Unlike what happens for Anderson disorder,
at high values of $\lambda$ a gap opens for values of $k_{x}$ around $k_{x}=0$
(although the bulk as a whole remains gapless). This is reflected in the
density of states, that drops around $E=0$ for higher disorder values.
Similarly to what was observed for Anderson disorder, the tilt of the energy
spectrum is preserved as disorder increases.

Figure 27:  Density of states evolution with a) Anderson disorder and b)
Aubry-André (AA) disorder, for $B_{y}=d$, $\alpha=0.2d$ and $\Delta_{s}=0.5d$
and with c) Anderson disorder and d) Aubry-André disorder for $B_{y}=4d$,
$\alpha=0.2d$, and $\Delta_{s}=0.3d$.

In Fig. 26, the clean system is in the regime where unidirectional MESs
appear. The values of the $p$-wave pairing and spin orbit term are kept
constant in relation to the case of Fig. 25, but the $s$-wave pairing term is
decreased from $\Delta_{s}=0.5d$ to $\Delta_{s}=0.3d$ and the magnetic field
is increased from $B_{y}=d$ to $B_{y}=4d$. The spectrum acquires a tilt in the
opposite direction if compared to the clean system in Fig. 25, as a result of
the increased magnetic field. With Anderson disorder, Fig. 26(a), The
unidirectional Majorana edge states are robust to small values of disorder
strength but as disorder increases the structure of the band is lost, as bulk
states fill the lower energy values. This differs from Fig. 25(a) where the
tilt of the spectrum is preserved even at higher values of disorder. As
disorder is increased, there is at first an increase in the value of the DOS
at zero energy, which then decreases for higher values of disorder. For
$\lambda>1.8$ the density of states becomes nearly constant in the considered
range of $E\in[-0.4,0.4]$.

In Fig. 26(b), when a certain value of Aubry-André disorder is reached,
”flipped” unidirectional states appear in the system. This is seen clearly in
Fig. 26(b) for the values of $\lambda=1.2$ and $\lambda=1.4$, as a band with
negative slope appears for values of $k_{x}$ around $k_{x}=0$. At
$\lambda=1.2$ there is a coexistence of unidirectional ”flipped” left-moving
edge modes (with negative slope) around $k_{x}=0$ and right-moving edge modes
(with positive slope) for higher (absolute) values of $k_{x}$. A backflow
current that balances the current on the edges is created on the bulk: extra
right or left moving modes will appear depending on the net current on the
edges.

## References

  * [1] C. W. J. Beenakker, Annu. Rev. Condens. Matter Phys. 4, 113 (2013).
  * [2] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * [3] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * [4] J. Alicea, Rep. Prog. Phys. 75, 076501 (2012).
  * [5] M. Sato and Y. Ando, Rep. Prog. Phys. 80, 076501 (2017).
  * [6] Y. Tanaka, Y. Mizuno, T. Yokoyama, K. Yada, and M. Sato, Phys. Rev. Lett. 105, 097002 (2010).
  * [7] M. Sato, Y. Tanaka, K. Yada, and T. Yokoyama, Phys. Rev. B 83, 224511 (2011).
  * [8] C. L. M. Wong, J. Liu, K. T. Law, and P. A. Lee, Phys. Rev. B 88, 060504 (2013).
  * [9] A. P. Schnyder and S. Ryu, Phys. Rev. B 84, 060504 (2011).
  * [10] A. P. Schnyder, P. M. R. Brydon, and C. Timm, Phys. Rev. B 85, 024522 (2012).
  * [11] G. E. Volovik, Jetp Lett. 107, 516 (2018).
  * [12] S. Peotta and P. Törmaä, Nature Communications 6, 8944 (2015).
  * [13] S. Bose, A. M. García-García, M. M. Ugeda, J. D. Urbina, C. H. Michaelis, I. Bruhega, and K. Kern, LABEL:https://www.nature.com/articles/nmat2768Nat. Materials 9, 550 (2010).
  * [14] M. A. N. Araújo, A. M. García-García, and P. D. Sacramento, Phys. Rev. B 84, 172502 (2011).
  * [15] I. S. Burmistrov, I. V. Gornyi, and A. D. Mirlin, Phys. Rev. Lett. 108, 017002 (2012).
  * [16] J. Mayoh and A. M. García-García, Phys. Rev. B 90, 134513 (2014).
  * [17] B. Fan and A. M. García-García, Phys. Rev. B 102, 184507 (2020).
  * [18] C. Carbillet, V. Cherkez, M. A. Skvortsov, M. V. Feigelḿan, F. Debontridder, L. B. Ioffe, V. S. Stolyarov, K. Ilin, M. Siegel, D. Roditchev, T. Cren, and C. Brun, Phys. Rev. B 102, 024504 (2020).
  * [19] A. Sĺebarski, M. Fijalkowski, P. Zajdel, M.M. Maśka, J. Deniszczyk, M. Zubko, O. Pavlosiuk, K. Sasmal, and M. B. Maple, Phys. Rev. B 102, 054514 (2020).
  * [20] A. Sĺebarski, M. Fijalkowski, M. M. Maśka, J. Deniszczyk, P. Zajdel, B. Trump, and A. Yakovenko, Phys. Rev. B 103, 155133 (2021).
  * [21] M. V. Feigelḿan, L. B. Ioffe, V. E. Kravtsov, and E. A. Yuzbashyan, Phys. Rev. Lett. 98, 027001 (2007).
  * [22] M. V. Feigelḿan, L. B. Ioffe, V. E. Kravtsov, and E. Cuevas, Ann. Phys. 325, 1390 (2010).
  * [23] J. Mayoh, and A. M. García-García, Phys. Rev. B 92, 174526 (2015).
  * [24] K. Zhao, H. Lin, X. Xiao, W. Huang, W. Yao, M. Yan, Y. Xing, Q. Zhang, Z.-X. Li, S. Hoshino, J. Wang, S. Zhou, L. Gu, M. S. Bahramy, H. Yao, N. Nagaosa, Q.-K. Xue, K. T. Law, X. Chen, and S.-H. Ji, Nat. Phys. 15, 904 (2019).
  * [25] C. Rubio-Verdú, A. M. García-García, H. Ryu, D.-J. Choi, J. Zaldívar, S. Tang, B. Fan, Z.-X. Shen, S.-K. Mo, J. I. Pascual, and M. M. Ugeda, Nano Letters 20, 5111 (2020).
  * [26] B. Fan and A. M. García-García, Phys. Rev. B 101, 104509 (2020).
  * [27] Z. Fan, G.-W. Chern, and S.-Z. Lin, Phys. Rev. Research 3, 023195 (2021).
  * [28] M. Stosiek, F. Evers, and I. S. Burmistrov, Phys. Rev. Research 3, L042016 (2021).
  * [29] X. Zhang and M. Foster, Phys. Rev. B 106, L180503 (2022).
  * [30] K.-E. Huhtinen, J. Herzog-Arbeitman, A. Chew, B. A. Bernevig, and P. Törmä, Phys. Rev. B 106, 014518 (2022).
  * [31] Y. Hu, S. M. L. Teicher, B. R. Ortiz, Y. Luo, S. Peng, L. Huai, J. Ma, N. C. Plumb, S. D. Wilson, J. He, and M. Shi, Science Bulletin 67, 495 (2022).
  * [32] K. Jiang, T. Wu, J.-X. Yin, Z. Wang, M. Z. Hasan, S. D. Wilson, X. Chen, and J. Hu, National Science Review nwac199 (2022).
  * [33] X. Wu, T. Schwemmer, T. Müller, A. Consiglio, G. Sangiovanni, D. Di Sante, Y. Iqbal, W. Hanke, A. P. Schnyder, M. M. Denner, M. H. Fischer, T. Neupert, and R. Thomale, Phys. Rev. Lett. 127, 177001 (2021).
  * [34] B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, and S. D. Wilson, Phys. Rev. Lett. 125, 247002 (2020).
  * [35] C. Mielke, D. Das, J.-X. Yin, H. Liu, R. Gupta, Y.-X. Jiang, M. Medarde, X. Wu, H. C. Lei, J. Chang, P. Dai, Q. Si, H. Miao, R. Thomale, T. Neupert, Y. Shi, R. Khasanov, M. Z. Hasan, H. Luetkens, and Z. Guguchia, Nature 602, 245 (2022).
  * [36] J. Sutradhar, S. Mukerjee, R. Pandit, and S. Banerjee, Phys. Rev. B 99, 224204 (2019).
  * [37] M. Saha, B. P. Venkatesh, and B. K. Agarwalla, Phys. Rev. B 105, 224204 (2022).
  * [38] L.-J. Lang, X. Cai, and S. Chen, Phys. Rev. Lett. 108, 220401 (2012).
  * [39] M. Tezuka and N. Kawakami, Phys. Rev. B 85, 140508 (2012).
  * [40] W. DeGottardi, D. Sen, and S. Vishveshwara, Phys. Rev. Lett. 110, 146404 (2013).
  * [41] X. Cai, L.-J. Lang, S. Chen, and Y. Wang, Phys. Rev. Lett. 110, 176403 (2013).
  * [42] J. Wang, X.-J. Liu, G. Xianlong, and H. Hu, Phys. Rev. B 93, 104504 (2016).
  * [43] X. Cai, J. Phys.: Condens. Matter 26, 155701 (2014).
  * [44] Y. E. Kraus, Y. Lahini, Z. Ringel, M. Verbin, and O. Zilberberg, Phys. Rev. Lett. 109, 106402 (2012).
  * [45] S. Nakajima, N. Takei, K. Sakuma, Y. Kuno, P. Marra, and Y. Takahashi, Nat. Phys. 17, 844 (2021).
  * [46] S. Longhi, Opt. Lett. 45, 4036 (2020).
  * [47] O. Zilberberg, Opt. Mater. Express 11, 1143 (2021).
  * [48] T. Liu and H. Guo, Phys. Letters A 382, 3287 (2018).
  * [49] M. I. N. Rosa, Y. Guo, and M. Ruzzene, Appl. Phys. Lett. 118, 131901 (2021)
  * [50] M. Verbin, O. Zilberberg, Y. E. Kraus, Y. Lahini, and Y. Silberberg,  Phys. Rev. Lett. 110, 076403 (2013).
  * [51] Y. E. Kraus and O. Zilberberg, Phys. Rev. Lett. 109, 116404 (2012).
  * [52] J. H. Pixley, J. H. Wilson, D. A. Huse, and S. Gopalakrishnan, Phys. Rev. Lett. 120, 207604 (2018).
  * [53] V. Goblot, A. Štrkalj, N. Pernet, J. L. Lado, C. Dorow, A. Lemaître, L. Le Gratiet, A. Harouri, I. Sagnes, S. Ravets, A. Amo, J. Bloch, and O. Zilberberg, Nat. Phys. 16, 832 (2020).
  * [54] T. Xiao, D. Xie, Z. Dong, T. Chen, W. Yi, and B. Yan, Science Bulletin 66, 2175 (2021).
  * [55] F. A. An, E. J. Meier, and B. Gadway, Phys. Rev. X 8, 031045 (2018).
  * [56] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, and M. Inguscio, Nature 453, 895 (2008).
  * [57] L. Guidoni, C. Triché, P. Verkerk, and G. Grynberg, Phys. Rev. Lett. 79, 3363 (1997).
  * [58] H. P. Lüschen, S. Scherg, T. Kohlert, M. Schreiber, P. Bordia, X. Li, S. Das Sarma, and I. Bloch, Phys. Rev. Lett. 120, 160404 (2018).
  * [59] Y. Lahini, R. Pugatch, F. Pozzi, M. Sorel, R. Morandotti, N. Davidson, and Y. Silberberg, Phys. Rev. Lett. 103, 013901 (2009).
  * [60] S. Aubry and G. André, Ann. Israel Phys. Soc 3, 18 (1980).
  * [61] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * [62] R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad. Sci. U.S.A. 108, 12233 (2011).
  * [63] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 86, 155449 (2012).
  * [64] D. Mao and T. Senthil, Phys. Rev. B 103, 115110 (2021).
  * [65] M. Gonçalves, H. Z. Olyaei, B. Amorim, R. Mondaini, P. Ribeiro, and E. V. Castro, 2D Mater. 9, 011001 (2022).
  * [66] J. Wang, X.-J. Liu, G. Xianlong, and H. Hu, Phys. Rev. B 93, 104504 (2016).
  * [67] J. Fraxanet, U. Bhattacharya, T. Grass, D. Rakshit, M. Lewenstein, and A. Dauphin, Phys. Rev. Research 3, 013148 (2021).
  * [68] X. Tong, Y.-M. Meng, X. Jiang, C. Lee, G. D. de M. Neto, and G. Xianlong, Phys. Rev. B 103, 104202 (2021).
  * [69] T. Lv, T.-C. Yi, L. Li, G. Sun, and W.-L. You, Phys. Rev. A 105, 013315 (2022).
  * [70] D. R. Hofstadter, Phys. Rev. B 14, 2239 (1976).
  * [71] S. Das Sarma, Song He, and X. C. Xie Phys. Rev. Lett. 61, 2144 (1988).
  * [72] S. Das Sarma, Song He, and X. C. Xie, Phys. Rev. B 41, 5544 (1990).
  * [73] J. Biddle and S. Das Sarma, Phys. Rev. Lett. 104, 070601 (2010).
  * [74] S. Ganesham, J. H. Pixley and S. Das Sarma, Phys. Rev. Lett. 114, 146601 (2015).
  * [75] S. Gopalakrishnan, Phys. Rev. B 96, 054202 (2017).
  * [76] A. Duthie, S. Roy, and D. E. Logan, Phys. Rev. B 103, L060201 (2021).
  * [77] T. Liu, X. Xia, S. Longhi, L. S.-Palencia, SciPost Phys. 12, 027 (2022).
  * [78] S. Ganeshan, K. Sun, and S. Das Sarma, Phys. Rev. Lett. 110, 180403 (2013).
  * [79] M. Gonçalves, B. Amorim, E. V. Castro and P. Ribeiro, SciPost Phys. 13, 046 (2022).
  * [80] A. P. Mackenzie and Y. Maeno, Rev. Mod. Phys. 75, 657 (2003).
  * [81] R. Joynt and L. Taillefer, Rev. Mod. Phys. 74, 235 (2002).
  * [82] K. Matano, M. Kriener, K. Segawa, Y. Ando, and G. Zheng, Nature Phys 12, 852 (2016).
  * [83] E. Bauer, G. Hilscher, H. Michor, Ch. Paul, E. W. Scheidt, A. Gribanov, Yu. Seropegin, H. Noël, M. Sigrist, and P. Rogl, Phys. Rev. Lett. 92, 027003 (2004).
  * [84] I. Sugitani, Y. Okuda, H. Shishido, T. Yamada, A. Thamizhavel, E. Yamamoto, T. D. Matsuda, Y. Haga, T. Takeuchi, R. Settai, and Y. Ōnuki, J. Phys. Soc. Jpn. 75, 043703 (2006).
  * [85] N. Kimura, K. Ito, K. Saitoh, Y. Umeda, H. Aoki, and T. Terashima, Phys. Rev. Lett. 95, 247004 (2005).
  * [86] S. Fujimoto, J. Phys. Soc. Jpn. 76, 051008 (2007).
  * [87] P. A. Frigeri, D. F. Agterberg, A. Koga, and M. Sigrist, Phys. Rev. Lett. 92, 097001 (2004).
  * [88] M. Sato and S. Fujimoto, Phys. Rev. B 79, 094504 (2009).
  * [89] A. Daido and Y. Yanase, Phys. Rev. B 95, 134507 (2017).
  * [90] A. V. Balatsky, I. Vekhter, and J.-X. Zhu, Rev. Mod. Phys. 78, 373 (2006).
  * [91] I. Reis, D. J. J. Marchand, and M. Franz, Phys. Rev. B 90, 085124 (2014).
  * [92] P. D. Sacramento, J. Phys.: Condens. Matter 27, 445702 (2015).
  * [93] K. Pöyhönen, I. Sahlberg, A. Westström, and T. Ojanen, Nat Commun 9, 2103 (2018).
  * [94] T. C̆adez̆, R. Mondaini, E. V. Castro, and P. D. Sacramento, Acta Phys. Pol. A 135, 1180 (2019).
  * [95] T. Devakul and D. A. Huse, Phys. Rev. B 96, 214201 (2017).
  * [96] J. T. Chalker and G. J. Daniell, Phys. Rev. Lett. 61, 593 (1988).
  * [97] J. T. Chalker, Physica A: Stat. Mech. Appl. 167, 253 (1990).
  * [98] B. Huang and W. V. Liu, Phys. Rev. B 100, 144202 (2019).
  * [99] J. H. Pixley, P. Goswami, and S. Das Sarma, Phys. Rev. Lett. 115, 076601 (2015).
  * [100] Y. Fu, J. H. Wilson and J. H. Pixley, Phys. Rev. B 104, L041106 (2021).
  * [101] M. Buchhold, S. Diehl, and A. Altland, Phys. Rev. B 98, 205134 (2018).
  * [102] M. Gonçalves, P. Ribeiro, E. V. Castro and M. A. N. Araújo, Phys. Rev. Lett. 124, 136405 (2020).
  * [103] Y. Fu, E. J. König, J. H. Wilson, Y.-Z. Chou, and J. H. Pixley, npj Quantum Materials (2020)5:71.
  * [104] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Phys. Rev. B 96, 245115 (2017).
  * [105] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * [106] R. Resta, J. Phys.: Condens. Matter 22, 123201 (2010).
  * [107] T. Fukui, Y. Hatsugai, and H. Suzuki, J. Phys. Soc. Jpn. 74, 1674 (2005).
  * [108] Y.-F. Zhang, Y.-Y. Yang, Y. Ju, L. Sheng, R. Shen, D.-N. Sheng, and D.-Y. Xing, Chinese Phys. B 22, 117312 (2013).
  * [109] Here we limit ourselves to the analysis of a system of size $41\times 41$, although draft results for a system of size $20\times 20$ followed the same behaviour as the one discussed.
  * [110] R. Queiroz and A. P. Schnyder, Phys. Rev. B 89, 054501 (2014).
  * [111] A. Mildenberger, F. Evers, A. D. Mirlin, and J. T. Chalker, Phys. Rev. B 75, 245321 (2007).
  * [112] N. Sedlmayr, J. M. Aguiar-Hualde, and C. Bena, Phys. Rev. B 91, 115415 (2015).
  * [113] K. Kobayashi, T. Ohtsuki, K.-I. Imura, and I. F. Herbut, Phys. Rev. Lett. 112, 016402 (2014).
  * [114] J. C. C. Cestari, A. Foerster, M. A. Gusmão, and M. Continentino, Phys. Rev. A 84, 055601 (2011).
  * [115] F. Wegner, Z Physik B 36, 209 (1980).
  * [116] F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).