Thermodynamic properties of disordered quantum spin ladders

In this paper, we study the thermodynamic properties of spin-$1/2$ antiferromagnetic Heisenberg ladders by means of the stochastic series expansion quantum Monte Carlo technique. This includes the thermal properties of the specific heat, uniform and staggered susceptibilities, spin gap, and structure factor. Our numerical simulations are probed over a large ensemble of random realizations in a wide range of disorder strengths $r$, from the clean ($r=0$) case up to the diluted ($r \rightarrow 1$) limit, and for selected choices of number of legs $L_y$ per site. Our results show some interesting phenomena, like the presence of crossing points in the temperature plane for both the specific heat and uniform susceptibility curves which appear to be universal in $r$, as well as a variable dependence of the spin gap in the amount of disorder upon increasing $L_y$.


Introduction
Understanding the effect of quenched randomness on zero-and finite-temperature properties of spin models is one of the most intriguing problems in theoretical physics.Recent model systems including quenched randomness that are under investigation include the spin-1/2 Heisenberg spin chains [1,2], the spin-1/2 J − Q model on the square lattice [3], and quantum spin chains with power-law long-range antiferromagnetic couplings [4].Low-dimensional spin systems have also been at the centre of intense investigation due to the development of advanced theoretical, computational and experimental methods that paved the way for new results [5,6,7,8,9].Among the wealth of low-dimensional systems, quantum spin systems with antiferromagnetic interactions share plentiful physical properties, even in one dimension.For example, Haldane's conjecture that antiferromagnetic spin chains with integer spin values exhibit a gapped spectrum has been backed up by solid theoretical [10], numerical [11], and experimental [12] studies.It is worth noting here that the spin-1/2 Heisenberg coupled chains with an even number of legs have a finite spin gap (∆) to the lowest triplet excitation.Also, some ladder systems have an exponentially decaying spin-spin correlation function and uniform susceptibility, manifesting the existence of a spin gap [13,14].The width of the spin gap for a two-leg ladder spin system with the isotropic coupling constant J can be roughly estimated as ∆ ∼ 0.5J via quantum Monte Carlo techniques [15].What is more, the value of the spin gap can be modified or even completely eliminated in various ways, among others, by the choice of a different lattice topology [16,17], the diffusion of disorder in the spin-spin coupling [18], or the application of an external magnetic field [19].
Different ladder models have been studied in the literature, such as spin-ladder systems with dimerization [20,21,22,23], zig-zag ladders [24,25], and mixed ladders [26,27,28].The presence of quenched randomness is expected, in most cases, to have a significant effect on the thermal and magnetic properties of the system, even at the low disorder limit.Weakly disordered anisotropic spin-1/2 ladders have been handled perturbatively [29] and critical properties of strongly disordered systems have been analyzed benefiting from the strong disorder (density matrix) renormalization-group method [16,30,31].
Modern quantum Monte Carlo techniques are powerful and robust tools for the study of disordered spinladder systems.For instance, the stochastic series expansion quantum Monte Carlo technique has already been successfully used to investigate the spin-1/2 Heisenberg ladders under the presence of quenched bond randomness [32], indicating that neighbouring bond energies change sensitively with the position of disorder in the spin-spin coupling term.Furthermore, in Ref. [33], some unusual and interesting effects of disorder on collective excitations have been reported with the calculation of the ground-state dynamic structure factor for a ladder system with bond disorder along the legs and rungs of the ladder.More recently, the thermodynamic properties of a two-leg quantum spin-ladder system including quenched bond randomness only along the rung direction of the system have been discussed [34], for varying values of the system parameters.
In this framework, we attempt to provide here an additional insight into some critical aspects of disordered quantum spin-ladder systems, not fully considered previously.Therefore, the goal of the present study is threefold: (i) Monitor the trend of the crossing points of the specific heat and uniform susceptibility curves, (ii) scan the response of the spin-gap parameter, and (iii) scrutinize the low-temperature behaviour of the staggered susceptibility and structure factor, while increasing the number of legs of the spin-ladder system and for different values of the disorder parameter.Whenever it is possible, we compare our results with the main outcomes of Ref. [34].Before going on further, we would like to underline that the existence of the crossing point phenomenon, observed also in the present work, is one of the most interesting and puzzling effects in strongly correlated electron systems.Up to now, and to the best of our knowledge, no unified mechanism for this phenomenon has been proposed.However, it is generally believed that crossing points occur in certain systems that are in the vicinity of a quantum or a second-order phase transition, or some physical systems showing magnetic instability, such that their properties are prominently sensitive to changes in thermodynamic variables, such as the magnetic field and the temperature.
The structure of the paper is as follows: In Section 2 we briefly describe the model and the implementation of the stochastic series expansion quantum Monte Carlo method.In Section 3 we provide the simulation details and the physical observables studied here.The numerical results are presented and analysed in Section 4. Finally, Section 5 contains a summary and a brief outlook.
2 Model and method 2.1 Antiferromagnetic spin-1/2 random-bond Heisenberg spin ladders We introduce the Hamiltonian of the quantum spinladder model in a general framework in order to align with the formulation of the stochastic series expansion technique for convenience.Namely, the following Hamiltonian describes spin models that include a number of bonds N b with site index i(b) and spin operators S i(b) interacting with coupling strength J b .For the present spinladder system with N sites, we consider all spin-spin interactions between nearest neighbours to be antiferromagnetic (J b > 0).Note that for the bonds connecting the spins along the legs of the spin ladder, we have J b = J.The bonds located at the rungs of the system are randomly selected from a uniform bimodal distribution of the form where (J 1 + J 2 )/2 = 1 and J 1 > J 2 > 0 following previous practice [34,35,36,37].The present study concerns only the case of 50%/50% weak/strong bonds following the traditional bond disorder implementation commonly used in the case of classical Ising ferromagnets with disorder [38].With the help of the usual control parameter r = (J 1 − J 2 )/2 that reflects the strength of the bond randomness on the rungs of the system we define J 1 = 1 + r (> 1) and J 2 = 1 − r (< 1).An example of the quenched bond disorder in the system is shown in Figure 1.For r = 0 the clean spin-ladder system is recovered, whereas the case r = 1 corresponds to the bond diluted case among the rungs of the system and is therefore beyond the scope of this study.

Stochastic series expansion quantum Monte Carlo method
Within the framework of the stochastic series expansion technique for the isotropic S = 1/2 Heisenberg J 1 J 2 J Fig. 1 An example of the quenched bond-disorder configuration on a quantum spin-ladder system.The couplings J 1 and J 2 along the rungs are depicted by solid blue and red dotted lines, respectively.The interactions J are represented by solid black lines.
antiferromagnetic model, it is possible that bond operators can be disentangled into diagonal (H 1,b ) and off-diagonal (H 2,b ) terms as follows Benefiting from equation ( 3), we can thus rewrite Hamiltonian ( 1) as We note here that the constant energy term is not necessary for the implementation of the algorithm (yet, it should be added when calculating the energy).
According to the stochastic series expansion quantum Monte Carlo protocol [39,40,41], the partition function of the system can be Taylor expanded with a chosen spin basis as follows The summation appearing above refers to all configurations α and all possible operator strings S L , including the unit operator H 0,0 to ensure that the length of the strings is fixed and also a unit bond coupling for convenient implementation; for this extra index, we have J 0 = 1.Here, n is the number of non-unit operators in the string, while n 2 denotes the number of off-diagonal operators.Finally, as usual, β is the reduced inverse temperature.Since the present system includes random bonds among the rungs, the non-zero weights are bonddependent for an allowed spin configuration.Hence, the weight W (α, S L ) can be expressed in the following manner 3 Numerical simulations and observables

Simulation details
Using the above scheme we generated numerical data for the quenched random-bond spin-ladder system for different system sizes at a wide temperature spectrum and for several values of the disorder strength parameter within the range r = {0 − 0.9}.As already mentioned above, our geometrical setup consists of L x × L y systems, where L x defines the linear dimension along the legs of the lattice and is fixed at L x = 256 during all simulations.We note in passing that we have performed some additional test simulations using L x values up to 512 to identify possible finite-size and crossover effects in the observed quantities.However, our analysis, not shown here for brevity, disclosed that the results are almost L x -independent so the considered L x = 256 dimension is enough for an accurate determination of the properties.L y now stands for the number of legs of the ladder system.In fact, one of the most interesting aspects of the physics of quantum Heisenberg ladders is that their low-energy properties are distinct when comparing even and odd L y values [42].It is wellknown that there is a spin gap for the case of even L y values, while ladder systems having odd L y values are gapless [39].In this framework, we only studied systems with L y = {2, 4, 6}.Closing, some numerical details: We performed simulations over 500 independent random realizations for each pair of (r, L y ) values.We applied boundary conditions which are periodic along the legs and free along the rungs of the spin ladder.In our protocol, the first 10 5 Monte Carlo steps were discarded during the thermalization process and numerical data were collected and analyzed during the following 5 × 10 5 Monte Carlo steps.We note in passing that the simulation time needed for a single realization on a node of a Dual Intel Xeon E5-2690 V4 processor was approximately 30 minutes for the largest system size (L x = 256, L y = 6) and the minimum temperature value that we have reached during the simulations.Finally, errors were estimated using the standard jackknife method [43,44].

Observables
The specific heat C of the system is straightforward to estimate with the help of the number of non-unit operators n in the operator sequence [39] Also, by constructing estimators from the Kubo integral, we gain access to static susceptibilities via [39] χ AB = Here the integrand denotes the ensemble average of an imaginary time-dependent product with operators A(τ ) = e τ H A(0) e −τ H .For the case of diagonal operators A and B with eigenvalues a(k) and b(k), respectively, the static susceptibility is retrieved by including eigenvalues from all propagated states, as follows [45,46] For the total magnetization M , equation ( 9) reduces to the uniform susceptibility χ u with a(k A low-temperature form of the uniform susceptibility of the spin-ladder system with an even number of rungs L y is given by [39] where ∆ is the well-known spin gap that can be extracted by linearization of equation (11) and making use of a least-squares fit -typically and as it will be shown below, ground-state values of ∆ are obtained while varying (r, L y ) parameters.
In addition to the physical observables introduced above, we also scrutinized in this work the thermal variations of the staggered susceptibility χ s and structure factor S. For the staggered magnetization M s , equation (9) with a(k) = b(k) = M s (k) provides the staggered susceptibility Finally, the staggered structure factor can be computed via where N = L x × L y the total number of spins on the system.

Results
We start the presentation of our results with Figure 2, showcasing the effects of disorder on the specific heat C of the system for the three selected values of L y considered in this work.As it is evident from all panels, the specific-heat maxima decrease with increasing values of the disorder-strength ratio r for all given values of L y .In fact, the broad maxima of the curves shift towards a higher temperature regime and become flatter as r gets larger.Interestingly, the considered system has a particular temperature point where the specific-heat curves intersect, taking the same value independent of r, but dependent on L y .The inset in all three panels highlights the fine sweeping around these points and the vertical dashed lines mark their positions in the temperature axis, denoted as (T /J) (C) cross .To obtain the crossing points for the specific-heat curves we used numerical data obtained for six [r, r + 0.4] pairs of the disorderstrength ratio, where r = {0, 0.1, 0.2, 0.3, 0.4, 0.5}.We quote the estimates (T /J) (C) cross = 0.946(3), 1.004(1), and 1.018(2) for L y = 2, 4, and L y = 6, respectively.Our analysis indicates that upon increasing L y the location of the crossing point shifts to higher temperatures.Note that the results obtained here for the case with L y = 2 agree nicely within error bars with those from the analysis of Ref. [34].
We would like to underline that for this class of models including randomness on the local bonds, the presence of crossing points has been identified and studied in the context of experimental [47,48,49,50] and theoretical [51,52,53] systems.What is more, their physical origin has been scrutinized in detail for lattice models and continuum systems [47], and a wealth of numerical outcomes has been obtained in particular for the case of the half-filled Hubbard model at all dimensions [48].Based on these studies one may safely conclude that the specific-heat values at these distinctive points are characterized by an almost universal value of ∼ 0.34/k B , whereas the corresponding temperatures are different at all dimensions.It is also worth stressing here that, for the present system, the rate of change in the specific-heat values as a function of the disorder parameter changes its sign at the intersection leading to a total entropy change of zero.Thereby, the crossing point of the specific heat may be also considered as an inflection point.
In analogy to Figure 2, we present in Figure 3 the uniform susceptibility (χ u )-curves.Several comments are in order at this point: 1.At high temperatures, χ u curves obtained for varying values of r are almost identical and the typical Curie behaviour is present in the system for all values of L y .2. In contrast to C, the broad maxima of χ u shift to a lower temperature regime with increasing r, also for all values of L y .3. Similarly to the observations recorded in Figure 2, there exist crossing points in the χ u curves; see the corresponding insets of Figure 3 for a zoom in of the criss-crossing area.To determine these points in χ u we followed the procedure outlined above for the case of the specific heat.We quote the estimates (T /J) We underline here that the absolute difference as L y increases so that these discrepancies can be possibly attributed to the presence of finite-size effects. 5. Finally, an exponentially decreasing shift is noticeable at the relatively low-temperature regimes for all L y and r values considered which clearly manifests the presence of a spin gap [42].values of r and L y at the low-temperature regime, i.e., (T /J) ≤ 10 −1 .The solid and dashed lines are fits of the form (11) within the temperature range 5 × 10 −2 ≤ (T /J) ≤ 10 −1 where a clear linear response is observed.These extrapolations allow us to retrieve the spin-gap values ∆ as (T /J) → 0. The relevant r-dependence of the spin-gap values for varying values of L y is sketched in Figure 5. Remarkably, while inspecting Figure 5 we deduce that when the amount of disorder among the rungs gets more robust, the spin-gap values tend to decrease, leading to an increment in the slopes of the relevant lines as shown in Figure 4.The most pronounced effect is found for the case L y = 2 where the spin gap shows an almost linear dependence in r.Additionally, our simulation results suggest that upon increasing L y the impact of the disorder on the spin gap becomes less significant and, for the case with L y = 6, almost irrelevant; see the red curve in Figure 5 where ∆ appears as almost constant.As it is well-known from previous studies [13,55], for a larger number of legs the drop in the uniform susceptibility sets in at smaller temperature values and is steeper.Hence, the spin gap is expected to decrease with increasing L y and the effect of randomness among the rungs of the system on the spin gap becomes trivial.We stress that our numerical estimates of ∆ for the clean (r = 0) cases are in good agreement also with previous works [39,56].
Along with the specific-heat and uniform susceptibility curves that show crossing points as detailed above, we have also considered in the last part of this section the staggered susceptibility (χ s (π, π)) and structure factor (S(π, π)) curves for the same set of system parameters.Our numerical results (not shown here for brevity) indicate that: (i) both quantities do not exhibit any instance of joining in the temperature plane for all studied values of the disorder-strength ratio r, and (ii) the structure-factor curves peak at a temperature well below that of the corresponding spin gap, an observation which is in full agreement with earlier studies [34,56], giving further credit to the results presented in this work.In this framework, we consider it more beneficial to present in Figure 6 the r-dependence of the low-temperature behaviour of χ s (π, π), panel (a), and S(π, π), panel (b), for L y = 2, 4, and 6.Note that the numerical data were collected at the minimum temperature value reached during our simulations, i.e., (T /J) = 5 × 10 −2 .The main aftermath from Figure 6 is that while the χ s (π, π) and S(π, π) curves are almost flat for a wide range of r values for the case L y = 2, they tend to decrease with increasing r for L y = 4, and this trend becomes even more prominent for L y = 6.Thus, in contrast to the case of the the spin-gap, see Figure 5, here the impact of disorder becomes much more striking with increasing L y .To clarify this point it is possible to say that due to the presence of a spin gap in the system, the static structure factor is expected to remain finite, as (T /J) → 0. With increasing L y for a chosen r value, the spin gap tends to decrease, and as a result of this, S(π, π) increases.Analogous conclusions have been reported for the case of spin-1/2 Heisenberg ladders yet under the presence of a different type of disorder [56].

Summary
To summarize, we investigated the effects of quenched bond randomness generated along the rungs of the antiferromagnetic spin-1/2 Heisenberg ladders using the stochastic series expansion quantum Monte Carlo technique.Numerical data were obtained for a wide spectrum of the disorder parameter r, from the clean case (r = 0) up to the bond-dilution limit (r → 1) and for three values of L y of the system.Our simulations showed that while the specific-heat maxima decrease with increasing r values for all studied cases of L y , the uniform susceptibility curves appear to be almost identical, typical of the expected Curie behaviour.Only in the low-temperature regime, an increase in the χ u values with increasing disorder was observed.Remarkably, for both C and χ u r-dependent curves the footprint of some characteristic crossing points in the temperature plane was unveiled, whose values are independent of r, however still moderately depend on L y .Notably, the numerical values of the specific heat at these special points may be considered as (nearly) universal for the spin-ladder system under study.An analogous occurrence cannot be excluded also for the uniform susceptibility, although there appears to be a slim dependence of the χ u values exactly at the crossing points on L y , possibly attributed to finite-size effects.Another important result emerging in our study refers to the spin gap, whose presence was revealed from the exponentially decreasing current of the uniform susceptibility curves at the low-temperature regime.From our analysis we deduce that the spin gap decreases with increasing disorder parameter, analogously to the case of decreasing rung coupling values in the clean system, and that upon increasing L y the impact of disorder becomes less significant.Finally, we also studied the low-temperature behaviour of the staggered susceptibility and structure factor.In contrast to the spin-gap case, the impact of disorder in this case becomes more salient with increasing L y .The results of this paper can be extended in several directions: -In terms of model systems, another fruitful platform would be the spin-1 quantum Heisenberg ladder, where the presence of the reported crossing points (and possible universality) has not been yet clarified.-In terms of the diffused randomness, it would be beneficial to consider also other types of disorder distributions, such as bond dilution or non-magnetic impurity atoms which could mimic some aspects of quantum magnetic materials.-Finally, the role of an additional external magnetic field on the spin ladder could provide insight to other physical phenomena, not touched upon the current work.

Fig. 2
Fig.2Specific-heat curves vs. temperature for a wide range of the disorder parameter r and three values of L y .The insets are an enlargement of the intersection area (see also the discussion in the main text).

Fig. 3
Fig.3Same as in Figure2, but for the uniform susceptibility curves.

Fig. 4
Fig. 4 Spin-gap fit lines of the form (11) for the clean case and all disorder ratios studied in this work for L y = 2 (a), L y = 4 (b), and L y = 6 (c).

Figure 4 Fig. 5
Figure 4 now displays the rescaled uniform susceptibility of the system, −T ln (T 1/2 χ u ), for all studied

Fig. 6
Fig.6Disorder-ratio dependence of the low-temperature behavior of the staggered susceptibility (a) and structure factor (b) for all L y values considered.Lines are a simple guide to the eye, while error bars appear smaller than the symbol sizes used.