Magnetization Process of the Spin-1/2 Triangular-Lattice Heisenberg Antiferromagnet with Next-Nearest-Neighbor Interactions -- Plateau or Nonplateau

An $S=1/2$ triangular-lattice Heisenberg antiferromagnet with next-nearest-neighbor interactions is investigated under a magnetic field by the numerical-diagonalization method. It is known that, in both cases of weak and strong next-nearest-neighbor interactions, this system reveals a magnetization plateau at one-third of the saturated magnetization. We examine the stability of this magnetization plateau when the amplitude of next-nearest-neighbor interactions is varied. We find that a nonplateau region appears between the plateau phases in the cases of weak and strong next-nearest-neighbor interactions.


Introduction
Frustration has attracted the attention of many researchers in condensed-matter physics because frustration becomes a source of nontrivial quantum phenomena in various systems of physics. In magnetic materials, frustrations occur, for example, when antiferromagnetic interactions of a system form a triangle. The triangular-lattice antiferromagnet is a typical case. Since Anderson 1 pointed out that this model is a possible candidate for the realization of a spin-liquid ground state owing to frustrations in the system, the S = 1/2 Heisenberg model in particular has been studied.  On the basis of extensive studies of the ground state of this system, many physicists believe that the symmetry-breaking state with the so-called 120-degree structure is realized in the ground state. On the other hand, a recent large-scale numerical study 23 suggested the absence of such breaking; this issue remains controversial even now.
One nontrivial phenomenon in the quantum Heisenberg antiferromagnet on the triangular lattice under magnetic fields is the magnetization plateau. This phenomenon appears at one-third of the saturated magnetization in the magnetization curve. It is known that the plateau on the triangular-lattice antiferromagnet appears by the so-called order-by-disorder mechanism at zero temperature even though the corresponding classical system does not show the plateau. 24 In an early stage of investigations of the plateau, the presence of the plateau was just a theoretical prediction. However, experimental realizations have been reported, for examples, for Ba 3 CoSb 2 O 9 in the S = 1/2 case 25 and for Ba 3 NiSb 2 O 9 in the S = 1 case. 26 On the other hand, the destabilization of this plateau was studied from the viewpoint of the effect of randomness in the system. 27,28 Under the above circumstances, we are faced with a question: What else can destabilize the plateau?
A possible candidate is additional interactions at a next-nearest-neighbor (NNN) pair. NNN interactions in the S = 1/2 triangular-lattice Heisenberg antiferromagnet have already been studied; 5,6,9,29,30 however, the studies were carried out for the system without an external magnetic field. The purpose of the present study is to clarify how the plateau behaves in the presence of NNN interactions at zero temperature. In particular, we focus our attention on whether or not a plateau of this height is present during the variation of NNN interactions. Our numerical-diagonalization study provides us with information on a new phase transition driven by NNN interactions. This paper is organized as follows. In the next section, the model studied here is introduced. The method is also explained. The third section is devoted to the presentation and discussion of our results. We first observe magnetization processes for various amplitudes of NNN interactions. In order to understand the behavior observed in magnetization processes, several analyses are carried out. In the final section, we present our conclusion.

Model Hamiltonian and Method
The Hamiltonian studied here is given by and the Zeeman term is given by Here, S i denotes the S = 1/2 spin operator at site i. In this study, we consider the case of an isotropic interaction in spin space. Site i is assumed to be the vertices of a triangular lattice composed of bonds of nearest-neighbor interactions J 1 illustrated by the black lines in Fig. 1  sublattices: A, B, and C. Note here that each sublattice also forms a triangular lattice illustrated by colored lines in Fig. 1. The colored lines represent NNN interactions whose amplitudes are all given by J 2 . We denote the ratio J 2 /J 1 by r. We consider that all the interactions are antiferromagnetic, namely, J 1 > 0 and J 2 > 0. Energies are measured in units of J 1 ; hereafter, we set J 1 = 1 and examine the case of r ≥ 0. Note here that for J 2 = 0, namely, r = 0, the present lattice is identical to the triangular lattice without NNN interactions and that for infinitely large J 2 , namely, r → ∞, the system reduces to three isolated triangular-lattice antiferromagnets. The finite-size clusters that we treat in the present study are depicted in Fig. 2. We examine the cases of N s =27 and 36 under the periodic boundary condition. In order to detect the magnetization plateau at one-third of the saturated magnetization in both limiting cases for a vanishing J 2 and an infinitely large J 2 , it is necessary for N s /9 to be an integer. In order to capture two-dimensionality well, additionally, the cluster shapes are assumed to be rhombic and to have an inner angle of π/3. The rhombic condition is satisfied not only for the triangular lattice of nearest-neighbor interactions but also for any of the triangular lattices of each sublattice.
We calculate the lowest energy of H 0 in the subspace belonging to j S z j = M by numerical diagonalizations based on the Lanczos algorithm and/or Householder algorithm. Our diagonalizations are carried out in the basis where the z-axis is taken as the quantized axis of each spin. The numerical-diagonalization calculations are unbiased; one can therefore obtain reliable information on the system. The energy is denoted by E r (N s , M ), where M takes an integer or a half odd integer up to the saturation value M sat (= N s S) for the N s -site system with the ratio r. The normalized magnetization is denoted by m = M/M sat . We focus our attention on the magnetization process at zero temperature. For given N s and r, we evaluate the magnetic field where the magnetization increases from M to M + 1 at the field Some of the Lanczos diagonalizations were carried out using an MPI-parallelized code that was originally developed in the study of Haldane gaps. 31 The usefulness of our program was confirmed in large-scale parallelized calculations. 21, 32-35

Results and Discussion
Now, we depict our numerical results of the magnetization process in the three cases of r = J 2 /J 1 = 0.1, 0.3, and 2.2 in Fig Since NNN interactions are small, significant differences are not observed between the two cases. Note here that the presence of the m = 1/3 plateau is clearly detected even when NNN interactions are switched on when the amplitude of the interaction is not large. Next, let us observe the case of r = 2.2 presented in Fig. 3(c). The result for N s = 27 shows that there are some finite-size steps with large widths and other steps with significantly smaller widths. The large-width steps appear once every three steps. This behavior should be compared with that in the case of an infinitely large J 2 in the inset of Fig. 3(c). In the limiting case, the system is reduced to three isolated triangular-lattice antiferromagnets of N s = 9, each of which shows finite-size steps owing to sublattice systems of N s = 9. Let us return to the main panel of Fig. 3(c); one can recognize that the behavior of the large-width steps for N s = 27 comes from the finite-size characteristics of sublattice systems. Therefore, the small-width steps originate from the interactions of J 1 , which is smaller than J 2 , as an effect of perturbation. A similar behavior of the large-width steps can clearly be observed in the result for N s = 36 at m = 1/3, 1/2, and 2/3, although the behaviors at m = 1/6 and 5/6 seem present but weaker. At least for m = 1/3, the large-width steps survive for both N s = 27 and 36. In the case of r = 0.3 in Fig. 3(b), on the other hand, there is a marked difference from the two cases of r = 0.1 and 2.2; one cannot find clear plateaulike behaviors at any height including m = 1/3.
Next, let us investigate the stability of the nonplateau behavior at m = 1/3 in Fig. 3(b) when r is changed. To know the stability, we start our analysis under the assumption of a nonplateau situation. Then, the energy per site in the thermodynamic limit ǫ(m) as a function of m is defined as If we assume that ǫ(m) is an analytic function of m, the spin excitation energy would become magnetization curve as follows: Minimizing the energy of the total Hamiltonian H 0 + H Zeeman yields the magnetization curve at zero temperature using h = ǫ ′ (m)/S. The field derivative of the magnetization is defined as χ mag ≡ dm/dh = S/ǫ ′′ (m). When the nonplateau behavior appears, the derivative should become a nonzero value, namely, χ mag = 0. In this case, Eq. (6) can be rewritten as This means that, if there is a region of r where N s ∆ Ns is almost constant with increasing N s , one recognizes that the region is a nonplateau. If a plateau exists, on the other hand, N s ∆ Ns would increase with N s . Let us apply this argument to the case of m = 1/3 and discuss our present numerical result of N s ∆ Ns shown in Fig. 4, which is expected to privide some valuable information concerning whether or not the plateau at m = 1/3 is present. From the criteria explained above, Fig. 4 shows that the plateau is present from r = 0 to r ∼ 0.15. The behavior markedly changes between r ∼ 0.15 and r ∼ 0.2. From r ∼ 0.3 to r ∼ 0.7, one cannot find a significant size dependence of N s ∆ Ns , which strongly suggests that the system is in the nonplateau region. The presence of this region is a primary result of the present study. At r ∼ 0.8, the behavior of N s ∆ Ns for N s = 27 shows a change in the r dependence; on the other hand, a corresponding change is not observed in the result for N s = 36. As a consequence of the change in N s = 27, from r ∼ 0.8 to r ∼ 1.5, N s ∆ Ns gradually decreases as N s increases; the characteristics of this region are unclear at present. At r ∼ 1.5, N s ∆ Ns for N s = 36 shows a continuous but marked change in its increase as a function of r. Above r ∼ 1.6, therefore, N s ∆ Ns for a given r clearly increases as N s increases, which strongly suggests that the system is in the plateau region. The presence of the nonplateau region necessarily indicates that there is a transition point at r = r c1 between the point with a plateau at r = 0 and the nonplateau region. The presence of the nonplateau region also indicates that there is a transition point at r = r c2 between the point with a plateau in the limit of r → ∞ and the nonplateau region.
In order to capture the boundaries r c1 and r c2 well for a given N s when the transition is continuous, it is useful to observe the second derivative of the ground-state energy with respect to the parameter of the model of interest. 36,37 Such a derivative is defined as χ ene ≡ −∂ 2 [E r (N s , M )/N s ]/∂r 2 . Numerically, we evaluate this quantity using First, let us discuss the behavior of the smaller-r peaks. It is notable that the position of the smaller-r peaks almost does not change with respect to N s . The position is in good agreement with the position where the transition between the small-r plateau region and the nonplateau region occurs, which is suggested from the examination of N s ∆ Ns in Fig. 4. It is reasonable to consider that the smaller-r peaks for each N s correspond to the transition. Therefore, the occurrence of the transition is evident at r ∼ 0.17, namely, at r c1 ∼ 0.17. Next, let us discuss the behavior of the larger-r peaks. Although the position of the larger-r peaks shows a significantly large change with respect to N s , both r ∼ 0.71 for N s = 27 and r ∼ 1.41 for N s = 36 are in agreement with the position where the behavior of N s ∆ Ns for each N s markedly changes. The observation of the larger-r peaks strongly suggests that the transition certainly occurs between the nonplateau region and the large-r plateau region, although it is difficult to precisely estimate its transition point. To summarize, therefore, our present analysis results suggest that the system shows a plateau at m = 1/3 for a small r, that the plateau disappears once at r = r c1 , and that the plateau opens again for r > ∼ r c2 . In particular, we successfully estimate r c1 ∼ 0.17. On the other hand, it is still difficult to estimate r c2 precisely. Precise estimation of r c2 should be carried out in future studies. The present model of the two limiting cases J 2 = 0 and J 2 → ∞ forms the up-up-down states at m = 1/3. These states are magnetized under a magnetic field and are both considered to be collinear. Such collinear states include components that are typical, namely, they have significantly large weights in eigenstates. Since we take the basis in our diagonalizations so that each element is expressed by either an up spin or a down spin, the behaviors of collinear states are captured by observing the weights of such elements. In the plateau region of small r, all spins in one subspace among A, B, and C subspaces are down; all spins in the remaining two subspaces are up. This situation is illustrated in Fig. 6(a).
In the plateau region of large r, on the other hand, in each of the A, B, and C subspaces, (1/9)N s spins are down and (2/9)N s spins are up. Within each sublattice, the positions of (1/9)N s down spins are located in the pattern illustrated in Fig. 6(a). Under such situations, there are two possible patterns of spin configurations of the whole system; the two patterns are illustrated in Figs. 6(b) and 6(c). In the configuration in Fig. 6(b), all down spins are located along a linear line. In the configuration in Fig. 6(c), on the other hand, an island is composed of three down spins; each island is located as far from each other as possible. Both patterns can produce a plateau at m = 1/3. Let us then examine the weights of the three patterns. The weights are evaluated as the sum of squared coefficients of spin configurations in the normalized ground state at m = 1/3 when the selected configurations are linked to a pattern illustrated in Figs. 6(a)-6(c) by a translational or rotational symmetry of the system. Let us observe the weights for the N s = 36 system that show no discontinuities in the results presented in Fig. 5; the results are depicted in Fig. 7. It is observed that the pattern in Fig. 6(a) has a large weight for r that is smaller than r ∼ 0.15. At r ∼ 0.2 and above, this pattern loses its weight. This behavior is in good agreement with the behavior in Fig. 4 and r c1 ∼ 0.17 obtained from Fig. 5. For r > 1, on the other hand, the weights gradually increase for the patterns in Figs. 6(b) and 6(c). At r ∼ 1.4 and above, both weights are significantly large; the weight for the pattern in Fig. 6(b) is larger than that in Fig. 6(c). This behavior suggests that the m = 1/3 plateau in the large-r region originates from the formation of the spin configuration of the pattern in Fig. 6(b). This suggestion should be confirmed by other calculations in future studies. In the present system, the above-mentioned states are magnetized under a magnetic field and are both considered to be collinear; the positions of up or down spins are different between the two cases, as described in the last paragraph. The present study reveals that an intermediate region without a plateau appears between the two different up-up-down states. A similar situation is known in the square-lattice Heisenberg antiferromagnet with NNN interactions, the so-called J 1 -J 2 model. [38][39][40][41][42][43][44][45][46][47][48][49][50][51][52] The properties of this system have long been investigated between the two phases, one is the ordinary Néel phase and the other is a phase of another collinear state when J 2 /J 1 is varied. A recent variational Monte-Carlo study 51 suggested that such an intermediate region is divided into more than one phase. From an analogy of the square-lattice J 1 -J 2 model, the obtained intermediate region in the present model (1) may be composed of complex phases. This issue should be investigated in future studies.
At m = 1/3, the magnetization plateau appears in various frustrated systems. In the kagome-lattice antiferromagnet, the plateau is accompanied by anomalous critical behavior with critical exponents different from the exponent δ = 1, which is typical of two-dimensional systems, just outside of the edges of the m = 1/3 state in the magnetization curve, 33,53 where the critical exponent δ is defined as |m − m c | ∼ |h − h c | 1/δ near the transition point h c . The Cairo-pentagon-lattice antiferromagnet 54, 55 and square-kagome lattice antiferromagnet 56, 57 reveal a magnetization jump at an edge of the m = 1/3 state. A similar jump also appears for the kagome-lattice antiferromagnet with a distortion. 58,59 On the other hand, the magnetization plateau of the triangular-lattice antiferromagnet without NNN interactions shows the typical exponent δ = 1 in the two-dimensional systems. 11 The critical behavior around the edges of the m = 1/3 state of the present system should be clarified in future studies in which NNN interactions are switched on.

Conclusion
We investigated the ground-state magnetization process of the spin-1/2 Heisenberg antiferromagnet on a triangular lattice with next-nearest-neighbor interactions by the numerical-diagonalization method. For small amplitudes of NNN interactions, the m = 1/3 plateau of this system survives; a plateau of the same height exists for large amplitudes of NNN interactions. We have found that, in an intermediate region, this plateau disappears. In particular, the boundary of the intermediate region on the smaller-J 2 side is found to be J 2 ∼ 0.17J 1 , which should be examined to obtain a precise estimate in the future. To precisely know where the boundary on the larger-J 2 side is, investigations of larger systems are required. Further study of phenomena due to NNN interactions in a triangular-lattice antiferromagnet would greatly contribute to our understanding of the frustration effect in quantum spin systems. Grant Numbers 16K05418, 16K05419, and 16H01080 (JPhysics). Nonhybrid thread-parallel calculations in numerical diagonalizations were based on TITPACK version 2 coded by H. Nishimori. In this research, we used the computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research projects (Project ID: hp170017, hp170028, hp170070, and hp170207). Some of the computations were performed using facilities of the Department of Simulation Science, National Institute for Fusion Science; Institute for Solid State Physics, The University of Tokyo; and Supercomputing Division, Information Technology Center, The University of Tokyo. This work was partly supported by the Strategic Programs for Innovative Research; the Ministry of Education, Culture, Sports, Science and Technology of Japan; and the Computational Materials Science Initiative, Japan.