Ground-State Phase Diagram of an Anisotropic S=1/2 Ladder with Different Leg Interactions

We explore the ground-state phase diagram of the $S=1/2$ two-leg ladder with different leg interactions. The $xy$ and $z$ components of the leg interactions between nearest-neighbor spins in the $a$ ($b$) leg are respectively denoted by $J_{{\rm l},a}$ and $\Delta_{\rm l} J_{{\rm l},a}$ ($J_{{\rm l},b}$ and $\Delta_{\rm l} J_{{\rm l},b}$). On the other hand, the $xy$ and $z$ components of the uniform rung interactions are respectively denoted by $\Gamma_{\rm r} J_{{\rm r}}$ and $J_{{\rm r}}$. In the above, $\Delta_{\rm l}$ and $\Gamma_{\rm r}$ are the $XXZ$-type anisotropy parameters for the leg and rung interactions, respectively. This system has a frustration when $J_{{\rm l},a} J_{{\rm l},b}<0$ irrespective of the sign of $J_{\rm r}$. The phase diagram on the $\Delta_{\rm l}$ ($|\Delta_{\rm l}| \leq 1.0$) versus $J_{{\rm l},b}$ ($-2.0\leq J_{{\rm l},b}\leq 3.0$) plane in the case where $J_{{\rm l},a}=0.2$, $J_{{\rm r}}=-1.0$, and $\Gamma_{\rm r} = 0.5$ is determined numerically. We employ the physical consideration, and the level spectroscopy and phenomenological renormalization-group analyses of the numerical date obtained by the exact diagonalization method. The resultant phase diagram contains the ferromagnetic, Haldane, N{\'e}el, nematic Tomonaga-Luttinger liquid (TLL), partial ferrimagnetic, and $XY1$ phases. Interestingly enough, the nematic TLL phase appears in the strong-rung unfrustrated region as well as in the strong-rung frustrated one. We perform the first-order perturbational calculations from the strong rung coupling limit to elucidate the characteristic features of the phase diagram. Furthermore, we make the density-matrix renormalization-group calculations for some physical quantities such as the energy gaps, the local magnetization, and the spin correlation functions to supplement the reliability of the phase diagram.


Introduction
Over the past several decades, a great deal of effort has been devoted to clarifying the frustration effect on the ground-state properties of low-dimensional quantum spin systems with competing interactions. It is now well known experimentally as well as theoretically and numerically that a reciprocal influence between the strong quantum fluctuation and the geometrical frustration in these systems induces various exotic quantum ground states such as the dimer state accompanying spontaneous translational symmetry breaking, the spinnematic Tomonaga-Luttinger liquid (TLL) state, and so on. The former state is realized in an S = 1/2 zigzag chain with antiferromagnetic nearest-neighbor (nn) and nextnearest-neighbor (nnn) interactions. [1][2][3][4][5][6] On the other hand, the latter state is realized in an anisotropic S = 1 chain, [7][8][9] in an S = 1 chain with bilinear and biquadratic interactions, 10) and in an S = 1/2 zigzag chain with ferromagnetic nn and antiferromagnetic nnn interactions under an external magnetic field. [11][12][13] In the case of an S = 1/2 two-leg ladder system, * E-mail address: tone0115@vivid.ocn.ne.jp the frustration is usually introduced by adding nnn leg and/or diagonal interactions to an original system with nn leg and rung interactions, as has been extensively investigated. [14][15][16] As another example, an S = 1/2 ladder with uniform nn leg and alternating rung interactions has also been studied by several authors. [17][18][19] We 19) have discussed the ground-state phase diagram in the frustrated case where rung interactions are ferromagnetically-antiferromagnetically alternating and have a common Ising-type anisotropy, while leg interactions are antiferromagnetically uniform and isotropic. Our results show that, when the leg interactions are relatively weak compared with the rung interactions, the incommensurate Haldane state as well as the commensurate one appears as the ground state in the whole range of the Ising-type anisotropy parameter. This appearance of the Haldane state in the case where the Ising character of rung interactions is strong is contrary to the ordinary situation, and is called the inversion phenomenon concerning the interaction anisotropy. [20][21][22][23] In the present paper, we explore the ground-state phase diagram of an anisotropic S = 1/2 ladder with different leg interactions, which is schematically sketched in S x j,a S x j+1,a + S y j,a S y j+1,a + ∆ l S z j,a S z j+1,a Here, S x j,ℓ , S y j,ℓ , and S z j,ℓ are, respectively, the x-, y-, and z-components of the S = 1/2 operator S j,ℓ at the (j, ℓ) site assigned by the jth rung and the ℓ(= a or b) leg; J l,a and J l,b denote, respectively, the magnitudes of the a leg and b leg interactions, while J r denotes that of the rung interaction; ∆ l and Γ r are, respectively, the parameters representing the XXZ-type anisotropies of the former and latter interactions; L is the total number of rungs, which is assumed to be even. It is emphasized that this system has frustration when J l,a J l,b < 0 irrespective of the sign of J r .
It is unfortunate that real materials for which the above H is a good model Hamiltonian have been neither synthesized nor found so far. We do not think, however, that it is physically unrealistic. Actually, Yamaguchi et al. 24,25) have recently demonstrated the modulation of magnetic interactions in spin ladder systems by using verdazyl-radical crystals. The flexibility of molecular arrangements in such organic-radical materials is expected to realize S = 1/2 ladder systems with different leg interactions.
Our Hamiltonian H contains five interaction parameters. Throughout the following discussions, we focus our attention upon the case where J l,a = 0.2, −2.0 ≤ J l,b ≤ 3.0, J r = −1.0, |∆ l | ≤ 1.0, and Γ r = 0.5, unless otherwise stated. Here, we choose |J r | as the unit of energy. It is noted that the anisotropies of the leg and rung interactions are, respectively, of the XY and Ising types. The motivation for treating this case is as follows. When the ferromagnetic rung interactions with Ising-type anisotropy are much stronger than both kinds of leg interactions (|J l,b | ≪ 1.0), a pair of S = 1/2 spins at each rung forms a bound state of two magnons with S z j,a +S z j,b = ±1. This may lead to the nematic TLL state, which accompanies two-magnon bound states, as the ground state at least in the frustrated region. Furthermore, the XY -type anisotropy of the leg interactions is expected to stabilize the nematic TLL state.
In this paper, we use the term "nematic TLL" for the TLL characterized not only by the formation of twomagnon bound pairs but also by the dominant nematic four-spin correlation function. In the nematic TLL state composed of the bosons of two-magnon bound pairs, 7,12) the nematic four-spin and longitudinal two-spin correlations are dual to each other, and the decay exponent of the former, η ++−− , and that of the latter, η zz ℓ , obey the relation η ++−− η zz ℓ = 1. [See Eqs. (38) and (37) for the definitions of η ++−− and η zz ℓ , respectively.] If η ++−− < 1, the nematic correlation is dominant, leading to the nematic TLL state defined above. On the other hand, the TLL state with the dominant longitudinal two-spin correlation is not allowed in the present system, since, if η zz ℓ < 1, the Umklapp scattering becomes relevant and the TLL state turns into the Néel state. The situation is distinct from that of zigzag chains under a magnetic field, in which the TLL state accompanied with two-magnon bound pairs and dominant incommensurate spin-densitywave correlation appears in wide regions of the magnetic phase diagrams. 12,13,26,27) We note, at this junction, that ground states of the ladder system governed by the Hamiltonian H have already been discussed in a few other cases. [28][29][30] We 29) have determined the ground-state phase diagram on the 1/Γ r (Γ r ≥ 1.0) versus J l,b plane in the case where J l,a = ±0.2, J r Γ r = −1.0, and ∆ l = 1.0. The obtained phase diagrams consist of the XY 1 phase, the triplet-dimer phase, the partial ferrimagnetic phase (which is called the noncollinear ferrimagnetic phase in the paper), and the Haldane phase. Interestingly, the direct-product tripletdimer state becomes to be the exact ground state 28,31) when J l,b = −J l,a and Γ r > ∼ 1.2. On the other hand, extending Tsukano and Takahashi's work, 28) Sekiguchi and Hida 30) have shown, in the case where J l,a J l,b < 0 and ∆ l = Γ r = 1.0, that the partial ferrimagnetic state, which is a spontaneously magnetized TLL state with incommensurate magnetic correlation, appears as the ground state over a wide parameter range.
The remainder of this paper is organized as follows. In Sect. 2, with the help of physical considerations we numerically determine, by using the exact diagonalization (ED) method, the ground-state phase diagram on the ∆ l versus J l,b plane in the case where J l,a = 0.2, J r = −1.0, and Γ r = 0.5. In Sect. 3, we show that the distinguishing features of the obtained phase diagram are well explained by performing perturbational calculations. In Sect. 4, in order to supplement the reliability of the phase diagram, we carry out density-matrix renormalization-group (DMRG) calculations 32,33) for some physical quantities such as the energy gaps, the local magnetization, and the spin correlation functions. The final section (Sect. 5) contains concluding remarks. In Appendix, the groundstate phase diagram on the Γ r versus J l,b plane in the case where J l,a = 0.2, J r = −1.0, and ∆ l = 1.0 is briefly discussed.

Numerical Determination of the Ground-State Phase Diagram
In this section we present the result for the groundstate phase diagram on the ∆ l versus J l,b plane in the case where J l,a = 0.2, J r = −1.0, and Γ r = 0.5, which has been determined by using a variety of numerical methods based on the ED calculation. First of all, we show the obtained phase diagram in Fig. 2. This phase diagram consists of six kinds of phases; these are the ferromagnetic, partial ferrimagnetic, XY 1, Haldane, Néel, and nematic TLL phases. It is noted that the nematic TLL phase appears as the ground state in the strongrung frustrated region (−1.0 < ∼ J l,b < 0.0), as is expected, and survives even in the strong-rung unfrustrated region (0.0 ≤ J l,b < ∼ 0.4).
Before explaining how to numerically determine the phase boundary lines, we introduce some related physical quantities. We denote, respectively, by E P 0 (L, M ) and E P 1 (L, M ), the lowest-and second-lowest-energy eigenvalues of the Hamiltonian H under the periodic boundary condition (PBC), S L+1,ℓ = S 1,ℓ , within the subspace characterized by L and the total magnetization  = −0.280 ± 0.002 at the L → ∞ limit, where the numerical error is estimated from the difference between the above two extrapolated results.
Here, M is a good quantum number with the eigenvalues of M = 0, ±1, · · · , ±L. Furthermore, the quantity E T 0 (L, M ) is the lowest-energy eigenvalue of H under the twisted boundary condition (TBC), S x L+1,ℓ = −S x 1,ℓ , S y L+1,ℓ = −S y 1,ℓ , and S z L+1,ℓ = S z 1,ℓ , within the subspace determined by L and M . Then, we define several energy differences as follows: It is noted that the values of the ground-state magnetization M g (L) are M g (L) = L, 0 < M g (L) < L, and M g (L) = 0 in the ferromagnetic, partial ferrimagnetic, and remaining four phases, respectively. Let us first discuss the phase boundary lines between the Néel and nematic TLL phases, which are depicted by the brown lines in Fig. 2. It should be emphasized that the nematic TLL state is the TLL state which accompanies two-magnon bound states and is equivalent to the XY 2 state, originally proposed by Schulz. 7) The corresponding phase transition is of the Berezinskii-Kosterlitz-Thouless (BKT) type 34,35) accompanying the spontaneous translational-symmetry breaking (STSB), and therefore, as is well known, the phase boundary line can be accurately estimated by the level spectroscopy (LS) method developed by Okamoto and Nomura. 3,5,6) In this method, the difficulty coming from the logarithmic-correction problem associated with the BKT transition is removed. According to this method, the critical value J of J l,b for a given value of ∆ l can be evaluated as follows. We numerically solve the equation to compute the finite-size critical value J  various values of L. Then, we extrapolate these finitesize results to the thermodynamic (L → ∞) limit. Practically, we have carried out the ED calculation to compute J 's as functions of ∆ l .
In the following discussion in this section, we mainly describe how to compute the finite-size critical values only, since the evaluation of the critical value at the thermodynamic limit has been done in the same way.
Secondly (see the green line in Fig. 2), we discuss the phase boundary line between the XY 1 and nematic TLL phases. The nematic TLL state accompanies two-magnon bound states, as mentioned above, while the XY 1 state does not. These lead to the following fact. In the groundstate magnetization curve for a given finite-size system under an external magnetic field H applied along the z-axis, the magnetization increases stepwisely with increasing | H|; the first step occurs from the M g (L) = 0 state to the M g (L) = 2 state in the nematic TLL phase, while it occurs from the M g (L) = 0 state to the M g (L) = 1 state in the XY 1 phase. From these arguments, it is easy to see that the finite-size critical value J , for a given value of ∆ l , for the phase transition between the XY 1 and nematic TLL phases is calculated by solving numerically. Figure 4 shows the determination of the critical value J (c;XY 1,nTLL) l,b at the thermodynamic limit for ∆ l = −0.2.
Thirdly (the blue line in Fig. 2), the phase transition between the XY 1 and Haldane phases is the BKT transition without the STSB. It is also well known that, in this case, Nomura and Kitazawa's LS method, 36) in which the difficulty coming from the logarithmic-correction problem is also removed, is very powerful for determining the phase boundary line. The finite-size critical value ∆ (c;XY 1,H) l (L) of ∆ l for a given value of J l,b is calculated by using the equation In Fig. 5 we show the determination of the critical value ∆ (c;XY 1,H) l at the thermodynamic limit for J l,b = 2.0. Fourthly (the umber lines in Fig. 2), the phase transitions between the ferromagnetic phase and one of the XY 1, Néel, and nematic TLL phases are of the firstorder, because these transitions are between the M g (L) = L and M g (L) = 0 states at the L → ∞ limit. Thus, it is apparent that the phase boundary lines are obtained from the equation As an example, we discuss here the phase transition between the ferromagnetic and nematic TLL phases in the case of J l,b = −0.7, denoting the finite-size critical value by ∆ (c;Fro,nTLL) l (L). In Fig. 6 the determination of the critical value ∆ (c;Fro,nTLL) l at the thermodynamic limit in this case is illustrated.
Fifthly (the magenta, red, and black lines in Fig. 2), according to the result of our calculations, the partial ferrimagnetic phase appears in a rather narrow region of ∆ l , which lies between the ferromagnetic-phase region and the XY 1-phase region, only when J l,b > ∼ 1.2. It shows further that the phase transition between the ferromagnetic and partial ferrimagnetic phases is between the M g (L) = L and M g (L) = L − 1 states, while that between the partial ferrimagnetic and XY 1 phases is between the M g (L) = 1 and M g (L) = 0 states or between the L − 1 ≥ M g (L) ≥ 2 and M g (L) = 0 states depending upon whether J l,b > ∼ 1.55 or 1.55 > ∼ J l,b > ∼ 1.2. Thus, it may be considered that the former transition is always of the second-order, while the latter transition is of the second-order when J l,b > ∼ 1.55 and of the first-order when , respectively. The phase boundary line between the ferromagnetic and partial ferrimagnetic phases can be analytically calculated by examining the instability of the one-magnon excitation from the ferromagnetic state. However, since this calculation is straightforward but too tedious, we show in Fig. 2 the numerical result obtained by solving It is noted that the obtained finite-size critical value is independent of L (the magenta line). This means that in the present case, the M g (L) = L − 1 state with the lowest energy is a commensurate state. The finite-size critical value ∆ (c;PFri,XY 1) l (L) of ∆ l , for a given value of J l,b , for the second-order phase transition between the partial ferrimagnetic and XY 1 phases is calculated by solving Figure 8 illustrates the determination of the corresponding critical value ∆ (c;PFri,XY 1) l at the thermodynamic limit in the case where J l,b = 2.0 (the red line). Generally speaking, it is difficult to accurately determine the critical value ∆ (c;PFri,XY 1) l of ∆ l at the thermodynamic limit, for a given value of J l,b , for the first-order phase transition between the partial ferrimagnetic and XY 1 phases by performing only the ED calculations for finite-L systems with up to L = 14 rungs. Fortunately, however, in the present case the corresponding finite-size critical value ∆ (c;PFri,XY 1) l (L) of ∆ l is almost independent of L, as can been seen from Fig. 7(b), and therefore we use ∆ (c;PFri,XY 1) l (14) for ∆ (c;PFri,XY 1) l (the black line). We consider that the numerical error due to this simplification is less than 0.001.
Lastly (the cyan line in Fig. 2), the phase transition between the Haldane and Néel phases is the 2D Ising-type transition. In this case, the critical value ∆ (c;H,N) l of ∆ l at the thermodynamic limit for a given value of J l,b can be evaluated by the phenomenological renormalizationgroup (PRG) method 37) in the following way. We cal- Then, we extrapolate the finite-size results to the L → ∞ limit to evaluate ∆ (c;H,N) l . In a practical extrapolation we fitted them to a quadratic function of (L + 1) −2 by using the least-squares method. In Fig. 9 we illustrate this procedure for the case of J l,b = 2.0.
Before closing this section, we make two more comments on the phase diagram shown in Fig. 2. First, there are two tetracritical points, at one of which the ferromagnetic-nematic TLL and Néel-nematic TLL phase boundary lines come in contact with each other, and at the other the XY 1-Haldane, Haldane-Néel, Néelnematic TLL, and nematic TLL-XY 1 phase boundary lines merge. As regards the latter tetracritical point, we have concluded that this is the case based on the physical consideration by referring to the result of the perturbational calculation discussed in Sect. 3 [see Fig. 10(b)]. We note that it is hard to draw this conclusion only from the present ED calculation for finite-L systems with up to L = 14 rungs. Secondly, the results of our ED calculation suggest that, when the value of J l,b is sufficiently small (at least when J l,b < ∼ −4.0), the XY 1 phase may appear below the Néel phase. The phase transition between these two phases is the BKT transition with the STSB, and the phase boundary line can be estimated from the equation 3,5,6) ∆E P 1 (L, 0) = ∆E P 0 (L, 1) instead of Eq. (5). We have tried to estimate this phase boundary line in a similar way. Although the finite-size critical points which satisfy Eq. (12) exist, the numerical error coming from the extrapolation L → ∞ is fairly large in this case. Thus, we believe that the XY 1 phase exists in this region, but it is difficult to draw the phase boundary line. In order to overcome this difficulty, it is indispensable to attempt more sophisticated analyses, which are beyond the scope of the present study.

Perturbational Calculations
In order to elucidate the characteristic features of the phase diagram shown in Fig. 2, we have performed the first-order perturbational calculation from the strongrung coupling limit. We begin with the two-spin problem of the jth rung; The eigenstates and the eigenenergies of H r are as follows: where | ↑ j,ℓ and | ↓ j,ℓ denote, respectively, the S z j,ℓ = 1 2 and S z j,ℓ = − 1 2 states. We take first the three states |ψ 1 , |ψ 2 , and |ψ 3 , and neglect the |ψ 4 state, because we suppose that J r (= −1) is ferromagnetic. In this restricted space, we interpret these three states as being the T z j = +1, T z j = 0, and T z j = −1 states, respectively, of the pseudospin operator T j with T j = 1. Comparing the matrix elements, we see that the original spin operators S µ j,ℓ (µ = x,y,z) can be expressed by within the restricted space. Thus, in the first-order perturbation theory, we can obtain the effective Hamiltonian where the D eff term comes from the energy differences E 1 −E 2 and E 3 −E 2 . The effective Hamiltonian (19) is a fundamental model for the anisotropic T = 1 chain, which has been investigated by several authors. [7][8][9] In particular, Chen et al. 9) have numerically determined the phase diagram of H eff on the ∆ eff versus D eff plane. Figure 10(a) schematically shows the above-mentioned phase diagram of the effective Hamiltonian H eff obtained by Chen et al. 9) We have chosen that J l,a = 0.2, J r = −1.0, and Γ r = 0.5, and therefore the parameters in H eff are given by From Fig. 10(a), H eff , and the above three equations, we obtain the phase diagram given in Fig. 10(b). This phase diagram has a point symmetry with respect to the point (∆ l , J l,b ) = (0.0, −0.2), which can be explained as follows.
If we make the unitary transformations, the effective Hamiltonian H eff is transformed into Thus, we see that the two points, (J eff , ∆ eff , D eff ) and (−J eff , −∆ eff , −D eff ), belong to the same phase. The transformation can be carried out by which explains the point symmetry of the phase diagram shown in Fig. 10(b). The above arguments concerning the point symmetry well explain the reason why the ferromagnetic, Néel, and nematic TLL phases appear in two places in Fig. 2(a). They suggest that the XY 1 state appears as the ground state on the lower outside of Fig. 2(a) (J l,b < −2.0). Although we believe the existence of the XY 1 phase, we are unfortunately not able to determine its phase boundary, as already noted in the last paragraph in Sect. 2. Also, the Haldane phase does not appear in the region of negative J l,b in Fig. 2(a), and further the partial ferrimagnetic phase found in Fig. 2(a) does not appear in Fig. 10(a). We think that the appearance of the partial ferrimagnetic phase is attributed to the strong frustration effect, which is not taken into account in the present first-order perturbational calculation.
Let us explain the correspondence between Fig. 2(a) and Fig. 10(a) a little bit more in detail. When J l,b changes, on the vertical line ∆ l = 1.0 in Fig. 2(a), from a large positive value to −0.2, the parameter J eff is always positive, ∆ eff = 1.0, and D eff changes from a small negative value to a large negative value. Thus, the corresponding point in Fig. 10(a) moves along the downwardarrowed blue solid line. On the other hand, when J l,b changes, again on the vertical line ∆ l = 1.0 in Fig. 2(a), from −0.2 to a large negative value, the parameter J eff is always negative, ∆ eff = 1.0, and D eff changes from a large positive value to a small positive value. Since the points (J eff , ∆ eff , D eff ) and (−J eff , −∆ eff , −D eff ) belong to the same phase as already explained, these changes in the parameters in H eff are equivalent to those in the case where J eff is positive, ∆ eff = −1.0, and D eff changes from a large negative value to a small negative value. Namely, the corresponding point in Fig. 10(a) moves along the upward-arrowed blue solid line. Although the two blue solid lines are separated from each other, we surmise that they are connected with each other by the blue dotted line in the higher-order perturbation theory at least near ∆ eff = 0.0. Another possibility is that they are connected through ∆ eff = ±∞. Thus, the phase of the original system H successively changes to the Néel, nematic TLL, and ferromagnetic phases with decreasing J l,b , which well explains the phases on the ∆ l = 1.0 line in Fig. 2(a). We can make a very similar discussion, for example, for the line ∆ l = −0.5 in Fig. 2(a) by use of the red dot-dashed lines in Fig. 10(a). Namely, the phase successively changes to the XY 1, ferromagnetic, nematic TLL, and Néel phases with decreasing J l,b in this case. Thus, the perturbational calculation together with the numerical result obtained by Chen et al. 9) qualitatively explains the numerically determined phase diagram of the present system (H).
Here, we discuss the validity of the perturbational calculation. At a glance, this perturbational calculation is thought to be valid for the strong-rung coupling case, |J l,b | ≪ |J r |. However, the bosonization theory from the weak-rung coupling limit shows that the rung coupling is relevant and is renormalized to the strong-rung coupling case, as is known in the simple ladder case. This fact strongly suggests that the present perturbational calculation is qualitatively valid even for the weak-rung coupling case.

DMRG Calculations
In this section, we present several results of the DMRG calculation, which supplement the reliability of the phase diagram shown in Fig. 2. In this calculation, we assume the open boundary condition (OBC) for the Hamiltonian H, where the sum over j is taken from j = 1 to j = L−1. We treat the finite-size systems with up to L = 96 rungs, L being assumed to be even. The number of block states in the DMRG calculation required to achieve the desired accuracy depends on the state targeted. In our calculation, we have kept up to 600 block states for the most severe case and checked that the results are accurate enough for our arguments below by monitoring the dependences of the data and the truncation error on the number of kept states. The where depending upon whether j(> 0) is odd or even, j 0 = L+1−j 2 or j 0 = L−j 2 , respectively, which means that the two rungs j 0 and j 0 + j are at (almost) equal distances from the center of the open ladder. In the above three equations, · · · L,M denotes the expectation value associated with the lowest-energy state within the subspace determined by L and M in the system under the OBC; when M = M g (L), this expectation value is simply the ground-state expectation value for the system with L rungs.
where α 1,ℓ , α 2,ℓ , and α are numerical constants. The decay exponents obey the relation η zz ℓ η ++−− = 1. We have performed the least-squares fitting of the nematic fourspin correlation function ω ++−−,O (j; 96, 0) to Eq. (38) using the data for 6 ≤ j ≤ 36. We have thereby obtained the estimate of the decay exponent, η ++−− ∼ 0.577, which is much smaller than 1 and supports the prediction of the dominant nematic correlation. For the longitudinal two-spin correlations, we have unfortunately not been able to achieve a precise estimation of the decay exponent due to rather large open-boundary effects. We have fitted ω zz,O ℓ (j; 96, 0) to Eq. (37) using the data for several ranges of j, namely, 3 ≤ j ≤ 20, 4 ≤ j ≤ 24, 5 ≤ j ≤ 30, and 6 ≤ j ≤ 36. Then, we have obtained the estimates   in the ranges of 1.70 ≤ η zz a ≤ 2.11 and 1.63 ≤ η zz b ≤ 2.31. Although these estimates show rather large variations, the values are consistent with the expectation η zz ℓ = 1/η ++−− , namely, η zz ℓ ∼ 1/0.577 = 1.73, and indicate that the longitudinal spin correlation decays much faster than the nematic correlation. The DMRG data thereby support the result of the ED analysis that the system at the point (∆ l = 0.5, J l,b = −0.5) belongs to the nematic TLL phase. It is also emphasized that η zz a and η zz b are almost equal to each other in spite of the fact that J l,a is positive, while J l,b is negative. This result reflects the fact that the nematic TLL is a one-component TLL, 11,12) which leads to the same decay exponents of the correlation functions within the a leg and within the b leg.
where α ′ ℓ and α ′ are constants. It is also expected that the longitudinal and transverse spin correlations are dual to each other, resulting in the relation η xx ℓ η zz ℓ = 1, while the exponent of the nematic correlation function is related to that of the transverse-spin correlation function as η ++−− = 4η xx ℓ . We have estimated η xx ℓ and η ++−− by the least-squares fitting of ω xx,O ℓ (j; 96, 0) and ω ++−−,O (j; 96, 0) to Eqs. (39) and (40), respectively, using the numerical data for 6 ≤ j ≤ 36. The results are η xx a ∼ 0.15, η xx b ∼ 0.16, and η ++−− ∼ 0.61, which support the expectation that η ++−− = 4η xx ℓ . For the longitudinal-spin correlation function, we have found that ω zz,O ℓ (j; 96, 0) is negative for all j and that the oscillating component is far from sizable, indicating that the second (uniform) term in Eq. (37) is much more dominant than the first (oscillating) term. This observation, with our estimate η xx ℓ ∼ 0.15 -0.16 above, is found to be consistent with the expectation that η xx ℓ η zz ℓ = 1, since η zz ℓ = 1/η xx ℓ ∼ 6.3 -6.7 is much larger than 2, the exponent of the uniform component. Figure 15 shows  Fig. 2. From this figure, we see that the strong Néel long-range order (LRO) certainly exists, supporting our conclusion obtained by the ED analysis that the system at this parameter point indeed belongs to the Néel phase. It is noted that we have also obtained a similar result at the point (∆ l = 0.75, J l,b = 1.0) which belongs to another Néel-phase region.

Concluding Remarks
We have numerically determined the ground-state phase diagram on the ∆ l versus J l,b plane of the S = 1/2 two-leg ladder with different leg interactions, which is governed by the Hamiltonian H of Eq. (1), confining ourselves to the case where J l,a = 0.2, J r = −1.0, and Γ r = 0.5. The resultant phase diagram, which is presented in Fig. 2, consists of the ferromagnetic, Haldane, Néel, nematic TLL, partial ferrimagnetic, and XY 1 phases. In addition, the phase diagram on the Γ r versus J l,b plane in the case where J l,a = 0.2, J r = −1.0, and ∆ l = 1.0, in which the nematic TLL phase also appears, has been obtained (see Fig. A·1 later). By estimating the phase boundary lines, we have successfully applied the LS 3, 5, 6, 36) or PRG 37) analysis of the numerical data obtained by the ED method to the transition between two of the Haldane, Néel, nematic TLL, and XY 1 phases, while we have compared the energies E P 0 (L, M )'s for the transition associated with the ferromagnetic and partial ferrimagnetic phases. In our opinion, one of the significant results of the present work is the fact that the nematic TLL phase appears in the strong-rung unfrustrated region as well as in the strong-rung frustrated region.
We have carried out first-order perturbational calculations from the strong-rung coupling limit, in which the present system (H) is mapped onto the anisotropic S = 1 chain system (H eff ), in order to outline the essential features of the phase diagram. Furthermore, performing the DMRG 32, 33) calculations, we have clarified the following. In the Haldane phase, there is a finite excitation gap (the Haldane gap) above the ground state, and the edge-localized spins appear. In the Néel phase, the magnetic LRO of the Néel type is found. In the nematic TLL phase, the two-magnon bound state characterizes this phase, and the nematic four-spin correlation function with a quasi-LRO is more dominant than the longitudinal two-spin correlation function with a quasi-LRO and the transverse two-spin correlation function with a short-range order. In the XY 1 phase, all of the three correlation functions exhibit a quasi-LRO, and among them the transverse two-spin correlation function is most dominant.
The recent result of the bosonization analysis for the present system, attempted by one of the authors (S. C. F.) suggests that the nematic TLL phase also appears in the weak-rung frustrated region. This prediction is highly plausible for the following reason. Let us consider a two-magnon excited state in the system of almost independent chains with (J l,a , ∆ l ) and (J l,b , ∆ l ) coupled via the weak rung interaction (J r , Γ r ). We assume that J l,a > 0.0, J l,b < 0.0, J r < 0.0, |∆ l | < 1.0, and 0.0 ≤ Γ r < 1.0. The ferromagnetic Ising-like rung interaction favors the formation of the two-magnon bound state in a rung, as mentioned in Sect. 1. On the other hand, if the onemagnon excitation energies in each chain are different, the leg interactions result in two magnons sitting in the chain with the smaller excitation energy, resulting in a scattered state of two magnons. Therefore, at the level of a crude approximation, one can expect that the two magnons form a bound state when |J r | is sufficiently larger than the difference between the one-magnon excitation energies in each chain, where (∆E) om (L; J ℓ , ∆ l ) is the one-magnon excitation energy in an independent chain with J l,ℓ and ∆ l . The L → ∞ asymptotic form of the one-magnon excitation energy is obtained from the conformal field theory 41) and the Bethe ansatz 42, 43) as Note that the chain with (J l,b , ∆ l ) is transformed into the chain with (|J l,b |, −∆ l ) by the π-rotation around the z axis for every other spin. Thus, solving the equation |∆(∆E) om | = 0 with Eqs. (42)-(45), we can determine the value∆ l of ∆ l at which the difference in the one-magnon excitation energies vanishes. The result is ∆ l = cos π J l,a /|J l,b | 1 + J l,a /|J l,b | , from which, when J l,a = 0.2, for example,∆ l = 0.0 for J l,b = −0.2 and∆ l = 0.158 for J l,b = −0.3. We thereby expect that, at least for a parameter region around the line ∆ l =∆ l , the inclusion of weak ferromagnetic, Ising-like rung interactions of J r leads to the formation of a twomagnon bound state. Here, we should mention that the above argument is not sufficient for one-dimensional systems with strong quantum fluctuations, and more sophisticated approaches such as the weak-coupling renormalization group method are required. We are now investigating the ground-state phase diagram in this weak-rung region, and the results will be reported elsewhere. Finally, we hope that the present work will stimulate future experimental studies on related subjects, which include the synthesis of a ladder system with different leg interactions. and P18H04330 (J-Physics). We also thank the Supercomputer Center, Institute for Solid State Physics, University of Tokyo and the Computer Room, Yukawa Institute for Theoretical Physics, Kyoto University for computational facilities.
Appendix Figure A·1 demonstrates the obtained ground-state phase diagram on the Γ r versus J l,b plane in the case where J l,a = 0.2, J r = −1.0, and ∆ l = 1.0. This consists of the Néel, Haldane, nematic TLL, XY 1, ferromagnetic, and partial ferrimagnetic phases, all of which are the same as those in Fig. 2. In this phase diagram, the phase boundary line between the Néel and XY 1 phases is drawn by the blue line, and the colors of the other lines have the same meanings as those in Fig. 2.
All of the phase boundary lines, except for that for the first-order phase transitions between the partial ferrimagnetic phase and one of the nematic TLL and XY 1 phases, have been numerically determined in a way similar to those discussed in Sect. 2; in this case, Okamoto and Nomura's LS method, 3,5,6) where we use Eq. (12), works very well for estimating the phase boundary line between the Néel and XY 1 phases. Obtaining the phase boundary line for the above first-order phase transitions, we have calculated M g (L) by carrying out DMRG calculations for the finite system with L = 36 rungs under the OBC. We consider that the result for the L = 36 system gives a fairly good approximation of the result at the L → ∞ limit.
The perturbation theory developed in Sect. 3 predicts the phase diagram shown in Fig. A·2 in the following way. Here, the parameter ∆ eff is given by ∆ eff = 1.0, since ∆ l = 1.0. Thus, the J l,b > −0.2 and J l,b < −0.2 regions in Fig. A·1 correspond, respectively, to the downwardand upward-arrowed blue solid lines in Fig. 10(a) Fig. 10(a) leads to the horizontal nematic TLL band in Fig. A·2. The predicted phase diagram given in Fig. A·2 qualitatively explains the numerical one given in Fig. A·1, although the XY 1 and partial ferrimagnetic phases do not appear in Fig. A·2. The appearance of these phases in Fig. A·1 is due to the frustration effect and is beyond the description of the first-order perturbation theory. We note that the existence of the XY 1 phase is considered to be an example of the inversion phenomenon concerning the interaction anisotropy 20-23) stated in Sect. 1.