Random free-fermion quantum spin chain with multispin interactions

We study the effects of quenched disorder in a class of quantum chains with ( p + 1)-multispin interactions exhibiting a free fermionic spectrum, paying special attention to the case p = 2. Depending if disorder couples to (i) all the couplings or just to (ii) some of them, we have two distinct physical scenarios. In case (i), we find that the transitions of the model are governed by a universal infinite-randomness critical point surrounded by quantum Griffiths phases similarly as happens to the random transverse-field Ising chain. In case (ii), we find that quenched disorder becomes an irrelevant perturbation: the clean critical behavior is stable and Griffiths phases are absent. Beyond the perturbative regime, disorder stabilizes a line of finite-randomness critical points (with nonuniversal critical exponents), that ends in a multicritical point of infinite-randomness type. In that case, quantum Griffiths phases also appear surrounding the finite-disorder transition point. We have characterized the correlation functions and the low-temperature thermodynamics of these chains. Our results are derived from a strong-disorder renormalization-group technique and from finite-size scaling analysis of the spectral gap computed exactly (up to sizes ∼ 10 7 ) via an efficient new numerical method recently introduced in the literature [Phys. Rev. B 104 , 174206 (2021)].


I. INTRODUCTION
The importance of studying one-dimensional quantum models is invaluable.Due to the peculiarities of the phase space in d = 1, many models can be precisely described and, thus, serve as important test beds for physical insights [1].Among those models, there is an important class that goes under the general name of free systems.These are systems (of volume V ) whose the exponentially large Hilbert space of dimension ∼ e V can be expressed as a combination of ∼ V quasienergies of a free-particle system.Despite the name, they exhibit nontrivial phenomena, such as phase transitions and zero-energy fractional edge modes, and are, thus, good starting points to understand many complex behaviors of more general systems [2][3][4].
Initially, the known free systems models were linked to a Jordan-Wigner transformation which maps interacting spin models into free fermionic particles, i.e., into models which are bilinear in the fermionic operators or fields.Later, it was realized that the existence of this transformation is not a necessary condition [5].It is possible to construct the eigenspectrum of the spin system from the free-fermion pseudoenergies obtained from the roots of a characteristic polynomial.This polynomial is constructed thanks to the infinite number of conserved charges of the model [6][7][8].Interestingly, these latter free systems, a priori, cannot be written in term of local bilinears of fermion operators.
The first models in this class are the Z N -symmetric free parafermionic models [5,[9][10][11] and the three-spin interacting Fendley model [6].Actually, it has been found that Fendley model belongs to a large family of spin chains with multispin interactions and Z N -symmetry [7,8].For N > 2, the spectrum is complex and has a free parafermionic form.Furthermore, it was shown that the free-particle pseudo-energies of some of the above models are also the pseudo-energies of a multispin U(1) symmetric XY model [12,13].Additional developments related to Fendley model can be found in Refs.[14][15][16][17].
Recently, it was shown that, in general, the roots of the characteristic polynomial yielding the free fermionic spectrum can be efficiently obtained numerically [12].As a consequence, the finite-size gaps of the spin system can be computed straightforwardly for quite large system sizes.Furthermore, the numerical cost for exactly computing the finitesize gap with machine precision (and, hence, the associated dynamical critical exponent) is minimum at criticality increasing only linearly with the chain size.While this may seems innocuous for translation-invariant systems where analytical results can often be obtained, it is of great applicability to quenched disordered systems where analytical results are scarce and exactly numerical results are plagued by numerical instabilities inherent to the extreme slow critical dynamics of these systems [18].This brings us to the topic of our research.What are the effects of quenched disorder on these generalized free fermionic systems?
Despite all the progress on characterizing the infiniterandomness criticality, we currently do not know what are the necessary conditions ensuring its appearance. 1 It is then desirable to further study other aspects which were not considered in the previous studies.Here, we consider interactions involving more than the conventional two-body interactions.We pay special attention to the Z 2 -symmetric free fermionic case in which the interactions involve three consecutive spins [6].The clean phase diagram has three critical lines separating three gapful phases which are related to each other by triality. 2The universality classes of the associated transitions are that of the transverse-field Ising chain, but that of the multicritical point (where these three critical lines meet) is yet to be determined.Currently, only its dynamical critical exponent z = 3  2 [6] and the specific-heat exponent α = 0 [8] are known.An inhomogeneous quantum Ising chain sharing the same energy spectrum of this multicritical point was introduced in Ref. 16.In this related Ising chaing the order-parameter exponent is β = 1  8 like the standard Ising chain.Similarly to the random transverse-field Ising chain and to the spin-1/2 XX chain with random couplings, we show in this paper that generic quenched disorder stabilizes quantum Griffiths phases surrounding the three transition lines.In addition, the corresponding universality classes of these transitions and that of the multicritical point are the same and are of infiniterandomness type (with disorder-independent tunneling exponent ψ = 1  2 ).Our results are based on a generalization of the strong-disorder renormalization-group (SDRG) method [74][75][76] and on exact numerical calculations of the spectral gap using the aforementioned method based on the characteristic polynomial that allows us to obtain exact numerical results for very large lattice sizes [12].After unveiling the structure of the SDRG flow, we generalize our results to the case of (p + 1)-multispin interactions (p = 1, 2, 3, . . . ) and arrive basically at the same conclusions.All the transitions are in the same universality class of infinite-randomness type with the same tunneling exponent ψ = 1 2 .In addition, we have also studied the case in which disorder couples only to one type of coupling constants (the others remaining homogeneous).In that case, some of the clean transition lines remain stable for weak disorder strength while a line of finite-disorder fixed points (with nonuniversal dynamical critical exponents z) appear for large disorder strengths, and terminates in an infinite-randomness multicritical point.
This paper is organized as follows.In Sec.II we define the model studied and review some key results important to our purposes.In Sec.III we overview the expected effects of quenched disorder in a heuristic way.Our arguments are mostly based on the effects caused by rare-regions in the nearcritical Griffiths phases.In Sec.IV we review the SDRG method for the standard 2-spin interacting case and generalize it to the 3-spin interacting case.In Sec.V we report our numerical study of the finite-size gap statistics of the model which are in agreement with the SDRG results.We present further discussions and concluding remarks in Sec.VI.Finally, the details of the renormalization-group flow are presented in Appendices A and B.

II. THE MODEL AND REVIEW OF KEY RESULTS
We consider the (p + 1)-multispin interacting quantum chains (p = 1, 2, 3, . . .), whose Hamiltonian, introduced in Refs.7, 8, is given by σ x,y i are Pauli matrices associated with spin-1 2 degrees of freedom at site i, and L is the total number of {h i } energy-density operators (which is also the total number of spins in the chain).The case p = 1 is equivalent (up to global degeneracies) to the inhomogeneous transverse-field Ising quantum chain.The local interaction operator h i involves min {p + 1, i} spins and satisfy the algebra (for i ̸ = j) h i , h j = 0, if |i − j| ≤ p, and [h i , h j ] = 0 otherwise.(3) In other words, h i and h j commute if they are farther than p lattice units (|i − j| > p), and anticommute otherwise.Evidently, from Eq. ( 1), h 2 i = 1.Finally, λ i is the local multispin energy coupling.In this work, we introduce quenched disorder by considering {λ i } as independent random variables.Their precise distribution will be defined later.Interestingly, it was shown [7,8] that the spectrum of (1) has the free fermionic form where s k = ±1, and the free fermionic pseudo-energies ε k = 1/ √ x k , with {x k } being the roots of the polynomial of degree L = L+p p+1 (with ⌊x⌋ being the integer part of x): whose coefficients C L (ℓ) are obtained from the recurrence relation with the initial condition P j (x) = 1 for j ≤ 0.
It is important to notice that the free fermionic character is guaranteed when open boundary conditions are applied.For other boundary conditions, very likely this is not true [7,8], and the solution of the model remains an open problem.
The first gap in the spectrum energy is ∆ = 2/ √ x max , where x max = max {x k } is the largest root of the polynomial (5).At and near criticality, it was recently shown [12] that x max can be efficiently computed for very large chains even though C L (ℓ) grows factorially with ℓ. (Indeed, there is no need to compute all C L (ℓ).)This is accomplished when one uses the Laguerre's upper bound (LB) for the roots of a polynomial where as the starting initial guess for the x max .
In the (critical) homogeneous case λ i = λ , it was shown [12] that, for any p, the quantity with 0 < ε < 1 being a constant.Thus, ∆ LB has the same finite-size scaling properties of the finitesize gap ∆ and, thus, can be used to obtain the dynamical critical exponent z, i.e., with z LB = z = p+1 2 [7,8].Numerically, this is convenient since ∆ LB can be efficiently computed.
For p = 1, it was shown [12] that, in the critical quenched disordered case ({λ i } being independent and identically distributed random variables), ∆ LB = (1 − ε L ) ∆ with ε L vanishing slowly as L → ∞.This provides a convenient tool to obtain the dynamical scaling relation.In this case, with universal tunneling exponent ψ = 1/2.By universal we mean that ψ does not depend on the particular distribution of the disorder variables. 3This result is expected since the model Hamiltonian (1) for p = 1, apart from global degeneracies, has the same spectrum as the transverse-field Ising chain.It is worthy noting that the activated dynamical scaling (10) has been confirmed by many analytical and numerical studies [66].In addition, ∆ LB can also be used to study the finite-size gap in the near critical Griffiths phase.In this phase, the system is gapless even though it is short-range correlated (finite spin-spin correlation length) [29].The finite-size gap ∆ and the associated LB estimate, ∆ LB , obeys the power-law scaling (9) with z being the off-critical (Griffiths) dynamical exponent which depends on the distance from criticality.As criticality is approached, the off-critical dynamical exponent increases and becomes infinite at the critical point, in accordance to the activated critical dynamical scaling (10).
The model (1) for p = 2 and nondisordered couplings was studied in Ref. 6 for couplings of period 3 [which is the natural period given by the algebra (3)], i.e., Exploring the triality of the model, the phase diagram was determined (see Fig. 1).We now want to discuss what we mean by triality in this model.In the bulk limit, the lattice translations , and the lattice reflection (i → L + 1 − i) do not change the algebra (3) and, therefore, the spectrum.Thus, Now, consider the case α = λ A /λ C < 1 fixed.If the transition happens at β = λ B /λ C ̸ = 1, the relation (12) implies the existence of another transition at β −1 for the same α.Since there is only a global Z 2 symmetry to be broken, then there is only a single phase transition.Therefore, the transition must happen at β c = β −1 c = 1.Successive applications of the relation (12) imply that the self-triality curves (red lines in Fig. 1) are the transition lines.Notice that this argument can be generalized to higher values of p where the self-(p + 1)-ality hyperplanes are the transition manyfolds.For p = 3, see Ref. 14.
The three different phases are characterized by the following expectation values: h A = |∑ i ⟨h 3i−2 ⟩|, h B = |∑ i ⟨h 3i−1 ⟩|, and h C = |∑ i ⟨h 3i ⟩|.For λ A,B < λ C , the system is in a phase where h C > h A,B .By symmetry or, more precisely, by triality, there are other two phases in which h A > h B,C and h B > h A,C which happen when λ C,B < λ A and λ A,C < λ B , respectively.There are three phase transition boundaries: As a final remark motivating this work, we recall that multibody interactions are not an uncommon feature in condensedmatter physics.It naturally arises in Mott insulators and, specifically, is a key ingredient for stabilizing chiral spin liquids [80,81].In addition to that, transitional metal oxides naturally have four-body interactions as described by the Kugel-Khomskii model [82].Multispin interactions also appear in quantum spin chains due to another reason: the mapping between different models.For instance, there is an equivalence of the Ising model with two-and three-spin interactions with the four-state Potts model [83].In addition, it is well known a mapping between the quantum Ashkin-Teller model in 1D, which has four-spin interactions, and the XXZ spin-1/2 chain, which has only conventional two-spin interactions [84].Finally, the multispin interactions should not be seen as the main ingredient for novel physics.The main ingredient is the algebra of the local Hamiltonian operators (3) that solely determines the free fermionic spectrum.The representation of this algebra in our work involves multispin interactions.But this is not a necessary condition.

III. OVERVIEW OF THE EFFECTS OF QUENCHED DISORDER
In this work, we study the effects of quenched disorder on the Hamiltonian (1), paying special attention to the first nontrivial case p = 2.We inquire how disorder on the coupling constants changes the clean phase diagram Fig. 1 as well as the universality classes of the transitions.
We build our arguments taking as the starting point the physical behavior of the clean system [6] revised in Sec.II.For this sake, we assume that the set of couplings {λ A,i } ≡ {λ 3i−2 }, {λ B,i } ≡ {λ 3i−1 }, and {λ C,i } ≡ {λ 3i } are random variables respectively distributed according to the probability distributions P A (λ ), P B (λ ) and P C (λ ).

A. The simpler case of a vanishing coupling λ A,i
When one of the couplings is vanishing, say λ A,i = 0, the effective algebra (3) is that of the model with p = 1 [see, also, Eq. ( 6)], which corresponds to the standard algebra of the transverse-field Ising chain.In that case, the phase diagram is that of the random transverse-field Ising chain which is well known [29,85].The transition takes place when the typical values (geometric mean) of the remaining couplings equal each other [86], i.e., the system is critical when δ = 0, where δ = δ i , with • • • denoting the disorder average and For δ > 0 (δ < 0), the system is in the B-(C-)phase.

Uncorrelated disorder
According to the Harris' criterion [19], uncorrelated quenched disorder is a relevant perturbation at δ = 0 (since dν = 1 < 2, with d = 1 being the number of spatial dimensions in which disorder is uncorrelated and ν = 1 being the correlation length exponent of the clean theory) and, thus, the universality class of the transition must change.As shown by Fisher [28,29], the universality class is of infinite-randomness type with activated dynamical scaling (10).
In addition, surrounding this exotic quantum critical point, there are Griffiths phases whose spectral gap vanishes and the spin-spin correlations are short-ranged.The off-critical dynamical scaling is a power-law (9) with effective dynamical exponent z ∝ δ −1 diverging as the system approaches criticality.

Appearance of locally correlated disorder
On the other hand for λ B,i = e ε λ C,i with ε being a constant, i.e., for perfectly correlation between the local random variables in lattices B and C, the Harris criterion has to be applied with caution.This is because the local distance from criticality δ i = ε is uniform throughout the chain.It turns out that disorder is an irrelevant perturbation [85].The clean critical behavior is stable up to some critical disorder strength, beyond which it changes to a finite-randomness critical behavior, where the dynamical critical scaling is a conventional power-law (9) but with a nonuniversal critical dynamical exponent z, i.e., it depends on the disorder strength [18,87].In addition, no Griffiths effects exist in this case.
In this work, we do not consider explicitly the case where λ B,i and λ C,i are locally correlated.However, we do consider the case in which both couplings are uniform (λ B,i = λ B and λ C,i = λ C ) and the third one (λ A,i ) is random.If the typical The multicritical point (in both cases) is also in the same infiniterandomness universality class.The blue dotted line in panel (a) is a transition line where the universality class is of finite-randomness type (1 < z < ∞).The phases have the same nature as the homogeneous case (see Fig. 1) and the shaded regions delimits the associated Griffiths phases where the spectrum gap vanishes.value of λ A,i is sufficiently small, we show that quenched disorder is perturbatively irrelevant in the renormalization-group (RG) sense.Thus, disorder can be simply ignored, and the transition at λ B = λ C is in the universality class of the clean system.As we increase the values of λ A,i 's beyond the perturbative regime, disorder induces randomness in the renormalized couplings λB,C .Interestingly, the induced disorder has a perfect correlation, namely, λB,i = λC,i .Thus, the longdistance physics of locally correlated random couplings naturally appears in this model.Consequently, a finite-randomness critical point governs the transition for sufficiently large λ A,i 's.
B. The boundary phases: the case of small λ A,i couplings What are the effects of small coupling λ A,i ?In the clean case, due to triality, small λ A cannot shift the location of the critical point λ B = λ C (see Fig. 1) (see Ref. 6 for another argument).In the disordered case, we numerically show (see Sec. V) that this remains true, i.e., the critical point remains at δ = ln λ B,i − ln λ C,i = 0 provided that ln λ A,i < ln λ B,i .Thus, the phase transition lines of the phase diagram, in the random system, are equal to those of the clean system with λ A , λ B , and λ C replaced by their typical values λ A,typ , λ B,typ , and λ C,typ as sketched in Fig. 2.
How about the universality classes of the transitions?Here, we consider two cases: two competing (strongest) couplings (a) not generating random mass, and (b) do generating random mass δ i , see Eq. (13).
In case (a) the two strongest couplings are homogeneous (λ B,i = λ B and λ C,i = λ C ) and λ A,i is typically much smaller than λ B,C , the transition is in the universality class of the clean system (Ising) as shown by the solid red line in Fig. 2(a).When approaching the multicritical point, however, the weak disordered couplings {λ A,i } become nonperturbative and a line of finite-randomness fixed points emerges (as previously mentioned).The resulting universality class of the transition [blue dotted line in Fig. 2(a)] has critical dynamical exponent z larger than the unity and diverges as the infinite-randomness multicritical point is reached.
In case (b) the two competing couplings generate a random mass (δ 2 − δ 2 ̸ = 0).This happens whenever either one or both couplings are independent random variables.The clean (Ising) universality class in unstable since the Harris criterion is violated.The resulting universality class is the one of the random transverse-field Ising chain with activated dynamical scaling (10), and the associated phase boundaries are the dashed lines in the phase diagram of Fig. 2.
C. The multicritical point: the case of strong λ A,i couplings What is the change of the clean multicritical point in the presence of disorder?We cannot apply the Harris criterion since the clean correlation length exponent ν is not known.
To answer this question, we develop an appropriate strongdisorder renormalization-group technique (see Sec. IV).Our results strongly indicate that quenched disorder is a relevant perturbation.In addition, we show that the resulting universality class is of infinite-randomness type with activated dynamics (10) and universal tunneling exponent ψ = 1 2 .We have also confirmed these results by numerically studying the finite-size gap of the system (see Sec. V).

D. The presence of quantum Griffiths phases
We now inquire about the off-critical properties.Are they affected by quenched disorder?The nature of the phases do not change since disorder on the coupling constants does not break neither a symmetry of the Hamiltonian nor a symmetry of the order-parameter field, i.e., disorder does not couple directly to the order-parameter field in the associated underlying field theory.Thus, the phase diagram has the same phases as the clean one.However, near the transition lines, random mass induces Griffiths phases.In those regions of the phase diagram (shaded areas in Fig. 2), the spectral gap vanishes and the correlations remain short-ranged.The finite-size gap scaling is of power-law type (9) with nonuniversal (i.e., disorderdependent) effective dynamical exponent z.
Here, the Griffiths phases can be understood through the lenses of the so-called rare regions (RRs): large and rare spatial regions in a phase locally different from the bulk.For definiteness, consider the case in which λ C,i = λ C (i.e., uniform) and λ A,i 's are distributed between λ A,min and λ A,max , with 0 < λ A,min < λ A,max .Evidently the typical value λ A,typ = exp(ln λ A ) is between λ A,min and λ A,max .To start, let us analyze the phase transition between the A-(h A > h B,C ) and Cphases (h C > h A,B ) where we can disregard the weak coupling λ B,i (say, for simplicity, that max{λ B,i } < min{λ A,min , λ C }).As previously discussed, the transition takes place when λ C = λ A,typ .When λ C ≫ λ A,max , the system is deep in the homogeneous C-phase.Its properties are just the one of the clean systems with the random couplings λ A(B),i replaced by their typical value λ A(B),typ .Importantly, the spectral gap is finite.
When λ A,typ < λ C < λ A,max , on the other hand, there are RRs where the local couplings λ A,i are typically greater than λ C .Being locally in the A phase, they endow the system a high A-phase susceptibility.The spin in the domain walls between the Aand C-phases can be arbitrarily weakly coupled and are responsible for the low-lying excitations closing the spectral gap.As neither the bulk nor the RRs are critical, the corresponding correlation length is finite.By duality, an analogous Griffiths phase appears when λ A,min < λ C < λ A,typ .
Interestingly, there is a simple quantitative argument providing the closing of the spectral gap in the Griffiths phase.Consider a RR of size L RR .The effective interaction between the domain wall spins are thus of order J DW ∼ e −L RR /ξ RR , where ξ RR is the corresponding correlation length in that particular RR.For simplicity, we will consider ξ RR = ξ to be RR-independent (a more precise treatment can be found in Ref. 20).The reason for J DW being exponentially small is because the RR itself does not harbor Goldstone modes since the symmetry of the Hamiltonian interactions is discrete.The system low-energy density of states ρ DOS is dominated by the excitations of the weakly coupled domain walls.Ignoring the even weaker coupling to other domain wall spins belonging to other RRs, then ρ DOS (ω) ∼ dL RR w RR (L RR ) δ (ω − J DW ).
Here, we simply sum of all possible RRs weighting their contribution by their existence probability w RR ∼ e −L RR /ℓ , with ℓ ∝ −1/ ln p, and p being the probability of λ A,i being greater than λ C .Notice that the probability of finding a RR decreases exponentially with its volume, and ℓ is a constant that depends on the distribution's details of the coupling constants.Consequently, one finds that ρ DOS ∼ ω −1+1/z , with dynamical Griffiths exponent z = ℓ/ξ .Notice the absence of a gap or a pseudogap in ρ DOS .Actually, there is a divergence in the low-energy density of states when z > 1 and ω → 0. We see that z diverges ∼ 1/δ when approaching the transition.
By triality, the resulting Griffiths phases are those sketched in Fig. 2(b) if the B-couplings are also randomly distributed between λ B,min and λ B,max (0 < λ B,min < λ B,max ).
Finally, near the multicritical point there are Griffiths phases with RRs locally belonging to, say, either A-phase or B-phase, while the bulk is in the C-phase.A similar feature also appeared in the Griffiths phase of quantum Ashikin-Teller chain [88,89].In those cases, the effective dynamical exponent z = max{z A , z B }, where z A(B) is the dynamical exponent provided by the Griffiths singularities of the A-(B-)RRs.

E. The absence of Griffiths phases
If, on the other hand, the B-couplings are also homogeneous (λ B,i = λ B ), the resulting transition between the Band C-phases is the one of the clean transverse-field Ising chain for sufficiently weak λ A,max (as previously discussed).In addition, there is no associated Griffiths phase since there are no RRs (λ A,max < λ B = λ C ) as sketched in Fig. 2(a).However, when approaching the multicritical point, RRs in the Aphase appear and enhance the low-energy density of states.As a result, the gap closes around the transition.At criticality, those A-RRs can even provide a larger dynamical exponent z.In that case, the clean critical point is replaced by a line of finite-disorder critical points [dotted boundary line in Fig. 2(a)].Finally, at the multicritical point, the approximation of weak A-couplings completely breaks down and the most general theory contains a random-mass term.Therefore, this critical point is of infinite-randomness type.

IV. THE STRONG-DISORDER RENORMALIZATION-GROUP METHOD
In this section, we develop a strong-disorder renormalization-group (SDRG) method suitable for studying the long-distance physics of the Hamiltonian (1) for p = 1 and 2 and with random coupling constants.It is an energy-based RG method where strongly coupled degrees of freedom are locally decimated out hierarchically.Namely, we search for the strongest coupled local degrees of freedom and freeze them in their local ground state.The couplings between the remaining degrees of freedom are renormalized perturbatively.This procedure becomes more and more accurate if the local energy scales become more and more disordered (broadly distributed).In that case, the perturbative renormalization procedure becomes more accurate after each RG decimation step.This method was originally devised to conventional spin-1/2 models [74][75][76] and later on generalized to many other models.For a review, see Refs.66, 67.
A. Case p = 1

The decimation procedure
To start, let us consider the case p = 1 in which the model Hamiltonian (1) simplifies to For simplicity, we disregard the boundary conditions.Although ( 14) and the random transverse-field Ising chain are distinct, they share (apart of global degeneracies) the same eigenenergies.
Following the SDRG philosophy, we search for the largest local energy scale Ω = max λ j , say |λ 2 |.We then project the Hamiltonian onto the low-energy sector of As a result, the projection can be interpreted as a replacement of spins σ 2 and σ 3 by an effective spin-1/2 degree of freedom σ which is defined by σ In second order of perturbation theory, we find that and h = σ x 1 σ z σ z 4 .For our purposes, the constant term is harmless and can be disregarded.
Notice that h is now a three-spin interaction.However, since σ appears only in h, it is a local gauge variable whose role is simply to double the degeneracy of the spectrum.The SDRG decimation procedure ( 15) can be straightforwardly generalized to operators involving an arbitrary number of "internal" degrees of freedom since the algebra (3) is preserved.
Alternatively, the additional degeneracy induced by the effective internal spin σ can be interpreted as if the renormalized chain is, actually, two decoupled new chains.In the first one, σ is fixed in state | +⟩ (the original spins 2 and 3 fixed in the ground state |+⟩ of H 0 ) and the corresponding renormalized Hamiltonian is simply H1 = − λ σ x 1 σ z 4 , while in the second one σ is fixed at | −⟩ (the original spins frozen in the state |−⟩) and H1 = + λ σ x 1 σ z 4 .These two chains have the same spectrum and the subsequent SDRG decimations are identical (apart from the signs of the renormalized couplings which, for our purposes, are not important).
Thus, h can be simplified back to a two-spin interaction at the expense of dealing with two "twins" renormalized chains.The only difference between them being the sign of the renormalized coupling constant λ .
To use this simplification and keep track of the degeneracies, we are then required to introduce the quantity g Ω .It measures the total number of gauge (extra) spin-1/2 degrees of freedom at the energy scale Ω.Clearly, after each decimation it renormalizes to with the initial condition g Ω 0 = 0.The total number of effective degrees of freedom in the chain N Ω renormalizes to with the initial condition Figure 3. Decimation scheme for the Hamiltonian (1) in the p = 1 case.In (a), the decimation is sketched in the real space with points and lines representing spin sites and coupling constants, respectively.In (b), the decimation is sketched in the Hamiltonian space with circles representing the local energy operators and the lines connecting anticommuting operators.
In Fig. 3(a) we sketch the decimation procedure ( 15)- (17).Regarding the local energy scales, the decimation procedure ( 15) is identical to that of the random spin-1/2 XX chain [33,90] and that of the random transverse-field Ising chain [29].This is not a surprise since the free-particle spectra of all these models are the same.
Finally, it is instructive to recast the decimation procedure in the Hamiltonian space as shown in Fig. 3(b).The jth circle represents the local energy operator h j .A line connecting different circles means that the sharing operators anticommute with each other.Disconnected operators act on different Hilbert spaces and, thus, trivially commute with each other.In the decimation procedure, h 2 and the "neighboring" operators are replaced by h which anticommutes with the neighboring operators h 0 and h 4 .The algebra structure is, thus, preserved along the SDRG flow.

The SDRG flow
Since the SDRG decimation rule (15) is the same as that for the spin-1/2 XX chain, the renormalization-group flow of the coupling constants is already known [29,33,90] with σ 2 x = x 2 − x 2 being the variance of x.For δ ≫ 1, the SDRG flow is towards a stable fixed point in which only the odd couplings are decimated.This implies that only the even couplings are renormalized and, thus, are much smaller than the odd ones.This corresponds to a phase in which |⟨h 2i−1 ⟩| > |⟨h 2i ⟩|.In the spin-1/2 XX chain, this corresponds to the odd-dimer phase where spin singlets are formed over the odd bonds, i.e., |⟨S 2i−1 • S 2i ⟩| > |⟨S 2i • S 2i+1 ⟩|.The correspondence of this phase in the transverse-field Ising chain is not so straightforward since the XX spin-1/2 chains maps into two independent random transverse-field Ising chains [33].In the first one, the odd couplings of our model play the role of the transverse fields of the Ising chain.In the second, these roles are exchanged.Thus, the phase |⟨h 2i−1 ⟩| > |⟨h 2i ⟩| corresponds to the paramagnetic (ferromagnetic) phase in the first (second) quantum Ising chain.
If 0 < δ ≪ 1, the system is in the associated Griffiths phase.Typically, |⟨h odd ⟩| > |⟨h even ⟩|, but there are some "defects" inside which |⟨h odd ⟩| < |⟨h even ⟩|.These defects form the rare regions discussed in Sec III.Surrounding a rare region, there are two (domain-wall) spins weakly coupled.As a result of their weak coupling, the typical and mean values of the finitesize gap vanish ∼ L −z [Eq.( 9)], which defines an off-critical dynamical (Griffiths) exponent z.As the critical point is approached, this exponent diverges as z ∼ |δ | −1 .
We now further discuss on the effects of the RRs in the Griffiths phase through the lenses of the SDRG method.Recall that, in the transverse-field Ising chain, the origin of the gapless modes in, say, the paramagnetic phase is due the RRs which are locally in the ferromagnetic phase and fluctuate between the two ferromagnetic states.As this is a coherent tunneling process involving many spins, the associated relaxation time increasing exponentially with the RR's volume.In the XX spin-1/2 chain, these RRs correspond to patches which are locally in the even-dimer phase while the bulk is in the odddimer phase.The two domain walls delimiting a RR are, in zeroth order of approximation, simply free (unpaired) spins.To lowest nonvanishing order in perturbation theory, these spins (say, at sites 1 and ℓ) actually interact via an effective coupling constant equal to λℓ = λ 1 λ 3 . . .λ ℓ−1 /λ 2 λ 4 . . .λ ℓ−2 [which could also be obtained by a successive iteration of Eq. ( 15)].Thus, the gap of a finite chain is simply the excitation energy of the weakest coupled spins in these domain walls: min{| λℓ |}.Notice that, as expected, λℓ vanishes exponentially with the RR size implying an exponentially large relaxation time.An analogous physical picture appears in our model.For simplicity, consider a compact RR in which the local even couplings λ 2 , λ 4 , . . ., λ ℓ−2 (ℓ even) are greater than the local odd couplings λ 1 , λ 3 , . . ., λ ℓ−1 .In that case, after decimating the even operators h 2i in that RR, an effective operator linking spins 1 and ℓ appear h1 and a longer correlation between spins 1 and ℓ develop, σ x 1 σ z ℓ ̸ = 0. Consequently, a low-energy mode arises with excitation energy of order λ .Evidently, by duality, there are analogous conventional and Griffiths phases for δ < 0.
At criticality δ = 0, the fixed point is universal in the sense that critical exponents do not depend on the details of the disorder distributions.For instance, the finite-size gap distribution for sufficiently large systems is [91] where ln λ even , Ω 0 is the maximum value of λ in the bare system, ∆ is the finite-size gap, and ψ = 1/2 is the universal tunneling exponent.The distribution P SDRG is L independent for L ≫ γ −1 D = π/8σ 2 ln λ , the inverse of the Lyapunov exponent [92] which plays the role of a clean-dirty crossover length [93,94].The relation between length and energy scales follow from the scaling variable in Eq. ( 20), from which follows the activated dynamical scaling ln ∆ ∼ −L ψ in Eq. ( 10).

Thermodynamics
The thermodynamic observables follow straightforwardly.For instance, the low-temperature entropy is simply S ∼ 1 L (N T + g T ) ln 2, which simply counts the total number of active (undecimated) spins at the energy scale Ω = T .The reasoning is the following [76].At low temperatures T ≪ Ω 0 , the distribution of effective coupling constants is singular and, thus, the majority of the couplings are much smaller than the maximum energy scale Ω = T .In sum, the active spins are essentially free.At criticality (δ = 0), [95,96] and Notice the residual zero-temperature entropy coming from the exponentially large ground-state degeneracy.
The specific heat follows from C = T ∂ S/∂ T .In the lowtemperature limit, which is similar to that of the random spin-1/2 XX chain.

Ground-state degeneracy and the spectrum of other models
It is interesting to further explore the connection between the p = 1 model, the spin-1/2 XX chain, and the transversefield Ising chain in view of the SDRG decimation procedure.
Within the SDRG framework, one can obtain the whole spectrum of the transverse-field Ising chain in the following way [97].When performing the decimating procedure, one can either search for the ground state (and then project onto ground state of the local Hamiltonian) or, alternatively, search for the excited states (and then project onto the excited state of the local Hamiltonian).If one projects onto the excited states, the effective coupling of field pics up a different sign but the magnitude is the same [see Eq. ( 15)].However, for decimation purposes, all remains the same as the sign is irrelevant.Performing the decimation for all possibilities of low and excited states, one constructs the entire spectrum of the Ising chain: 2 L/2 states in total.(Recall that the associated Ising chain has half of the sites, but the same number of operators in the Hamiltonian.)Notice that this is equivalent to consider two "twins" renormalized chains as previously outlined.
Thus, all ground states of the p = 1 model are equivalent to the entire spectrum of the transverse-field Ising chain.
The connection to the XX model is alike.Here, the local Hamiltonian has 3 energy levels: one corresponds to the spin-0 singlet state, another to the zero-magnetization spin-1 triplet state, and the remaining one is doubly degenerate corresponding to the ±1-magnetization spin-1 states.If one disregards the doubly degenerate ±1-magnetization states, the decimation procedure recovers that of the Ising chain.Thus, the many ground-states of the p = 1 model corresponds to a small fraction of the states in the spectrum of the XX spin-1/2 chain.
We now inquire about the spectrum of the p = 1 chain.When projecting the system onto to local excited states, the only difference with respect to projecting onto the local ground state is a sign picked up by the effective energy scales and a flip of one of the spins.Thus, the entire spectrum can be easily related to the states in the ground-state manifold.There will be 2 L 2 states each of which is 2 L 2 degenerate.Precisely, all states can be represented by where φ j is either Here, { j 1 , j 2 } is the pair of spin sites decimated together in the jth decimation which is the same pair for all states.The precise state φ j (if φ j,+ or φ j,− and with the + or − sign) depends on the sign of coupling constant and on whether the projection was made into the local ground or excited states.All of these, are determined by the history of decimation procedure.

Spin-spin correlations
In the SDRG approach, any of the 2 L/2 ground states of H 0 is a simply product state as specified in (23).In that case, the spins in the jth pair are strongly correlated and dominates the average value of the spin-spin correlation.The probability that a pair of length r = j 2 − j 1 is formed along the SDRG flow is proportional to r −2 for sufficiently large r [33,96].Hence, the mean correlation function decays only algebraically, with universal exponent η = 2. (Here, • • • denotes the disorder average.)The typical value of the correlation, on the other hand, decays stretched exponentially fast [33], The distribution of the values of the correlation function is also known.For that, we refer the reader to Refs.18 and 94.
The off-critical (δ ̸ = 0) correlations are also known [29,33,85].They decay exponentially faster ∼ e −r/ξ with a diverging correlation length where ν = 2 for the mean correlations and ν = 1 for the typical correlations.Interestingly, all those results apply for any state as well provided that only the magnitude of the correlations are concerned.

The usual SDRG method
We now derive the SDRG decimation rules for the Hamiltonian (1) with p = 2.
Following the usual SDRG receipt, we treat ) as a perturbation to H 0 = −λ 5 h 5 .(We are assuming that the energy cutoff Ω = |λ 5 |.) Up to second order in perturbation theory, the renormalized Hamiltonian is H1 = P 0 H 2 1 P 0 / (−2 |λ 5 |), where P 0 is the projector onto the ground state of H 0 .The direct terms in the square of H 1 are unimportant constants since h 2 i = 1.The cross terms are proportional to the anticommutators h i , h j .Thus, many of those vanish because of the algebra (3).The surviving ones are those which commute with each other.Thus, H1 simplifies to H1 = − λC hC − λAC hAC − λA hA where the operators are The decimation procedure ( 27) and ( 28) is represented in the Hamiltonian space in Fig. 4(a).Unlike the p = 1 case [see Fig. 3(b)], the structure of the anticommuting algebra (3) of the original Hamiltonian is not preserved and new operators (involving σ y ) appear.

The block SDRG method
The arising of the new operator hAC in (27) hinders the practical implementation of the method.First, it prevents us from projecting the renormalized system onto the ground states of the new effective spin operators σ and τ as in the p = 1 case.Second, and more importantly, it requires a generalization of the SDRG procedure to take into account these new operators.While this is possible and cumbersome, we found, surprisingly, a quite simpler route guided by the algebra of the energy-density operators (3).We generalize the "usual" SDRG approach above reported to (in the lack of a better therminology) a "block" SDRG approach.In the latter, we consider a larger unperturbed Hamiltonian (and thus, a larger Hilbert space) when performing the decimation procedure (see details in Appendix A).The size of the block is the maximum number of operators which anticommute among themselves.
When decimating the largest local energy scale Ω = |λ 5 |, instead of considering H 0 = −λ 5 h 5 (which is a B-type operator), we consider a larger block involving the Aand C-type "nearest-neighbor" operators, i.e., H 0 = −λ 4 h 4 −λ 5 h 5 −λ 6 h 6 .This is the largest block which encompass h 5 and still have only two energy levels (as h 5 ), 4 i.e., the eigenenergies of H 0 are ± λ 2 4 + λ 2 5 + λ 2 6 .Then, we project the H − H 0 on the ground-state subspace of H 0 .The degeneracy of the ground state is 2 4 and, thus, can be spanned by four effective spin-1/2 degrees of freedom σ a , σ b , σ c , and σ d .In practice, we have to project λ 2 h 2 , λ 3 h 3 , λ 7 h 7 and λ 8 h 8 .The result is that, in the regime |λ 5 | ≫ λ 4,6 , the renormalized operators are where Surprisingly, the renormalized operators are different in character from those in Eq. (27).Interestingly, they preserve the algebra structure (3) of the original system as depicted in Fig. 4(b) if we identify h2 → h 2 , h3 → hC , h7 → hA , and h8 → h 8 .Furthermore, the "hybrid" operator hAC [generated from the usual SDRG approach, see Fig. 4(a)] is not generated in the block SDRG approach [see Fig. 4(b)].Instead, there are only "pure" operators except for the BCtype operator in (29) and the AB-type operator in (31).This is very convenient because they can be neglected at strongdisorder fixed points (corresponding to situations near and at the dashed transitions in Fig. 2).The reasoning is the following.Since the effective disorder at and near the transition is very large ( where The decimation procedure (34) and ( 35) is schematically depicted in Fig. 4(b).By symmetry, an analogous decimation procedure is obtained in the case |λ 4 | ≪ |λ 6 |, with the exchanges λ 2 ⇋ λ 8 and λA ⇋ λC , which is obtained after a convenient change in the definition of the effective operators σ a , σ b , σ c , and σ d .The decimation procedure (34) and (35) [see Fig. 4(b)] is very convenient.It preserves the algebra structure (3) and the operators of the original Hamiltonian hi = σ x i σ z i+1 σ z i+2 .This procedure is a straightforward generalization of the decimation procedure of the p = 1 case in the following sense.For p = 1, a decimation of an A-type coupling implies the renormalization of the neighboring B-type couplings and vice versa [see Eq. ( 15)].For p = 2, due to triality, the decimation of a Btype operator implies the renormalization of the neighboring Aand C-type couplings [see Eq. ( 35)].

On the equivalence between the usual and the block SDRG approaches
In the regime |λ 5 | ≫ |λ 4,3 |, there should be no difference between the usual and block SDRG approaches as the ground state of the block H 0 = −λ 4 h 4 − λ 5 h 5 − λ 6 h 6 is simply that of the local H 0 = −λ 5 h 5 furnished with trivial degeneracies.It is somewhat surprising that seemingly fundamentally different decimation procedures arise from these approaches.Thus, we inquire whether these two decimating procedures are really different.
In Ref. 39 the SDRG method was carried out in the antiferromagnetic Heisenberg model in a zigzag and two-leg ladder geometries.In both cases, it was shown that, in the early stages of the SDRG flow, further neighbors interactions arise, just like what was found in the usual SDRG decimation [see Fig. 4(a)].However, in the final stages of the SDRG flow, the renormalized geometry of the system always converged to a chain geometry.(Evidently, one had to disregard exceedingly small couplings that were generated along the flow.)This may be a general feature of quasi-one-dimensional critical or nearcritical systems.The low-energy long-wavelength effective theory is that of a chain.Technically, couplings of "long" operators that connect exceedingly distant spins arise after many renormalizations and, thus, are typically much smaller than those of "shorter" operators.
In sum, the longer range character of the new hybrid operator hAC in the usual SDRG approach may be an irrelevant "operator" which vanishes in the latter stages of the SDRG flow.In the block SDRG approach, the new hybrid operators hAB and hBC could be clearly determined as irrelevant ones near and at criticality, where the flow is towards strong disorder, a self-consistent assumption that we have yet to prove.Therefore, it is plausible that hAC in the usual SDRG approach can be neglected by the same reasons.In that case, both the usual and block SDRG approaches are equivalent near and at criticality.This is a fundamental feature of the renormalizationgroup philosophy.Details on defining the coarse-grained operators should not matter in the long-wavelength regime.

SDRG flow corresponding to conventional phases
Having derivied the SDRG decimation rules, we now analyze the flow of the coupling constants.
Let us start by analyzing the simpler case corresponding to flow towards the gapped phases.In this case, one of the coupling constants, say, of C-type, is always greater than the others.Precisely, min {λ 3i } > max {λ 3i−1 , λ 3i−2 }.In that case, notice that only the bare h 3i operators are decimated.Therefore, the corresponding fixed point is noncritical and, thus, represents a phase: the h C > h A,B conventional phase in Fig. 2. Notice that both the usual and the block SDRG can be applied here since only the original h 3i operators are decimated.Thus, in the SDRG framework, the spectrum gap is simply ∆ = 2 min {λ 3i }.
Evidently, by triality, there are additional two phases in which h A > h B,C and h B > h A,C as shown in Fig. 2. Finally, we call attention to the fact that, in the SDRG framework, these phases exist because a decimation of h i , renormalizes the neighboring interactions h j with | j − i| ≤ p = 2. Thus, for a generic value of p in the Hamiltonian (1), we expect, at least, p + 1 different phases.First, we analyze the case in which at least two of the three types of couplings are random, i.e., the SDRG flow associated with the transitions and Griffiths phases shown in Fig. 2(b) away from the multicritical point.
Let us consider the case in which max {λ 3i−2 } < min {λ 3i−1 , λ 3i }.In this case, the SDRG flow is dictated only by the competition between {λ 3i } and {λ 3i−1 } as the A-type couplings will never be decimated.Thus, we can simply drop the λ A 's in the decimation procedure (34) and ( 35) as the distribution P A renormalizes to an extremely singular one.In that case, the zigzag-chain geometry of the renormalized system becomes that of a single chain.Therefore, the flow is simply that of the case p = 1.
To analyze it, we simply generalize the distance from criticality Eq. ( 18) to From here, all energy-related results from the p = 1 case follows straightforwardly.The transition is governed by an infinite-randomness critical point and happens for δ BC = 0.The gap distribution is that of Eq. ( 19) with the scaling variable (20) redefined as This redefinition is due to the fact that the total number of decimations in the p = 1 case is L/2 while it is L/3 for the p = 2 case.Thus, L/2 in the p = 1 case translates to L/3 in the p = 2 case.In addition, there are associated Griffiths phases for |δ BC | ≪ 1.The off-critical dynamical exponent z diverges as z ∼ |δ BC | −1 when it approaches criticality.The nature of the low-energy modes are also associated with domain wall spins surrounding a rare region, just like for p = 1.
By triality, there are other two boundary transitions for δ AB = 0 and δ AC = 0.These boundaries and the Griffiths phases are shown in Fig. 2 as dashed lines and shaded regions, respectively.The SDRG results here reported, thus, put the heuristic arguments of Sec.III in solid grounds.
Now we analyze the case in which the two competing (strongest) couplings are uniform, i.e., λ B,i = λ B , λ C,i = λ C , and λ A,i random with λ A,typ < λ B,C .This corresponds to the region in the phase diagram of Fig. 2(a) surrounding the clean and finite-randomness transition but sufficiently far from the multicritical point.
When max {λ A,i } < λ B,C , the SDRG method does not need to be applied since weak λ A coupling is irrelevant.The transition at λ B = λ C is in the clean Ising universality class, and the spectral gap vanishes only at criticality.When λ A,typ < λ B < λ C < max {λ A,i }, the system in the Griffiths C-phase with high A-susceptibility as in the generic case discussed above.(Analogously for λ A,typ < λ C < λ B < max {λ A,i }.) Finally, we discuss the interesting case when λ A,typ < λ B = λ C < max {λ A,i }.The system is globally critical between Band C-phases but there are rare regions locally in the A-phase.What are their effects?Applying the block-SDRG decimation procedure, we simply decimate A-operators in the first stages of the flow.After decimating all A-operators with λ B = λ C < λ A,i < max {λ A,i }, we can simply ignore the remaining A-operators since they would be surrounded by locally stronger Band C-operators.The effective chain is, thus, that with Band C-operators only.However, the effective coupling constants are not uniform.Interestingly, notice that the effective couplings obey the condition λB,i = λC,i .The precise condition for the absence of random distances from criticality (random mass) discussed in Sec.III A 2. When the A-rare regions are small, the effective couplings λB,C,i are weakly renormalized and their effective disorder (variance of their distribution) is weak.As a result, the clean critical behavior is stable.However, when approaching the multicritical point we expect the renormalization of λB,C,i to become more and more relevant.This means that their effective disorder increases to the point where the clean critical behavior is destabilized.The resulting phase transition is thus of finite-randomness type with nonuniversal critical exponent z.This result should be valid up to the multicritical point where z reaches its maximum value.Interesting, z is formally infinity in the other transition lines meeting at the multicritical point.
We end this section by noticing that the above results straightforwardly generalizes to any value of p.

The SDRG flow at the multicritical point
Finally, we deal with the case in which all couplings compete.
At first glance, one could think in analyzing the flow equations for the distributions P A , P B , and P C analytically.The flow equation for P A is (see Appendix B) where P X = P X (λ ; Ω), P X (Ω) = P X (Ω; Ω), and (39) By triality, the flow equations for P B and P C follow from exchanging the labels accordingly.The last term on the r.h.s. of (38) implements the renormalization of the A-type coupling [ λA in Eq. ( 35), and we are considering only the magnitude of the coupling constants] when a Bor C-type coupling is decimated.The remaining terms ensure that P A remains normalized when the cutoff energy scale Ω is changed.As shown in App.B, the critical point (where P A = P B = P C ) of the flow equation ( 38) is of infinite-randomness type [see Eq. (B5)] and is related to a different universal tunneling exponent ψ = 2 3 which is greater than that ψ = 1 2 of the case p = 1.A larger tunneling exponent means larger effective disorder as the typical value of the finite-size gap is even smaller [see Eq. ( 37)].
This sounds intuitive since two coupling constants are renormalized in each decimation instead of one renormalization per decimation as happens in the p = 1 case.
Although this gives us a prescription to find a new universality class beyond that of the permutation-symmetric universality class (where ψ = 1/N, with N being an integer [37,40,46,47,98]), the analytical analysis of this problem is not correct.The flow equation (38) assumes no correlation between the three types of couplings.However, when, say, an B-type coupling is decimated the renormalized Aand C-type couplings are added as neighbors in the renormalized system [see Fig. 4(b)].Being neighbors diminishes their chance of being further renormalized (which would increase the effective disorder strength by producing even more singular couplings) when compared with the case described by Eq. (38) where they are inserted in the chain in an uncorrelated fashion.
As treating the correlated case analytically is not simple, we then proceed our study numerically.We implement the block SDRG rules (34) and (35) for a system where all the coupling constants are independent random variable and identically distributed according to where Ω 0 is the initial energy cutoff and D 0 parametrizes the disorder strength of the bare system.
We first study the moments of the distribution of the renormalized couplings along the SDRG flow.As usual in the SDRG approach, it is convenient to define the logarithmic couplings ζ i = ln (Ω/λ i ) and the logarithmic energy cutoff Γ = ln (Ω 0 /Ω).We then study the mean value ζ and stan- as a function of Γ.Our results for D 0 = 1/2 and D 0 = 1 and system size of L = 3 15 spins is shown in Fig. 5(a).Clearly, ζ ≈ σ ζ ≈ Γ + const in the Γ → ∞ limit.This implies an infinite-disorder fixed-point critical distribution.Thus, the SDRG method here employed is justified.
Next, we study the infinite-disorder fixed-point distribution.For comparison, that distribution for p = 1 is well known [33,96,99].It is with  Energy and length scales The relation between energy and length scales along the SDRG flow is shown in Fig. 5(a).Simply the length scale ρ −1 ∼ Γ 1/ψ with universal tunneling exponent ψ = 1/2.We also expect the finite-size gap ∆ to obey the activated dynamical scaling (37).This is confirmed in our numerics as shown in Fig. 5(c).Both the mean value and the width of the distribu-tion of ln (∆) behave similarly ∼ L ψ .Here, the finite-size gap is obtained by decimating the entire chain.The last decimated coupling constant λfinal provides the finite-size gap ∆ = 2 λfinal .

Thermodynamics and correlations
As the fixed-point distribution is of infinite-randomness type, the critical thermal entropy and specific heat behave similarly as in the p = 1 case discussed in Sec.IV A 3. The only difference is the residual entropy which modifies Eq. ( 21) to the reasoning being that the ground state have degeneracy 2 2L/3 .The specific heat follows straightforwardly and recovers (22).
In the SDRG framework, the ground-state spin correlation C i, j,k = σ x i σ z j σ z k is ±1 if the spins i < j < k are decimated together, and weakly vanishing (we set to zero) otherwise.To obtain the average value of C i, j,k , we then build a normalized histogram |C(r 1 , r 2 )| in log-scale.After each SDRG decimation (decimation of spins i, j, and k), we add a unity to the bin (r 1 , r 2 ) if ln ( j − i) is between ln r 1 − 1 2 and ln r 1 + 1 2 , and ln (k − j) is between ln r 2 − 1 2 and ln r 2 + 1 2 .This procedure is repeated until the entire chain is decimated and averaged over many disorder realizations.The histogram |C(r 1 , r 2 )| is then a proxy to average value of C i, j,k , i.e., We plot in Fig. 6(a) the mean value of the spin correlation C (r 1 , r 2 ) as a function of the internal distances r 1 and r 2 .Clearly, the correlation is maximum when r 1 = r 2 , as expected from symmetry.
In Fig. 6(c), we plot the average correlation C (r 1 , r 2 ) for various values of r 2 .Interestingly, our numerical data indicates that C (r 1 , r 2 ) ∼ r −(ϕ−φ ) 1 for fixed r 2 and r 1 ≫ r 2 .Analyzing these exponents for other sample sizes, we estimate their error to be of order 10%.Unfortunately, we do not have an analytical derivation for these exponents.
How about the typical value of these correlations?Typically, the triad of spins in (42) are not decimated together and, thus, develop quite weak correlations.As argued by Fisher [33] and shown in many numerical works [18,94,100] for the case p = 1, this weak correlations are of order of the typical value of the coupling constants involved, C typ (r) ∼ J typ (r).Very plausible, this is also the case for any p with the caveat that more than one coupling constant is involved.Thus, C typ ∼ J typ (r 1 )J typ (r 2 ) . . .J typ (r p ).Thus, from the dynamical scaling (10) we then conclude that the typical value of the correlations decays stretched exponentially fast with the spin separations.Roughly, we expect where ℓ(r 1 , r 2 ) is a function which equals max {r 1 , r 2 } for r 1 ≫ r 2 or r 2 ≫ r 1 and ≈ cr when r 1 ≈ r 2 ≈ r (with c being disorder-independent constant of order unity), γ D is a disorder-dependent Lyapunov exponent [18,94] related to the clean-dirty crossover length [93], and A is a constant of order unity.

Generalization to p > 2
The usual SDRG method of Subsec.IV B 1 can be easily generalized to any p.The key point is that new effective operators h generated under the SDRG flow either commute or anticommutate and are of two-level type h2 = 1.In this case, the only type of decimation that can take place is of second order in perturbation theory, and, thus, the renormalized coupling constants have the multiplicative structure of Eq. ( 28).As a consequence, the relation between length and timescales at a phase transition can only be of activated type Eq. ( 10), meaning that, at sufficiently strong disorder where the SDRG is applicable, the phase transition is governed by a infiniterandomness fixed point.
The analysis of the block SDRG flow of Subsec.IV B 1 allowed us to conclude that this infinite-randomness critical point has tunneling exponent ψ = 1 2 and very different mean and typical values of the correlation function for the cases p = 1 and 2.Although the generalization of the block SDRG method to other values of p is not straightforward (as it involves many nontrivial projections, see Appendix A), it is plausible to conjecture that the phase transitions for uncorrelated disorder are governed by an analogous infiniterandomness critical point for any p.The tunneling exponent is universal and equals ψ = 1 2 .It not only governs the dynamical scaling but also the typical correlations.On the other hand, the exponent for the mean correlation function is universal, i.e., disorder independent, but does depend on the value of p.
A detailed and quantitative investigation of the critical behavior and of the associated Griffiths phases is out of the scope of the present work and is left for future research.

V. FINITE-SIZE GAP STATISTICS
In this section, we present our numerically exact results on the spectral gap of the model Hamiltonian (1) for p = 2.This quantity is computed from the roots of the polynomial as described in Sec.II and compared with the SDRG predictions described in Sec.IV B.

A. Technical details
We refer the reader to Ref. 12 for useful numerical methods for calculating the roots of the polynomial ( 5)- (6).We emphasize that our method allow us to evaluate the finite-size gap in lattice sizes up to ∼ 10 7 sites in the neighborhood of the transition lines with a numerical cost that grows linearly with the system size competing with the SDRG numerical cost.For systems that large, typically, the polynomials (5) have coefficients spanning in 500 orders of magnitude.Therefore, the entire evaluation of these coefficients can only be achieved using high numerical precision (500 digits).However, as shown in Ref. 12, we only need the last 100 coefficients of the polynomial to obtain the finite-size gap with quadruple standard numerical precision (32 digits).When applying this method to our Hamiltonian, we found useful to keep the last coefficient of (5) always of order unity.This is accomplished by factorizing these coefficients when iterating the recursion relation (6).With that, the entire procedure can be accomplished using only standard routines in FORTRAN with quadruple precision.

B. Off-critical finite-size gap I
We start analyzing the off-critical Griffiths phases sketched in the phase diagram of Fig. 2(b).For such, we study three sets of chains: In all cases, λ B is distributed uniformly in the interval [0.5, 1.5], and λ C is a constant (running from 1.2 down to 0.9) which serves as a tuning parameter.The running λ C passes through the critical point λ * = λ B,typ = e 3 ln 3−2 ln 2 2 −1 ≈ 0.95578 [see Eq. ( 36)] but does not include it since the analysis of the critical system is reported in other sections.Notice that set III of chains is critical for λ C < λ * C .For that reason, we consider only λ C > λ * C in this section.In case I, we explore the interplay between the two strongest couplings λ B and λ C while λ A is a much weaker coupling.The Rare Regions are of B-type inside a bulk in the C phase.In case II, the value of λ A is increased such that some Rare Regions of A-type also appear.Finally, in case III, both Rare Regions of Aand B-type appear equally.At the final point λ C = λ * C , the multicritical point is reached.In Fig. 7(a) we show the relation between ln ∆ (with ∆ being the system finite-size gap) and the system size L for λ C = 1.From this, we can obtain the off-critical dynamical exponent z by fitting ∆ ∼ L −z to the data.Repeating this procedure for other values of λ C , we determine how z diverges as the critical point is approached, see Fig. 7(b).It diverges as z ≈ 1/ (2δ ), with δ being the distance from criticality as defined in (36).We emphasize that this numerically exact result is in agreement with the SDRG analytical predictions in the z → ∞ limit.

C. Off-critical finite-size gap II
We now study the case in which only one of the couplings is disordered, say, λ A uniformly distributed in the interval [0, 0.2].The corresponding phase diagram is shown in  The coupling constants λ A,i are uniformly distributed in the interval [0, 0.2] (set I), [0.2, 1.2] (set II) and [0.5, 1.5] (set III).In all cases, λ B,i are uniformly distributed in the interval [0.5, 1.5], and λ C is uniform and plays the role of a running parameter throughout the Griffiths phases.(a) The typical finite-size gap ln ∆ as a function of the system size L for λ C = 1.The data is averaged over N samples = 6 000 disorder realizations and the error bars are about the size of the data points.Linear fits for the data are provided in the legends from which the dynamical exponent z is obtained.The best fits are restricted to system sizes L > 10 6 .(b) The dynamical exponent z as a function of the distance from criticality Eq. (36).Lines are guide to the eyes.Inset: same data in the main panel in log-log scale.The dashed line is the SDRG asymptotic behavior (δ ≪ 1) z = 1/(2δ ).
Fig. 2(a).Here, we focus on the off-critical region near the transition between the Band C-phases.Notice that there are no Griffiths phases surrounding the transition for sufficiently small λ A,typ .Thus, the gap is vanishing only at the transition point λ B = λ C .Closer to the multicritical point, Griffiths phases appear.
We plot in Fig. 8 the typical value of the system gap ∆ as a function of the distance from criticality δ .Since λ C and λ B are homogeneous, the definition ( 18) is not useful because of the vanishing denominator.Here, we simply use δ = ln (λ C − λ B ). Also, we fix λ C = 1 > max {λ A } and use λ B as a running parameter.Notice that ∆ diminishes as in the clean system: ∆ ∼ δ .For small δ , ∆ saturates due to finite-size effects, meaning that the correlation length is greater than L. We recall that our method is not optimal for gapful systems.This is because larger the gap larger is the number of coefficients required in the characteristic polynomial [12].For that reason, we studied chains of "small" sizes up to 10 267 sites.
Having shown the absence of Griffiths phases far from the multicritical point (λ C > max {λ A }), we now show their existence otherwise.We thus study chains in which λ C = 1, λ A is uniformly distributed in the interval 0, eλ A,typ , and λ B is the tuning parameter.We report that the finite-size gap vanishes as ∆ ∼ L −z [as in Fig. 7(a)], with the value of the offcritical dynamical exponent z increasing up to a finite value as the transition is approached, see Fig. 9. Notice that far from the critical point λ B = 1, the off-critical dynamical exponent is practically insensitive to the distance from criticality δ = 1 − λ B as the low-energy behavior is dominated by Rare Regions which are locally in the A phase.In Fig. 10 we plot the typical value of the finitesize gap [and the corresponding Laguerre bound ∆ LB , see Eqs. ( 7)- (10)] for those critical chains as a function of the system size L. As expected from the activated dynamical scaling (10), the finite-size gap vanishes stretched exponentially fast with universal (i.e., disorder-independent) tunneling exponent ψ = 1/2, compatible with the asymptotic behavior (L ≫ 1) of our data.We call attention to the fact that ∆ and ∆ LB become virtual indistinguishable for L ≫ 1.This feature was already noticed for the case p = 1 [12].We conjecture that finite-size values of ∆ LB → ∆ for L → ∞ in the case of infinite-randomness criticality.As argued in Ref. 12, this is because the largest root of the characteristic polynomial separates from the other ones as L increases.Thus, only the few largest coefficients of the polynomial are necessary to accurately compute it (see Sec. V A).

E. Critical finite-size gap statistics II
Here, we investigate the intriguing critical line in which the two major couplings are uniform, i.e., the horizontal boundary line in the phase diagram Fig. 2(a).Thus, we take λ B = λ C = 1 and λ A uniformly distributed between 0 and eλ A,typ .The typical value of the A couplings, λ A,typ , is used as a tuning parameter and reaches the multicritical point when λ A,typ = 1.
Our results are shown in Fig. (11).Clearly, sufficiently far from the multicritical point (λ A,typ ≤ λ * A,typ ∼ 0.6) the critical point is in the clean Ising universality class z = 1 [solid line in Fig. 2(a)].Approaching the multicritical point (λ A,typ > λ * A,typ ), a line of finite-disorder fixed points [dotted line in Fig. 2(a)] is tuned with varying critical dynamical exponent z that diverges as the multicritical point is approached (see inset).For comparison, we show the line z = (2γ) −1 as predicted by the SDRG method when γ = 1 − λ A,typ → 0. This seems to be the trend for γ 0.01.We cannot rule out that we are plagued by finite-size effect for smaller γ.

F. multicritical finite-size gap statistics
Here, we finally investigate the multicritical point.We find that in both phase diagrams of Fig. 2, the multicritical point is of infinite-randomness type with universal tunneling exponent ψ = 1/2.
In Fig. 12 we plot the distribution P of the finite-size gap ∆ properly rescaled according to Eq. (37) for various system lengths and disorder parameters.Clearly, the scaling variable η is sufficient to produce that data collapse for the system sizes used and this gives us further confidence on the activated dynamical scaling (10) with universal (disorderindependent) tunneling exponent ψ = 1/2.We notice, in addition, that the probability distribution P is not universal and Figure 12.Finite-size gap distribution P for different system sizes L. The finite-size gap ∆ is rescaled according to the scaling variable η in Eq. (37).The coupling constants are either homogeneous equal to e −1 or uniformly distributed in the interval [0, 1] as indicated in the legends.The width σ 0 is, respectively, 1/3, 2/3, and 1 in panels (a), (b), and (c).The normalized histograms were built using 105 to 10 6 disorder realizations.The dashed line is the SDRG prediction (19) for the transverse-field Ising chain (p = 1).is weakly dependent on the disorder details. 5 Furthermore, P is quite different from the SDRG prediction Eq. ( 19) for the transverse-field Ising chain (p = 1), even though they are governed by the same infinite-disorder fixed point.Finally, we report (not shown for clarity) that this distribution depends on the modularity of the lattice size L, i.e., P (η) depends on L mod 3.This is not a surprise since a similar difference also appears in the clean system.There, the finite-size gap amplitudes also depend on the modularity, i.e., ∆ ∼ aL −z with a ≡ a (L mod 3).

VI. CONCLUSIONS
We have studied the effects of quenched disorder on a family of free fermionic models (1) with (p + 1)-multispin interaction, paying special attention to the case p = 2 corresponding to the random version of the three-spin interacting Fendley model.
When all the coupling constants are generically disordered, the clean phase transitions (including the multicritical points) are destabilized.The replacing transition is of infiniterandomness character where the critical dynamical scaling is activated (10) with universal (disorder-and p-independent) tunneling exponent ψ = 1/2.This exponent governs the lowtemperature singular thermodynamic behavior.The average value of the correlation functions are also universal and decay algebraically with the spin separations.We do not have an analytical theory for those exponents and a detailed study is left for future research.The typical value of the correlations, on the other hand, behaves quite differently.It decays stretched exponentially fast with the spin-spin separation.In addition, strong Griffiths singularities surround the transitions.Although the system is noncritical with shortrange correlations, the spectral gap vanishes.The singularities are related to the slow dynamics of the domain walls surrounding the so-called rare regions.This phenomena is very similar to the domain-wall-induced (or rare-regions-induced) Griffiths singularities in the dimerized XXZ spin-1/2 chain and in the random transverse-field Ising model.
When only one (of few) type of coupling constants are disordered, the scenario is more involved.When the disordered coupling is weak, the singular critical behavior of the clean system is stable and there is no surrounding Griffiths phases.Upon increasing the magnitude of the disordered coupling, the universality of the transition changes and is of finiterandomness character.The critical scaling is power-law conventional (9) but with nonuniversal (disorder-dependent) dynamical critical exponent z.The increasing of the magnitude of the disordered couplings yields to another effect.It nucleates rare regions which contributes to off-critical Griffiths singularities.Interestingly, the finite-randomness criticality extends up to the multicritical point.
Although we have explicitly worked out those results for the cases of p = 1 and 2, from symmetry grounds we expect them to be valid for all p in the family model (1) for sufficiently large disorder.Evidently, we cannot exclude other scenarios appearing when p is very large and disorder is weak.
Finally, the agreement between our numerically exact results obtained for quite large lattice sizes (up to L ∼ 10 7 sites) with the "block" SDRG is remarkable.We stress that the SDRG methed is not an exact one.It took more than a decade after its conception to realize that it can provide asymptotically exact critical exponents and other universal quantities.
This was accomplished comparing the SDRG with exact diagonalization of a few models such as the transverse-field Ising chain and the XX spin-1/2 chain and for moderate lattice sizes (L ∼ 10 3 ) and the Heisenberg chain (L ∼ 200).Here, we have the rare opportunity to show that this is also true for a different family of models and for quite large chains.
We believe that our method for evaluating the finite-size gaps for large system sizes with minimal numerical cost (∼ L) will be a useful tool to study the effects of disorder in the phase transition of other systems.
Say that Ω = λ 2 A,2 + λ 2 B,2 .In this case, we treat H 0 = −λ A,2 h A,2 − λ B,2 h B,2 exactly and project the operators H 1 = −λ B,1 h B,1 − λ A,3 h A,3 onto its ground-state subspace.More precisely, we are interested in the projections of the operators σ z 3 and σ x 5 .The ground-state subspace of H 0 is spanned by |s , which, at first glance, does not seem to recover the usual SDRG result (15).This is not the case.One can correct the signs of the couplings by appropriately redefining the effective spin-1/2 degrees of freedom σ and τ.
In sum, the block SDRG approach here derived is similar to the usual SDRG derived in Sec.IV A. Qualitatively, they are the same in the sense that 2 original operators are removed from the system in each decimation step and that the renormalized couplings are smaller than the original ones (both approaches are self-consistent).The only difference is quantitative: the effective couplings are renormalized differently.However, this difference vanishes in the strong-disorder limit.

Case p = 2
The Hamiltonian of interest is  As in the simpler SDRG approach (Sec.IV B) the algebra (3) is not preserved.However, it is almost preserved.The only operators changing the algebra structure are the operators hB,1 and hB,3 .They anticommutate with the nearest-and next-nearest-neighbor operators and commute with the farther-neighbor operators, preserving the algebra (3).However, they do not commute with each other.Instead, λB,1 hB,1 , λB,3 hB,3 ∝ λ B,1 λ A,2 λ C,2 λ B,3 /E ) can be neglected.The reason is the following.Consider for instance the terms inside parenthesis in Eq. (A5), clearly, the term originating the AB-type operators (the first term) can be viewed as a small tilt to the molecular field of the B-type operator (the second term).In the regime |λ A,2 | ≫ |λ C,2 |, this small transverse molecular field can be neglected.This possibility of neglecting a hybrid operator in detrimental of a "pure" one does not appear in the usual SDRG.Maybe because the local Hilbert space (that of H 0 ) is not large enough to accommodate more than one possibility of renormalization.

Appendix B: Simplified SDRG flow
In this section, we consider the simplified version of the SDRG decimation procedure for the p = 2 case.As stated in Sec.IV B, the first decimation is such that five operators are removed h 1,2,...,5 and three new ones are inserted h1,2,3 in the system (see Fig. 4).If, for some reason, one could neglect h2 , the algebra structure would not change after decimation.This allows for a simplification of the problem.The new operator h1 ( h3 ) corresponds to a renormalization of the couplings on the sites 3i − 2 (3i − 1).We then can write an equation for the transformation of the coupling constant distributions.The transformation for the distribution of the couplings 3i − 2 when the cutoff Ω diminished to Ω − dΩ is is a functional quantifying the change in P Y when a X-type coupling constant is decimated.
Here, P X (Ω, Ω) P Y (λ 1 ) P Y (λ 4 ) dΩdλ 1 dλ 4 is the probability of having the decimation of an X-type coupling which involves the neighboring Y -type couplings λ 1 and λ 4 .We need to sum over all possibilities for the values of these neighboring couplings.The first two deltas correspond to the removal of these Y -type couplings.The last one corresponds to the addition of the renormalized coupling λ = λ 1 λ 4 Ω .The normalization constant is important to keep the distribution P A normalized after the cutoff is reduced.Integrating both sides of Eq. (B1) from λ = 0 to Ω − dΩ, the l.h.s. is simply N .Up to linear order in dΩ, the r.h.s is 1 − (P A (Ω, Ω) + P B (Ω, Ω) + P C (Ω, Ω)) dΩ.This is the expected result if one counts that a decimation of type B removes a net fraction of P B (Ω, Ω) dΩ couplings of type A. In addition, a decimation of A type removes a fraction of P A (Ω, Ω) dΩ couplings of type A. The beta function for the distribution P A simplifies to where P X (Ω) = P X (Ω, Ω), P X = P X (λ , Ω), P B∪C (Ω) = P B (Ω) + P C (Ω), and
and λ B < λ C = λ A .All of them are in the 2D Ising universality class, and, thus, the dynamical exponent is z = 1.The multicritical point λ A = λ B = λ C is, on the other hand, in a different universality class where z = 3 2 [6] and the specific-heat exponent α = 0 8.The purpose of present work is the study of the quenched disorder effects in the phase transitions of the model Hamiltonian (1) for p ≥ 2.

Figure 2 .
Figure 2. The phase diagram of the Hamiltonian (1) for p = 2.In panel (a), {λ 3i−2 } = λ A,i is a set of independent random variables while the remaining couplings λ 3i−1 = λ B and λ 3i = λ C are homogeneous.In panel (b), at least two of the coupling constants are independent random variables (see text).The solid red line is a transition line in the Ising universality class of the clean system [z = 1, see Eq. (9)].The red dashed lines are transitions in the infinite-randomness universality class [ψ = 1/2, see Eq. (10)].The multicritical point (in both cases) is also in the same infiniterandomness universality class.The blue dotted line in panel (a) is a transition line where the universality class is of finite-randomness type (1 < z < ∞).The phases have the same nature as the homogeneous case (see Fig.1) and the shaded regions delimits the associated Griffiths phases where the spectrum gap vanishes.

Figure 4 .
Figure 4. Strong-disorder RG decimation scheme for the Hamiltonian (1) in the p = 2 case.The decimation is depicted in the Hamiltonian space analogous to Fig.3(b)where circles represent the local energy operators and lines connect operators which anticommute.In (a), the SDRG method is implemented in its simpler form (usual SDRG) where a newly generated operator hAB disrupts the algebra (3).In (b), the SDRG method is implemented in a slightly more general fashion (block SDRG) which preserves the algebra (3) in the renormalized chain .

= σ x d σ z 9 σ z 10 .
which we show a posteriori), very likely either |λ 4 | ≫ |λ 6 | or the |λ 4 | ≪ |λ 6 |.In the former case, sin θ ≈ 0 in Eq. (30), and, thus, the hybrid type operators can be neglected.The renormalized B-type operators then simplify to h2 = −sign (λ 5 ) σ x Taking sin θ = 1 − cos θ = 0 and noticing that the new effective spin-1/2 degrees of freedom σ a and σ b appear only in combinations that commute with each other ( σ z a σ z b in h2 and σ y a σ y b in h3 ), we can project the resulting renormalized Hamiltonian in the common eigenstates of σ z a σ z b and σ y a σ y b .As a result, we have four "twins" renormalized systems.They are

5 .
SDRG flow corresponding to conventional Griffiths phases and phase transitions between two phases.

Fig. 5 (
b), we show the distribution of log couplings at different stages of the SDRG flow for the case D 0 = 1.(We find statistically identical result for D 0 = 1/2 which is not shown for clarity.)After the initial stages of the SDRG flow, the fixed point is reached.Here, the different stages of the SDRG flow is parametrized by the density ρ = N Ω /L, where N Ω is the number of active spins at the cutoff energy scale Ω. Recall that, after each decimation step of the block SDRG procedure, 3 spins are removed.Clearly, the fixed point distribution differs from the p = 1 case Eq. (40) (dashed line y = e −x ), but only slightly.We attribute this small difference to the correlations arising among renormalized couplings under the flow.

Figure 5 .
Figure 5. (a) The mean value of ζ , the standard deviation σ ζ and the density of active spins ρ as a function of the logarithmic cutoff energy Γ = ln Ω 0 /Ω (Ω = max{λ i } ).(b) Snapshots of the coupling constant distributions at different densities ρ along the SDRG flow.(c) The typical value of the finite-size gap as a function of the system size L. We have considered chains of up to L = 3 15 spins, with bare disorder strengths D 0 = 1/2 and 1.Here, ψ = 1 2 .The error bars are about the size of the symbol sizes.

5 C = r 2 -Figure 6 .
Figure 6.Average value of the spin correlations (42).(a) The average value of the spin correlation |C(r 1 , r 2 )| [see Eq. (42)] as a function of the internal spin-spin separations r 1 and r 2 .(b) The average values of the diagonal |C(r, r)| (open circles), the off-diagonal |C(r, 2r)| (stars), the integrated ∑ r 2 |C(r, r 2 )| = ∑ r 1 |C(r 1 , r)| (open squares), and the cluster-size ∑ r 1 ,r 2 δ r,r 1 +r 2 |C (r 1 , r 2 )| (closed diamonds) correlations.(c) The average correlation |C(r 1 , r 2 )| as a function of r 1 for different values of r 2 .The system size is L = 3 15 averaged over 1 500 disorder realizations which is sufficient to yield error bars of the size of the symbols.Solid lines are simple guide to the eyes.

Figure 7 .
Figure 7. Finite-size gap analysis of the Griffiths phases in Fig. 2(b).The coupling constants λ A,i are uniformly distributed in the interval [0, 0.2] (set I), [0.2, 1.2] (set II) and [0.5, 1.5] (set III).In all cases, λ B,i are uniformly distributed in the interval [0.5, 1.5], and λ C is uniform and plays the role of a running parameter throughout the Griffiths phases.(a) The typical finite-size gap ln ∆ as a function of the system size L for λ C = 1.The data is averaged over N samples = 6 000 disorder realizations and the error bars are about the size of the data points.Linear fits for the data are provided in the legends from which the dynamical exponent z is obtained.The best fits are restricted to system sizes L > 10 6 .(b) The dynamical exponent z as a function of the distance from criticality Eq.(36).Lines are guide to the eyes.Inset: same data in the main panel in log-log scale.The dashed line is the SDRG asymptotic behavior (δ ≪ 1) z = 1/(2δ ).

Figure 8 . 1 N samples = 6000 Figure 9 .
Figure 8.The typical value of the system gap ∆ as a function of the distance from criticality δ = ln (1 − λ B ) for chains of different sizes L. The coupling constants are such that λ A is uniformly distributed in the interval [0, 0.2], and λ B and λ C are homogeneous with λ C = 1 and λ B being the tuning parameter.The data is average over N samples = 1 500 disorder realizations.The error bars are about the symbol sizes.Solid lines are simple guide to the eyes.The dashed line is the clean behavior in the thermodynamic limit: ∆ ∼ δ φ ∆ , with φ ∆ = 1.

Figure 10 . 2 − 1 ≈
Figure 10.Finite-size gap ∆ (and the corresponding Laguerre bound ∆ LB ) as a function of the system size for chains A and B (see text).The typical value is compatible with activated dynamical scaling ln ∆ ∼ −L ψ , with universal (disorder-independent) tunneling expoente ψ = 1 2 , see dashed line and Eq.(37).The data is averaged over N samples = 10 3 (4 × 10 4 ) disorder realizations for chain A (B), and the error bars are about the symbol sizes.

λ B = λ C = 1 NFigure 11 .
Figure 11.The typical value of the finite-size gap ∆ as a function of the system size L for chains with homogeneous coupling constants λ B = λ C = 1 and {λ A } uniformly distributed in the interval 0, e × λ A,typ for various values of the tuning parameter (from top to bottom) λ A,typ = 0.6, 0.8, 0.9, 0.95, 0.99, and 0.999.The data is averaged over N samples = 10 4 disorder realizations and the error bars are about the symbol sizes.The solid lines are best linear fits to ln ∆ ∼ −z ln L (constrained to L > e 13 ) from which the critical dynamical exponent is obtained (respectively, z = 1.00, 2.08, 4.56, 9.76,  19.48, and 49.82)  and shown in the inset as a function of the distance from the multicritical point γ = 1 − λ A,typ .