Infinite randomness with continuously varying critical exponents in the random XYZ spin chain

We study the antiferromagnetic XYZ spin chain with quenched bond randomness, focusing on a critical line between localized Ising magnetic phases. A previous calculation using the spectrum-bifurcation renormalization group, and assuming marginal many-body localization, proposed that critical indices vary continuously. In this work we solve the low-energy physics using an unbiased numerically exact tensor network method named the"rigorous renormalization group."We find a line of fixed points consistent with infinite-randomness phenomenology, with indeed continuously varying critical exponents for average spin correlations. A self-consistent Hartree-Fock-type treatment of the $z$ couplings as interactions added to the free-fermion random XY model captures much of the important physics including the varying exponents; we provide an understanding of this as a result of local correlation induced between the mean-field couplings. We solve the problem of the locally-correlated XY spin chain with arbitrary degree of correlation and provide analytical strong-disorder renormalization group proofs of continuously varying exponents based on an associated classical random walk problem. This is also an example of a line of fixed points with continuously varying exponents in the equivalent disordered free-fermion chain. We argue that this line of fixed points also controls an extended region of the critical interacting XYZ spin chain.


I. INTRODUCTION
In many situations, phases of many-body quantum systems are stable under weak static, or "quenched," disorder in the presence of a gap, and the disorder average of certain quantities can be calculated in a related clean system via either the replica trick or supersymmetry arguments for non-interacting models [1,2]. However, these methods are not suitable for relevant disorder, or disorder along with interactions, which together produce a rich variety of behaviors. In contrast, real-space thinking should be suitable for directly accounting for spatial inhomogeneity. Interestingly, strong disorder causes certain classes of disordered systems to become tractable on long scales, making real-space renormalization group (RG) approaches amenable to analytical treatments controlled by the flow to infinite randomness. In this work we investigate a modern application of real-space RG to a random XYZ spin chain [3,4], where we use exact numerics to perform unbiased exploration and validation, and also use the strong-disorder renormalization group (SDRG) to demonstrate and characterize such fixed points using the language of random walks.
The original development of a real-space RG appropriate for strong-disorder physics in one dimension (1d) is due to Ma, Dasgupta, and Hu [5]. The feature distinguishing SDRG from, e.g., spin blocking, is that effective degrees of freedom are explicitly associated with an energy scale rather than with a spatial grouping. In this way the disorder realization determines the pattern of integrating out fluctuations.
Such an approach is now understood to be wellmotivated by the idea of an infinite-randomness fixed point (IRFP), a stable solution of the SDRG equations discovered by Fisher in Refs. [3,6,7] at which effective disorder strength grows with the scale without bound, and SDRG predictions become asymptotically exact. In an IRFP, disorder dominates the low-energy physics and physical observables are not self-averaging; average behaviors are instead often determined by rare regions within a disorder realization. Interestingly, although such fixed points lack conformal symmetry, the phenomenology can resemble that of CFT fixed points: for instance, the scaling of average entanglement follows the conformal form with an effective central charge which in some cases is related to the central charge of the clean theory (but does not obey the same rules under RG) [8][9][10].
Since its introduction, the SDRG has been specialized to a variety of classical and quantum systems, and the original scheme has seen many generalizations; see recent reviews [11]. For example, applications in twodimensional (2d) random models also yield IRFPs in these settings [12][13][14][15][16][17][18]. In another direction, SDRG methods were extended to treat all eigenstates of a quantum Hamiltonian [19][20][21][22], in order to assess the possibility of many-body localization (MBL) of excited states. (There are by now multiple reviews of MBL, for instance see Refs. [23,24].) The many-body extended SDRG procedures do not perform an iterative targeting of the lowenergy space, but instead tabulate emergent conservation laws corresponding to the local integrals of motion of an MBL phase; nevertheless, the equations are formally quite similar to the original picture implementing a more traditional RG.
One of the extended many-body SDRG procedures, the "spectrum bifurcation renormalization group" (SBRG) developed in Ref. [21] for Hamiltonians comprising Pauli strings, was applied to the random XYZ spin chain by Slagle et al. [4]. There, along a phase boundary between localized Ising antiferromagnets (proposed to be MBL), disorder-and energy-averaged Edwards-Anderson spin correlations were found to decay as power laws with continuously varying critical exponents. Average entanglement entropy scaling also exhibited a stable effective central charge. The phase transition was conjectured to be "marginal MBL," meaning that eigenstates do not thermalize but exhibit a logarithmic violation of the area law. However, it has recently been argued that such marginal MBL Hamiltonians are perturbatively unstable to ergodicity at finite energy density due to resonances [25,26]. As is true of all excited-state SDRG schemes, Refs. [4,21] rely on MBL for validity, and these recent arguments call this assumption into question.
In the present work we investigate the SBRG findings using unbiased numerics for the ground state and lowenergy excited states. We emphasize that our focus is entirely on low-energy properties, and we will not have anything to say about MBL physics at arbitrary energy density. However, we find the possibility of continuously varying power laws in IRFPs already very interesting and worth further study. The random XYZ chain-while suspected to support infinite-randomness phenomenology in Fisher's original work, Ref. [3]-has eluded understanding due to the lack of a closed-form SDRG solution, and developing a stronger grasp of such instances would constitute an important advance.
Strongly disordered models pose an especially difficult challenge for unbiased numerics, and have long been recognized as among the only 1d models to be resistant to standard methods, chiefly the density matrix renormalization group (DMRG). We apply a relatively new tensor network numerical method named the rigorous renormalization group (RRG) to this problem, as it has already been shown to be effective in the related random XY model [27]. Our goal for the unbiased tensor network computations is to test the findings of Ref. [4], and better understand the disordered fixed points associated with the critical line.
As a brief overview of our results, the data found by RRG are in support of both infinite-randomness physics as well as continuously varying critical indices for disorder-averaged correlations. These conclusions are based on direct measurements in MPS, along with scaling of low-energy spectral gaps, which we solve for in the various symmetry sectors of the model up to systems of length N = 80 spins. Our findings are in general agreement with the SBRG results, namely, that critical indices controlling decay of correlations, as well as long-range mutual information, vary along the critical line, while the "central charge" is fixed. We additionally study the critical exponent ψ, which characterizes IRFP dynamics through the relationship log(1/E) ∼ L ψ between energy scale and length, and find that its value is close to, but may be varying away from, the free-fermion fixed point with ψ = 1 2 . These numerical results for the critical line are captured reasonably well by a self-consistent Hartree-Fock mean-field that treats J z couplings as interactions added to the free-fermion XY chain [throughout, J x,y,z j refer to terms in the XYZ chain as in Eq. (1)]; the Hartree-Fock also apparently produces continuously varying exponents. This finding motivates study of a "locallycorrelated" XY chain with correlations only between terms on the same link of the lattice. The locally correlated model again exhibits similar behavior, and its SDRG structure has an advantageous mathematical connection to the theory of random walks. Within this setting we write rigorous bounds fully determining the critical exponent for power-law decay of a certain average spin correlation function. This exponent indeed varies continuously, proving that the free-fermion critical line of the locally-correlated model is marginal, and is described by a line of IRFPs. This result resolves a question posed by Fisher [3], as illustrated in Fig. 1. In this figure, we parameterize correlations between J x j and J y j by a generic parameter δ varying between δ = 0 (completely uncorrelated or XY model) and δ = 1 (completely correlated or XX model) [for a specific example, see Eq. (19)]; deviation of δ from 1 can also be viewed as introducing random anisotropy to the XX model.
Returning to the interacting model, based on the above understanding of the noninteracting case and the RRG numerical data, we conjecture that at least in the neighborhood of the free-fermion model, interactions are irrelevant and the local correlations generated in the SDRG drive the interacting theory to the line of noninteracting IRFPs at long distances. This scenario is presented in Fig. 2 and represents our conjectured explanation for the continuously varying critical exponents in the XYZ chain.  Fig. 4 of Ref. [3]. In this work we prove the line of fixed points along the exactly marginal direction δ, which describes the degree of correlation between bond terms J x and J y , in the notation of Eq. (1). (Note that δ = 1 corresponds to σ 2 a = 0 in Fisher's notation.) The average anisotropy a is as defined in Ref. [3]; in the present work we consider only the line a = 0.
The outline of this paper is as follows. In Sec. II we present the XYZ spin model and summarize the history of its SDRG, along with explicitly developing the RG  Figure 2. We propose the following schematic flows for the XYZ antiferromagnet, where δ is the degree of correlation between local J x and J y couplings [as defined in Eq. (1)] and J z is the bandwidth of the J z distribution, with statistical isotropy corresponding toJ z = 1. The line of fixed points at J z = 0 is the same as in Fig. 1, andJ z is argued to be perturbatively irrelevant. We conjecture that anyJ z < 1 is irrelevant at δ = 0, but through generation of finite δ flows to the line of non-interacting IRFPs. The methods we employ cannot access the statistically isotropic XYZC or U(1)-symmetric XXZC fixed points, but XXZC was previously described by Fisher [3]. The flows on the dashed line between XYZC and XXZC lie on a manifold separating the basins for XY and ZAF, which is not well described by this slice through parameter space. We avoid any specific conjecture on this matter but remark that it is an interesting topic for further study.
rules in the many-body language. In Sec. III we perform an unbiased study of the ground state using RRG. In Sec. IV, based on our numerical results, we develop both a Hartree-Fock mean-field theory and the free-fermion locally correlated effective model. In Sec. V, we use a picture of the SDRG procedure in terms of random walks to prove continuously varying critical exponents in the locally correlated effective model. In Sec. VI we conjecture a possible long-distance fate of the RG flow for the critical XYZ spin chain, and finally in Sec. VII we discuss the implications of all of these results taken together.

A. Spin chain Hamiltonian
As our most general model we consider the antiferromagnetic XYZ spin chain with quenched randomness in all couplings; that is, The couplings J α j > 0, α = x, y, z, are independent. This model generically has a Z 2 × Z 2 global symmetry, with generators given by the Ising-type operators g x = N j=1 σ x j and g y = N j=1 σ y j . In particular, local field terms are excluded by this symmetry. This model also respects time reversal on the spins, which we implement as g y K, where K is complex conjugation in the z basis.
We impose the same functional form on the disorder distributions of J x j , J y j , and J z j (though delay specification until Sec. III), with bandwidths specified by a set of parametersJ x ,J y ,J z > 0. If the value of any one of these is larger than the other two, the ground state of the model displays Ising antiferromagnetic (AFM) order. As we are considering strong disorder, we anticipate that these phases are localized. If two bandwidths are equal and of the largest magnitude, the model lies on a boundary between localized phases with distinct types of magnetic order; we will primarily consider this case. If all three disorder bandwidths are equal, the model has a statistical S 3 permutation symmetry and sits at a tricritical point in the phase diagram [3,4].
Many exact results are known for phases of the Hamiltonian Eq. (1) in certain limits, and we provide a brief recap here. The SDRG was in fact originally introduced by Ma, Dasgupta, and Hu in order to study the random Heisenberg antiferromagnet with SU(2) symmetry [5], achieved in the present notation by fixing J x j = J y j = J z j for all bonds j. These works argued for the asymptotic development of a power-law singularity in the distribution of couplings and computed leading contributions to critical indices, which vary slowly along the flow. Fisher [3] generalized this analysis to account for anisotropy and performed a thorough study of the resulting phase diagram. The SDRG rules for the random XX model (J x j = J y j and J z j = 0 for all j), which breaks the SU(2) spin rotation symmetry to a U(1) subgroup, are very similar to those of the isotropic model, and in particular both realize random-singlet (RS) phases [28]. In the ground state the microscopic spins are paired up into singlet states at arbitrarily long scales. Correlations between the spins in a singlet are of order unity, and are strongly suppressed with the rest of the system. Thus typical spin correlations are short-ranged, whereas the average correlations are dominated by rare paired spins. This is one hallmark of an IRFP: that a distribution which is broad on a logarithmic scale leads to exponential separation between typical and averaged properties of the state. From the density of paired spins one finds that average spin correlations exhibit power-law decay, scaling as r −2 for separation r. This defines the XX fixed point exponents η x = η y = η z = 2. The characteristic energy scale of the singlets in the RS phase follows where ψ = 1 2 . As a consequence for the density of states, the dynamical exponent is formally infinite.
The random XY chain (i.e., independent J x j and J y j but withJ x =J y ,J z = 0), in contrast, does not realize the RS phase. With the mean in-plane anisotropyJ x − J y serving as the quantum control parameter, Fisher [3] computed the critical exponents ν = 2 and β = 3 − √ 5 for the transition separating Ising x-and y-AFM phases. This was accomplished through a lattice duality mapping to two decoupled copies of the random transverse-field Ising model (RTFIM), whose SDRG equations are also well-studied [6,7,29]. Translating the RTFIM results to the present XY chain, at the phase transition the critical exponent for the decay of x and y components of spin correlations is η is the golden ratio.
Starting from the opposite limit of the XX model, with J x j = J y j for all j, it was also found by Fisher [3] that weak random in-plane anisotropy, which moves along the phase transition toward the XY point, is a marginal perturbation. It was not clear whether this is the case along the entire phase boundary, and we will in fact be led to take up this question in some detail in Sec. V.
The set of exponents for disorder-averaged spin correlations can be completed using the mapping of the XX and XY models to free fermions [30]. For the anisotropic model with S 2 permutation symmetry, η z = 4. In a chain with open boundaries, consideration of the form of the surface magnetization leads to the scaling of the end-toend spin correlations η e x = η e z = 1 for the XX model and η e x = 1, η e z = 2 for the XY model. Focusing on a different type of spin chain, Damle and Huse [31] studied permutation-symmetric multicritical points arising from effective low-energy theories of partially dimerized spin-S models with SU(2) symmetry. They performed a fixed-point analysis of the SDRG equations for degrees of freedom localized at the boundaries between distinct domains of n = 2S + 1 different types of local order (i.e., topological phases distinguished by the properties of edge modes localized near the ends of open chains). Their primary result is a generalization of the n = 2 random-singlet criticality to a countably infinite set of IRFPs with critical exponents ψ = 1 n and ν = 2n √ 4n+1−1 . The permutation symmetry refers to the interchange of distributions for the different types of order, which mediate effective couplings between the domain walls. While the permutation-symmetric tricritical point atJ x =J y =J z in our model shares the statistical symmetry of these theories for n = 3, its microscopic details are dissimilar and it is not clear a priori whether this category of universality applies. Indeed, our estimates of the exponent ψ at the XYZ tricritical point in Sec. III B 3 appear to rule out the applicability of the Damle-Huse universality in this case.

B. Majorana representation
Aspects of this problem become more evident in the language of fermions, for which we use the Jordan-Wigner transformation. Equation (1) maps to a spinless p-wave superconductor with density-density interactions: which has position-dependent hopping t j = J x j + J y j and pairing potential ∆ j = J x j − J y j . Following the idea of Kitaev [32], Motrunich et al. [33], it is enlightening to introduce two species of Majorana fermion, The η j and ζ j are Hermitian, and normalized so that (η j ) 2 = (ζ j ) 2 = 1. In terms of these operators the Hamiltonian is written The symmetry group of the problem is somewhat more expressive in the Majorana language. In the following we specialize to even system sizes N ∈ 2Z. The generators of the global symmetry translate to The symmetries measure fermion parity on two disjoint sets partitioning the Majorana orbitals. The Hamiltonian Eq. (5) takes the form of separate "imaginary random hopping" problems (see Ref. [33]) on these two chains of Majoranas of length N , which we denote X = {ζ 1 , η 2 , ζ 3 , . . . , η N } and Y = {η 1 , ζ 2 , η 3 , . . . , ζ N }. On each chain the coefficients of the Majorana hopping termswhich are fermion parity measurements on adjacent orbitals within a chain-alternate between iJ x j and −iJ y j . There are also inter-chain coupling terms with coefficients −J z j . A single "rung" term iη j ζ j is odd under the parity symmetries, and H instead includes the double-rung interactions −η j ζ j η j+1 ζ j+1 .
The anti-unitary symmetry K (i.e., complex conjugation in the σ z basis) acts on the Majoranas as {i, η j , ζ j } → {−i, η j , −ζ j }. This symmetry prohibits nonzero expectation values of the form iη j η k or iζ j ζ k , even when these orbitals belong to the same Majorana chain.
Constraining J z j = 0 for all j, the resulting Hamiltonian H xy ≡ H[J x ,J y ,J z = 0] is quadratic and can be solved for any particular disorder realization by diagonalization of the auxiliary Bogoliubov-de Gennes (BdG) matrix in the particle-hole basis. The mapping to the Majoranas in Eq. (4) transforms the BdG matrix into a particular form decoupling the two Majorana chains X and Y. This further simplifies the solution for the singleparticle eigenstates to diagonalization of a pair of N × N tridiagonal matrices.
As we are considering boundaries between Ising ordered phases, the natural observables are the corresponding magnetic order parameters σ α , α = x, y, z. Written in terms of fermion operators, the spin correlation func- From Wick's theorem, in the ground state of any specific disorder realization C x (j, j + r) and C y (j, j + r) can be computed as Pfaffians of antisymmetric 2r × 2r matrices, and the calculation further simplifies due to the separation into two Majorana chains. We focus on this case and consider the angle brackets · as denoting expectation values measured in the ground state, although the expressions Eqs. (8)-(10) apply more generally. We will be discussing disorder-averaged correlations C α (j, j + r) and when this is clear we will drop the overline. In the following we work exclusively along the line with statistical symmetry between J x j and J y j and will often collectively refer to C x,y (j, j + r), as C ⊥ (j, j + r).
C. Strong-disorder renormalization group

Decoupled Majorana chains
Examining the Hamiltonian on Majorana chains X and Y also clarifies the form of the analytic SDRG. In the decoupled model H xy , the RG proceeds independently on each of the chains, which are endowed with parity conservation. The SDRG for a single such chain was developed explicitly in the single-particle spectrum language by Motrunich et al. [33] and in the many-body Hamiltonian language by Monthus [22]. We review the result here, specialized to our case, in the many-body language, which naturally extends to the interacting problem [22]. For now we consider only a single Majorana chain, and relabel the orbitals as γ n , n = 1, . . . , N . The Hamiltonian acting on this chain is H M = N −1 n=1 ih n γ n γ n+1 . Suppose that the largest energy scale is set by H 0 = ih k γ k γ k+1 for some k ∈ [1, N − 1]. H 0 measures fermion parity on the two orbitals, with eigenvalues ±h k associated with the two parity states; denote the splitting by Ω = 2h k . Accordingly, this term is diagonalized by the complex fermion mode f † 0 = 1 2 (γ k + iγ k+1 ), which has projectors π + = f 0 f † 0 and π − = 1−π + = f † 0 f 0 into the even and odd parity sectors, respectively. In terms of the projectors we have H 0 = (Ω/2)(π + − π − ).
The rest of the terms in H M ≡ H 0 +V can be treated as a perturbation if the nearby couplings are much smaller than the local gap |Ω|. Although this condition may not be satisfied initially, the validity of the assumption improves during the RG flow because the SDRG generates an effective disorder distribution with increasingly broad logarithm. The rest of the Hamiltonian can be divided into diagonal and off-diagonal components with respect to H 0 ; specifically, Note that V od contains only a constant number of local terms. We denote the small scale of these terms relative to H 0 by the parameter . The effective Hamiltonian with emergent good quantum number f † 0 f 0 is found by a Schrieffer-Wolff transformation eliminating V od up to O( 2 ) [34][35][36][37]. That is, H M = e iS H M e −iS , where the Hermitian generator of the rotation can be expanded in powers of as S = S [1] + S [2] + · · · . The conditions on the rotation are that S [1] is off-diagonal and satisfies V od = [H 0 , iS [1] ], and S [2] eliminates off-diagonal terms at O( 2 ) (but we will not need to write it explicitly). A suitable generator is iS [1] the final line being Eq. (17) of Ref. [22]. The off-diagonal terms are those which share an odd number of Majoranas with H 0 and thus anticommute. Consequently V od = ih k−1 γ k−1 γ k + ih k+1 γ k+1 γ k+2 and Finally the rotated Hamiltonian is This result includes a renormalization of the strength of the H 0 term which increases the magnitude of the splitting, in addition to a new term ih k−1 γ k−1 γ k+2 . By projecting into the low-energy sector of H 0 (which depends on the sign of h k ), the Majoranas γ k and γ k+1 are frozen into one of the definite parity states of the complex fermion mode, and thereby decoupled, or "decimated," from the effective Hamiltonian. The single effective coupling h k−1 replaces three hopping terms in H M . Because the new term maintains the imaginary random-hopping form, the SDRG is closed in this model space and can be iterated, with the flow acting on the disorder distribution of the couplings {h n }. During the RG flow, some of the terms involved in decimations will be themselves renormalized couplings from prior steps; they can be made to fit the present format by re-indexing the chain after every step to remove the decimated Majorana orbitals. In addition, the specific form of the renormalized coupling h k−1 permits a framing of the SDRG in terms of a classical random walk; this approach will be developed in detail in Sec. V. The many-body Hilbert space is therefore decomposed into a tensor product of non-interacting complex fermions in definite parity states. Returning to the XY model viewed as two decoupled Majorana chains and running the above procedure independently on each of the chains, one can deduce from the signs of the couplings in Eq. (5) that the ground state is even under g x and g y if N mod 4 = 0 and odd under g x and g y if N mod 4 = 2. The ground state spin correlations in an eigenstate of the Hamiltonian can also be understood from this picture; see Sec. II D.
As a technical remark, one way to deal with the signs of the couplings in Eq. (5)-needed to deduce g x and g y quantum numbers as well as the signs of the correlation functions-is to perform a gauge transformation of the Majorana fermions as η j = s j η j , where s j = 1 if j = 4n + 1 or 4n + 2 and s j = −1 if j = 4n + 3 or 4n + 4, while ζ j = s j (−1) j+1 ζ j . The Hamiltonian written in terms of the primed Majoranas takes the form j iJ x j ζ j η j+1 +iJ y j η j ζ j+1 , i.e., all Majorana hopping amplitudes are positive in the convention where the Majoranas are written in the same order as they appear on the chain: ih nm γ n γ m with n < m has h nm > 0. This property is preserved under the SDRG, which simplifies analysis of the signs. For example, for Majoranas γ n , γ m with n < m decimated as a pair we then have iγ n γ m = −1 at the zeroth order in the SDRG, and using the non-crossing property of the pairs in each Majorana chain fixes the signs of correlations in Eqs. (8)-(10) to be (−1) j−k . To avoid confusion, in formulas we keep using the original Majoranas as in Eq. (5).

Majorana problem with inter-chain interaction terms
In the presence of interactions coupling the two Majorana chains, it is necessary to consider the full Hamilto-nian Eq. (5). In the notation of the present section we Because all of the terms in H are measurements of fermion parity, the general framework from the previous section-in particular Eq. (16)-still applies. Now there are two cases: the largest energy scale can be set by one of either the hopping terms {h M n } or the interactions {K n }. While one can in principle consider both cases following Ref. [22], for our purposes we will study only the hoppingdominated case. Suppose that The components appearing in each off-diagonal block of the Hamiltonian are The effect of the interactions in perturbation theory is simply to modify the couplings into operators which we refer to as "interacting couplings:" h X k±1 → h X ,int k±1 . This is a reasonable shorthand because the interacting couplings commute with each other and all fermion operators appearing in the formula. Then from the result Eq. (19), Projecting into the low-energy sector sets iγ X k γ X k+1 → −sgn(h X k ) and again decouples the Majorana operators γ X k and γ X k+1 from the rest of the system, decimating them by creating a complex fermion mode with definite parity. As in the non-interacting case, the magnitude of the splitting is increased by renormalization of H 0 , and a new hopping term h X k−1 is added to the X chain. However, the leading-order effect of the interactions, at O( ), arises from V d , where the "degradation" of the term As a result, correlations develop between the hopping terms on the same bond. This aspect of the perturbation will constitute the basis of a mean-field study of the interacting system, presented in Sec. IV.
The effective Hamiltonian also includes renormalized couplings h Y k−1 and h Y k+1 , as well as new four-fermion terms which change the structure of the lattice graph, and a six-fermion term. The appearance of these terms breaking the form of H, as well as the generation of correlations between terms, are an indication that the RG flow cannot be tracked exactly in the interacting model. However, if the interaction terms already tend to be weak compared to the hopping, the higher-order terms generated by this process will accordingly be weaker still. This is the situation, at least initially, in the random XYZ model with smallJ z ; however there is no guarantee at this point that the relative strengths of the different types of couplings are maintained asymptotically. We will return to this question more systematically in Sec. VI, after we understand the non-interacting problem with correlated Majorana hopping amplitudes in the two chains in Sec. V.

D. XY model spin correlations in SDRG
From the controlled SDRG for the random XY model one can deduce that average correlations in the ground state follow power laws-although typical correlations are short-ranged-and even calculate the exponents. One also obtains a more qualitative picture of the behavior of the spin correlation functions.
Expanding Eq. (10) in the ground state at distance r, Other terms vanish due to symmetry. One sees immediately that C z (j, j + r) = 0 if r is even. For odd r, C z (j, j + r) assumes a large value if and only if the sites j and j + r were decimated together on both Majorana chains, in which case both expectation values iη j ζ j+r and iζ j η j+r have approximately unit magnitude and opposite sign, so the sign of C z is negative. Otherwise if this decimation did not occur in one or both Majorana chains the contribution is suppressed, arising only from higher-order terms in the perturbation theory. Consider the correlations averaged over sites j as well as over disorder realizations, which average we denote C z (r). Nearly all terms will be vanishingly small, with rare terms of roughly unit magnitude occurring with some density; these dominate the average. It is a result of Ref. [3] for the RS phase that at sufficiently large separation the likelihood of such a decimation scales as r −2 ; thus for two independent Majorana chains η z = 4. The transverse correlations Eqs. (8) and (9), summarized as C ⊥ (j, j +r), are the expectation values of strings of 2r Majoranas. Such operators are evaluated as the sum of r-fold products of expectation values of symmetryallowed bilinear contractions, with signs arising from the signature of each permutation. A term in the sum has a large value if and only if it contracts all Majoranas with their decimation partners in the SDRG. This will be the case for exactly one term if all decimations of the Majoranas appearing in the string expectation value are "internal;" that is, if all decimation partners are also included. If any Majoranas were decimated with orbitals which do not appear in the string, the expectation value will be small. We again define C ⊥ (r) as the average over sites and disorder realizations.
If on both chains X and Y the sites j and j + r are decimation partners, then as described above, this pair contributes a large value to C z (r). The pair also necessarily contributes a large value to C ⊥ (r), as pairing the extremal Majorana orbitals in a string implies that all decimations are internal to the string. Thus, the critical exponent η ⊥ lower-bounds η z . As reviewed earlier, for the random XY model η ⊥ = 3 − √ 5 ≈ 0.764; the bound is saturated in the XX model where η ⊥ = η z = 2 [3].
Finally, the SDRG picture also tells us about the endto-end spin correlations in the XX and XY models. The expectation value C z (1, N ) ≡ C z (N ) obtains large contributions if on both Majorana chains the end sites 1 and N are paired in the SDRG. While such occurrences in the two chains are perfectly matched in the XX model and have probability 1/N or η e z = 1, in the XY model the occurrences are independent, giving η e z = 2. On the other hand, the expectation value C ⊥ (1, N ) ≡ C ⊥ e (N ) includes all Majorana orbitals on one chain, and all but those at sites 1 and N on the other. This string has a large expectation value if all of these Majoranas are paired internally, which is to say that the two excluded Majoranas are decimated together. As this is occurs on a single chain only, it has the same probability in both the random XX and XY models. Indeed, η e ⊥ = 1 in both cases [30].

A. "Rigorous RG" numerical method
The standard numerical technique for equilibrium states of many-body quantum systems in 1d is the density matrix renormalization group (DMRG) [38,39], which has been remarkably effective in conjunction with matrix product state (MPS) representations of low-energy wavefunctions [40,41]. Over nearly 30 years, DMRG has seen enormous practical success in a wide range of models of physical interest. However, for some time its effectiveness was not well explained: even as MPS attained a rigorous footing with the proof of the area law of entanglement in 1d [42][43][44], the existence of an efficient algorithm for eigenstates given an area-law Hamiltonian remained unclear. It was not until the work of Landau et al. [45] in 2015 that a polynomial-time algorithm was developed for ground states of gapped models, proving that an efficient method is possible in principle.
However, the algorithm exhibited in Ref. [45] bears little resemblance in its particulars to DMRG, and a similar proof for the DMRG algorithm appears to be challenging; in fact, it is known that popular multi-site variants can be NP-hard in the worst case [46]. As a practical matter, in systems with strong disorder DMRG is susceptible to spurious convergence to excited states, an outcome which cannot be readily diagnosed [47]. This is fundamentally a consequence of performing an iterated local optimization over MPS parameters. The rigorous algorithm is distinguished by a reliance on an approximate ground state projector (AGSP), an operator derived from the Hamiltonian, which was introduced by Arad et al. [48]. The role of the AGSP is to provide global information, ensuring that intermediate states can be efficiently represented and directing the algorithm along a computationally tractable route to the ground state.
AGSP-based methods were later generalized to lowenergy excited states in models with slightly relaxed conditions on the density of states [49]. Based on this work, in collaboration with Vidick we introduced the rigorous renormalization group (RRG), a numerical implementation for low-energy states of local Hamiltonians in one dimension [27]. While the implemented method differs slightly from the proof construction and does not strictly satisfy the conditions of the guarantee-whose parameters are not known a priori regardless-it inherits the intuitive benefits of the AGSP and has been seen to be effective in practice for nontrivial low-energy spectra like those of strongly disordered systems, or in the presence of nearly degenerate manifolds [27,50], where DMRG may be unreliable.
In the following sections, we perform a numerical study of the lineJ z ∈ [0, 1],J x =J y = 1, in the phase diagram of Eq. (1), using RRG. Our objective is primarily to verify by unbiased numerics the observation of continuously varying critical exponents in the SBRG study of Slagle et al. [4], and then to shed additional light on the nature of the low-energy theory. (Here we focus solely on the ground state properties and low-energy physics, rather than the question of MBL.) For concreteness, we use the disorder distribution described in Eqs. (3) and (4) of Ref. [4], namely, We use a milder disorder strength Γ = 2, as compared to Γ = 4 for the previous work [4]. Both choices lead to strong disorder physics and the specific value should have little effect on the universal low-energy physics for large enough systems. However, we find that the logarithm of the distribution of the energy gaps depends significantly on Γ, with smaller values tending to lead to larger gaps; this eases the challenge to the numerics which in any case are limited by double-precision floating-point errors on the order of 10 −16 . In RRG we are capable of accurately resolving energy scales down to log 10 (Ω/ ) ∼ −12, and validate our results against the free-fermion solution at the soluble pointJ z = 0.
To construct the AGSP for RRG we use a Trotter approximation to a thermal operator e −βH . The output of the RRG algorithm is a subspace of constant dimension approximating the low-energy states of the model. We use an implementation based on ITensor [51], in which we explicitly realize the Z 2 × Z 2 symmetry and solve for the lowest two eigenstates in each of the four symmetry sectors [52]. In each case the MPSs generated by RRG are then further optimized using DMRG in order to minimize the overlap with high-energy states. The RRG "hyperparameters" s and D (see Ref. [27] for details) are chosen so that for the majority of disorder realizations DMRG can optimize the RRG output in a small number of sweeps. For approximately the most challenging 1% of realizations, DMRG requires many sweeps to converge. In these instances we repeat the calculation, increasing  Table I. RRG hyperparameters are shown for values ofJ z studied numerically. As described in the text, we optimize the output of RRG using DMRG, and for finiteJ z take as a measure of accuracy the number of sweeps required for convergence. These values of s and D were chosen in order to accurately converge approximately 99% of disorder realizations on N = 80 spins. For the small fraction of more difficult realizations which are not solved by the hyperparameters above we repeat the algorithm with increased values, finding that convergence is achieved this way.
the RRG hyperparameters, and find that the improved RRG states are easily converged by DMRG. From comparison with exact free-fermion results forJ z = 0 obtained by numerical matrix diagonalization, we find that if RRG produces states which are successfully converged by DMRG and the excitation gap is larger than the target threshold 10 −12 , the ground state energy and gap are numerically exact in 99.5% of realizations. As we will show in the following section, atJ z > 0 the finite-size gaps tend to be larger than those atJ z = 0 and should be easier for RRG; thus we believe our results are even more reliable for these points.

Critical spin correlations
We measure spin correlations in the RRG ground state of H[J x = 1,J y = 1,J z ] withJ z ranging from 0 to 1 and microscopic disorder strength Γ = 2 throughout. Bulk correlations in an open chain of length N are measured for r ≤ N 2 including only sites j, j + r ∈ { N 4 , . . . , 3N 4 }, in order to distinguish the power law from the end-to-end correlations closer to the boundaries. We show disorderaveraged correlations data measured in chains of length N = 80 sites in Fig. 3, which includes slices at values ofJ z moving along the phase boundary from the freefermion model to the tricritical point. Already the raw data clearly shows power laws with varying exponents for both C ⊥ and C z in the bulk.
End-to-end spin correlations are measured only between the single pair of sites 1 and N for each disorder realization, and exhibit correspondingly larger statistical fluctuations. In addition, reproducing C z e (N ) correlations presents a singular challenge for the RRG algorithm. As discussed in Sec. II D, in the SDRG the likelihood of a nonzero value of σ z 1 σ z N at the XY freefermion point is the square of the probability of an endto-end singlet in a spin chain of length N in the RS phase. That is, the distribution is broad on a logarithmic scale, with the average being dominated by a very small tail. More importantly, the disorder realizations located in the tail-of outsize importance in the average-are those on  Open circles indicate C ⊥ (r) data, while filled circles mark C z (r). The disorder averages for each value ofJ z include 1500 realizations. In the spatial average we include only the middle half of the spin chain-that is, only sites in { N 4 , . . . , 3N 4 }-in order to separate the bulk correlations from the ends, which exhibit different scaling laws. See Fig. 5 for the critical power law decay exponents extracted from this data. In order to measure the power laws we show the absolute value of the correlations, which originally have a staggered sign pattern (−1) r . In addition, only odd r are shown for C z data because the values for even r, though demonstrating a similar power law, are much smaller (atJ z = 0 they are identically 0, see Sec. II D). which sites 1 and N were decimated together on both Majorana chains, which correlate with the smallest excitation gaps in the low-energy spectrum and are the most difficult realizations for the method to solve accurately. We show disorder-averaged end-to-end correlations as a function of N in chains up to N = 80 in Fig. 4. One sees that the C ⊥ e correlations depend weakly onJ z and have close slopes on the log-log plot, suggesting similar power law exponents. On the other hand, the C z e correlations depend strongly onJ z and despite evident statistical scatter appear to have varying slopes.
Our unbiased numerical results for the bulk correlations are in broad agreement with the finding of Slagle et al. [4] of critical exponents governing the decay of spin correlations that vary continuously withJ z . In contrast to the previous approach, we perform direct measurements in optimized MPS for the ground state. We show the extracted power law exponents for the bulk and end-to-end correlations in Fig. 5 as a function ofJ z . As expected, the C ⊥ and C z exponents approach each other at the tricritical (permutation-symmetric) point J z = 1, where we estimate the bulk critical index to be η ⊥ = η z ≈ 1.48.  Both bulk and end-to-end exponents are included, with known results for the bulk correlations in the free-fermion model atJ z = 0 indicated by red stars, and results for end-to-end correlations by yellow diamonds. An increase in statistical noise is evident in the end-to-end correlations as compared to the bulk. The reason that these computations, particularly C z e (N ), are more difficult, is discussed in the text.

Entanglement structure
We also study measures of entanglement in the RRG ground states for varyingJ z . The average bipartite entanglement entropy of a connected subsystem of length adjacent to the system boundary is known to scale ac-cording to the conformal field theory result S b ( ) =c 6 ln , with a universal constantc. In some cases the "effective central charge"c is apparently related to the central charge of the clean model [8]; for example, in the critical phase of a single Majorana chainc = ln 2 2 = c ln 2, where c = 1 2 is the central charge of a clean Majorana fermion chain. Accordingly, the XY fixed point hasc = ln 2, being equivalent to two decoupled critical random Majorana chains. From finite-size scaling of the disorder-averaged half-system bipartite entanglement entropy S b (N/2) we find with fair precision thatc is stable at this value for any interaction strengthJ z along the critical line, in agreement with Ref. [4]. The subsystems A and B considered in this case are single spins separated by a distance r, and the average is taken over sites in the middle half of the chain. Also shown is the effective central chargec, found from finite-size scaling of the half-chain entanglement entropy. Whilec appears to be insensitive to the coupling between the two Majorana chains, the LRMI exponent varies continuously.
We also measure long-range mutual information (LRMI) between disconnected regions; the formula for this entropic quantity in terms of the entanglement entropy of a subsystem is I(A : B) = S(A) + S(B) − S(A ∪ B). We will take A and B to be single spins separated by a distance r; Ref. [4] found that up to appropriate rescaling, the lengths of the subsystems do not affect the asymptotic behavior. The disorder-averaged LRMI we denote I(r), and this quantity will decay no faster than the slowest observable. That is, in the symmetric ground state of an ordered phase I(r) will be long-ranged; in a phase without order one expects exponential decay; and at a critical point the exponent ρ, I(r) ∼ r −ρ , lowerbounds the power-law decay exponent of any local observable. We show disorder-averaged LRMI data in the upper panel of Fig. 6. The critical exponent ρ varies continuously withJ z , as is the case with the other critical indices measured, and is very close to the exponent η ⊥ , suggesting that the correlations of the order parameters for the adjacent phases saturate the lower bound everywhere along the boundary. Our RRG results for ρ as well as the effective central chargec are shown in the lower panel of Fig. 6. AtJ z = 1 we estimate ρ ≈ 1.73, which is somewhat larger than the estimates of η ⊥,z but is in general agreement and is also similar to the SBRG estimates in Ref. [4].

Scaling of excitation gap
Because RRG produces not only the ground state but a constant number of low-energy states, it is possible in principle to study spectral properties as well. We focus first on the simplest of these, the energy gap to the lowest excitation in a finite system. From the SDRG for the free-fermion point one observes that this excitation consists of flipping the parity of the complex fermion associated with the lowest-energy (i.e., the last decimated) pairing on either Majorana chain. As we consider chains with lengths that are multiples of 4, the ground state is found in the (g x , g y ) = (+1, +1) sector of the global (Z 2 ) 2 symmetry and the first excited state will be found in either the (+1, −1) or (−1, +1) sector.
The distribution of excitation gaps is known exactly via the mapping to two decoupled copies of the RTFIM, where the universal form of the gap distribution is known from the work of Fisher and Young [29]. The gap in the random XY model is the minimum of two independent random variables sampled from the distribution of Ref. [29]. In Fig. 7 we show histograms of the (logarithmic) excitation gaps for the random XYZ model with varyingJ z for chains of length N = 80. The exact distribution for theJ z = 0 point is indicated with a dotted line.
Indicated on Fig. 7 by vertical lines and the labels MJ z are the medians of the histograms; these are provided as a characterization of the distributions that is not overly sensitive to the tails, where the energy gaps can be close to the numerical threshold. While the precise tails are not accessible, it is rare for RRG to make an error which would move a disorder realization out of the tail into the bulk of the distribution. Thus, the median provides an accurate summary of the gap distribution although the mean cannot be reliably estimated. In Fig. 8 the scaling with chain length of the median of the gap distribution is shown with varyingJ z . This allows an estimate of  The medians include long tails that are not shown, as they contain energy gaps too small to be accurately measured by the RRG algorithm; however the estimate of the median is not sensitive to these uncertainties. The trace for each value ofJ z includes 1500 disorder realizations.
the exponent ψ controlling the length-energy relationship Eq. (2), which takes the value ψ = 1 2 at the freefermion point. The RRG scaling data suggest that there may be a systematic drift in ψ asJ z is varied toward the permutation-symmetric pointJ z = 1, however it is difficult to exclude the possibility of a stable ψ with a long crossover aroundJ z = 1. In either case, this result does not support the n = 3 Damle-Huse universality for this tricritical point.

Symmetry properties of low-energy states
As described in Sec. III B 3, in the non-interacting model H xy , the symmetry properties of the ground and low-lying states can be deduced from the single-particle excitations used to build the many-body states. For convenience we relabel the Z 2 × Z 2 symmetry sectors (always working on systems with N ∈ 4Z): denote the freefermion ground state sector (g x , g y ) = (+1, +1) as 0; the sector (−1, −1) as 1; (+1, −1) as 2; and (−1, +1) as 3.
The value of the critical exponent ψ extracted from finite-size scaling of excitation gaps in RRG is shown. The upper panel shows the finite-size scaling of the medians MJz (shown in Fig. 7 for N = 80), with each data point including 1500 disorder realizations. The lower panel shows the extracted power law exponents for both the first gap, denoted E1 − E0 (found from the data shown in the upper panel) as well as the second and third energy gaps. At the free-fermion pointJ z = 0, ψ = 1 2 , and the systematic deviation from the exact value is likely due to finite-size corrections. At this point the first and third energy gaps are very often identical, both being associated with the lowest-energy decimation on one chain. Away from this point, this is no longer necessarily the case and a drift in ψ withJ z is visible in the E1 − E0 curve. excitation of the two lowest energy single-particle states. With a logarithmically broad disorder distribution, as at an IRFP, the third excited state is very likely to be of the latter type; thus we expect that for sufficiently long N , the four lowest-energy states of H xy will most often come from the sectors {0, 2, 3, 1} or {0, 2, 2, 0}, or their Z stat At the tricritical point this picture cannot apply, as the S 3 counterparts of the free-fermion-allowed configurations (these include, e.g., {0, 1, 2, 3} and {0, 1, 1, 0}) must also occur and with equal likelihood; thus we study the critical line by tabulating occurrences of free-fermiondisallowed low-energy configurations in disorder realizations with finiteJ z . We classify the various configura-   Figure 9. Sampled estimates of the likelihood of the various symmetry patterns of low-energy states are shown as a function ofJ z . The lower panel shows the same data as the upper, zoomed in on the bottom of the y-axis. The free-fermionallowed Types 1, 2, and 3 are defined above and drawn with solid lines, and the free-fermion-disallowed Types 1 * , 2 * , and 3 * consist of all other partners under the action of the S3 statistical symmetry, and are drawn with dashed lines. Here we provide summary data which is averaged over system sizes N = 32, 48, 64, 80, with 6000 total disorder realizations for each value ofJ z . (In Fig. 10 we study the dependence on N .) AtJ z = 0 we assume that only Types 1, 2, and 3 are present and include eigenstate permutations of the exact symmetry pattern for very small splittings < 10 −12 ; nevertheless there is still a low rate of "Other" instances.
tions as described in the table in Fig. 9, and their likelihood in our sample of disorder realizations is plotted. Note that in this plot we have averaged over all system sizes, in order to provide an initial summary of the typical behavior (we will study the scaling behavior with N later).
For H xy the dominant pattern is Type 1, with a substantial minority of Type 2 and very few of Type 3. The S 3 counterparts, which are forbidden in the picture of decoupled Majorana chains, are labeled Types 1 * , 2 * , and 3 * . The category "Other" includes all low-energy configurations not matching any of the types already described. There is a very small, though finite, fraction of such instances; however these are nearly entirely associated with very small excitation gaps. As already described, in such situations with very small splitting RRG cannot systematically identify the lowest-energy state or the exact sequence of excitations, so the precise order of symmetry sectors is not reproduced. AtJ z = 0 we are able to "interpret" many such cases by assuming that the energy-permuted free-fermion-allowed symmetry pattern is the correct one, though away from this point a corrected type cannot be uniquely determined. (AtJ z = 0 some low-energy patterns found by RRG cannot be interpreted as one of the free-fermion-allowed configurations, and these are the realizations classified as "Other" at this point.) Moving away fromJ z = 0, the Types 1 * , 2 * , and 3 * occur with increasing probability. We find that Type 2 decreases more quickly for smallJ z than Type 1, which is in line with our understanding, developed in Sec. IV, of the interaction as introducing correlations between the Majorana chains (such correlations make it less likely that the two lowest-energy single-particle states occur in the same Majorana chain). The rate of "Other" instances is very low and decreasing with increasingJ z , suggesting that these remain attributable to errors due to small energy gaps, and the only new types of symmetry pattern appearing at low energy are those related to the free-fermion-allowed types by S 3 . As one expects from the definitions of each type, the frequency of Types 1 * , 2, and 3 * , are roughly twice those of Types 1, 2 * , and 3, respectively, atJ z = 1. Here the S 3 partners Types 1 + 1 * describe roughly 91% of disorder realizations, with Types 2 + 2 * and 3 + 3 * describing roughly 4.5% each.  Figure 10. The ratio of the combined likelihood of the freefermion-disallowed Types 1 * + 2 * + 3 * to the combined likelihood of Types 1+2+3 is shown as a function ofJ z , separately for system sizes N = 32, 48, 64, 80. Each data point includes 1500 disorder realizations. For intermediateJ z ∈ (0, 1), there is a consistent trend toward lower probabilities as the system size increases from N = 32 to 64, meaning that the low-energy symmetry patterns of longer systems are more likely to be free-fermion-like. The quantity p(1 * +2 * +3 * ) is very similar for system sizes N = 64 and 80 at all values ofJ z , with the difference being within the apparent statistical scatter. At the tricritical pointJ z = 1 the predominant scaling behavior is reversed, and the quantity appears to be converging toward its long-distance fixed value from below with increasing system size N .
From the above general picture of the low-energy states we learn that the critical line is characterized by the increasing probability of the free-fermion-disallowed sym-metry partners Types 1 * , 2 * , and 3 * with increasing interaction strengthJ z . The dependence of these probabilities on system size provides a hint about the RG relevance or irrelevance of the interaction. In Fig. 10, we show the ratio of the combined likelihood of Types 1 * + 2 * + 3 * to that of Types 1 + 2 + 3 as a function ofJ z for each system size separately [53]. While these data suffer from poorer statistics than those of Fig. 9, there is a trend for allJ z ∈ (0, 1) toward lower probabilities with increasing N , meaning that at longer scales the disorder realizations appear more free-fermion-like. The system sizes N = 64 and 80 are quite similar by this measure, and the differences between these values are smaller than the apparent statistical noise. In contrast, the dependence on system size is opposite at the tricritical pointJ z = 1, as the likelihoods converge to their asymptotic value from below with increasing length scale. In Secs. VI and VII we make a conjecture consistent with this observation, that the interactions may in fact be irrelevant but the SDRG generates a marginal perturbation (corresponding to the local correlation of renormalized terms, see Sec. IV) which ultimately takes the system to a line of free-fermion fixed points with variable exponents.

IV. MEAN FIELD THEORY OF INTERACTION
Turning onJ z > 0 introduces four-fermion interaction terms to the quadratic Hamiltonian H xy . These terms couple the Majorana chains X and Y in such a way that the ground state is no longer analytically tractable under SDRG, which generates multi-fermion terms in the effective Hamiltonian that proliferate with increasing RG scale. However, as mentioned in Sec. II C 2, if at some point in the RG the interaction terms are typically weaker than the hopping terms then the effective higher-order descendants will be even weaker. One might hope, then, that by beginning with a bandwidthJ z J x ,J y the strength of these terms may be suppressed at all scales, leading to only a minimal effect on the criticality.
Based on this understanding, we consider the mean field theory by "expanding" the interaction into fermion bilinear terms. In the Majorana language, the meanfield structure is particularly transparent; here the only symmetry-allowed bilinear terms act internally on the chains. For J z j 1, This can also be seen in terms of the original spins, where the mean field theory takes the form The effect of the allowed terms is to renormalize the existing couplings in the following way: With expectation values · understood to be evaluated in the ground state of the mean field Hamiltonian with parameters (J x j ) mf , (J y j ) mf , the above represent selfconsistency equations (i.e., minimization equations in the variational perspective of the mean field theory). Because the Majorana chains remain decoupled, the meanfield theory can be solved in the analytic SDRG, at least in principle, by accounting for the distributions of effective J x j and J y j couplings no longer being independent. In the following subsections we numerically investigate the universal behavior of this mean-field theory, and provide exact results from the analytic SDRG in Sec. V.
A. Self-consistent Hartree-Fock treatment of interaction terms 10 0 10 1 r 10 −5 Figure 11. Bulk correlations data from the self-consistent Hartree-Fock mean-field theory are shown with varying band-widthJ z , up to separation r = 64 in chains of length N = 128.
Filled markers indicate C z (r) data, and open C ⊥ (r). The disorder averages for each value ofJ z are taken over 25000 realizations and include only the middle half of the spin chain, as described in the caption to Fig. 3. These simpler free-fermion calculations are cheaper to perform, and accordingly exhibit better statistics than those of Figs. 3-10.
We first perform a self-consistent numerical study of the interaction term in the quadratic mean-field theory by directly implementing Eqs. (34) and (35) in the BdG Hamiltonian, iteratively solving the ground state of the Hamiltonian and updating the mean-field couplings until reaching convergence. The bulk correlations data in the thus determined mean field ground state are shown in Fig. 11, end-to-end correlations in Fig. 12, and a summary of the critical exponents in Fig. 13.
The key finding of the mean field treatment is that the power law exponents in all correlation functions do evolve  Figure 13.
Critical exponents are shown for the selfconsistent Hartree-Fock mean-field theory with varying interaction strengthJ z ∈ [0, 1], extracted from the correlations data in Figs. 11 and 12. Both bulk and end-to-end exponents are included, with known results for the bulk correlations in the free-fermion model atJ z = 0 indicated by red stars, and results for the end-to-end correlations by yellow diamonds. The pointJ z = 1 in this model does not feature any special symmetry.
withJ z in a similar way to those of the interacting model. This not necessarily expected since, e.g., in a clean XXZ model the mean field, while capturing some short-range energetics, cannot capture varying power laws in the fully interacting theory. By understanding the features in the mean field responsible for capturing the varying power laws in the random XYZ chain, in the following sections we will be led to a plausible scenario for the physics of this system.
While the mean field theory is reasonably accurate for J z ≤ 0.4, it is evident from Fig. 11 that the magnitudes of the mean-field correlation functions aroundJ z = 1 do not approach their actual values. At the tricritical point of the interacting model the statistical S 3 symmetry of the Hamiltonian leads to the equivalence of the averages C ⊥ and C z ; as the mean field lacks this symmetry, it is not surprising that the distinction persists. Moreover, there is nothing special aboutJ z = 1 in the mean-field model; note also that this specific mean field does not allow any symmetry breaking, and we see that the best it can do upon increasingJ z is to approach the XX chain, which is a poor approximation forJ z 1.
Nevertheless, buoyed by the success of the mean field at smallJ z , we now examine more closely the effective parameters (J x,y j ) mf . As the interaction strength is increased, the J x j and J y j terms tend to become more similar. We can clearly see how this happens in the spin formulation of the self-consistent mean field of Eqs. (34) and (35): a large bare AFM J x j > 0 will tend to correlate σ x j and σ x j+1 strongly antiferromagnetically (achieving σ x j σ x j+1 ≈ −1 if this is the dominant coupling), and in the presence of AFM J z j > 0 this will lead to an increase of the effective AFM J y j coupling, and vice versa. However it is not clear what sort of model the full self-consistent mean field treatment actually constitutes, as the iterated nature of the solution could lead to long-range correlations effects among the couplings. In the following section we propose a more straightforward model intended to broadly capture the features of this self-consistent Hartree-Fock mean field. We will see that the ultra-short-range correlations among J x j and J y j identified above can already explain continuously varying power laws.
B. Numerical study of random XY chain with locally correlated couplings

Definition of locally-correlated XY model
The rules Eqs. (34) and (35) for the mean-field couplings modify bonds on one Majorana chain based on expectation values across the same bond on the other chain. As a result, recalling that J z j > 0 for all j, the terms on a given bond-which at the mean-field level are strengthened by the interactions-develop correlations among themselves. Terms on separate bonds also get correlated in less obvious ways, since the mean field ground state is influenced by all bonds, but we will proceed by ignoring such longer-range correlations among the couplings. We refer to such an effective model as having "local correlations," in order to distinguish from spatial correlations between terms on separated bonds. One can mimic the behavior of the mean field theory and explore the effects of such correlations using the following parameterization of the couplings: for A j , B j independent random variables and δ ∈ [0, 1], let Tuning δ from 0 to 1 interpolates between fully independent couplings and the perfectly correlated case with U(1) symmetry. That is, the parameterization runs along the line between the random XY and XX spin chains. As mentioned in Sec. II A, Fisher [3] found that weak random anisotropy is marginal around the XX point, which is in the RS phase. However it was not resolved whether this perturbation is truly marginal, or perhaps instead marginally relevant or irrelevant. The mean-field numerical results in this section provide an investigation into this question, a topic which will be discussed in more detail within the analytic SDRG in Sec. V. It is not immediately clear to what extent the locally-correlated free fermion effective model defined in Eqs. (36) and (37) shares the qualitative features of the XYZ model, or indeed the self-consistent mean field theory. We investigate this by repeating the measurements of bulk and end-to-end spin correlations in chains of similar length to the previous studies, now varying the coupling correlation parameter δ. Figures 14, 15, and 16 demonstrate that these critical indices do vary continuously in a similar way to the interacting case. Our observation that this mean-field approach indeed exhibits  Figure 16.

Exact diagonalization study of locally correlated Majorana chains
Critical exponents governing spin correlations in the locally-correlated XY model with varying correlation parameter δ are shown, extracted from data shown in Figs. 14 and 15. Both bulk and end-to-end exponents are shown, with known results for the bulk correlations in the uncorrelated XY model at δ = 0 indicated by red stars and known end-to-end critical spin exponents by yellow diamonds. Known critical exponents for the U(1)-symmetric XX model at δ = 1 are similarly indicated; in this case η e ⊥ = η e z = 1 and η ⊥ = ηz = 2. The discrepancy in η ⊥ is likely a result of a long crossover, as the disorder distribution of Eqs. (36) and (37) is somewhat weaker than Eq. (31) for the same value Γ = 2. many of the qualitative features of the original case suggests that at least for smallJ z , the primary effect of the interactions is to correlate the coefficients of the hopping terms on the two Majorana chains. However, we emphasize that although the η z and η ⊥ converge to similar values at the XX point δ = 1 and the tricritical XYZ pointJ z = 1, the reasons for this are not necessarily the same. The mean field should not be taken too seri-ously as a picture of the interacting phase away from the perturbative regime.

V. LOCALLY CORRELATED XY MODEL IN THE RANDOM WALK FORMALISM
Some types of disordered quantum Hamiltonian can be uniquely associated with a classical random walk (RW). An alternative picture of the SDRG viewed through this connection is useful for understanding the properties of IRFP phases. The RW formulation has previously been applied to both the RTFIM [54,55] and AFM quantum spin chains [30,56]. In this section we first review the RW for a single Majorana chain based on the SDRG procedure of Sec. II C 1. While all results for correlation functions in this case are known from Fisher's analytic solutions for flows approaching the RS fixed point, we demonstrate how to obtain some power law exponents from different arguments, which will generalize to the locally correlated XY chain where we do not have analytic flows. We first obtain rigorous bounds in the continuum limit on the asymptotic scaling of the Majorana pairing probability (which determines the correlations of the z component of spin in the random XX and XY chains) based on RW survival probability, a connection which had previously been noted in Ref. [30]. We then consider the problem of two locally correlated RWs, one for each Majorana chain, following the effective model developed in Sec. IV B. This system turns out to correspond to an anisotropic two-dimensional RW. We again rigorously bound the likelihood of decimation using the RW survival probability, where we find that the power law exponent varies continuously with the local correlation parameter. As a result, we are able to prove a specific form for continuously varying critical exponents of spin correlations in the locally-correlated effective model.

A. RW formulation of SDRG for the Majorana chain
Returning to the notation of Sec. II C 1, define the logarithm of the energy associated with each bond in the Majorana chain Hamiltonian H M as u n = ln(J/|h n |), n = 1, . . . , N −1. HereJ is a bare bandwidth for the coupling terms, meant to evoke the parameters of the Hamiltonian Eq. (1). From Eq. (5) one sees that ifJ x =J y , in each Majorana chain of the random XY model the hopping terms are identically distributed. Note that the signs of h n are not important for the discussion of probabilities of site pairings below, and are only needed to fix sign factors for the spin correlation functions, as discussed at the end of Sec. II C 1. We consider the specific disorder distribution Eq. (31) withJ x =J y =J = 1. Then the distribution of log-energies is exponential, with distribution parameter Γ: which has mean u = Γ and variance Var(u) = Γ 2 . The Majorana model H M on N sites is associated with a 1d RW m, a Markov chain with state variables (x n , σ n ), n = 1, . . . , N , where x n ∈ R is a cumulative log-energy defined below and σ n = (−1) n−1 is an internal Z 2 variable determining the sign of the next step to be taken [57]. The discrete RW time n matches the spatial index of the quantum chain. A given disorder realization {h j } 1≤j<N corresponds to a RW step sequence {σ j u j } 1≤j<N : that is, the state of m at time n = 1, . . . , N is In the following we will sometimes leave the σ n state variable implicit, and refer to x n as m[n]. Let Prob(x, σ, n) be the distribution of m[n], which is governed by the master equation (40) We now consider the behavior under the SDRG of a RW m associated with a Majorana chain H M . The largest local energy scale |h k |, for some k, corresponds to the smallest log-energy u k . The effect of the Shreiffer-Wolff transformation up to second order is to eliminate the following hopping terms: and to introduce the renormalized bond term (There is also a shift of the leading energy scale, but this will not be important here.) For the RW the new step is In this way the SDRG transformation corresponds to a sequential "smoothing" of the RW, in which the global step of smallest magnitude and its neighbors are removed, and replaced by a treble step directly connecting x k−1 and x k+2 . For an illustration, the reader is referred to Fig. 8 in App. B of the arXiv version of Ref. [56], or Fig. 1 of Ref. [11]. We define an inversion operation I acting on a RW m of length N as That is, I flips the spatial and time coordinates of m.
(The constant shifts the starting point of Im to 0.) We also define reflection R a of the spatial coordinate about the line x = a: We will make extensive use of a "gluing" operation ⊕ which joins two RWs at their endpoints. For RWs m 1,2 with lengths N 1,2 , then, n = 1, . . . , N 1 + N 2 , (46) That is, the combined RW m 1 ⊕ m 2 first performs the N 1 − 1 steps of m 1 , followed by the N 2 − 1 steps of m 2 . It is assumed that the first step of m 2 has opposite σ state variable as compared to the last step of m 1 ; this is required on the spin chain, where m 2 begins and m 1 ends on the same sublattice.
Using the above definitions a precise statement can be made about the decimation of a site n = k, which we suppose without loss of generality to be a local minimum. For k to have decimation partner k > k in the SDRG, with k − k = r, a RW m must admit a decomposition where m ext,L has length k, m ext,R has length N − k + 1, (For a pictorial description, see also App. B of the arXiv version of Ref. [56].) These conditions relate the likelihood of a decimation pairing sites k and k to the survival probability of the "interior" and "exterior" partial RWs on the fully bounded interval (0, ∆). The physical interest of this quantity follows from the strong correlations shared by sites paired in the SDRG; in particular, the scaling of the decimation probability determines average spin correlations, as described in Sec. II D.
Note that the writing of Eq. (47) is chosen so that the exterior RWs m ext,L and m ext,R have identical structure to the interior RW m int . That is, all walks evolve forward in time starting at step 1 with the first step being positive. Implicit in this is the assumption that the inversion and reflection operations used result in identical probabilities for the RWs because the microscopic distributions for u n are identical for n even and odd.
Focusing on asymptotic scaling (i.e., n, r 1), we describe the RW in continuous time, passing from n → t. The central limit theorem specifies that a sum of random variables approaches a Gaussian distribution for sufficiently large n, provided only that the moments of the constituent distributions are bounded. The variance of the continuum distribution is Var(x) = Var(u)t. The effect of the internal state variable σ can be accounted for by noting that sites which decimate together necessarily inhabit distinct sublattices. This means that one additional σ = +1 step is always taken. The mean of the probability distribution, then, is the expectation value of this step: x ≡ x 0 = u [58]. The asymptotic density in free space we denote by Now the continuum limit of Eq. (40) is the diffusion equation [59] with diffusion constant D = Var(u)/2. Eq. (48) is the Green's function of Eq. (49) on x ∈ R with initial condition G(x, t = 0) = δ(x − x 0 ). This illustrates that the continuum limit of the RW can be treated as a diffusing particle initially localized at x = x 0 . Accordingly, in the following sections we use the language of the diffusion problem, referring to the counterparts of discrete RWs associated with particular Majorana Hamiltonians as "paths," "histories," or "trajectories." We also sometimes write the initial condition explicitly, as G(x, t; x 0 ). Finally, we will use the notation defined in this section for the discrete case, e.g., I, R a , and ⊕, to also refer to the counterparts of these operations in the continuum.

B. Rigorous bounds on critical exponents in the Majorana chain from RW survival
The diffusion equation on the fully bounded interval (0, ∆), i.e., with absorbing boundary conditions at x = 0 and x = ∆, can be solved straightforwardly by harmonic expansion. From the time-dependent solution one can directly calculate the scaling of the asymptotic decimation probability and reproduce Fisher's detailed results in Refs. [3,7]. However, in Sec. V D the fully bounded geometry for two locally correlated Majorana chains becomes too complicated to solve this way. Instead we employ a different approach by proving upper and lower bounds with the same power-law scaling, based on the survival probability in a semi-infinite domain. A similar method will work also for the locally correlated effective model with an arbitrary degree of correlation.
First consider the survival probability of a RW in the semi-infinite interval at time t > 0. As in the free case Eq. (48), the initial condition on the constrained density G(x, t) is G(x, t = 0) = δ(x − x 0 ), but an absorbing boundary is present at x = 0, restricting the solution domain to x ∈ (0, ∞) and terminating trajectories that reach x = 0. The boundary condition G(x = 0, t) = 0 is accounted for by placing an "image charge" at x = −x 0 and superposing the distributions: G(x, t) = G free (x, t; x 0 ) − G free (x, t; −x 0 ). We generally work in a "scaling limit," where assuming in the last line x 0 √ Dt. This approximation is valid at late times in integrals over the spatial coordinate, as the exponential factor strongly mitigates the error introduced, and allows us to extract leading power-law behaviors. The survival probability in the semi-infinite geometry in the scaling limit is 1. End-to-end decimation probability for a single finite Majorana chain In order to support end-to-end decimation between sites 1 and N , the RW m[n = N ] associated with a finite Majorana chain of length N need only satisfy Condition 1 of the previous section, with r = N . In the continuum limit for the RW (N → L), the likelihood that the left end t = 0 is involved in the final decimation is given by the survival probability S(t = L) ∼ 1/ √ L; however Condition 1 additionally requires that its decimation partner be the right end t = L. Applying I to m, one sees that the requirement to reach a maximum at t = L takes the same form as the absorbing boundary condition x = 0 near t = 0. Thus a naive estimate of the end-to-end decimation probability p e (L) is the independent survival of the two ends, or S(L) 2 ∼ 1/L. Although these events are not actually independent, we will show that the naive estimate turns out to give the correct scaling. Some intuition for this is that surviving histories tend to be located increasingly far away from the absorbing boundary [60]: consequently, the "special" low-probability behavior is confined to the neighborhood of the ends, while the middle of the RW can be allowed to be nearly typical. A precise statement of these schematic remarks is that we are able to determine the scaling of p e (L) by considering two independent "half-RWs" m 1,2 of length t = L/2, constructing RWs of length L which satisfy Condition 1 as m = m 1 ⊕ Im 2 .
To be more concrete, we first give a rigorous upper bound on the end-to-end decimation probability p e (L). Any RW m can be decomposed as m = m 1 ⊕ Im 2 , that is, into two independent "half-RWs" running up to time t = L/2, one running over times t ∈ [0, L/2], and the other over t ∈ [L/2, L], with the two RWs properly glued at their respective time t = L/2. It may be the case that m 1 and m 2 never reach the absorbing boundary, and thus each is considered a surviving RW in the semi-infinite geometry. Any RW instance of length L producing an end-to-end pairing in the SDRG, i.e., satisfying Condition 1, indeed decomposes in this way, with only one absorbing boundary in each case. The converse statement is not true, because when such two surviving trajectories are joined, we cannot guarantee that the full RW satisfies Condition 1. Thus, the desired probability p e (L) ≤ S(L/2) 2 ∼ 1/L.
To prove a lower bound on p e we construct a subset of all paths satisfying Condition 1 by considering certain m 1 and m 2 , each of length t = L/2, which when glued together as m 1 ⊕ Im 2 satisfy the criterion. Again, in the present case we can solve the problem with two absorbing boundaries, but we want to demonstrate how to extract the behavior using the semi-infinite solution, where the geometry is simpler, as this will be the only option for the locally correlated model. Specify constants α and β, 0 < α < β ≤ 2α, and define a target window x ∈ [α √ Dt, β √ Dt] for a time t > 0. In the problem with one absorbing boundary at x = 0, the fraction of surviving trajectories contained in the target window at t is That is, a constant fraction p w (α, β) of the surviving density of RWs at time t is located within the target window. The above calculation Eq. (53) leads to an overcounting of valid paths which can be glued to satisfy Condition 1, because it includes "dangerous" histories which take an excursion to large x values before returning to the target window at time t. Half-RWs m 1 and m 2 constrained in this way and glued as m 1 ⊕ Im 2 may cross the eventual decimation log-energy scale ∆ prematurely, which would spoil the lower bound. To account for the dangerous cases, we exclude those histories which ever cross x = β √ Dt and then return to the target window. The way we achieve the exclusion is the following. Suppose that a history m[t ], t ∈ [0, t], performs q crossings of the line x = β √ Dt at times {t 1 , t 2 , . . . , t q } before returning to the target window at t = t. Immediately after t q , the history must travel downwards and remain below x = β √ Dt until t = t. We apply the following transformation: (54) where as indicated by the subscripts m t ≤tq describes the RW up to time t = t q and m t >tq the section t ∈ (t q , t]. T does not change the earlier partial RW but reflects the later about the line . Moreover, the likelihood of the trajectory is unaffected by T. Now every dangerous path with q ≥ 1 crossings can be identified with a transformed partner terminating in the shadow window and having the same probability. Thus the density in the shadow window at time t upper bounds the contribution to the density in the target window arising from dangerous histories. (The upper bound is not saturated, because a trajectory included in the shadow window could deviate above x = 2β √ Dt for some t ∈ (t q , t], and this RW would have no T −1 counterpart due to the absorbing boundary at x = 0.) From the previous calculation, the fraction of the surviving density contained in the shadow window is p sw (α, β) = e −β 2 /4 − e −(2β−α) 2 /4 . Consequently a lower bound on the density of valid surviving histories in the target window at time t is given by There is an extended region of (α, β) for which the coefficient is positive; for example, p corr w (α = 2, β = 4) ≈ 0.33. Figure 17. A dangerous trajectory contributing to the counting pw of the density in the target window, colored in blue, is illustrated. The shadow window used to eliminate these trajectories is also shown, colored in orange. The particular history m shown has q = 4 crossings of the upper limit of the target window and the reflected partial path R β √ Dt m t >tq , terminating in the shadow window, is shown in green. Because the diffusion is unbiased, both m and the transformed Tm path have the same probability, and as any such dangerous trajectory has a counterpart under the transformation, the density in the shadow window upper-bounds the associated contribution to the density in the target window. Now take t = L/2. Two RWs m 1 and m 2 fulfilling the criteria above are suitable for constructing a RW of length L which satisfies Condition 1 as m = m 1 ⊕ Im 2 . The result is a trajectory of length L reaching a maximum at t = L (assured by taking β ≤ 2α) without crossing x = 0. Not all RWs of length L which support end-toend decimation in the SDRG can be constructed this way, only those with m[L/2] lying in the target window and m[t ≤ L/2] below the upper limit of the target window, but every RW coming from this construction evidently satisfies Condition 1. Thus this probability is a lower bound on p e (L) ≥ [p corr w (α, β)S(L/2)] 2 ∼ 1/L. Together with the upper bound, this establishes the scaling of end-to-end decimation probability p e (L)-and thus the power law for end-to-end correlations in a single random Majorana chain-as 1/L.

Bulk decimation probability in a single Majorana chain
Guaranteeing decimation away from the edges of a Majorana chain requires satisfying both Conditions 1 and 2. To find the probability p b (r) of decimation at scale r in the bulk-i.e., that two fixed sites separated by r are decimated as a pair-we decorate interior RWs m int by gluing exterior RWs to the left and right. We showed that the probability of such an m int is p e (L = r) ∼ 1/r, so we need only find suitable exterior RWs satisfying Condition 2 (while bearing in mind conditions involving both interior and exterior RWs).
For the probabilities associated with the exterior walks, we are interested in the likelihood ω(x; A) that a RW with spatial coordinate x starting from x = x ≥ 0 eventually reaches a value x = A before being absorbed at the domain boundary x = 0. We require the consistency condition ω(x; A) = ω(x − dx; A) , where the average is taken over sufficiently small displacements dx, and dx = 0, (dx) 2 = 0 (reflective of the microscopic step distribution) [60,61]. Taylor expanding leads to Laplace's equation ∇ 2 ω = 0 which, together with the boundary conditions ω(0) = 0 and ω(A) = 1, has solution ω(x; A) = x/A.
A lower bound on p b (r) is now straightforward based on m int as defined in Sec. V B 1, coming from a subset of all RWs of length L = r supporting end-to-end decimation. Any such m int is constructed from two glued half-RWs, each terminating at t = r/2 inside of a target window x ∈ [α Dr/2, β Dr/2]; thus the total and maximum deviation at t = r is bounded above by ∆(r) = β √ 2Dr. Given m int , the probability of a suitable exterior RW m ext,L or m ext,R is greater than or equal to ω(x 0 ; ∆(r)); writing a full RW satisfying all conditions, we find that p b (r) ≥ [p corr w (α, β)S(r/2)] 2 ω(x 0 ; ∆(r)) 2 ∼ r −2 .
In the same spirit as the upper bound on end-to-end decimation probability, consider m int = m 1 ⊕ Im 2 ; that is, decomposed as two half-RWs surviving until t = r/2, with final spatial deviations ∆ 1 and ∆ 2 and likelihoods G(∆ 1 , r/2; x 0 ) and G(∆ 2 , r/2; x 0 ), respectively. All RWs with end-to-end decimation are of this form. Now incorporating the probability of exterior RWs which must reach a height ∆ 1 + ∆ 2 , the likelihood of the full RW provides an upper bound on the probability of bulk decimation: Making use of ω(x 0 ; ∆ 1 + ∆ 2 ) 2 ≤ 1 2 ω(x 0 ; ∆ 1 )ω(x 0 ; ∆ 2 ) the integrals factorize, and we find Again these upper and lower bounds exhibit the same scaling, proving that p b (r) ∼ r −2 for a single Majorana chain, in agreement with known results (see the XX case in Sec. II D).

C. Locally-correlated Majorana chains as a two-dimensional RW
To make statements about locally correlated Majorana chains requires dealing simultaneously with two RWs (returning for the moment to the discrete formulation) m x [n] and m y [n], associated respectively with the X and Y Majorana hopping chains. In the general case, the steps taken by each at time n are not independent, being instead drawn from a joint distribution µ(u, v). If the full state of the system is specified by variables (x n , y n , n), the master equation for the probability distribution Prob(x, y, n) is (60) This is however just the master equation for a RW in two dimensions (2d). In the natural 2d vector notation with x = (x, y) and u = (u, v) , The continuum limit of the master equation Eq. (61) is determined by the details of the microscopic distribution µ, and does not in general reduce to the simple Laplacian. As a remedy we begin by transforming the problem into isotropic diffusion.
Let µ be centered, with covariance matrix [62] where corr(u, v) = cov(u, v)/σ 2 ≡ δ ∈ [0, 1], with fixed σ 2 = Var(u) = Var(v). (The value of δ here is related to, but not necessarily the same as, the bare δ defined in Sec. IV B. δ > 0 implies positive correlation between u and v, as observed in the mean field for the AFM spin chain.) The continuum limit of evolution driven by µ is anisotropic diffusion along the eigenvectors of Σ,ê ± = 1 √ 2 (1, ±1) , with diffusion coefficients D ± = σ 2 2 (1 ± δ). The 2d RW evolves by isotropic diffusion under a linear transformation of the plane W : W performs a rotation about the origin by π/4, followed by a δ-dependent anisotropic rescaling. There is a divergence at δ = 1, where Σ is rank-deficient; this reflects the fundamentally one-dimensional nature of the perfectly correlated case. We will refer to the (x, y) coordinates of the original problem as the "physical geometry," and the image (x,ỹ) of W as the "solution geometry," where the governing equation is isotropic diffusion, now with D. Rigorous bounds on critical exponents in the locally correlated model

End-to-end decimation probability for two locally correlated finite Majorana chains
Investigating end-to-end decimation directly in the exact solution for the fully bounded geometry would necessitate solving Eq. (64) in a parallelogram. A harmonic decomposition is not possible here, and as far as we are aware the solution requires a prohibitively complicated Schwarz-Christoffel conformal transformation usually performed numerically [63]. Nevertheless, analytic results for two Majorana chains with arbitrary local correlations are possible by utilizing the connection to the survival probability in the simpler semi-infinite geometry.
As was the case for the single Majorana chain, we employ a semi-infinite domain, now bounded by the lines x = 0 and y = 0. The origin is evidently fixed by W, and the boundaries map to the linesỹ = ±λ 2x , wherex lies in theê − direction andỹ inê + . These boundaries delimit an absorbing wedge geometry with opening angle Θ given by cos Θ = −δ, which runs from Θ = π/2 at δ = 0 to Θ = π at δ = 1. In terms of the wedge halfangle θ ≡ Θ/2, the domain boundaries areỹ = ±(cot θ)x. For easy reference, we collect some relationships between these geometric parameters: The Green's function in the infinite wedge can be found from the free-space distribution by the method of images for opening angles Θ = π/m, with m a positive integer. This entails 2m − 1 image charges with alternating sign, arranged symmetrically around the wedge apex. However this approach is of limited use, as we need Θ ∈ [ π 2 , π), and instead we will use the Green's function known for arbitrary opening angle from an alternative solution. In polar coordinates, with the wedge apex at radius ρ = 0 and solution domain bounded by absorbing walls G(ρ, φ = 0, t) = G(ρ, φ = Θ, t) = 0 (i.e., the angle φ is defined relative to one of the absorbing boundaries), we have [64] G(ρ, φ, t; ρ 0 , φ 0 ) = e −(ρ 2 +ρ 2 where ν = π/Θ and I lν is a modified Bessel function of the first kind: In the physical geometry the initial condition is (x 0 , y 0 ) = ( u , v ), where u = v is again the result of each 1d RW taking one additional positive step according to the discrete microscopic distribution. In the solution geometry this point maps to ρ 0ê+ , where ρ 0 = √ 2λ u . In polar coordinates the source point is (ρ 0 , φ 0 = θ). Consequently, in Eq. (67) the factor sin(lνφ 0 ) vanishes for even l and for odd l is equal to a sign (−1) (l−1)/2 . As in the 1d case, we work in the scaling regime at late times t, where we are able to extract the leading power-law behavior. Again, spatial integrals are regulated by the exponential factor, which decays fast enough to suppress errors arising at large ρ. Because ν ∈ (1, 2] the leading behavior requires only the l = 1, m = 0 term in the double sum, and sets e −ρ 2 0 /4Dt → 1. The survival probability is determined from the Green's function by integration over the wedge. Explicitly, in the scaling limit The survival exponent depends on the opening angle as This result for a RW in a 2d wedge is in fact well known [60,61,65]. As Θ is a function of the correlation coefficient δ, continuously varying behavior of this type is in agreement with the numerical observations in Sec. IV B. Specifically, again relying on the naive assumption that the two ends of the chain decimate independently, the likelihood of this pairing scales as [S(L)] 2 ∼ L −π/Θ , which matches the known end-to-end scaling exponents η e z = 2 for the uncorrelated model at δ = 0 and η e z = 1 for δ = 1.
Our strategy for rigorously bounding the probability of end-to-end decimation occurring on both chains using the infinite wedge results is analogous to that of Sec. V B. From the Green's function we establish that at late times a constant fraction of surviving RWs are suitable for subsequent gluing to contribute to this probability, being found in a specified target window, using a shadow window to exclude dangerous trajectories. By gluing the ends of two RWs at time t = L/2 we establish bounds on the power law. We will use the notation of the previous section, namely I and ⊕, to refer to the generalizations of the relevant transformations to 2d.
In particular, we can write an upper bound immediately. Any 2d RW of length (duration) L corresponding to two locally correlated Majorana chains can be decomposed into half-chains of length L/2 as m = m 1 ⊕ Im 2 , as in the 1d case. m 1 and m 2 may be valid surviving trajectories in their semi-infinite wedge, and some will produce end-to-end decimations on both physical Majorana chains described by the 2d RW m. Trajectories that do not decompose in this way into surviving half-chains will not satisfy Condition 1. Because not every pair of surviving m 1 and m 2 will do so either, the probability is upper-bounded as p e (L) ≤ S(L/2) 2 ∼ L −π/Θ . Now in order to prove a lower bound on p e (L), let α and β be positive constants, α < β ≤ 2α, and define the target window for a 2d RW at time t to be the square In the physical geometry the window is a square; however, when mapped to the solution geometry the window becomes a parallelogram. The corners {a, b, c, d} map to as illustrated in Fig. 18. Treating this exact shape in the polar coordinates of Eq. (67) is complicated; instead we define an integration volume that is a subset of the target window, with the same t scaling, but which leads to a simpler bound. Consider the midpoints of the edges of the target window in the solution geometry, which we denote {ẽ,f ,g,h}. They describe the four corners of a rectangle, symmetric about the line φ = θ, with edges in the directionsê − andê + (see Fig. 18). We define an integration domain bounded by radial values ρ + (of pointsf andh) and ρ − (ofẽ andg), and the angular deviation ψ of pointsf andh from the midline φ = θ.
The proof that this "sector" geometry is indeed a subvolume of the target domain for any opening angle Θ < π can be seen by drawing a picture. The specific integration bounds can be found straightforwardly from Eq. (73), but the crucial property is their scaling with t. Define the radial limits as ρ ± = C ± (α, β, δ) √ Dt; the angular integration half-width ψ = ψ(α, β, δ) turns out to be purely geometric, with no t dependence. Again extracting the leading behavior for late times t, the fraction of surviving paths whose position at time t is in the integration window is where I(α, β, δ) = So p 2d w is indeed a constant, determined only by the correlation coefficient δ and the constants α and β.
As was the case for the 1d RW, the calculation above includes a "dangerous" contribution which should be subtracted in order to lower-bound the decimation probability by subsequent gluing of half-chains m 1 and m 2 . Again we upper-bound this contribution by calculating the fraction in a shadow window. We consider those paths to be dangerous which ever cross the lines x = β √ Dt or y = β √ Dt in the physical space before returning to the target window at time t. In the solution geometry these lines map to We define the boundary for dangerous trajectories piecewise as (see Fig. 18) Suppose a trajectory with time parameter t makes q crossings of D at times {t 1 , . . . , t q } at various points {(ρ 1 , φ 1 ), . . . , (ρ q , φ q )} before returning to the target window at time t = t. After its last crossing at (ρ q , φ q ), it must stay within the allowed region for times (t q , t]. We transform the trajectory by reflecting the partial RW for times t ∈ (t q , t] about the component of D that was crossed at t = t q , either D R if φ q ∈ (0, θ] or D L if φ q ∈ (θ, Θ). This is the counterpart in 2d to the 1d transformation T. Because the step distribution in the solution geometry is isotropic, the transformed path has the same probability as the dangerous original. (The reflection must be performed in the solution geometry, and does not commute with W.) The shadow window in this case has two components, which are disconnected for Θ < 2π 3 but overlap for Θ >  Figure 18. The solution geometry is illustrated for the 2d RW problem in the wedge with opening angle Θ, found from the correlation coefficient by cos Θ = −δ. The exact target window is colored in blue, and the sector defining the easier integration subregion for the target in yellow. The two components of the shadow window are found by reflecting the exact target window across the lines DL and DR and are colored in orange, with the easier bounding shadow integration region, which necessarily covers these areas, in green.
andb about D R , andã andc about D L . The coordinates of the points reflected about D R are with similar forms forã L andc L . The four-sided figures described by the exact shadow window are evidently complicated. As with the target window, we bound the area using a sector which scales in the same way, however in this case an upper bound is required. The upper limit ρ sw + is the radial coordinate of pointsc L andb R , and the lower limit ρ sw − is that shared by the cornersb andc. The angular half-width is the maximum of the angular half-widths of pointsc andã R ; this depends on the specific value of Θ. Again we find integration limits ρ sw ± = C sw ± (α, β, δ) √ Dt, and ψ sw = ψ sw (α, β, δ).

Bulk decimation probability in two locally correlated Majorana chains
Once again we can extend the result for end-to-end decimation p e -requiring that both Majorana chains satisfy Condition 1-to the bulk likelihood p b (r) (for two fixed spins separated by r) by considering also Condition 2. We first write a lower bound on the bulk pair decimation probability by identifying exterior RWs which are guaranteed to satisfy Condition 2 when properly adjoined to an interior RW of the type used for the lower bound on p e in the previous section. Specifically, we restrict to exterior RWs with endpoints at time t ≡ r (for concreteness, but any constant multiple of r would do as well) within a particular sector (specified below) in the solution geometry. In the physical geometry, ∆(r) = β √ 2Dr is an upper bound on the total deviation of each of the 1d RWs m x and m y described by the 2d interior RW m int .
One way to guarantee the bulk decimation is to require that each of the physical 1d RWs described by each of the exterior 2d RWs m ext,L and m ext,R survive, and exceed ∆(r) at t = r. A point (ρ, φ) in the solution geometry corresponds to in the physical geometry. Employing angular integration limits φ ∈ (θ − ψ, θ + ψ), where ψ can be chosen to be the same value used for m int , sufficient radial limits for our purposes are ρ ext − = ∆(r) √ sin Θ/ sin(θ−ψ) and ρ ext + → ∞ (noticing that sin(θ − ψ) ≤ sin(θ + ψ) for all ψ ∈ [0, θ]). From the calculation of the previous section there is a constant probability κ(α, β, δ) that any surviving RW lies in a window bounded by ρ ∈ [ρ ext − , ρ ext + ] and φ ∈ [θ − ψ, θ + ψ] at t = r. Such a RW has deviation at least ∆(r) in the physical x and y coordinates and thus as either m ext,L or m ext,R is suitable for satisfying Condition 2 for bulk decimation when properly adjoined to m int as constructed previously; thus p b (r) ≥ p e (r)[κS(t = r)] 2 ∼ r −2π/Θ . Similar to the case of a single Majorana chain, for an upper bound we make use of the probability ω(ρ, φ; A) of a RW with spatial coordinates (ρ , φ ) reaching radius ρ = A in the wedge given a starting point (ρ, φ). This probability follows Laplace's equation ∇ 2 ω = 0, now with boundary conditions ω(ρ, φ = 0) = ω(ρ, φ = Θ) = 0, ω(ρ = A, φ) = 1. Assuming a separable solution ω(ρ, φ) = R(ρ)T (φ), we find that for the angular coordinate the solutions are T n (φ) = sin(nνφ), n = 1, 2, 3, . . . , where as before ν = π/Θ. For the radial coordinate which has solutions of the form R n (ρ) = ρ ±nν . Determining the constants from the boundary conditions, Along the relevant line φ = θ, the probability simplifies to In order to write an upper bound on the bulk decimation probability, we consider a full RW satisfying both Conditions assembled from an m int = m 1 ⊕ Im 2 , where each of m 1 and m 2 must survive until t ≡ r/2, along with exterior RWs m ext,L and m ext,R which must reach a particular radial coordinate (determined from m int as specified below) without being absorbed. Suppose that m 1 and m 2 terminate at coordinates (ρ 1 , φ 1 ) and (ρ 2 , φ 2 ), which define the deviations of the physical RWs (∆ x,1 , ∆ y,1 ), and (∆ x,2 , ∆ y,2 ) according to Eq. (84). The full deviation of the interior walk m int in the physical coordinates is (∆ x , ∆ y ) = (∆ x,1 + ∆ x,2 , ∆ y,1 + ∆ y,2 ) and the physical 1d RWs described by m ext,L and m ext,R must exceed the corresponding ∆ x or ∆ y before being absorbed. For this to be the case it is necessary, but not sufficient, that the 2d exterior RWs each survive in the wedge until reaching radial coordinate A ≡ √ sin Θ min(∆ x , ∆ y ) in the solution geometry. Defining for m 1 and m 2 similar A 1 ≡ √ sin Θ min(∆ x,1 , ∆ y,1 ) and A 2 ≡ √ sin Θ min(∆ x,2 , ∆ y,2 ), we note that A ≥ A 1 , A 2 . The probability of finding two such m ext,L and m ext,R given the terminating locations of m 1 and m 2 is Then, integrating over the distribution of the interior half-chain coordinates, p b (r) = ρ 1 dρ 1 dφ 1 G(ρ 1 , φ 1 , r/2; ρ 0 , θ)× ρ 2 dρ 2 dφ 2 G(ρ 2 , φ 2 , r/2; ρ 0 , θ)× We restrict to the right half-wedge, as the integrand is symmetric about φ = θ. The angular integral is which converges for Θ > π/2, equivalently δ > 0. (The exponent we are bounding is known at δ = 0, and follows from the result of Sec. V B 2.) Combining the upper and lower bounds, we prove that p b (r) ∼ r −2ν , and the bulk correlations exponent for two locally correlated Majorana chains with parameter δ is E. Numerical SDRG study The final results of this section, Eqs. (83) and (95), are in qualitative agreement with the quantum simulations of Sec. IV B for relatively short Majorana chains, and are consistent with previously-known results at the points δ = 0, 1, where the locally-correlated model describes the random uncorrelated XY and perfectly correlated XX IRFPs. For further verification we implement the SDRG update Eq. (19) directly for two Majorana chains with locally-correlated terms, and are able to access larger system sizes. This also allows us to study the bulk C ⊥ (r) power laws, which are not analytically tractable in the mapping to RWs used in the preceding subsections.
The numerically extracted exponents are shown in Fig. 19. The bare correlation coefficient δ may become slightly renormalized from the lattice scale definition in Eqs. (36)-(37) compared to the meaning in the continuum 2d RW treatment in Sec. V D, but these simulations are in good agreement with the analytic forms for η e z (δ) and η z (δ) we obtained. While we have precise analytical knowledge only of the critical exponents η e z and η z , we observe that η ⊥ also varies continuously. In contrast, η e ⊥ = 1 for any value of δ, by the argument presented in Sec. II D.  Figure 19. Numerical SDRG data are shown for two locallycorrelated Majorana chains, with the end-to-end and bulk decimation probability exponents-equivalent to η e z and ηz, respectively, in the quantum model-compared to the analytic forms Eqs. (83) and (95) (dashed lines). Also shown are critical exponents η e ⊥ and η ⊥ measured in the numerical SDRG, as well as red stars indicating known values of bulk correlations exponents at δ = 0 and 1, and yellow diamonds indicating known values of end-to-end correlations exponents. The end-to-end correlations data were taken from 1 000 000 disorder realizations each for system sizes up to N = 128, and the bulk correlations data were taken from 100 000 disorder realizations at system size N = 256, utilizing the middle half of each of the two Majorana chains.

VI. FIXED POINTS FOR THE INTERACTING MODEL
In Sec. V we performed a study of the behavior of critical exponents under a varying degree of local correlations in a random free-fermion model. Despite the lack of tractable SDRG flow equations, we showed that the local correlation controlled by δ is a marginal perturbation which tunes along a line of IRFPs. In the present section we advance the perspective that this line of non-interacting fixed points in fact also controls the long-distance behavior of the interacting model for small J z strength below the transition to the z-AFM phase.
To do so requires a study of the SDRG at intermediate stages, taking into account more general terms produced by the interactions. Equation (29) describes the result of an initial decimation, but eventually descendant terms will be frequent and must also be taken into account. We change our conventions here from those of Sec. II C 2 for convenience: namely, we denote the Majorana chains by I, II rather than X , Y; and by a gauge transformation (described at the end of Sec. II C 1) we set the signs of h I n , h II n > 0 for all n = 1, . . . , N − 1, and K n ≡ K n,n < 0. In order to capture the effect of iterated decimations, we observe that in Eq. (29) descendants of the form K n,m (iγ I n γ I n+1 )(iγ II m γ II m+1 ) are produced, which generalize the K n of Eq. (22). We enlarge the space of couplings to include all such terms, with initial distribution K n,m = 0, n = m. If the average K ≡ | K n,m | can be considered to be a small parameter (for weak interactions K < | h n |), the higher-fermion term in Eq. (29) appears at order O(K 2 ) and can thus be neglected. We will demonstrate that the space of couplings including all K n,m is closed under RG flow up to O(K), and that the structure of the signs is preserved. Furthermore, we will show that the strength of the K terms decreases in some sense relative to the h terms, suggesting that interactions are irrelevant, at least in the neighborhood of the free-fermion fixed point.
Following the approach of Sec. II C, denote the largest term as H 0 = ih I k γ I k γ I k+1 and associate with the eigenstates of this term a complex fermion f † 0 = 1 2 (γ I k + iγ I k+1 ) with projectors π + = f 0 f † 0 and π − = f † 0 f 0 into the evenand odd-parity sectors, or the high-and low-energy eigenstates, respectively, of H 0 . The off-diagonal terms in the Schrieffer-Wolff treatment share exactly one Majorana operator with H 0 : Separating V od into symmetry sectors, we find that We make use of the "interacting couplings" notation used also in Sec. II C 2 to connect with the non-interacting case, but here it is not evident that these couplingswhich are really operators-all commute. Nevertheless, a suitably generalized version of Eq. (19) implements the Schrieffer-Wolff transformation: The effective terms in the first line of Eq. (103) (and the first term of the second line) are h-type, with positive coefficients in the low-energy sector of H 0 where iγ I k γ I k+1 = −1. Conversely, the remaining terms in the second line are K-type (recalling that γ I k−1 and γ I k+2 become adjacent after the decimation of γ I k and γ I k+1 ), and have coefficients with negative signs. One sees that the signs of the initial distributions, namely h I,II n > 0 and K n,m < 0, are maintained during the RG flow, and it is evident from Eq. (103) that these types of terms are closed under the SDRG up to O(K).
As a measure of the evolution of the relative strength of K terms to h terms under this RG step, we compare the renormalized K eff k−1,m to the geometric mean of the proximate h terms h I,eff k−1 and h II m : This suggests that if the h terms are dominant initially, they will be even more so during the SDRG and will asymptotically constitute the entirety of the decimations.
The diagonal terms which contain both decimated Majoranas are Upon decimation, setting iγ I k γ I k+1 = −1 in the ground state gives O(K) contributions to the Majorana hoppings in the other chain, h II,eff Given the opposite signs of the h and K couplings, this increases the overall strength of the remaining Majorana hoppings. This is the local SDRG analog of the "mean field" of Eqs. (34) and (35) where the J z interactions renormalize the J x and J y couplings by strengthening and correlating them, as was already noted in Sec. II C 2 and discussed in Sec. IV. Here we note that including these renormalizations of the h couplings only improves our arguments for the persistence of the dominance of these couplings over the K couplings.
The terms omitted from Eq. (103) at O(K 2 ) are the following: The first terms in each line are corrections to the groundstate energy and the strength of the renormalized bond coupling on chain I [which again preserves the sign struc-ture and strengthens this hopping compared to the leading contribution in Eq. (103)]. Along with these, fourfermion terms within chain II and six-fermion inter-chain terms appear at O(K 2 ). The former are expected to be ultimately irrelevant, based on previous studies of a single Majorana chain realized in the quantum Ising model [3]. However these four-fermion terms and the six-fermion terms will produce yet more complicated descendants in subsequent RG steps, and there will also be "degradation" processes leading to fewer-fermion terms, including renormalization of the two-fermion terms, similar to the discussion after Eq. (105) [22]. In this case we must rely on the perturbative argument to justify dropping them, viewing them as irrelevant other than feeding into strictly marginal correlations among the effective Majorana hoppings in the two chains.
Together with the understanding of the locally correlated XY model in the previous section, this leads us to propose the following picture for the critical XYZ chain along the line separating the x-AFM and y-AFM phases. For smallJ z , this critical line is actually controlled by the line of free Majorana fixed points with locally correlated hoppings characterized in Sec. V. The effect of the interactions J z in the original model with no correlations among the couplings (δ = 0) is to develop such correlations among the renormalized J x and J y couplings under RG while the J z couplings flow to zero. The ultimate degree of such correlations (i.e., the fully renormalized parameter δ eff ) then determines the long-distance power laws in the average spin correlation functions. We further conjecture that this persists for allJ z <J z crit = 1 below the transition to the z-AFM phase. While we do not have perturbative control close to this transition, any alternative would require yet another transition belowJ z crit which we did not observe and consider to be less natural. Note that in this scenario the transition to the z-AFM phase is controlled by a different non-free-fermion fixed point, and we do not have access to this S 3 -symmetric fixed point in the present study. We will further discuss the above conjecture, its corollaries and possible tests, as well as open questions in the concluding section.

VII. DISCUSSION
In this paper, motivated by the observations of Slagle et al. [4], we have performed a study of the low-energy properties of the random XYZ model using unbiased numerics. We focus on the line separating the x-AFM and y-AFM phases, which exhibits statistical symmetry between J x and J y couplings. At all points allowing comparison our results are in general agreement with the previous findings of Ref. [4] which used SBRG and presumed critical MBL physics at arbitrary energy density. Our results strongly suggest that-regardless of the behavior of highly excited states-there is quantum critical behavior in the ground state and the critical line is described by IRFPs with continuously varying critical exponents in the disorder-averaged correlation functions. Perhaps surprisingly, a Hartree-Fock mean-field theory treating the J z interaction terms as perturbations around the ran-dom XY (free-fermion) fixed point yielded results that are qualitatively rather consistent with the full interacting model at small to moderate J z couplings, including continuously varying power laws. This is in contrast to the clean case, where the mean field model is not qualitatively accurate due to divergences in the perturbation theory [66].
The locally correlated XY effective model, introduced with the idea of distilling the essential feature of the mean field theory, again exhibited continuously varying critical exponents, which we were able to establish numerically in larger sizes than for the XYZ chain. Because of the particular free-fermion form of this effective model, we were able to treat it in the SDRG using the random walk formulation in two dimensions. By making use of a connection between survival probability and the structure of decimation in the RG, we showed analytically that critical exponents for end-to-end and bulk C z spin correlations vary continuously as the coupling correlation parameter δ is tuned, and we also observed varying exponents in the bulk C ⊥ correlations by running the SDRG numerically. This result singles out and proves one of the scenarios of Fisher [3] that random anisotropy is strictly marginal along the critical line connecting the random XX and random XY fixed points; that is, there is a line of fixed points connecting the XX and XY IRFPs as sketched in Fig. 1.
Motivated by the successful understanding of the locally correlated XY model, we revisited the SDRG for the full interacting XYZ chain in the regime of small interactions and proposed a scenario where these interactions are irrelevant, but during the initial flows they generate effective correlations between the local J x and J y couplings (i.e., Majorana hopping amplitudes on the two chains). Such flows are sketched in Fig. 2. These local correlations in the free-fermion couplings then lead to non-universal power laws in the average spin correlations: this is our story for the continuously varying criticality in the XYZ spin chain.
We note that continuously varying critical exponents were previously observed in IRFPs associated with correlated disorder by Rieger and Iglói [67], however in a qualitatively different setting than ours. Specifically, disordered fixed points perturbed by the introduction of long-range correlations ∼ r −a to the disorder in the random transverse-field Ising chain exhibit critical indices varying continuously with a for a < 1. Their setting has only one Majorana chain and the correlated disorder is within the chain. Also, in their case the ψ exponent varies continuously, which reflects a different character of the corresponding "random walker" imprinted by the long-range correlations in the disorder.
Non-universal exponents at IRFPs were also observed in cases with very broad (singular) distributions of random couplings [68,69]. This again occurs already in a single chain and has varying exponent ψ, and the variation can be traced directly to the singularity in the probability distribution of the microscopic couplings, while the exponents are universal for non-singular probability distributions.
The XYZ chain studied here is different from the above examples with varying exponents in that there are no long-range correlations or singular distributions input into the microscopic disorder. In this way the continuously varying exponents are intrinsic to this system rather than imprinted extrinsically. What is important in the XYZ chain is that we have two simultaneously critical Majorana chains whose couplings become locally correlated. This insight may be useful when looking for other IRFPs with intrinsic continuously varying critical indices.
We conclude by returning to the discussion of the proposed scenario for the fully interacting XYZ chain. This scenario is based on the conjecture that the four-fermion and higher terms are irrelevant other than feeding into correlations between the Majorana hoppings. While this is plausibly justified for small interactions in Sec. VI, we have not fully proved it and the status for intermediate interactions is less certain. In this respect, it would be useful to carry out a systematic numerical SDRG study of the fully interacting problem (e.g., using the scheme of Monthus [22]) keeping track of all generated interactions as well as allowing decimations of the interaction terms when they happen to be the strongest. If our scenario is correct, we should see the interaction terms progressively decreasing relative to the Majorana hoppings. One should be able to perform such a study also directly in the spin variables using the SBRG approach of Slagle et al. [4] projected onto the ground state branch, e.g., as used in Ref. [70] in a different problem. Employing the insights gained here, it should be helpful to interpret various Pauli string terms generated under the SBRG as either Majorana hoppings or specific multi-fermion interactions. The SBRG can also be indispensable for studying the putative S 3 -symmetric fixed point describing the transition to the z-AFM phase, as a possible new IRFP that is not tractable with available analytical tools.
Thinking about a broader phase diagram, our work suggests that it could be fruitful to add another parameter "axis" and study the XYZ chain with locally correlated J x and J y couplings in the bare model (analogous to parameter δ in the correlated XY model), in addition to the interactions J z . Figure 2 shows this parameter space, and constitutes a mild abuse inasmuch as it serves as both a phase diagram and a picture of RG flows, the latter of which occur in space not captured by just the two parameters. In the space shown, the bare δ = 0 corresponds to the present XYZ chain, with the transition from the critical phase to the z-AFM phase at the S 3 symmetric point, marked XYZC in Fig. 2. On the other hand, δ = 1 corresponds to the XXZ chain studied in the original work by Fisher [3]. ForJ z below some threshold value, the XXZ spin chain is critical and controlled by the free-fermion XX point, while for largerJ z it undergoes a transition to the z-AFM phase. Fisher concluded that this transition is controlled by the so-called XXZC fixed point which is essentially random singlet-like, also marked in Fig. 2. An interesting question is the nature of the transition to the z-AFM phase driven by theJ z coupling as we vary the disorder correlation parameter from δ = 1 (XXZC fixed point) to the statistically isotropic XYZC fixed point. This line is marked with a question mark in Fig. 2, and one possibility is that it is also described by a line of fixed points, but we cannot at present exclude other scenarios. We leave these questions for future investigations, noting that the possibility of novel IRFPs is quite tantalizing and worth further exploration.