Random singlets and permutation symmetry in the disordered spin-2 Heisenberg chain: A tensor network renormalization group study

We use a tensor network renormalization group method to study random $S=2$ antiferromagnetic Heisenberg chains with alternating bond strength distributions. In the absence of randomness, bond alternation induces two quantum critical points between the $S=2$ Haldane phase, a partially dimerized phase and a fully dimerized phase, depending on the strength of dimerization. These three phases, called ($\sigma$,$4-\sigma$)=(2,2), (3,1) and (4,0) phases, are valence-bond solid (VBS) states characterized by $\sigma$ valence bonds forming across even links and $4-\sigma$ valence bonds on odd links. Here we study the effects of bond randomness on the ground states of the dimerized spin chain, calculating disorder-averaged twist order parameters and spin correlations. We classify the types of random VBS phases depending on the strength of bond randomness $R$ and dimerization $D$ using the twist order parameter, which has a negative/positive sign for a VBS phase with odd/even $\sigma$. Our results demonstrate the existence of a multicritical point in the intermediate disorder regime with finite dimerization, where (2,2), (3,1) and (4,0) phases meet. This multicritical point is at the junction of three phase boundaries in the $R$-$D$ plane: the (2,2)-(3,1) and (3,1)-(4,0) boundaries that extend to zero randomness, and the (2,2)-(4,0) phase boundary that connects another multicritical point in the undimerized limit. The undimerized multicritical point separates a gapless Haldane phase and an infinite-randomness critical line with the diverging dynamic critical exponent in the large $R$ limit at $D=0$. Furthermore, we identify the (3,1)-(4,0) phase boundary as an infinite-randomness critical line even at small $R$, and find the signature of infinite randomness at the (2,2)-(3,1) phase boundary only in the vicinity of the multicritical point.


I. INTRODUCTION
The ground state properties of antiferromagnetic Heisenberg spin chains have attracted a lot of attention for many decades, in particular after Haldane's conjecture [1,2] that half-integer and integer spin chains are distinct from each other.Half-integer spin chains with Heisenberg interactions have ground-state properties generically similar to the exactly solvable spin-1/2 chain [3], which has a gapless excitation spectrum [4] and power-law decaying spin correlations.Integer spin chains, on the other hand, have gapped ground states with exponentially decaying spin correlations [1,2] and a hidden topological order that is characterized by a nonlocal string correlation function [5].Qualitative differences are also found between half-integer and integer spin chains when the coupling constants alternate between two different values [6,7].The ground state is dimerized by infinitesimal alternation for half-integer spins, while the Haldane phase in an integer spin chain changes to a dimer state via a phase transition only with sufficiently strong bond alternation.
Another theme of high interest in condensed matter physics is the interplay between disorder, interactions and quantum fluctuations.The ground-state properties of low-dimensional quantum systems are often modified dramatically by introducing quenched disorder (i.e.timeindependent disorder).Remarkably, there are properties resulting from quenched disorder that are universal for a broad class of quantum spin chains, independently of whether the spin is integer or half-integer.The socalled random-singlet (RS) phase [8] is one of the disorder induced phases that possesses such universal properties.This phase describes the ground state of a spin-1/2 Heisenberg chain with any amount of disorder in couplings, and is also the ground-state phase of the spin-1 chain in the strong disorder limit where the excitation spectrum becomes gapless and the string topological order vanishes [9][10][11][12][13].
The RS phase was first found in the ground state of the random spin-1/2 chain by using the strong-disorder renormalization group (SDRG) method [8,[14][15][16].The SDRG method for the spin-1/2 chain consists of iteratively locking the strongest coupled spin pair into a singlet state, which decouples from the rest of the chain after effective couplings are generated among the remaining spins.This SDRG scheme ultimately flows toward an RS fixed point [8] that asymptotically represents the system's ground state, in which pairs of strongly entangled spins form singlets over all length scales, mostly short ranged but occasionally very long ranged.Those long-ranged singlets are rare; however, they dominate the arXiv:2309.04249v2[cond-mat.str-el]15 Dec 2023 average spin-spin correlations that decay asymptotically with distance L as an inverse-square form L −2 .By contrast, the correlations between typical pairs of widely separated spins are very weak and decay exponentially with the square root of their distance.Furthermore, the characteristic energy scale, ϵ, and length scale, L, of the singlets in the RS phase follow: with ψ = 1/2.This energy-length scaling is very different from the standard dynamic scaling, ϵ ∼ L −z , and implies that the dynamic exponent diverges z → ∞.The RS fixed point is one example of an infinite-randomness fixed point [17], which is characterized by extremely broad distributions of physical properties, even on a logarithmic scale, leading to the distinction between average and typical behaviors.
The SDRG method has been extended to higher spin chains [9][10][11][18][19][20].In particular, Damle and Huse applied an extended SDRG scheme to disordered spin-S Heisenberg chains with arbitrary S and obtained a class of infinite-randomness fixed points, called permutation symmetric fixed points [19].In the valence-bond (VB) picture for spin S > 1/2, each spin S is replaced with 2S virtual spin-1/2 variables, and valence bonds (singlets) are pairwise created between spin-1/2 variables that belong to different spin-S sites.For a spin-S chain, there are 2S + 1 distinct valence-bond solid (VBS) domains, denoted by (σ, 2S − σ) with σ ∈ {0, 1, • • • , 2S}; a VBS domain of type (σ, 2S − σ) (type σ) consists of σ spin-1/2 singlets over each even bond and 2S − σ spin-1/2 singlets over each odd bond (Fig 1).Among these generalized VBS states, the Haldane phase in an integer spin-S chain is associated with the symmetric (S, S) domain.A dimerized state (σ, 2S − σ) with σ ̸ = S can be realized in chains with alternating bond strength distributions.In the extended SDRG scheme, the degrees of freedom are effective spins of magnitude localized at the boundaries between distinct domains of type σ and σ ′ .These domain-wall spins interact with spins in neighboring domain walls through effective couplings that can be antiferromagnetic or ferromagnetic, depending on the types of environment domains.The renormalization of strong effective couplings between domain-wall spins under the RG leads to reconfigurations of associated domains.Using this domain-wall picture, Damle and Huse have predicted a series of infiniterandomness multicritical points, P n , that result from a competition between n domains.At a P n multicritical point, n domains appear with equal probability; hence the multicritical point is called permutation symmetric critical point.For a P n multicritical point, we have for the energy-length scaling exponent given in Eq. ( 1); this is a generalization of the RS fixed point with n = 2, where two domains occur with equal probability under the action of the RG.The Griffiths singularity (also known as the Griffiths-McCoy singularity for quantum systems) [21][22][23] is another interesting phenomenon arising from the interplay between quantum fluctuations and randomness.This phenomenon is characterized by singular low-energy behavior of various observable, such as the susceptibility and the specific heat, in an off-critical phase.Griffiths effects occur when there are rare but arbitrarily large spatial regions that are locally in a phase B due to disorder fluctuations, while the system is overall in a phase A. Quantum fluctuations enhance the effects of the rare regions, leading to a power-law behavior of the density of low-energy excitations such that [24] ρ(ϵ) ∼ ϵ −1+1/z (3) for a one-dimensional system, with a non-universal continuously varying exponent z that also describes the length scale and the energy scale through ϵ ∼ L −z .This power-law density of states is responsible for power-law singularities of certain observables at low temperatures.
Using the SDRG analysis [9,18,25], the generalized VBS states in dimerized spin-S chains, including dimer phases for any S and the Haldane phase in an integer-S chain, in the presence of sufficiently strong disorder have been identified as Griffiths phases, where the energy gaps are filled in yet the topological order (Haldane order and dimer order) persists.Considerable numerical efforts have been devoted to examine the theoretical predictions about the groundstate phases of random spin S > 1/2 chains, but so far have been mainly focused on the S = 1 chain [12,13,[26][27][28].Here we first use a tensor network strong-disorder renormalization group (tSDRG) method [29] to study the zero-temperature phases of the random spin-2 chain with alternating bond strength distributions, which to our knowledge have not been previously investigated numerically.
The paper is organized as follows.In Sec.II we define the model and summarize some known properties of the ground-state phases.In Sec.III we outline the SDRG and tSDRG methods.In Sec.IV we provide our numerical results for the ground-state phases depending on randomness and dimerization, focusing on a VBS order parameter and end-to-end correlations.We conclude in Sec.V with a summary and discussion.

II. THE MODEL
The model we consider is the spin-2 antiferromagnetic Heisenberg chain, described by the Hamiltonian: where S i is the spin-2 operator at site i, and J i are random coupling constants given by where the parameter D, with |D| ≤ 1, measures the strength of bond alternation (dimerization), and K i are random positive variables with the following power-law distribution: where R > 0, being the standard deviation of ln(K), parameterizes the strength of the randomness.This powerlaw distribution of bond randomness at D = 0 has been widely used in previous numerical studies for disordered systems [12,13,27,30,31].
According to the Damle-Huse domain-wall picture [19], random-singlet RS S phases out of spin S = 1/2, 1, 2 variables and multicritical points P n with possible n = 3, 4, 5 can occur when disorder comes into play.Since the undimerized spin-1/2 chain with any weak randomness is in the RS 1/2 phase, the critical domain walls at |D c,1 | and |D c,2 | with effective spin-1/2 are expected to evolve into an RS 1/2 state for R > 0. On the other hand, the RS 2 phase, arising from a competition between the (4, 0) domain and the (0, 4) domain, occurs only in the strong disorder limit at D = 0.In the RS 2 phase, spin-2 singlets connect sites on different sublattices over arbitrarily long distances in a random fashion, completely analogous to the RS 1/2 state.
For the spin-2 chain, there are three possible arrangements of random-singlet RS S phases and multicritical points P n in the R-D plane between weak and very strong disorder, as shown in Fig. 2, among which the phase diagram with the maximally symmetric multicritical point P 5 may occur only with some additional fine-tuned parameter in the model [19], or can be realized only with higher symmetry groups [35][36][37].All multicritical points P n belong to a set of infinite-randomness fixed points with the critical exponent ψ n given in Eq. ( 2), and the correlation-length exponent given by [9,10,19] In particular, the random-singlet RS S phases for all S is the special case with n = 2, where we have ψ = 1/2 and ν = 2.The off-critical regions, including the Haldane phase associated with the (2, 2) domain and the dimerized phases, have excitation gaps and exhibit topological order in the absence of randomness.The topological order in the Haldane phase can be detected by a nonzero limiting value of the generalized string order parameter [5,38] O z j,k = − S z j exp(iθ in the |j − k| → ∞ limit, where S z j is the z component of the spin-2 operator at site j, and an angle parameter with θ = π/2 is mostly suitable for the spin-2 case since O z j,k becomes a smooth function of the distance |j − k| [39], similar to θ = π for the spin-1 case [5].The topological order in the dimer phases is enforced by the Hamiltonian, and its sign reflects whether the Hamiltonian favors singlet pairs to be formed over even or odd bonds.Both the Haldane topological order and dimer order can survive in the presence of randomness even when the energy gaps close up .Such a gapless region with topological order is the Griffiths phase.
The main goal of our numerical study is to identify which diagram in Fig. 2 corresponds to the phase diagram of our spin-2 random chain.

A. SDRG
The Hamiltonian of a disordered quantum many-body system has an intrinsic separation in energy scales that Three possible phase diagrams of spin-2 chains in the R-D plane, predicted in Ref. [19].
allows us to find its ground state using the SDRG technique.The essential idea of the SDRG method, introduced in Ref. 14 and 15, is to find the largest term in the Hamiltonian successively and put the subsystem associated with this term into its ground state.The couplings between this subsystem and the rest of the system are treated by perturbation theory and effective couplings across the subsystem are generated.For example, for the random spin-1/2 Heisenberg antiferromagnetic chain, an effective coupling is generated between spins k − 1 and k + 2 with strength J = 1 2 when J k is the strongest coupling that locks spins k and k + 1 into a singlet at a certain step of RG.The new energy scale is then the strength of the strongest remaining coupling: Ω = max{ J}.By repeating this RG procedure, we gradually lower the energy scale and reduce the number of degrees of freedom in the system.The RG flow equation describing the evolution of the probability distribution under the RG process has been solved by D. S. Fisher in Ref. [8] for the spin-1/2 chain.The multiplicative relation in Eq. ( 9) suggests that it is more convenient to measure bond strengths on a logarithmic scale.In terms of logarithmic variables defined as Γ = − log(Ω) and ζ = log(Ω/ J), the fixed point distribution for the spin-1/2 chain corresponds to For higher spin chains, renormalized couplings in the conventional SDRG method may become stronger than the decimated couplings when the randomness is not sufficiently strong, which makes perturbation theory invalid.Thus, several extended SDRG methods, based on effective S = 1/2 models with both antiferromagnetic and ferromagnetic couplings, have been proposed for higher spin chains [9-11, 19, 40].Using a domain-wall model, Damle and Huse have extended Fisher's RG analysis to arbitrarily high spin.In the domain-wall model, one defines ρ σ to be the probability that a specific domain is of type σ and W σσ ′ to be the probability of a domain of type σ followed by one of type σ ′ , which are required to formulate the RG flow equations.The fixed point solution, which controls the multicritical point P n , is found to be given by the bond strength distribution: and indicating the domain permutation symmetry.From the fixed point solution, one can deduce the energy-length scaling relation − ln ϵ ∼ L ψn with ψ n given in Eq. ( 2).

B. Tensor networks and SDRG
Here we use a tree tensor-network generalization of the SDRG, referred to as tSDRG, to study the spin-2 random chain.This generalized SDRG scheme, proposed by Goldsborough and Römer [29], formulates the RG procedure as a tree tensor network and refines the perturbative approximation of SDRG by including higher energy states at each RG step, along the lines of a previous SDRG extension [41].The tSDRG method has been applied in studies of the quantum Ising chain [42] and spin-1 chains [13], where accurate results that are compatible with the results obtained by non-approximate quantum Monte Carlo calculations [43] and the density matrix renormalization group [12,28] are achieved.
The starting point of tSDRG is to express a onedimensional Hamiltonian of L sites as a sum of twosite Hamiltonian in terms of matrix product operators (MPOs) [44]: where an MPO W [i] at site i is a matrix of operators.Specifically for the model considered in this work, the two-site Hamiltonian reads where S + and S − are the ladder operators.For a chain The three-leg tensor V (triangles), built from the χ ′ lowest energy eigenvectors, truncates a two-site tensor into a renormalized site.The blue-shaded squares are local tensors of the MPO.The vertical legs denote physical indices; the horizontal legs denote virtual indices; (b) Isometric property with open boundary conditions (OBC), we have the following explicit form of the W tensors: for a site in the bulk, i.e. i ̸ = 1, L, and for two edge sites i = 1 and L. For a chain with periodic boundary conditions (PBC), the MPO tensors for i = 1 • • • L are all bulk tensors as given in Eq. ( 15), where the coupling J L links between two end sites L and 1.An important observation is that the two-site Hamiltonian H i,i+1 is encoded in the local matrix product It is easy to verify that 5,1 .The essential information required for the tSDRG procedure is the list of MPO W [i] and the list of two-site Hamiltonian H i,i+1 .
Similar to the conventional SDRG, in each RG iteration one selects a pair of adjacent sites to be renormalized, depending on the local energy spectrum.Here the selection is based on the largest energy gap, rather than the strongest coupling (corresponding to the lowest gap) [41,45,46].For each two-site Hamiltonian H i,i+1 we identify the energy gap ∆ i,i+1 , which is measured as the difference between the highest energy of the χ ′ -lowest eigenstates that would be kept and the energy of the (χ ′ + 1)-th eigenstate.We set a bond dimension χ as the upper bound of the number χ ′ to control the accuracy of the calculation.The actual number χ ′ is adjusted to keep full SU(2) multiplets.By increasing the bond dimension χ, we can obtain more accurate results but at the cost of computational resources and time.
After obtaining all local gaps described above we identify the two sites with largest gap.These two sites will be merged into a renormalized site as follows.We first express the k-th eigenstate |Ψ k ⟩ of H i,i+1 in terms of the local two-site product basis where |s i ⟩|s i+1 ⟩ is the local product basis.We then construct a projector that projects the two-site space to the renormalized site space with dimension χ ′ , In terms of matrix notation one can write V as and its tensor network diagram is given in Fig. 3.Note that the tensor V has the isometric property V † V = 1.To obtain the renormalized MPO W [i,i+1] associated with the renormalized site, one uses V to renormalize each element of the product as follows The two-site Hamiltonians that contain the renormalized site can be decoded as follows: This completes one iteration of tSDRG.We have now an updated list of MPOs and two-site Hamiltonians.Conceptually, one can re-label the sites so that the renormalized MPO and two-site Hamiltonian are labeled as W [i]  and H i,i+1 respectively.In practice the relabeling is not necessary.
To obtain the ground state of the system, we should repeat the tSDRG iteration until the system contains only two renormalized sites.At this stage we diagonalize the top Hamiltonian H top to obtain its ground state |Ψ top 1 ⟩.From |Ψ top 1 ⟩ and the projector V at each iteration, one A tree tensor network for the ground-state expectation value of the end-to-end correlator which acts only on two edge sites.Since V † V = 1, only those isometric tensors (orange triangles) that affect the two edge sites are considered.
can generate an inhomogeneous binary tree tensor network as sketched in Fig. 4(a).The expectation value of a product of local operators can be obtained by contracting the operators with the tree and its conjugate as sketched in Fig. 4(b).The contraction can be evaluated efficiently thanks to the property that V † V = 1.

IV. NUMERICAL RESULTS
In this section we use the tSDRG method to explore ground-state phases of the random spin-2 chain with alternating bond strength distributions as defined in Eqs. ( 5) and ( 6), focusing on two observables: the VBS order parameter based on a unitary operator appearing in the Lieb-Schultz-Mattis theorem [47] and the end-toend correlation function.We first show the R-D phase diagram, extracted from our numerical data, in Fig. 5 before presenting detailed results.We restrict ourselves to the D ≥ 0 region because the results for D < 0 can be obtained simply via the parity symmetry.The phase diagram that we have identified largely agrees with the type in Fig. 2(c).

A. VBS order parameter
We consider a unitary operator, called the twist operator, defined for a chain with L spins as operator was first introduced in the Lieb-Schultz-Mattis theorem [47][48][49], which states that for the ground state, Ψ GS , of a half-integer spin chain, one has in the limit L → ∞, indicating a gapless excitation spectrum.Furthermore, the ground-state expectation value of this operator has been found to be capable of detecting and characterizing VBS order [34].For a VBS state of type-σ, the asymptotic form of the expectation value is given by [34] that is, it is positive (negative) in the L → ∞ limit if σ is even (odd).Using the properties of z L , groundstate phase diagrams of dimerized spin-S chains with S = 1/2, 1, 3/2, and 2 in the absence of randomness were determined in Ref. [34].Remarkably, this order parameter is applicable also for strongly disordered systems, as demonstrated in a quantum Monte Carlo study for the random spin-1 Heisenberg chain [27].
Here we calculate z L using the tSDRG method for the random spin-2 chain with periodic boundary conditions.We explore the behavior of the disorder-averaged order parameter, z L , for a wide range of randomness and dimerization, parametrized by R and D. We have considered system sizes up to L = 512 and more than 5,000 random coupling samples to obtain the disorder average.Figure 6 shows the disorder-averaged order parameter for different system sizes at R = 1.5, 1.0, and 0.5 as a function of D. Here, in the cases with R = 0.5 and R = 1.0, one can clearly observe that the order parameter changes its sign at certain values of D; there are two such zerocrossings for each L that can be identified as two phase transition points between different random VBS states.The sign of z L indicates that the domains from left (small D) to right (large D) are successively in (2, 2), (3, 1) and (4, 0)-dominant phases.By varying R, the two domain walls, located at the zero-crossings of z L , form critical lines that connect the two critical points in the clean limit, i.e.D c,1 ≈ 0.18 and D c,2 ≈ 0.55 at R = 0.The distance between these two critical lines decreases as R increases.For the case with R = 1.5, one can observe the minimum value of z L for L = 512 (the largest system size that we consider here) is about zero; we can than identify the associated (R p , D p ) with this minimum as the junction of the two critical lines for the finite chain.In Fig. 5, the point at (R p , D p ) is denoted by P > 3 .The results described above for the disorder-averaged twist order parameter suggest that there is a multicritical point, at which three phases (2, 2), (3, 1) and (4, 0) meet, in the region of D ̸ = 0.This implies that the diagram in Fig. 2(c) corresponds to the R-D phase diagram of our model.To determine the location of this multicritical point in the limit of L → ∞, we find the (R p (L), D p (L)) point at which the disorder-averaged twist order parameter reaches its minimum at z = 0, and then estimate the critical values for R p and D p from an extrapolation to L → ∞ using and with ν = 2.3 for a multicritical point P 3 (see Eq. ( 7)).By doing so, we obtain R p (∞) ≈ 1.15 and D p (∞) ≈ 0.04, as shown in Fig. 7.Here we round to two decimal places to point up D p (∞) ̸ = 0, which is the chief characteristic of the phase diagram in Fig. 2(c).We note that, for a control parameter λ, the deviation of a finite-size pseudocritical point λ c (L) from the true critical point λ c (∞) in the limit of L → ∞ is often parameterized as ln(λ c (L)) − ln(λ c (∞)) for an infinite randomness fixed point [8,[50][51][52].For the bond strength distributions given in Eq. ( 5) and Eq. ( 6), the distance from the critical point defined in Eq. ( 28) is thus consistent with such logarithmic parametrization; the distance given in Eq. ( 29) is also suitable because the values of D we consider here are so small that the approximation ln(1 + D) ≈ D is valid.Now we turn to the undimerized region with D = 0. Figure 8 shows the dependence of disorder-averaged twist order parameters on R.Here the twist order parameter for each system size L is positive before converging to z L = 0 in the large R region, consistent with the scenario in which the system's ground state changes from a (2, 2) phase to an RS phase when the randomness exceeds a critical value; this critical point at D = 0 is also a multicritical point P 3 at which three phases (2, 2), (4, 0) and (0, 4) meet (see Fig, 2(c)).Since the order parameter z L does not change its sign for the (2, 2)-RS phase transition nor for a nearby (2, 2)-(4, 0) transition, it is difficult to determine the transition point accurately by z L .Nevertheless, the results shown in Fig. 8 (and also in Fig. 6 (a)) suggest that the critical R value at D = 0 for the largest size L = 512 is about R ≈ 1.5 (or slightly higher).Thus the multicritical point at D = 0 (denoted as P 0 3 in Fig. 5) and the one at D > 0 ( P > 3 in Fig. 5) have critical R values that are not far apart.As seen in Fig. 6(c) and Fig. 8, there are intersection points developing at some nonzero values of z L .In Ref. [27], such intersection points of z L , instead of zerocrossings, were used to identify the multicritical point and the RS critical line for the random spin-1 chain.Here, by exploring a wide range of parameters (R, D) and accessing larger system sizes, we have found that the z L curve crossings appear only in the region of small R, or the crossing points tend toward z L = 0. Thus, the zero order parameter, z L = 0, turns out to be a more reasonable indicator for a transition point even in disordered systems.

B. End-to-end correlations
In this subsection we investigate the end-to-end correlations in an open chain, which considers correlations between two end spins.This quantity is defined as for the ground state |Ψ GS ⟩ of a chain with L spins and open boundary conditions.
Here we first summarize some previous results for the end-to-end correlations.The end-to-end correlation of an open chain is closely related to the energy gap [53,54].In the RS phase the end-to-end correlations are very broadly distributed, with the typical behavior [8], where ψ = 1/2.On the other hand, the average endto-end correlation function decays as a power of L at 10. Average end-to-end correlation versus system size for the undimerized case (D = 0) with various R in a log-log plot.All lines are fits to the form aL −η .From the slopes of the fitting lines, corresponding to the exponents η1, we estimate R ≈ 1.2 for the multicritical P3 point.For R ≥ 1.5, the slope approaches η1 = 1, indicating that the system enters an RS phase; here it is the RS2 phase.
criticality [53], with η 1 = 1 in the RS phase.The behavior of the average end-to-end correlation in Eq. ( 32) was first derived for the infinite-randomness fixed point of the random transversefield Ising spin chain [53] and is also valid for the RS phase, as verified numerically in Ref. [12,13].At the multicritical point P 3 , typical end-to-end correlations go like Eq. ( 31) but with ψ = 1/3, according to SDRG analytical results.There have been so far no theoretical conjectures about the exponent in Eq. ( 32) for the average end-to-end correlations at P 3 ; nevertheless, previous numerical results [12,13] for the spin-1 chain estimated η 1 ≈ 0.69 − 0.7 for the P 3 multicritical point, which is expected also for the P 3 point in the spin-2 chain.Away from an infinite-randomness critical point, the distribution of end-to-end correlations has a power-law tail that behaves as [54] with a finite dynamic exponent z; in a Griffiths phase, the singular low-energy behavior of various observables is characterized by a large and continuously variable dynamic exponent z > 1.
We first examine the behavior of end-to-end correlations at the P 3 points.For the dimerized cases, we consider open chains with odd numbers of spins to balance the numbers of strong and weak couplings.Figure 9(a) shows the average of the correlations at R = 1.15 and D = 0.04 (the point P > 3 in Fig. 5), where the location of the multicritical point P 3 is according to finitesize scaling of the twist order parameters discussed in Sec.IV A.Here we estimate η 1 ≈ 0.69 from the results for the average correlation as a function of the chain length L, in good agreement with previous numerical results [12,13].Also the distributions of the logarithmic correlations (shown in Fig. 9(b)), which become broader with increasing size, can collapse onto each other by using the scaled variable with ψ = 1/3 [Fig.9(c)], consistent with the theoretical prediction [19].
For the undimerized case D = 0, we show the average correlations plotted against L for various values of R in Fig. 10; all data here decay as a power law: C 1 (L) ∼ L −η1 .From the slopes of the fitting lines, corresponding to the exponents η 1 , we estimate R ≈ 1.2 for the multicritical P 3 point (denoted by P 0 3 in Fig. 5) by comparing the value of η 1 with previous numerical results for the P 3 point in the random spin-1 chain [12,13].For stronger randomness, such as R = 1.5 and R = 1.6, the slopes approach η 1 = 1, which is the theoretical value for an RS phase; in this case, it is the RS 2 phase.Now we turn to the correlations on the critical lines at D > 0 that are the (2,2)-(3,1) and (3,1)-(4,0) boundaries, determined using the zero-crossings of the z L order parameter for the largest system size L = 512.Specifically, we focus on the four points denoted by a, b, c and d in Fig. 5. Figure 11 and Figure 12 show the averages and the distributions of end-to-end correlations at four such finite-size critical points.For the critical points, (R = 1.0,D = 0.34) and (R = 0.5, D = 0.48), at the boundary between the (3, 1) and (4, 0) phases, our numerical results are found to be in good agreement with the analytic predictions for an RS phase: the average end-to-end correlations decay as C 1 (L) ∼ 1/L, as shown in Fig. 11(c) and (d), and the distributions of the correlations, shown in Fig. 12(c) and (d), are extremely broad and become broader with increasing L. Furthermore, the distributions for different system sizes collapse well onto each other by using the scaled variable defined in Eq. (34) with the theoretical value ψ = 1/2 for the RS phase, as shown in Fig. 13(c)and (d).
On the other hand, the results for the (2,2)-(3,1) phase boundary are not fully compatible with the scenario of an RS phase.The mean end-to-end correlations C 1 (L) at the (2,2)-(3,1) phase boundary, as shown in Fig. 11(a) and (b) for the two critical points (R = 1.0,D = 0.12) and (R = 0.5, D = 0.14), appear to decay much slower than 1/L, contrary to the behavior in an RS phase.In particular, the distribution of the correlations in the region of small R, such as at (R = 0.5, D = 0.14) (see Fig. 12(b)), does not broaden with increasing size, implying that the phase is not associated with an infiniterandomness critical point.Assuming the scaling form given in Eq. ( 33), we estimate the dynamic exponent z ≈ 0.41 from the slope of the power-law tail of the distribution, which gives the scaling plot shown in Fig. 13(b).With stronger disorder at (R = 1.0,D = 0.12), the distribution becomes broader with increasing size (see Fig. 12(a)), showing the signature of infinite randomness.We have used the scaled variable in Eq. ( 34) and ψ = 1/3 to achieve good data collapse, as shown in Fig. 13(a).
Strong finite size effects and the close distance from the multicritical point at R = 1.15 may lead to the discrepancy between the exponent ψ found here and the theoretical predicted value ψ = 1/2 for an RS 1/2 phase.The small dynamic exponent z < 1 estimated for (R = 0.5, D = 0.14) does not lead to divergence of the local susceptibility [12,24], indicating there is no pronounced disorder-induced singular behavior.Another way, suggested in Ref. [12], to identify a nonsingular region with z < 1 or a singular region with z > 1 is via the so-called inverse average of the end-to-end correlation, defined as where C −1 1 is the average of the inverses.The inverse average of C 1 is finite, C inv 1 > 0 if z < 1 , while it is zero if z > 1.As shown in Fig. 14, the inverse average C inv 1 at R = 0.5 is finite for a wide range of D before it converges to zero at D ≈ 0.4 for larger L, while the inverse average vanishes at R = 1 independent of D. The results for the inverse average C inv 1 also imply that there is no random-singlet phase at R = 0.5 in the region with small dimerization D < 0.4.This large nonsingular region also gives rise to the weak signature of infinite randomness in the larger R region at the (2,2)-(3,1) phase boundary.

V. SUMMARY AND DISCUSSION
Using a tensor-network SDRG method, we have explored the ground-state phases of the random spin-2 antiferromagnetic chain with alternating bond strength distributions.We have calculated the twist order parameter, defined as the ground-state expectation value of the unitary operator in Eq. ( 25), to classify the types of random VBS phases depending on the strength of bond randomness R and the dimerization D. For a disorder-free VBS phase (σ, 4 − σ) in a clean system, the twist order parameter is positive if σ is even and negative if σ is odd [34].In a random VBS domain, there is nonzero residual VBS order (dimerization) that can be detected by the disorder average of this order parameter, as we have demonstrated in this paper.Therefore, the zerocrossing of the disorder-average twist order parameter can serve to determine the phase transition point between different random VBS states, in the same way as for clean systems [34].Our results largely agree with the phase diagram sketched in Fig. 2(c).There is a multicritical point in the intermediate disorder regime with finite dimerization, where (2,2), (3,1) and (4,0) three phases meet.The (2,2)-(3,1) phase boundary and the (3,1)-(4,0) boundary extend to R = 0 and are predicted to be in the RS 1/2 phase for any R > 0 [19].However, from the results for end-to-end correlations we see no signs of an RS phase at the (2,2)-(3,1) boundary with small R and have instead found a large nonsingular region, characterized by z < 1, in the small R regime.Such a nonsingular region with z < 1 in the weak disorder limit has previously also been identified in numerical studies of the random undimerized (D = 0) spin-1 chain [12,13] and spin-3/2 chain [55]; these studies used the same power-law distribution of bond randomness, which is also identical to the distribution we consider here for the undimerized case.The nonsingular behavior in the weak disorder and weak dimerization limit and the absence of the RS 1/2 phase there can have more than one source.Compared to the Haldane gap value of ≈ 0.41J in the spin-1 case with nearest-neighbor coupling J, the finite gap in the clean spin-2 chain is much smaller, just about ≈ 0.09J.Even so, this small energy gap appears to have considerable impact on the ground-state properties of the (2,2) phase at weak disorder.Furthermore, we did not set a constraint to the strength difference between odd and even bonds while using the bond strength distribution (given in Eq.( 5)) in our calculations; this may produce unsharp dimerization for small systems, especially in the small D limit, which in turn leads to a weak critical signature or the absence of the RS 1/2 phase at the (2,2)-(3,1) boundary.In this regard we present in Fig. 15 the average energy gap of the chain with PBC at R = 0.5, obtained from the lowest-lying excitation of the top Hamiltonian H top , versus the dimerization D.Here we observe a clear (local) minimum in the average energy gap at D = 0.48 for all system sizes, with the (3, 1)-(4, 0) transition point indicated by the right dashed line -estimated from the zero crossing of the twist order parameter.On the other hand, in the small D region (corresponding to the nonsingular region), the curve for the largest size L = 512 develops a less clear minimum around the estimated (2,2)-(3,1) transition point (indicated by the left dashed line), while the curves for smaller sizes are flat in this region, showing strong finite-size effects.
The twist order parameter is a useful quantity, also for disordered systems, to determine the phase transition point between different random VBS states based on changes of the sign according to the valence-bond config-uration.However, the phase diagram of the spin-2 chain considered here has a (2,2)-(4,0) phase boundary where the twist order parameter does not change the sign, which makes it difficult to detect this phase transition.It would be desirable to find an order parameter that can accurately determine the (2,2)-(4,0) phase boundary which connects two P 3 multicritical points.
Recently, there has been an increasing interest in properties of higher-spin materials both from the theoretical and experimental perspective, especially in the context of Kitaev models and quantum spin liquids [56][57][58].Like the spin-2 chain that we consider here, the quasi-onedimensional versions of these higher-spin materials exhibit rich phase diagrams [59,60].Concerning disorder effects on the low-temperature phases, the tSDRG method and its variants [61][62][63][64][65][66] are certainly the most promising numerical tools for studying large-scale systems with accuracy.

FIG. 4 .
FIG. 4. (a)The tSDRG algorithm as a tree tensor network for a chain of 8 sites with periodic boundary conditions.The squares indicate the W -tensors in the MPO representation of the system's Hamiltonian, the triangles are the V -tensors used to truncate the two-site operators, and the circles represent the top Hamiltonian encoded in the final two-site tensor.The RG iteration proceeds upwards in the vertical dimension.The part below the W-tensors is the conjugate of the upper part.(b) A tree tensor network for the ground-state expectation value of the end-to-end correlator which acts only on two edge sites.Since V † V = 1, only those isometric tensors (orange triangles) that affect the two edge sites are considered.

FIG. 5 .
FIG.5.The phase diagram in the R-D plane with D ≥ 0, extracted from our numerical results (unfilled circles) for L = 512 and previous numerical data[34] at R = 0.The points labeled by P 0 3 and P > 3 are three-phase multicritical points for L = 512 at D = 0 and D > 0, respectively.The points labeled by P 0 3 and P > 3 are the multicritical points estimated for L → ∞.The four points denoted by a, b, c and d are finitesize critical points whose properties are further discussed in Sec.IV A and Sec.IV B.

FIG. 6 .
FIG. 6. Disorder-averaged twist order parameters zL for different system sizes at R = 1.5, R = 1.0 and R = 0.5, plotted versus the dimerization D. The red labels P > 3 , a, b, c and d indicate the points where zL = 0 for L = 512 (green data); these points are also shown in Fig. 5.

FIG. 7 .FIG. 8 .
FIG.7.Finite-size scaling of the location (Rp, Dp) at which the twist order parameter reaches its minimum at zL = 0.The multicritical point in the thermodynamic limit is estimated as Rp ≈ 1.15 and Dp ≈ 0.04 from an extrapolation to L → ∞, using Eqs.(28) and (29) with ν = 2.3.

48 FIG. 12 .
FIG. 12.The distributions of the end-to-end correlations for different sizes at the same critical points as investigated in Fig. 11.The distributions for (R = 1.0,D = 0.12) [Subfigure (a)] located at the (2,2)-(3,1) boundary with sufficiently strong disorder and the distributions at the (3,1)-(4,0) boundary [(c) and (d)] are broad and become broader with increasing L. The correlations for (R = 0.5, D = 0.14) [(b)] at the (2,2)-(3,1) boundary with smaller R, on the other hand, are not very broadly distributed, and the data curves for different sizes are very similar, implying the critical phase is not of infinite-randomness type.

2 FIG. 13 . 1 RFIG. 14 .
FIG. 13.Scaling plots of the distributions in Fig. 12.The correlations in (c) and (d), at the the (3,1)-(4,0) phase boundary, are rescaled as ln C1/L ψ with ψ = 1/2, and the data in (a) for a critical point at the (2,2)-(3,1) phase boundary are plotted in terms of the same rescaled variable but with ψ = 1/3 to achieve good data collapse.The correlations in (b), for a critical point at the (2,2)-(3,1) phase boundary with smaller R and D, are rescaled by assuming P (C1) ∼ C −1+1/z 1 with z = 0.41, which is estimated from the slopes of the tails at small values of C1.

FIG. 15 .
FIG. 15.The average gap at R = 0.5 as a function of the strength of dimerization D. The gap for each sample is obtained from the lowest-lying energy gap of the top Hamiltonian Htop in the tree tensor network.The two dashed lines indicate the (2, 2)-(3, 1) boundary (left line) and the (3, 1)-(4, 0) boundary (right line), based on the zero crossings of zL for L = 512.