Dynamic structure factor of disordered quantum spin ladders
DDynamic structure factor of disordered quantum spin ladders
Max H¨ormann, ∗ Paul Wunderlich, † and K. P. Schmidt ‡ Institute for Theoretical Physics, FAU Erlangen-N¨urnberg, Germany (Dated: June 6, 2018)We investigate the impact of quenched disorder on the dynamical correlation functions of two-leg quantum spin ladders. Perturbative continuous unitary transformations with the help of whitegraphs and bond-operator mean-field theory are used to calculate the one- and two-triplon con-tribution of the zero-temperature dynamical structure factor. Disorder results in huge effects onquasi-particles as well as composite bound states due to localization. This leads to intriguing quan-tum structures in dynamical correlation functions well observable in spectroscopic experiments.
Disorder is an inevitable ingredient of any condensedmatter. On the one hand disorder can change or evendestroy the physical behaviour of the associated cleansystems [1–3] or, on the other hand, it can induce fun-damentally new physics. This is especially true for cor-related quantum materials where the interplay of disor-der and quantum fluctuations can result in technologi-cal challenges or exotic phases of quantum matter likemany-body localization [4–7]. One important aspect inthe collective behaviour of correlated quantum matter isthe formation of quasi-particles and their role in quan-tum critical behaviour. While many studies have inves-tigated the static and thermodynamic properties of suchsystems in the presence of disorder [8–16], the fate ofquasi-particles under disorder is only rarely studied. Ex-perimentally, however, increasingly improving resolutionin spectroscopy like inelastic neutron or light scatteringas well as intentional doping to control disorder in quan-tum materials [17–23], demands theoretical predictionsfor dynamical correlation functions of correlated quan-tum matter in the presence of disorder.An outstanding arena, both experimentally and the-oretically, to investigate the influence of quenched dis-order on quasi-particle excitations are low-dimensionalquantum magnets, which host a variety of interestingquantum phases and elementary particles in the cleancase. Theoretically, the calculation of dynamical corre-lation functions for disordered quantum magnets is chal-lenging [24–26]. Only recently [25], the effect of disorderon single quasi-particles has been studied within bond-operator mean-field (MF) theory for the bilayer squarelattice Heisenberg model, but the fate of quasi-particlesunder disorder beyond MF theory is largely unexplored.An particularly promising system to advance are anti-ferromagnetic two-leg quantum spin ladders (QSLs) [28],which have been investigated successfully over the lastdecades in a number of experimental realizations [29–35] and by various theoretical tools [36–41]. Clean QSLshave non-ordered ground states and gapped triplon ex-citations [37, 42]. Inelastic neutron and light scatteringallow to access one-triplon dispersions but also the forma-tion of two-triplon bound states and continua reflectingthe presence of large quantum fluctuations. Furthermore,
FIG. 1. (a) QSL with bimodal disorder. The two differentrung ( J ⊥ (1 , ) and leg ( J (cid:107) (1 , ) couplings are shown as red andblue bonds. The unit cell ν is indicated by the rectangularbox. (b) Spectrum ω of the clean QSL for J (cid:107) /J ⊥ = 1 / k including one-triplon dispersion(black line), two-triplon continuum (grey) as well as spin-onetwo-triplon bound state (red line) calculated by pCUTs. it is possible experimentally to intentially dope QSL com-pounds so that quenched disorder is induced into the sys-tem [43]. It is therefore important to understand the in-fluence of disorder on the spectral signatures of triplonquasi-particles in QSLs. This is exactly the punchline ofthis letter. We calculate the one- and two-triplon contri-bution of the dynamic structure factor (DSF) at T = 0of two-leg QSLs in the presence of quenched disorder. Itis demonstrated that disorder has huge effects on collec-tive excitations yielding intriguing quantum structuresdirectly relevant for spectroscopic experiments. Disordered QSL —
The Hamiltonian of the disorderedQSL for a fixed disorder configuration { J } is given by H ( { J } ) = (cid:88) ν (cid:32) J ⊥ ν (cid:126)S ν, · (cid:126)S ν, + (cid:88) n =1 J (cid:107) ν,n (cid:126)S ν,n · (cid:126)S ν +1 ,n (cid:33) , (1)where the sum runs over all rungs and n = 1 , a r X i v : . [ c ond - m a t . s t r- e l ] J un the legs of the QSL (see Fig. 1a). The disorder config-uration { J } given by the antiferromagnetic J ⊥ ν and J (cid:107) ν,n depends on the type of quenched disorder. Here we focuson bimodal disorder, i.e. the rung and leg exchanges cantake either the value J κ with probability p or J κ withprobability 1 − p for κ ∈ {⊥ , (cid:107)} . However, our techni-cal treatment is more general and allows to consider anyform of quenched bond disorder.Clean QSLs with J ⊥ ν ≡ J ⊥ and J (cid:107) ν,n ≡ J (cid:107) havean unordered ground state and gapped triplon exci-tations for all J (cid:107) /J ⊥ , which are adiabatically con-nected to the isolated rung limit J (cid:107) = 0. In thislimit the ground state is a product state of singlets | s (cid:105) =( | ↑↓(cid:105) − | ↓↑(cid:105) ) / √ | t +1 (cid:105) = | ↑↑(cid:105) , | t +0 (cid:105) =( | ↑↓(cid:105) + | ↓↑(cid:105) ) / √
2, and | t − (cid:105) = | ↓↓(cid:105) .The low-energy spectrum of the clean QSL obtained bypCUTs is shown for J (cid:107) /J ⊥ = 1 / k = π . Two-triplon energies are ei-ther part of the two-triplon continuum or correspond totwo-triplon (anti)-bound states, whose energies dependon the total spin. We focus on the spin-one sector, whichis relevant for the DSF, where a two-triplon bound stateexists for a large range of momenta in the clean QSL dueto the attractive two-triplon interaction [38, 39].We then reformulate (1) in terms of triplet creationand annihilation operators t ( † ) ν,α with t † ν,α | s (cid:105) ≡ | t α (cid:105) and α ∈ {± , } on rung ν . Setting the average rung ex-change ¯ J ⊥ ≡ ( J ⊥ + J ⊥ ) / ≡ J ⊥± from it, allows to express (1) as H ( { J } ) = E + Q + (cid:88) n = − ˆ T n ( { J } ) , (2)where E = − N r / N r the number of rungs, thecounting operator Q = (cid:80) ν,α ˆ n ν,α with ˆ n ν,α = t † ν,α t ν,α ,and the ˆ T n with [ ˆ T n , Q ] = n ˆ T n change the triplet num-ber by n . The ˆ T n depend explicitly on ∆ ¯ J ⊥± as wellas J (cid:107) , and therefore on the disorder configuration { J } .Here ˆ T ± correspond to pair creation and annihilationprocesses, ˆ T contains triplet hopping as well as quar-tic triplet-triplet interactions, and ˆ T ± represent decayprocesses of one triplet into two or vice versa. Note thatˆ T ± = 0 holds for the clean case where the QSL possessesan exact reflection symmetry about the centerline givingrise to a conserved parity quantum number ± S ± ( k, ω ) ≡ lim N dc →∞ N dc S ± ( k, ω, { J } ) (3)with momentum k , frequency ω , number of disorder con- figurations N dc , and S ± ( k, ω, { J } ) ≡ − π Im (cid:104) |O †± ω − H + i0 + O ± | (cid:105) (4)where O ± ( k ) ≡ (cid:80) ν e i kν ( S zν, ± S zν, ) / (2 √N r ). The index ± reduces to the parity quantum number for the cleanQSL. pCUT — We perform pCUTs as our main approach aswell as bond-operator MF theory (for the latter see [27]).The major target of a pCUT is to unitarily transform (2),order by order in J (cid:107) ν,n and ∆ ¯ J ⊥± , to an effective Hamilto-nian H eff which conserves the number of triplons so that[ H eff , Q ] = 0 holds. As a consequence, the complicatedquantum many-body system is mapped to an effectivefew-body problem which is easier to solve. A pCUT ap-plication has a model-independent step, which expresses H eff in a sum of operator product sequences of the ˆ T n operators with exactly known coefficients. The most ef-ficient way of performing the second model-dependentstep, which amounts to normal order H eff , is a full-graphdecomposition using the linked-cluster theorem. Here theonly graphs are ladders segments so that the calculationfor the clean QSL is very simple. In contrast, in thepresence of bimodal disorder, there are 2 N r +1 differentgraphs for a ladder segment of N r rungs and a linked-cluster expansion becomes unefficient since N dc is large.At this point white graphs [46] are essential, since it al-lows to specify { J } only after the calculations on thegraphs. In practice, one determines the most generallinked contribution of a graph by allowing for a differentexchange coupling on every nearest-neighbor link of thegraph. The resulting multi-variable series can then beembedded on any specific { J } .We calculated H eff in the one- and two-triplon sectorup to order 8 and we determined the corresponding effec-tive observables O eff ± up to order 7. The convergence ofthe bare series is similar to the clean QSL where it givesquantitative results up to J (cid:107) /J ⊥ (cid:46) . N r = 100 and the DSF for a fixed dis-order configuration is obtained using a finite broadeningΓ = 0 .
01. Averaging over N dc = 1000 disorder config-urations then gives the averaged DSF (3). All relevantaspects discussed below are well converged despite theuse of finite N r and N dc [27], which is reasonable due tothe finite localization and correlation length of the QSL. One-triplon contribution —
For the clean QSL, thissector can be expressed as S − ( k, ω ) = a ( k ) δ ( ω − ω k )with the one-triplon dispersion ω k and the one-triplonspectral weights a ( k ) while S + ( k, ω ) = 0 due to the par-ity symmetry. The dispersion ω k is cosine shaped withgap at k = π and the spectral weight a ( k ) is monotonicin k with maximum at k = π and minimum at k = 0.Representative pCUT results for the DSF in the pres-ence of maximal pure rung or leg disorder are shown in FIG. 2. One-triplon contribution to the anti-symmetric DSF S − ( k, ω ) for pure rung (left) or leg (middle) disorder and to thesymmetric DSF S + ( k, ω ) for pure leg disorder (right). Left : bimodal rung disorder with p = 0 . J ⊥ = 1 . J ⊥ = 0 .
6. The constant leg exchange is J (cid:107) = 0 .
4. White lines represent one-triplon dispersions for a clean QSL with J (cid:107) = 0 . J ⊥ = 1 . J ⊥ = 0 . (cid:15) ( L )1 of open ladder segments of length L with J ⊥ = 0 . Middle/Right : bimodal leg disorder with p = 0 . J (cid:107) = 0 . J (cid:107) = 0 .
4. The black dashed line represents the one-triplon dispersion with mean hopping amplitudes.
Fig. 2 (see [27] for qualitatively the same MF results).These two cases behave differently with respect to reflec-tion about the centerline. For rung disorder the asso-ciated parity is still a conserved quantity and therefore S + ( k, ω ) = 0 holds for the one-triplon sector. In con-trast, since the leg disorder breaks the reflection sym-metry, a finite one-triplon contribution to S + ( k, ω ) ex-ists. However, this contribution is orders of magnitudesmaller and the dominant contribution to the DSF is still S − ( k, ω ). We find that the DSF of the disordered QSLsis fundamentally different to the one for the clean coun-terpart. These differences do not originate from largechanges of the spectral weights. Indeed, deviations to aclean QSL with constant mean exchange couplings aresmall [27]. Certainly, the most important effect of dis-order on the DSF is the disorder-induced localization ofeigenfunctions and their corresponding shape changes inmomentum space. Furthermore, the density of states ischanged by the bimodal disorder [27].There are fundamental differences between leg andrung disorder, although eigenfunctions are localized inboth cases. This is understood in leading-order pertur-bation theory: rung disorder leads to triplons hopping ina disordered chemical potential while leg disorder yieldsa disordered nearest-neighbor hopping and therefore to amomentum-dependent disorder potential. Thus localiza-tion is stronger for rung compared to leg disorder.We find two separated energy regions for rung disordercorresponding to J ⊥ and J ⊥ . This can be seen by com-paring to the one-triplon dispersions of a clean QSL withthese J ⊥ ν (white lines in Fig. 2). In first order this followsfrom Gerschgorin’s circle theorem which states that halfof the eigenvalues are bigger and half of them smallerthan one in the thermodynamic limit for that disorder. The spectral weights are higher in the low-energy re-gion due to localization. Considering an eigenstate in thelower (higher) region, it will be localized around stateswith dominantly low (high) rung values. For these eigen-states the average ratio of leg and rung coupling is thenlarger at lower energies. Hence the spectral weight islarger in the low-energy region. Most importantly, weobserve a fragmentation of the one-triplon DSF into en-ergies carrying maximal spectral weights close to k = π .These energies with largest intensity are located wellabove the one-triplon gap of the uniform QSL with thelower rung value. These energies can be understood byconsidering leading-order degenerate perturbation theoryfor small J (cid:107) / ∆ ¯ J ⊥ [27]. In this limit one has a fragmen-tation of the disordered QSL into decoupled ladder seg-ments with constant J ⊥ = 0 . L . In each opensegment one can easily solve the full one-triplon problemand determine the lowest one-triplon energy (cid:15) ( L )1 . For k = π , we find that the one-triplon DSF has indeed localintensity maxima at these energies (see solid white linesin Fig. 2). The approximate quantization of the one-triplon DSF is therefore directly linked to the discreteenergies (cid:15) ( L )1 arising from the strong rung disorder. Notethat similar quantum structures occur also for the high-energy region of the one-triplon DSF which can be againtraced back to a fragmentation of the disordered QSLinto almost decoupled ladder segments. For leg disorderone has only a single energy region which is well describedwith a one-triplon dispersion using mean hopping ampli-tudes (black line in Fig. 2). The spectral weight is largestat k = π and diminishes towards k = 0. However, the in-tensities are even smaller at small k as one expects fromthe spectral weight, since it is distributed over a particu-larly broad range of energies at k = 0. This comes from FIG. 3. Two-triplon contribution to S + ( k, ω ) for pure rung (left) or leg (middle) disorder and to S − ( k, ω ) for pure leg disorder(right). Left : bimodal rung disorder with p = 0 . J ⊥ = 1 . J ⊥ = 0 .
6. The constant leg exchange is J (cid:107) = 0 .
4. The three dashed lines represent the dispersion of the two-triplon bound state for the clean QSL using J ⊥ = 1 . J ⊥ = 1 .
0, and J ⊥ = 0 . Middle/Right : bimodal leg disorder with p = 0 . J (cid:107) = 0 . J (cid:107) = 0 .
4. Thedashed lines represent the dispersion of the two-triplon bound state for the clean QSL with J (cid:107) (white) and J (cid:107) (red). disorder in the local hopping amplitude (like for rung dis-order) which appears in subleading orders leading to ananisotropy between k = 0 and k = π [27]. Further, thefragmentation at k = 0 is of similar origin as for rungdisorder found for all momenta. Finally, the shape of theone-triplon contribution to S + ( k, ω ) can be fully tracedback to the different type of observable. Although thespectral weight vanishes exactly at k = 0, it is distributedalmost uniformly between k = π/ k = π [27]. Thesame is true for the intensities at fixed energy. For thelatter one has to recapitulate that the local observable inreal space only injects a finite weight for single triplonsif an asymmetry of leg couplings is present in the localsurrounding of the observable. The probability of suchan asymmetric and connected region of length L , where O + ( k ) is active, is exponentially descreasing with L . Itfollows that the projection of eigenstates on these regionsis almost flat in momentum space. Two-triplon contribution —
For the clean QSL, it iscontained solely in S + ( k, ω ) while S − ( k, ω ) = 0. It com-prises a continuum and a two-triplon bound state (seeFig. 1b). By far most of the spectral weight is carried bythis bound state, which therefore dominates the DSF.The two-triplon contributions obtained by pCUTs withmaximal pure rung or leg disorder are shown in Fig. 3.Note that bond-operator MF does not yield satisfactoryresults in this sector, since the attractive triplon interac-tion is neglected completely [27]. Localization also occursfor two-triplon states so that eigenstates have almost alltheir weight on a finite part of the two-triplon space inthe position basis [27]. Further, the states carrying mostof the weight will be two-triplon states localized on a fi-nite region of the lattice. For rung disorder, there are3 possible values for the sum of two local hopping am-plitudes resulting in 3 distinct structures in the DSF. The spectral weight of these structures decreases fromlower to higher energy, which can be understood fromthe decreasing average ratio of leg and rung couplingsas described for the one-triplon case above. Each regioncontains a bound state and a continuum, although onlythe bound state carries significant spectral weight (seeFig. 3 left). They gain a finite lifetime due to the disor-der. The same is true for leg disorder. One observes twobound-state structures, which fuse at k = π . The twostructures, one more dispersive than the other, can belinked to the bound states of clean QSLs taking the largeror lower leg coupling (see Fig. 3 middle). Interestingly,the maximal weight at k = π is located at a lower energycompared to the bound states of the clean QSLs. This islikely caused by scattering events of lower-energy boundstates with k < π yielding a final momentum k = π . Fi-nally, the two-triplon contribution to S − ( k, ω ) (see Fig. 3right), which is induced by the leg disorder, has a muchlarger weight compared to the one-triplon contributionto S + ( k, ω ) although the shape is similar. Notably, themaximum intensity is at the same energy as for S + ( π, ω )and is again associated with bound states. Conclusions —
The pCUT method is an efficient toolto investigate the fate of quasi-particles under quencheddisorder, which we exemplified by calculating the DSF ofdisordered QSLs. Disorder-induced quantum structureslike the quantization of collective triplon excitations inthe DSF emerge which are of direct relevance for spectro-scopic experiments. We therefore suspect that our find-ings trigger experimental as well as further theoreticalinvestigations on the fate of quasi-particles in disorderedcorrelated quantum matter.We thank Bruce Normand and Christian R¨uegg forfruitful discussions. MH and KPS acknowledge financialsupport from DFG project SCHM 2511/10-1. upplementary Materials to ”Dynamic structure factor of disordered quantum spinladders”.
Max H¨ormann, Paul Wunderlich and K. P. SchmidtAt first a complementary mean-field approach for the problem is discussed. Then the construction of Hamiltonianand observables with pCUTs is explained in more detail. The way the lifetime of a momentum state changes withdisorder and how this depends on the momentum itself is examined by considering self-energy calculations. For rungdisorder we found an approximate decoupling of the ladder into smaller ladder segments containing only one of thetwo rung values in the main body of the paper. This behaviour is explained with first-order degenerate perturbationtheory here. For all disorder setups of the main body of the paper the spectral weights, the density of states andthe inverse participation ratio is shown. Finally we check the convergence of the pCUT calculations with respect tonumber of rungs N r , number of disorder configurations N dc , and the order of the perturbative expansion. MEAN-FIELD APPROACH
In the framework of the mean-field approach [25, 47, 48] all cubic and quartic terms of the full Hamiltonian, thatwas introduced in the main body of the paper in terms of triplet operators, are discarded at first. Afterwards thehardcore constraint s † ν s ν + (cid:80) α t † ν,α t ν,α = 1 is taken into account by Lagrange multipliers µ ν and furthermore thesinglet expectation value (cid:104) s ν (cid:105) = s ν is introduced. The mean-field Hamiltonian reads (cid:101) H mf ( { J } ) = (cid:101) E mf0 + (cid:88) ν,α (cid:40)(cid:18) J ⊥ ν − µ ν (cid:19) t † ν,α t ν,α + 14 (cid:88) n =1 J (cid:107) ν,n s ν s ν +1 (cid:16) t ν,α t ν +1 ,α + t ν,α t † ν +1 ,α + h.c. (cid:17)(cid:41) (S1)with (cid:101) E mf0 = (cid:80) ν ( − J ⊥ ν s ν − µ ν s ν + µ ν ). To obtain the parameters µ ν and s ν it is necessary to solve self-consistentequations (cid:28) ∂ H mf ( { J } ) ∂µ ν (cid:29) = 0 (cid:28) ∂ H mf ( { J } ) ∂s ν (cid:29) = 0 , (S2)where (cid:104)·(cid:105) denotes the vacuum expectation value. If one tries to do so for a disordered ladder, a 2 N r systemof equations would remain which would be analytically and numerically demanding. Alternatively, self-consistentequations for ladders without disorder, that were also considered within the context of mean-field theory, are used tocalculate the energy properties and then the dynamic structure factor.For such clean ladders the parameters µ ν and s ν are uniform and the hardcore-constraint restricts to a global one.Starting from (S1) with µ ν = µ and s ν = s a Fourier transformation was done to transform the Hamiltonian intomomentum space and a Bogoliubov transformation for diagonalizing the Hamiltonian subsequently. Eventually self-consistent equations can be obtained analogue to (S2) which can be solved numerically, so all energy properties andthe dynamic structure factor of ladders without disorder can be obtained.These self-consistent equations for ladders without disorder are also used for ladders with disorder. To calculate µ ν and s ν , respectively, J ⊥ ν and an averaged value for J (cid:107) ν = (cid:80) n =1 ( J (cid:107) ν − ,n + J (cid:107) ν,n ) was used. Consider, however, if onewould use global parameters µ ν = µ and s ν = s for a disordered ladder, the outcome of this approach would be muchworse than using local parameters because the local fluctuations of the disorder are not sufficiently taken into account.On that basis a Bogoliubov transformation of (S1) can be conducted in a real-space setting H mf ( { J } ) = E mf0 + (cid:88) ν,α ω mf ν γ † ν,α γ ν,α (S3)with the one-particle energies of the Bogoliubov quasi-particles ω mf and E mf0 = (cid:80) ν ( − J ⊥ ν s ν − µ ν s ν + µ ν − J ⊥ ν + ω mf ν ). Starting from this result, ω mf and the corresponding eigenstates can be used to calculate the dynamicstructure factors S ± ( k, ω ) introduced in the paper.The mean-field results for the same bimodal leg and rung disorder setups as in the main body of the paper are shownin Fig. S1. In the one-triplon sector ( S − ( k, ω )), they exhibit a good qualitative accordance to the results obtainedby pCUT calculations. Nevertheless some deviations in shape and in spectral weight appear by construction, butwe see for coupling constants up to J (cid:107) ν,n /J ⊥ ν = 0 . S + ( k, ω ) that contains two-triploncontributions was also calculated using mean-field theory. As expected, only the two-triplon continuum can beobtained and the contribution of bound states is missing because of the neglection of triplon interactions. Note thatthe mean-field approach can be applied to arbitrary disorder configurations. FIG. S1. Uper panel: one-triplon contribution to the dynamic structure factor S − ( k, ω ) using the mean-field approach forbimodal rung disorder with p = 0 . J ⊥ = 1 . J ⊥ = 0 . J (cid:107) = 0 . left ) and bimodal leg disorderwith p = 0 . J (cid:107) = 0 . J (cid:107) = 0 . J ⊥ = 1 ( right ). White lines ( left ) represent one-triplon dispersionsfor a clean QSL with J (cid:107) = 0 . J ⊥ = 1 . J ⊥ = 0 . right ) representsthe one-triplon dispersion with mean hopping amplitudes. Lower panel: two-triplon contribution to the dynamic structurefactor S + ( k, ω ) using the mean-field approach for the same disorder setups as in the upper panel. The three dashed lines ( left )represent the dispersion of the two-triplon bound state for the clean QSL using J ⊥ = 1 . J ⊥ = 1 .
0, and J ⊥ = 0 .
6. In the legdisorder plot the dashed lines represent the dispersion of the two-triplon bound state for the clean QSL with J (cid:107) (white) and J (cid:107) (red). All lines were obtained by pCUTs and are the same as the ones shown in the main body of the paper. The calculationswere broadened with a Lorentzian of Γ = 0 . CONSTRUCTION OF HAMILTONIAN AND OBSERVABLE WITH PCUT
Rung and leg couplings are randomly chosen for a finite lattice of size N r according to a fixed disorder distribution.To construct the Hamiltonian one needs all matrix elements between one- and two-triplon states for that realization ofthe disorder. For a given order of the perturbation theory only linked subclusters with a length of at most that ordercontribute to these matrix elements. For all such subclusters of that particular realization the hopping amplitudes oftriplons or the interactions between two triplons are calculated. Afterwards these amplitudes are properly assignedto the corresponding matrix elements between one- or two-triplon states on the finite lattice. The remaining task isto diagonalize the obtained matrix.For the observable the treatment is similar. At first the local observables O ± ( ν ) and their action on the groundstate is calculated for every rung ν of the disorder realization. This action is not local anymore due to the pCUTtransformation and the observable can create one- or two-triplon states at most distanced by the perturbation orderfrom ν . The amplitudes for the triplon states that are created by the local observables are again created for allsubclusters of the lattice. The local actions have to be Fourier transformed to obtain the global action O ± ( k ). It isconvenient to calculate the action of O ± ( k ) in position space. For that one has to take care to properly account forthe phases caused by the non-local triplon actions of the local triplet observables. MOMENTUM DEPENDENCE OF MOMENTUM STATE LIFETIMES
In first-order Born approximation the mean self-energy isΣ( k, w ) Born = (cid:90) (cid:10) | V ( k, k (cid:48) | (cid:11) G ( k (cid:48) , w )dk (cid:48) . (S4)Here V ( k, k (cid:48) ) = H ( k, k (cid:48) ) for k (cid:54) = k (cid:48) , H is the one-triplon effective Hamiltonian and V ( k, k ) = H ( k, k ) − (cid:104)H ( k, k ) (cid:105) .The retarded unperturbed Green’s function is that of the mean hopping amplitudes, i.e. G ( k (cid:48) , w ) = 1 w − (cid:104)H ( k (cid:48) , k (cid:48) ) (cid:105) + i0 + . (S5)Note that in the thermodynamic limit H ( k (cid:48) , k (cid:48) ) = (cid:104)H ( k (cid:48) , k (cid:48) ) (cid:105) . For w = (cid:104)H ( k (cid:48) , k (cid:48) ) (cid:105) ! = (cid:104)H ( k, k ) (cid:105) the Dirac identityyields us two contributions to the imaginary part of the self-energy (assuming that k = 0 , π are the only max-ima/minima) Im (cid:0) Σ( k, (cid:104)H ( k, k ) (cid:105) ) Born (cid:1) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ddk (cid:48) (cid:104)H ( k (cid:48) , k (cid:48) ) (cid:105) | k (cid:48) = k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:0) Var( H ( k, k )) + (cid:10) |H ( k, π − k ) | (cid:11)(cid:1) . (S6)Thus the imaginary part of the self-energy at ( k, (cid:104)H ( k, k ) (cid:105) ) and the approximate broadening of the momentum state | k (cid:105) scales proportionally to Var( H ( k, k )) + (cid:10) |H (2 π − k, k ) | (cid:11) . (S7)Although | ddk (cid:48) (cid:104)H ( k (cid:48) ,k (cid:48) ) (cid:105)| k (cid:48) = k | = 0 the broadening at k = 0 , π is assumed to scale similarly.For increasing (cid:68) J (cid:107) ν (cid:69) there are big differences in the momentum lifetime at k = 0 and k = π . This can be seen in themain body of the paper for the leg disorder case in the one-triplon sector (Fig. 2, middle), where momentum stateshad a significantly longer lifetime at k = π than at k = 0. In Fig. S2 the scaling factor (S7) for the momentum statelifetimes described above is shown for the same leg disorder setup. The important values are those that coincide withthe mean dispersion (cid:104)H ( k (cid:48) , k (cid:48) ) (cid:105) depicted as dashed line in the figure. Clearly, one recognizes that at k = π momentumstates will have a significantly longer lifetime than at k = 0. In the following the matrix elements of H ,ν,ν (cid:48) are denotedas a ν,ν (cid:48) . As cause for the different behaviour at the mean band edges correlations between the local hopping term a ν, and other hopping terms were identified. To estimate the importance of the different correlations between hoppingelements on the lifetime of momentum states, the following relation proved useful.For correlated but stationary disorder, i.e. (cid:68) J (cid:107) ν J (cid:107) ν + δ (cid:69) (cid:54) = 0, (cid:10) J ⊥ ν J ⊥ ν + δ (cid:11) (cid:54) = 0 and (cid:68) J (cid:107) ν J ⊥ ν + δ (cid:69) (cid:54) = 0 but independent of ν ,the mean absolute square of the Hamiltonian in the momentum basis can be given using the cross-correlation theorem.Starting point is H transformed into the momentum basis. H ,k (cid:48) ,k = 1 N r (cid:42)(cid:88) ν (cid:48) e − i ν (cid:48) k (cid:48) | ν (cid:48) (cid:105) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) H (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88) ν e − i νk | ν (cid:105) (cid:43) = 1 N r δ k,k (cid:48) (cid:88) ν (cid:32) a ν, + 2 (cid:88) t =1 a ν,t cos( tk ) (cid:33) + 1 N r (1 − δ k,k (cid:48) ) (cid:88) ν e i ν ( − k + k (cid:48) ) (cid:32) a ν, + (cid:88) t =1 a ν,t ( e i tk (cid:48) + e − i tk ) (cid:33) . (S8)All that has to be known to use the cross-correlation theorem is the covariance structure of the hopping amplituderandom processes. One obtains N r (cid:10) |H ,k,k (cid:48) | (cid:11) = 2 (cid:88) d (cid:88) t,t (cid:48) =1 t>t (cid:48) Cov ( a ν,t a ν + d,t (cid:48) ) (cid:18) cos (( t (cid:48) − t ) k + d ( k − k (cid:48) ))+ cos (( t − t (cid:48) ) k (cid:48) + d ( k − k (cid:48) )) + cos ( t (cid:48) k + tk (cid:48) + d ( k − k (cid:48) )) + cos ( − t (cid:48) k (cid:48) − tk + d ( k − k (cid:48) )) (cid:19) + 2 (cid:88) d (cid:88) t =1 Cov ( a ν, a ν + d,t ) (cos ( tk + d ( k − k (cid:48) )) + cos ( − tk (cid:48) + d ( k − k (cid:48) )))+ (cid:88) d ≥ (cid:88) t =0 (2(1 − δ d, ) cos ( d ( k − k (cid:48) )) + δ d, ) · ( δ t, Cov ( a ν, a ν + d, ) + (1 − δ t, )Cov ( a ν,t a ν + d,t ) 2(1 + cos ( t ( k + k (cid:48) ))) . (S9) FIG. S2. (cid:104)|H ( k, k (cid:48) ( ω )) |(cid:105) + (cid:10) |H (2 π − k, k (cid:48) ( ω )) | (cid:11) is shown as a function of k and ω for the leg disorder setup of the main body ofthe paper. Note that ω ( k (cid:48) ) = (cid:104)H ( k (cid:48) , k (cid:48) ) (cid:105) . Further at the position of the dashed curve Var( H ( k, k (cid:48) ( ω )))+ (cid:10) |H (2 π − k, k (cid:48) ( ω )) | (cid:11) is shown. The dashed line depicts the dispersion of (cid:104)H ( k, k ) (cid:105) . With this formula, which is also valid in the thermodynamic limit, it is easy to estimate the importance of correlationsbetween hopping elements for the lifetime of momentum states.
LEADING-ORDER EFFECTS OF RUNG DISORDER
For strong rung disorder one can do degenerate perturbation theory around the unperturbed Hamiltonian of onlylocal triplon interactions. This yields in first order two Hamiltonians without interactions between rungs with differentrung couplings. The low-energy physics is described by the Hamiltonian with small rung values and the high-energyphysics by the one with the large values. An important issue of this perturbative picture are the lengths of chainsof consectutive low or high rung values. The mean amount of such chains of length L is N r (1 − p ) p L if p is theprobability for the specific rung value and N r = ∞ . For a model with only nearest-neighbour hopping all the physicsin first-order degenerate perturbation theory is happening in those finite ladder segments. When hopping to furtherneighbours is taken into account or when higher orders in the degenerate perturbation theory are not neglected thosefinite chains interact with each other. Fig. S3 shows the dynamic structure factor at k = π for the rung disordersetup of the main body of the paper and in the one-triplon sector (blue curve). In red the result obtained withfirst-order degenerate perturbation theory in the two different rung values applied on the one-triplon Hamiltonianused to calculate the blue curve is shown. The vertical black dashed lines show the smallest energies of consecutive J ⊥ = 0 . L with open boundary conditions. For L = 7 , , , , (3) these black lines match well withboth the maxima of the blue and red curves. This indicates that the discrete nature of the length of these laddersegments is the main reason for the appearance of multiple maxima in Fig. S3 and the approximate quantization ofthe DSF that was seen in the main body of the paper.The intensities of the red curves are biggest at L = 1 ,
2. This is what one would find assuming that the weight at k = π is proportional to L and neglecting hopping further than to next neighbours. Then the intensities are ∝ L N r (1 − p ) p l when only the smallest energy (cid:15) ( L )1 contributes to the DSF at k = π . There are several reasons why the intensities ofthe exact case differ. There is a shift towards lower energies caused by the interactions between several finite laddersegments. The blue curve rather consists of peaks on a plateau of constant intensity. It is caused by the splitting ofenergies induced by these additional interactions. Because more than L rungs can carry the weight of eigenfunctionsdue to additional interactions neglected in first-order degenerate perturbation theory the weight at k = π can scalestronger than with L . This is true for all L . However, because for a segment of length L = 1 the eigenstate is alocalized rung and thus non-dispersive, there is no preference in the L = 1-sector for k = π . In contrast, this isdifferent for larger L and is likely another reason why the blue curve shows stronger intensities there. FIG. S3. In blue the dynamic structure factor at k = π is depicted for the rung disorder setup of the main body of thepaper with parameters N r = 100, N dc = 1000 and Γ = 0 .
01. The red curve shows results by applying first-order degenerateperturbation theory on the one-triplon Hamiltonian used to calculate the blue curve. The vertical black lines indicate thelowest energy (cid:15) ( L )1 of a ladder segment of L rungs with value J ⊥ = 0 . SPECTRAL WEIGHTS
Fig. S4 shows the spectral weights S ± ( k ) := (cid:82) S ± ( k, ω )d ω obtained by pCUTs for all disorder setups discussed inthe main body of the paper. Interestingly, the spectral weights of the corresponding models with constant mean rungand leg couplings show only minor differences to the disodered one in the one-triplon sector. In the two-triplon sectorthe shape is similar but the scale is smaller.The observables S + ( k ) ( S − ( k )) in the one- respectively two-triplon sector that only arise due to the broken reflectionsymmetry around the centerline of the ladder have a very wide maximum around k = π showing only small k -dependence there. For k = 0 the weight of these observables is zero. The symmetric observable is zero at k = 0because it commutes with the Hamiltonian. The antisymmetric observable in the two-triplon sector also goes tozero at k = 0. For that, note that total spin-one two-triplon states have odd parity with respect to inversion of thelattice around their center of mass. At the same time O − (0) and the singlet reference state have even parity for thatsymmetry and since these parities are preserved by the pCUT O − (0) shows no weigth in that sector.One can conclude that the huge differences in the dynamic structure factor between the disordered and the non-disordered case are dominantly caused by changes in the dynamics and not by the spectral weights. FIG. S4. The spectral weights S ± ( k ) are shown as solid lines for all sorts of disorder considered in the paper (denoted as S ± ( J ⊥ ν , J (cid:107) ν )). The spectral weights with constant mean rung and leg values are plotted as dashed lines for the same disordersetups (denoted as S ± ( (cid:104) J ⊥ ν , J (cid:107) ν ) (cid:105) ). DENSITY OF STATES
The density of states (DOS) was calculated numerically for all disorder configurations { J } using a bin width of0 . N dc = 1000 samples of size N r = 100. It is normalized such that its integral over energy is1. The result is shown in the upper panel of Fig. S5. The DOS has large influence on the dynamic structure factorfor rung disorder in the one-triplon sector. For that disorder many regions with vanishing density of states occur.Also in the non-vanishing energy regions the DOS shows strong oscillatory behaviour for one-triplon rung disorder.For leg disorder the DOS is mainly smooth. At higher energies we see oscillatory behaviour again that correspondsto eigenfunctions centered around k = 0 where local hopping amplitudes play a larger role.In the two-triplon sector rung disorder does not lead anymore to regions of vanishing DOS. Still, we find an oscillatorybehaviour especially at larger energies. In contrast, for leg disorder differences to the non-disordered case becomesmall and we see an almost-smooth DOS. INVERSE PARTICIPATION RATIO
The inverse participation ratio (
IP R ) IP R = (cid:88) ν |(cid:104) n | ν (cid:105)| (S10)with | ν (cid:105) denoting position states is a simple and intuitive measure for the localization length of a normalized eigen-function | n (cid:105) . Suppose | n (cid:105) is a plane wave eigenstate of a periodic one-dimensional chain with one atom in the unitcell. Then the IP R will be 1 / N r , N r being the length of the chain. For all other extended states the IPR will bebigger but always remain ∝ / N r [49]. For a perfectly localized eigenstate with | n (cid:105) = | ν (cid:105) the value of the IP R isone. An exponentially localized state |(cid:104) n | ν (cid:105)| ∝ exp (cid:18) − | ν − ν | ξ (cid:19) (S11)has an IP R of (cid:80) ν exp (cid:16) − | ν − ν | ξ (cid:17)(cid:16)(cid:80) ν exp (cid:16) − | ν − ν | ξ (cid:17)(cid:17) = 2 − exp ( − ξ ) − (cid:18) − exp ( − ξ ) − (cid:19) ≈ ξ . (S12)The approximation used in the last step is already quite good for localization lengths ξ ≥
1. With that the localizationlength can be estimated as ξ ≈ / (4 IP R ) . The IP R ranges between 1 for a perfectly localized state and 0 for extendedstates as the system size N r → ∞ .For two-particle position states | ν, ν + δ (cid:105) and their corresponding eigenfunctions | n (cid:105) in position space we use ananalogous quantity to IP R : IP R = (cid:88) ν,δ |(cid:104) n | ν, ν + δ (cid:105)| (S13)Clearly IP R goes to zero with system size for extended states as well and only remains finite if the eigenfunctionsare localized on a finite subset of two-particle position states.The inverse participation ratio and its two-particle generalization are plotted for N r = 100 and N dc = 1000 (solidlines) and N r = 50 and N dc = 2000 (circles) in the lower panel of Fig. S5. In both cases a bin width of 0 . N r = 100 to N r = 50. Theconclusion is that the energy states fit with almost all their weight already on 50 rungs and that there should behardly any finite-size effects anymore for systems of 100 rungs or more.For leg disorder the IP R , increases towards higher energies. These have biggest weight on momentum k = 0. Higherorder effects lead to stronger broadening at k = 0 and this explains why the IP R , grows towards higher energies.The IP R in the two-triplon sector shows essentially the same features as the IP R in the one-triplon sector. Westress that it remains finite so that also two-triplon states have almost all their weight on a finite part of the two-triplon space in the position basis. One main difference to the IP R is that it varies less strongly for the bimodalrung distribution.The average IP R , is bigger in the rung disorder case. FINITE-SIZE EFFECTS, CONVERGENCE OF PERTURBATIVE EXPANSION AND STATISTICALERROR
As the dynamic structure factor could only be calculated for finite systems possible finite-size effects have to beexamined. Further the convergence of the perturbative expansion has to be checked. Due to the finite number ofsamples N dc a statistical error arises which has to be quantified.In Fig. S6 the dynamic structure factor at k = π is shown for all disorder setups of the main body of the paper. The FIG. S5. Upper panel: the inverse participation ratios
IP R , are plotted for N r = 100 and N dc = 1000 (solid lines) and N r = 50 and N dc = 2000 (circles) for all disorder setups used in the main body of the paper. Lower panel: the density of states(DOS), which is normalized with respect to energy, is shown for the same disorder setups and N r = 100 and N dc = 1000. Thebin width used to calculate both quantities was 0 . blue line shows results for N dc = 100 and N r = 100. The yellow dotted curve was calculated with the same disordersamples but with order 7 in the perturbative expansion of matrix elements of the Hamiltonian and order 6 for theobservables only. For leg disorder hardly any differences can be made out. Thus the perturbative expansion is wellconverged there. For rung disorder at the peaks differences can be seen which are due to a small offset in the energyof the peaks between the order 8 and 7 calculation. Anyway these offsets are quite small and the convergence ofperturbation theory is also expected to be good for rung disorder.Fig. S6 further shows the difference of the calculation with N dc = 200 and N r = 50 and N dc = 100 and N r = 100 byusing the same disorder samples again (red curve). For that each sample of the N r = 100 calculation was cut intotwo samples of N r = 50. The deviations found are quite small and biggest in the one-triplon leg disorder setups.The calculation can only be done for a finite number of samples and a finite width of energy bins or by folding witha Lorentzian of finite width Γ. Assuming that the value in each bin is a random variable and that the realizations ofdisorder are independent from sample to sample the behaviour of its standard deviation with the number of samples N dc , the number of rungs in that samples N r and broadening Γ of the Lorentzian is ∝ / √N dc N r Γ.Fig. S7 shows the relative standard deviation σ [ S ± ( π, ω )] / (cid:104) S ± ( π, ω ) (cid:105) for all disorder setups used in the main bodyof the paper and for parameters N r = 100 and N dc = 1000 (blue circles). Data points where the mean values weresmaller than 1 /
10 of the maximum value were excluded. One can see that the statistical error is smaller than 5% forthose values. ∗ [email protected] † [email protected] ‡ [email protected][1] T. Vojta, J. Phys. A , R143 (2006).[2] T. Vojta, J. Low Temp. Phys. , 299 (2010).[3] R.B. Griffith, Phys. Rev. Lett. , 17 (1969).[4] D. M. Basko, I. L. Aleiner and B. L. Altshuler, Ann. Phys. , 1126 (2006).[5] V. Oganesyan and D. A. Huse, Phys. Rev. B , 155111 (2007). [6] M. Znidaric, T. Prosen and P. Prelovsek, Phys. Rev. B , 064426 (2008).[7] A. Pal and D. A. Huse, Phys. Rev. B , 174411 (2010).[8] A. Furusaki, M. Sigrist, P.A. Lee, K. Tanaka, and N. Nagaosa, Phys. Rev. Lett. , 2622 (1994).[9] E. Westerberg, A. Furusaki, M. Sigrist, and P.A. Lee, Phys. Rev. Lett. , 4302 (1995).[10] N. Nagaosa, A. Furusaki, M. Sigrist, and H. Fukuyama, J. Phys. Soc. Jpn. , 3724 (1996).[11] T. Miyazaki, M. Troyer, M. Ogata, K. Ueda, and D. Yoshioka, J. Phys. Soc. Jpn. , 2580 (1997).[12] M. Greven and R.J. Birgeneau, Phys. Rev. Lett. , 1945 (1998).[13] R. M´elin, Y.-C. Lin, P. Lajko, H. Rieger, and F. Igloi, Phys. Rev. B , 104415 (2002).[14] J.A. Hoyos and E. Miranda, Phys. Rev. B , 214411 (2004).[15] K. Trinh, S. Haas, R. Yu, T. Roscilde, Phys. Rev. B , 035134 (2012).[16] K. Trinh and S. Haas, Phys. Rev. B , 075137 (2013).[17] H. Manaka, A.V. Kolomiets, and T. Goto, Phys. Rev. Lett. , 077204 (2008).[18] H. Manaka, H.A. Katori, A.V. Kolomiets, and T. Goto, Phys. Rev. B , 092401 (2009).[19] Tao Hong, A. Zheludev, H. Manaka, and L.-P. Regnault, Phys. Rev. B , 060410 (2010).[20] M.B. Stone, A. Podlesnyak, G. Ehlers, A. Hug, E.C. Samulon, M.C. Shapiro, and I.R. Fisher, J. Phys. Condens. Matter , 416003 (2011).[21] D. H¨uvonen, S. Zhao, M. Mansson, T. Yankova, E. Ressouche, C. Niedermayer, M. Laver, S.N. Gvasaliya, and A. Zheludev,Phys. Rev. B , 100410 (2012).[22] D. H¨uvonen, S. Zhao, G. Ehlers, M. Mansson, S.N. Gvasaliya, and A. Zheludev, Phys. Rev. B , 214408 (2012).[23] B. Nafradi, T. Keller, H. Manaka, U. Stuhr, A. Zheludev, and B. Keimer, Phys. Rev. B (R) , 020408 (2013).[24] O. Motrunich, K. Damle, D.A. Huse, Physical Review B , 134424 (2001).[25] M. Vojta, Phys. Rev. Lett. , 097202 (2013).[26] Y.-R. Shu, M. Dupont, D.-X. Yao, S. Capponi, and A.W. Sandvik, Phys. Rev. B , 104424 (2018).FIG. S6. In blue the dynamic structure factor at k = π is shown for N dc = 100 samples of N r = 100. The red dashed linesdisplay the difference of these values and the ones obtained with N dc = 200 samples of N r = 50 but with the same disordersamples. Yellow dotted lines show the dynamic structure factor for the same disorder samples and N r = 100 but with order 7 inthe perturbative expansion of the matrix elements of the Hamiltonian and order 6 for the observables only. The same disordersetups as in the main body of the paper were chosen and the same Lorentzian curve for broadening was used (Γ = 0 . S − , (b) leg disorder S − , (c) leg disorder S + . Two-triplon sector: (d) rung disorder S + ,(e) leg disorder S + , (f) leg disorder S − . FIG. S7. In blue the relative standard deviation σ [ S ± ( π, ω )] / (cid:104) S ± ( π, ω ) (cid:105) of systems with N r = 100 and N dc = 1000 is shown.The same disorder setups as in the main body of the paper were chosen and the same Lorentzian curve for broadening wasused (Γ = 0 . /
10 of the maximum value were excluded. One-triplonsector: (a) rung disorder S − , (b) leg disorder S − , (c) leg disorder S + . Two-triplon sector: (d) rung disorder S + , (e) leg disorder S + , (f) leg disorder S − .[27] The Supplementary Materials contain detailed information on the improved bond-operator mean-field theory and cor-responding results for the DSF, the momentum dependence of momentum state lifetimes, a detailed discussion of theleading-order effects for rung disorder, pCUT results for the spectral weight, density of states, and inverse participationratio, as well as the convergence properties of the pCUT approach.[28] E. Dagotto and T. M. Rice, Science , 618 (1996).[29] M. Windt, M. Gr¨uninger, T. Nunner, C. Knetter, K.P. Schmidt, G.S. Uhrig, T. Kopp, A. Freimuth, U. Ammerahl, B.B¨uchner, and A. Revcolevschi, Phys. Rev. Lett. , 127002 (2001).[30] S. Notbohm, P. Ribeiro, B. Lake, D.A. Tennant, K.P. Schmidt, G.S. Uhrig, C. Hess, R. Klingeler, G. Behr, B. B¨uchner,M. Reehuis, R.I. Bewle, C.D. Frost, P.Manuel, and R.S. Eccleston, Phys. Rev. Lett. , 027403 (2007).[31] Ch. R¨uegg, K. Kiefer, B. Thielemann, D. F. McMorrow, V. Zapf, B. Normand, M. B. Zvonarev, P. Bouillot, C. Kollath,T. Giamarchi, S. Capponi, D. Poilblanc, D. Biner, K. W. Kr¨amer, Phys. Rev. Lett. , 247202 (2008).[32] T. Hong, Y. H. Kim, C. Hotta, Y. Takano, G. Tremelling, M. M. Turnbull, C. P. Landee, H.-J. Kang, N. B. Christensen,K. Lefmann, K. P. Schmidt, G. S. Uhrig, C. Broholm, Phys. Rev. Lett. , 137207 (2010).[33] S. Ward, P. Bouillot, H. Ryll, H. Kiefer, K.W. Kr¨amer, Ch. R¨uegg, C. Kollath, and T. Giamarchi, J. Phys. Condens.Matter , 014004 (2013).[34] D. Schmidiger, P. Bouillot, T. Guidi, R. Bewley, C. Kollath, T. Giamarchi, and A. Zheludev, Phys. Rev. Lett. , 107202(2013).[35] S. Ward, P. Bouillot, C. Kollath, T. Giamarchi, K.P. Schmidt, B. Normand, K.W. Kr¨amer, D. Biner, R. Bewley, T. Guidi,M. Boehm, D.F. McMorrow, Ch. R¨uegg, Phys. Rev. Lett. , 177202 (2017).[36] T. Barnes, E. Dagotto, J. Riera, and E. S. Swanson, Phys. Rev. B , 3196 (1993).[37] D. G. Shelton, A. A. Nersesyan, and A. M. Tsvelik, Phys. Rev. B , 8521 (1996)[38] S. Trebst, H. Monien, C.J. Hamer, Z. Weihong, and R.R.P. Singh, Phys. Rev. Lett. , 4373 (2000).[39] C. Knetter, K.P. Schmidt, M. Grueninger, and G.S. Uhrig, Phys. Rev. Lett. , 167204 (2001).[40] K. P. Schmidt, C. Knetter, und G. S. Uhrig, Eur. Phys. Lett. , 877 (2001).[41] A. Lauchli, G. Schmidt, M. Troyer, Phys. Rev. B , 100409 (2003).[42] K. P. Schmidt und G. S. Uhrig, Phys. Rev. Lett. , 227204 (2003). [43] Simon Ward, PhD thesis entitled Spin ladder physics and the effect of random bond disorder , University College London(2014).[44] C. Knetter und G. S. Uhrig, Eur. Phys. J. B , 209 (2000).[45] C. Knetter, K. P. Schmidt, und G. S. Uhrig, J. Phys. A , 7889 (2003).[46] K. Coester and K.P. Schmidt, Physical Review E , 022118 (2015).[47] S. Gopalan, T.M. Rice, M. Sigrist, Phys. Rev. B , 8901 (1994).[48] B. Normand, Ch. R¨uegg, Phys. Rev. B , 054415 (2011).[49] D. J. Thouless, Physics Reports13