High-order coupled cluster method calculations of spontaneous symmetry breaking in the spin-half one-dimensional J 1 – J 2 model

In this article we present new formalism for high-order coupled cluster method (CCM) calculations for “generalized” ground-state expectation values and the excited states of quantum magnetic systems with spin quantum number s > 1 2 . We use high-order CCM to demonstrate spontaneous symmetry breaking in the spin-half J1–J2 model for the linear chain using the coupled cluster method (CCM). We show that we are able to reproduce exactly the dimerized ground (ket) state at the Majumdar-Ghosh point (J2/J1 = 1


Introduction
The formation of dimer-and plaquette-ordered singlet ground states (so-called valence-bond crystal (VBC) states) is an interesting and important topic in quantum spin systems.Often, the formation of enhanced dimer or plaquette correlations is driven by frustration, which can increase quantum fluctuations and which may result in such gapped rotationally-invariant quantum paramagnetic states .Usually, VBC states are complicated quantum many-body states, see, e. g., the Heisenberg antiferromagnet on the star lattice [10,24,25].However, for certain systems the VBC states are simple exact product eigenstates of the underlying Heisenberg interaction Hamiltonian.Examples for the appearance of such exact VBC product eigenstates are the spin-half J 1 -J 2 model on the linear chain [1][2][3][4][5][6][7][8] at the point J 2 /J 1 = 1 2 (the so-called Majumdar-Ghosh point) and the Shastry-Sutherland model [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23].Furthermore, it is often useful to distinguish between VBC phases that have the same translational symmetry as the spin Hamiltonian and those that spontaneously break the symmetry of the underlying spin lattice.Examples of the former case are the Shastry-Sutherland model and the Heisenberg antiferromagnet on the star lattice, whereas the J 1 -J 2 model on the linear chain is an example of spontaneous symmetry breaking.
Another mechanism for the formation of non-magnetic dimer or plaquette VBC ground states that does not involve frustration is the competition between non-equivalent antiferromagnetic nearest-neighbor bonds.This may lead to the formation of local singlets of two (or four) coupled spins if the strengths of the non-equivalent bonds differ sufficiently [10,[26][27][28][29][30][31][32][33][34][35][36][37][38].By contrast to frustration, which yields competition in quantum as well as in classical spin systems, this type of competition is present only in quantum systems.The symmetry of the ground state follows the symmetry of the Hamiltonian in such cases.
In this article we focus on the application of the coupled cluster method (CCM) [23,35,37, to the spin-half, one-dimensional J 1 -J 2 model.The Hamiltonian for this spin-half model has nearest-neighbor bonds of strength J 1 and next-nearest-neighbor bonds of strength J 2 is given by where the index i runs over sites on the lattice counting, ρ 1 runs over all nearest-neighbors to i, and ρ 2 runs over all next-nearest-neighbors to i.Here we use a Néel model state in which nearestneighbor spins on the linear chain are antiparallel.We rotate the spin coordinates of the 'up' spins so that notationally they become 'down' spins in these locally defined axes.The relevant Hamiltonian in these rotated coordinates is then given by (2) Henceforth we put J 1 = 1 and consider J 2 0.
The ground-state properties of this system have been studied using methods such as exact diagonalizations [2,7], DMRG [3][4][5]55], CCM [8,52,54,78,85], and field-theoretical approaches [5] (see [5,6] for a general review).Note that previous CCM studies of the model considering only independent-spin product model states that conserve the lattice symmetry are reported in [52,79].At J 2 /J 1 = 0 we have the unfrustrated Heisenberg antiferromagnet, where the exact solution is provided by the Bethe Ansatz.The ground state is gapless and the spin-spin correlation function s i • s j decays slowly to zero according to a power-law, i. e. no true Néel-like long-range order is observed.In the region J 2 /J 1 > 0 the nearest-neighbor (J 1 ) and next-nearest-neighbor interactions (J 2 ) compete, thus leading to frustration.At J 2 /J 1 = 0.2411 (1) the model exhibits a transition to a two-fold degenerate gapped dimerized phase with an exponential decay of the correlation function s i • s j [2,3,5,6].This state breaks the translational lattice symmetry.At the Majumdar-Ghosh point J 2 /J 1 = 1  2 there are two degenerate simple exact dimer-singlet product ground states corresponding to the dimerized product state for the Hamiltonian of equation ( 2) [1].
We note again that we treat the J 1 -J 2 model in the regions J 1 > 0 and 0 J 2 0.6 using the CCM in this article.We show how the dimerized ground-state at J 2 /J 1 = 1 2 is exactly reproduced by the CCM ket-state wave function.We consider the CCM correlation coefficients and show that the onset of the dimerized phase is indicated by a bifurcation of nearest-neighbour correlation coefficients.Finally, we demonstrate that the excited-state energy gap is predicted well by the CCM using a localized approximation scheme.New high-order CCM formalism for the ground and excited states for spin quantum number s 1  2 is presented in the appendices.

Method
The details of the practical application of high-order coupled cluster method (CCM) formalism to lattice quantum spin systems are given in [60,64,69,74] and also in the appendices to this article.The high-order CCCM code is freely available online [73].However, it is important for the reader who is not interested in the technical detail to note at this point that the exact ket and bra groundstate energy eigenvectors, |Ψ and Ψ|, of a general many-body system described by a Hamiltonian H, are given by Furthermore, the ket and bra states are parametrized within the single-reference CCM as follows: One of the most important features of the CCM is that one uses a single model or reference state |Φ that is normalized.This, in turn, leads to a normalization condition for the ground-state bra and ket wave functions ( Ψ|Ψ ≡ Φ|Φ = 1).The model state is required to have the property of being a cyclic vector with respect to two well-defined Abelian subalgebras of multi-configurational creation operators {C + I } and their Hermitian-adjoint destruction counterparts {C − I ≡ (C + I ) † }.The CCM formalism is exact in the limit of inclusion of all possible multi-spin cluster correlations within S and S, although this is usually impossible to achieve practically.It is therefore necessary to utilize various approximation schemes within S and S.Here we use the localized LSUBm scheme, in which all multi-spin correlations over distinct localities on the lattice defined by m or fewer contiguous sites are retained.Another important feature of the method is that the bra and ket states are not always explicitly constrained to be Hermitian conjugates when we make such approximations, although the important Helmann-Feynman theorem is always preserved.We remark that the CCM provides results in the infinite-lattice limit N → ∞ from the outset.
The manner in which high-order CCM may be solved has been discussed extensively elsewhere.For example, we refer the interested reader to [64,69] for additional extensive discussions of highorder CCM for the ground and excited states.Recently, however, considerable extensions to the basic high-order formalism for quantum spin systems have also been made.Hence, the formalism and computational algorithms used to derive and solve CCM ket-state and bra-state equations to high orders of approximation are also discussed in detail in Appendices A-C.Furthermore, we present new high-order formalism for "generalized" expectation values using the CCM in Appendix D. The manner in which these high-order techniques may be parallelized in order to reach extremely high orders of approximation for the ground state is considered in Appendix E. We remark that the excited state is constructed by solving the Schrödinger equation for the excitedstate, namely, H|Ψ e = E e |Ψ e , and where |Ψ e = X e e S |Φ .X e is called the excited-state operator and it may be used (as here) to form an excited state of different symmetry to the ground state.The basic CCM formalism for the excited state [64] is presented in Appendix F. Again, high-order CCM may be carried out for the excited state in an analogous manner as for the ground state.New high-order CCM formalism for the excited state for s 1  2 is given in Appendix G.The manner in which the excited-state equations are derived and solved is given in Appendix H. Finally, this problem may be solved using direct iteration of the CCM equations using a "shifted" power iteration approach and this approach is presented in Appendix I.
As mentioned previously, we use a Néel model state in which nearest-neighbour spins are antiparallel.Furthermore, the ground state lies in the subspace s z T ≡ i s z i = 0, whereas the excited state has s z T ≡ i s z i = +1.We use a "doubled" unit cell including two neighboring sites for this spin-half system on the linear chain at points (0,0,0) and (1,0,0) and a single Bravais vector (2,0,0) T to take into account the symmetry breaking.Thus, there are two distinct types of twospin nearest-neighbor ket-state correlation coefficients, namely, those connecting the sites inside the unit cell and those connecting different unit cells.These coefficients are denoted as S a 2 and S b 2 .It is straightforward to prove [85] that the ground state at J 2 /J 1 = 1 2 is obtained exactly by setting S a 2 = 1 and all other coefficients equal to zero.Starting from J 2 /J 1 = 1 2 we are able to track this exact solution at J 2 /J 1 = 1 2 within a certain LSUBm approximation for other values of J 2 /J 1 .

Results
The results for the nearest-neighbor ket-state correlation coefficients in LSUB14 approximation are presented in figure 1.Clearly, we see that the exact dimerized product-state solution for the ket ground state is obtained within LSUB14 level of approximation (and, indeed, at all LSUBm approximations with m 2) at J 2 /J 1 = 1 2 , i. e. S a 2 = 1 and all other coefficients equal to zero.Moving away from J 2 /J 1 = 1 2 we still find a CCM ground state that breaks the lattice symmetry.However, this dimerized state deviates from the simple product state, i. e. S a 2 = 1 and other non-zero coefficients S I occur.Furthermore, the solution (i.e. S a 2 = S b 2 ) having full translational symmetry is the only solution below a critical point J 2 /J 1 | c1 (< 1  2 ).Henceforth, we shall refer to this solution below J 2 /J 1 | c1 as the  "usual ('Néel-type') solution" because previous CCM calculations [52,54] for the J 1 -J 2 model have considered the non-symmetry breaking case only.For larger values of J 2 /J 1 a CCM termination point is observed at 2 ), shown by the boxes in figure 1.At this point, the real solution of the CCM dimerized solution is terminated.These CCM results indicate that a dimerized phase exists over a finite range of J 2 /J 1 , which is in agreement with the known results, see e. g. [3,5,6].Qualitatively similar results are observed at other levels of LSUBm approximation for the ket-state correlation coefficients as a function of J 2 /J 1 .The results for J 2 /J 1 | c1 and J 2 /J 1 | t are shown in table 1.It is obvious that the CCM critical point J 2 /J 1 | c1 becomes smaller (i.e., becomes closer to the true critical point J 2 /J 1 = 0.2411(1) [3,6]) with higher orders m of LSUBm approximation.However, the critical point J 2 /J 1 | c1 is still significantly too high even at the LSUB14 level of approximation.This subject is considered in more detail when we discuss the CCM results for the excitation energy gap below.On the side of J 2 /J 1 > 1 2 the existence of termination points can be related to the appearance of incommensurate spiral spin correlations at J 2 /J 1 > 0.6 [5,7,55,78] that are not taken into account in the model state used here.
The nearest-neighbor bra-state correlation coefficient at the LSUBm level of approximation at J 2 /J 1 = 1 2 has Sa 2 = 1/4 with m 4.This is shown in figure 2 for the LSUB14 level of approximation.We find that the bra-state solution for the nearest-neighbor correlation coefficients is close to 1/4 over the range However, we find that the nearestneighbor correlation coefficient diverges as J 2 /J 1 → J 2 /J 1 | c1 and this is also shown in figure 2. Again, the usual ('Néel-type') solution ( Sa The upper CCM termination point at J 2 /J 1 | t is also shown in figure 2 by the boxes on the right-hand side of the figure.Once more, qualitatively similar results are observed at other levels of LSUBm approximation for the bra-state correlation coefficients as a function of J 2 /J 1 .We now consider the ground-state energy of this system in the dimerized regime, and results of LSUBm approximation are shown in figure 3 as a function of J 2 /J 1 .We observe that the LSUBm converge rapidly with increasing level of approximation.Furthermore, we see that CCM results compare well to those results of exact diagonalizations in the regime 0 We study the details of the ground state energy for the new dimer solution and the usual ('Néel-type') solution for the LSUB14 level of approximation in figure 4. Again, we remark that our solution is an exact ground eigenstate at J 2 /J 1 = 1 2 .The exact ground-state energy of E g /N = −0.375J 1 is obtained at the point J 2 /J 1 = 1 2 at all levels of approximation, as expected.We also see that ground-state energy of the usual ('Néel-type') CCM solution in which S a 2 = S b 2 at the LSUB14 level of approximation actually lies below this exact solution.This indicates (i) that the usual ('Néeltype') CCM solution is a relatively poor choice at this point; and, (ii) that the CCM ground-state energy does not fulfill the variational principle [62].Furthermore, we see that CCM dimer solution  The results for the sublattice magnetization M of this model are presented graphically in figure 5. Since the one-dimensional J 1 -J 2 model does not possess Néel long-range order for any value of J 1 , J 2 0 the true value is M = 0.As is known from previous CCM calculations [52,64,78,79], the sublattice magnetization is nonzero (but small) using the usual Néel model state.However, the correct result M = 0 can be obtained [78,79] by extrapolating the 'raw' LSUBm data to m → ∞.Indeed we see that the CCM LSUBm values for M are non-negligible for Néel model state in the region J 2 /J 1 < J 2 /J 1 | c1 .It is also clear that M decreases with the level of approximation m approaching the true value M = 0 and that increasing the strength J 2 of the frustration weakens the magnetic order.More interestingly, we find that the sublattice magnetization behaves discontinuously at J 2 /J 1 | c1 by tracking the lattice symmetry-breaking dimerized solution, and then remains near to zero at LSUB10 to LSUB14 levels of approximation across the entire range (At the lower LSUB8 level of approximation the results for the sublattice magnetization differ from zero by a small amount in a small region above J 2 /J 1 > 1 2 .)On the other hand, by tracking the usual ('Néel-type') solution it is found (see [85]) that M changes continuously with J 2 and it is larger than for the dimerized solution for J 2 /J 1 | c1 < J 2 /J 1 < 1 2 .This behavior of M is another indication that the dimerized CCM solution describes the true physics of the model much better than the usual Néel solution.We note finally that the CCM sublattice magnetization is exactly zero at the Majumdar-Ghosh point J 2 /J 1 = 1 2 at all levels of LSUBm approximation using the dimerized product state.We consider the results for the excitation energy gap, , at J 2 /J 1 = 0 and J 2 /J 1 = 1 2 in figure 6.We find that the values for at J 2 /J 1 = 0 are given by 0.56040, 0.38346, 0.29025, 0.23310, 0.19458, 0.16689 at the LSUBm levels of approximation with m = {4, 6, 8, 10, 12, 14}.These results agree with previous results up to the LSUB12 level of approximation for the spin-half linear chain Heisenberg model [64].Furthermore, we find that an quadratic extrapolation rule at J 2 /J 1 = 0.0 gives −0.0036, also shown in figure 6.This result is in reasonable agreement with the exact result that is known is to be zero from the Bethe Ansatz at this point.By contrast, we find values for at J 2 /J 1 = 1 2 of 0.35250, 0.34170, 0.30548, 0.28732, 0.27559, and 0.26760 at the LSUBm levels of approximation with m = {4, 6, 8, 10, 12, 14}.These results are also shown in figure 6.We also see from this figure that the simple extrapolation of these results in the limit m → ∞ using a quadratic function gives a value for the gap of 0.2310.This result is in agreement with results of exact diagonalizations that predict a gap of 0.234 [86].We conclude that the CCM provides good results for the ground and the excited states at the point J 2 /J 1 = 1 2 for the new "symmetry breaking solution".
Results for the excitation energy gap using the LSUBm approximation presented as a function of J 2 /J 1 are given in figure 7. We see that the gap for all levels of LSUBm approximation is non-zero and positive in the region J 2 /J 1 0.6.We use the quadratic extrapolation rule used in figure 6 that provided reasonable results at J 2 /J 1 = 0 and J 2 /J 1 = 1 2 in order to extrapolate our results at other values of J 2 /J 1 .The results for the extrapolation of LSUBm data with m = {6, 8, 10, 12, 14} are given in figure 7. (Extrapolations in the region where J 2 /J 1 ≈ J 2 /J 1 | c1 are difficult to achieve accurately and so are not attempted here.)These results show that the gap is negative but close (i.e., | | < 0.003) for a finite region above J 2 /J 1 = 0. Indeed, the gap becomes zero at J 2 /J 1 = 0.238 and then becomes positive and finite for increasing values of J 2 /J 1 above this point.
This is an interesting result because we have obtained a gap for the extrapolated results for the "usual Néel" solution, i. e., without symmetry breaking.Indeed, this is somewhat disconcerting at first glance because one expects from the Lieb-Schulz-Mattis theorem that the ground state is either gapless and non-degenerate or that it is gapped and degenerate for half-integer spin Heisenberg chains.Hence, one might expect that the gap should appear at that point where we find symmetry broken solution, namely, for values of J 2 /J 1 greater than J 2 /J 1 | c1 .This is clearly not the case in figure 7.However, we remember that we are starting from a Néel model state and so we are explicitly breaking the symmetry.Furthermore, we note that we are using the localized (but systematic) LSUBm approximation scheme in order to include quantum effects.Hence, the gap is finite when determined at a given level of LSUBm approximation even at J 2 /J 1 = 0, as noted above.As was also found to be the case for the square and cubic lattice Heisenberg antiferromagnets [64], the correct solution at J 2 /J 1 = 0 of a gap of zero (i.e., the continuous symmetry breaking solution) is recovered only when we extrapolate in the limit m → ∞.Hence, this analysis does not necessarily exclude the possibility of a non-zero extrapolated gap opening up below J 2 /J 1 | c1 .Indeed, our result of J 2 /J 1 = 0.238 is very close to the "true" value J 2 /J 1 = 0.2411(1) from [3,6].
The presence of a gap would be consistent with an "effective system" with residual anisotropy in the z-direction, namely, on the s z s z terms in the Hamiltonian.The J 1 -J 2 model is known [3] to be a special case of the more general model with anisotropy ∆ on the s z s z terms in the Hamiltonian.Furthermore, the line ∆ = 1 (for J 2 /J 1 < 0.2411) is known [3] to be a boundary between a Luttinger spin liquid (gapless) for ∆ < 1 and a Néel ordered regime (with a finite gap) for ∆ > 1.If the "effective system" does indeed have a "residual anisotropy" then we would expect the value of the boundary to change as well.Indeed, from the results of [64], we would expect its value in terms of J 2 /J 1 to increase.This is exactly what we find; the spontaneous symmetry breaking solution occurs at values of J 2 /J 1 | c1 , which is well above 0.2411 (as shown in table 1).Presumably, as we carry out higher and higher orders of approximation then J 2 /J 1 | c1 will tend to J 2 /J 1 = 0.2411(1) [2,3,5,6].
We note that the results for the excitation energy gap in the region J 2 /J 1 > J 2 /J 1 | c1 are very highly converged.The LSUB14 result for the "usual Néel solution" is also given in figure 7. We see that these results lie much lower than the equivalent LSUB14 result for the dimer solution for J 2 /J 1 > J 2 /J 1 | c1 .Indeed, this usual Néel solution for LSUB14 provides poor results in the dimerized region.Furthermore, we see that the gap only becomes larger in the region J 2 /J 1 > 0.4, as predicted by DMRG results [5].However, the peak in the energy gap is predicted by the CCM to occur at J 2 /J 1 ≈ 0.52, whereas DMRG results [5] clearly indicate that this peak occurs at J 2 /J 1 ≈ 0.6.This probably indicates that our model state is a poor choice for J 2 /J 1 > 1 2 .

Conclusions
We have shown in this article that we can form dimer VBC ground states using the CCM by considering the spin-half J 1 -J 2 model for the linear chain.We showed that we are able to exactly reproduce the dimerized ground state at J 2 /J 1 = 0.5.Interestingly, a spontaneous symmetry-breaking dimerized CCM solution is observed for J 2 /J 1 < 1 2 , which only becomes equal to the usual ('Néel-type') solution that conserves the lattice symmetry at a CCM critical point J 2 /J 1 | c1 .Results for the bra state correlation coefficients diverged at this point as well.We took this to indicate the onset of the dimerized ground-state phase that breaks the translational lattice symmetry.We also found that termination points occurred in the CCM equations for J 2 /J 1 > 1 2 , and we took this to be an indication of another critical point in this region.Results for the ground-state energy for the dimerized CCM solution were found to agree extremely well with the results of exact diagonalizations for N = 28 and N = 32 chains in the dimerized regime.The effects of the bifurcation point at given levels of LSUBm approximation were also seen in the results for the sublattice magnetization that demonstrates a discontinuity at J 2 /J 1 | c1 .
Results for the excitation energy gap of this model were considered.We found that extrapolated results for the excitation energy gap at J 2 /J 1 = 0.0 and J 2 /J 1 = 1 2 of −0.0036 and 0.2310, respectively, were in reasonable agreement with the known results at these points.Interestingly, results for the extrapolated excitation energy gap indicated that a gap opens up at J 2 /J 1 = 0.238, which is in excellent agreement with the known result of J 2 /J 1 = 0.2411(1) from other studies [2,3,5,6].However, this result was found for the non-symmetry breaking solution to our equations.Superficially, this result seems to violate the Lieb-Shultz-Mattis theorem, namely, that the ground state is either gapless and non-degenerate or that it is gapped and degenerate for half-integer spin Heisenberg chains.Hence, one might expect a gap to occur only when we reach the spontaneous symmetry breaking solution to our equations, namely, for J 2 /J 1 J 2 /J 1 | c1 .However, we know already from [64] for the linear-chain, and square-and cubic-lattice antiferromagnets that LSUBm results yield a gap at a given level of LSUBm approximation.We start from a model state with spins in the z-direction (explicitly breaking the symmetry) and the zero gap result (i.e., for continuous symmetry breaking) is only ever reached when we extrapolate the LSUBm results.This "allows" a gap to form for J 2 /J 1 > 0.238 even for the non-spontaneous symmetry breaking solution to our CCM equations.We seem to have a new intermediate region from 0.24 < J 2 /J 1 < J 2 /J 1 | c1 , although this is simply an artefact and it will disappear as we increase the approximation level.However, we are clearly seeing a radical change in the nature of the ground state near to J 2 /J 1 ≈ 0.24 that is in good agreement with known results [2,3,5,6].Results for the excitation energy gap are also found to become large at J 2 /J 1 ≈ 0.4 in agreement with DMRG [5], and they show that the model state is probably a poor choice for J 2 /J 1 > 1 2 .Finally, new high-order CCM formalism for the ground state (including "generalized" expectation values) as well as for the excited state for spin quantum number s 1 2 is presented in the appendices.

A. CCM ground-state formalism
We begin the description of the CCM high-order formalism by noting again that we wish to solve the Schrödinger equation of equation (3).We remark again that the bra and ket states are given by equation ( 4).It may be proven from equations ( 3) and ( 4) in a straightforward manner that the ket-and bra-state equations are thus given by Φ|C − I e −S He S |Φ = 0, ∀I = 0, (5) The index I refers to a particular choice of cluster from the set of (N F ) fundamental clusters that are distinct under the symmetries of the crystallographic lattice and the Hamiltonian and for a given approximation scheme at a given level of approximation.We note that these equations are equivalent to the minimization of the expectation value of H = Ψ|H|Ψ with respect to the CCM bra-and ket-state correlation coefficients { SI , S I }.We note that equation ( 5) is equivalent to δ H/δ SI = 0, whereas equation ( 6) is equivalent to δ H/δS I = 0. Furthermore, we note that equation ( 5) leads directly to simple form for the ground-state energy given by The full set {S I , SI } provides a complete description of the ground state.For instance, an arbitrary operator A will have a ground-state expectation value given as The similarity transform of A is given by,

B. High-order ground-state operators and commutations
We begin the treatment of high-order CCM by introducing the ket-state correlation operator given, as usual, by However, it is an important point to note that each of the indices {i 1 , i 2 , • • • , i l } runs over all lattice sites.Furthermore, we assume that there are (l!) orderings of these indices (even for s > 1 2 ), although we never need to work out these factors explicitly in practice.The index I corresponds to one of the choices of {i 1 , • • • , i l } for the fundamental set of configurations.We may now write a set of high-order CCM ket-state operators, given by The indices k, m, n, and p depend on those sums in the Hamiltonian or of another given operator.We note that s ± = s x ±is y , [s z , s ± ] = ±s ± , and [s − , s + ] = −2s z .Hence, the following commutation relations may be proven: We may now write the similarity-transformed expressions of the single-spin operators s α ; α ≡ {+, −, z}, as e −S s + k e S ≡ s+ We see that there is a repeated index in G kk in the similarity transformed version of s − .Clearly, this term contributes only for systems with spin quantum number s > 1 2 .

C. Deriving and solving the CCM ground-state equations
We now wish to determine and solve the CCM ket-state equations, where the I-th such equation is given by (Note that we assume that Φ|C − I C + I |Φ = 1 in the above equation).Specific terms in the Hamiltonian are now explicitly written in terms of the high-order CCM operators as: (Note that s − |Φ = 0 is implicitly assumed in equation ( 15) above.)We now "pattern-match" the C − i operators to those relevant terms in the Hamiltonian from equation (15) above in order to form the CCM equations E I = 0 of equation ( 14) at a given level of approximation.These coupled non-linear equations are then solved readily by using, e.g, the Newton-Raphson method.However, these are solved via direct iteration for larger values of the approximation level because the cost of storing the Jacobian in local memory for the Newton-Raphson (or other) method becomes prohibitive.This may be also parallelized to achieve very high orders of approximation and this is discussed below.
We now define the following new set of CCM bra-state correlation coefficients given by x I ≡ S I and xI ≡ N B /N (l!)ν I SI and we assume again that Φ|C − I C + I |Φ = 1.Note that N B is the number of Bravais lattice sites.Note also that for a given cluster I then ν I is a symmetry factor which is dependent purely on the point-group symmetries (and not the translational symmetries) of the crystallographic lattice and that l is the number of spin operators.We note that the factors ν I , N , N B , and (l!) never need to be explicitly determined.The CCM bra-state operator may thus be rewritten as such that we have a particularly simple form for H, given by where x0 = 1.We note that the E 0 is defined by E 0 = 1 N Φ|e −S He S |Φ (and, thus, E 0 = 1 N E g ) and that E I is the I-th CCM ket-state equation defined by equation (14).The CCM ket-state equations are easily re-derived by taking the partial derivative of H/N with respect to xI , where We now take the partial derivative of H/N with respect to x I such that the bra-state equations take on a particularly simple form, given by The equation ẼI = 0 is easily solved computationally via LU decomposition for low to medium orders of approximation or via direct iteration (which may be parallelized, as discussed below) for even higher orders of approximation.This may be carried out once the CCM ket-state equations have been determined and solved.The numerical values of the coefficients {x I } may thus be obtained.We note that this approach greatly simplifies the task of determining the bra-state equations because we infer the bra-state equations directly from those of the ket-state equations via equation (19).Thus, we never need to evaluate equation ( 6) explicitly.

D. Generalized ground-state expectation values
The expectation value of a "generalized" spin operator that we shall call A may be treated in an analogous manner to that of the expectation value of the Hamiltonian, given by H.We write: and with C − 0 = 1.The similarity transform of A is defined by equation ( 8) and once again this process results in terms such as those shown in equation (15) may again be employed.However, unlike the J 1 -J 2 Hamiltonian of equation (1), we do not constrain k and m in the two-body terms to be only nearest-neighbors or next-nearest-neighbors.The expectation value of the generalized (spin) operator may again be written in a particularly simple form as: where x0 = 1 also and A 0 = 1 N Φ|e −S Ae S |Φ .The same code used to find ground-state equations may be used to find the generalized expectation values.Again we note that the index I in equation (21) runs from zero to N F .Again, we note that factors such as N B or ν I etc. do not need to be determined explicitly because they cancel because of the definition of {x I } given above.

E. Direct iteration of the ground-state equations and parallelization
The parallelization of the ground-state CCM problem for very high-order CCM is achieved by solving the ket-and bra-state equations (i.e., E I = 0 and ẼI = 0, respectively) via direct iteration.For the case of the ket-state equations this is slightly more complicated because there are non-linear terms with respect the ket-state correlation coefficients {x I }.We rearrange the ketstate equations such that the linear terms for x I for the i th CCM ket-state equation are on the left of the new equation and all other terms are on the right.The right-hand side of this new equation is denoted by E I after dividing through by the factor on the left-hand-side for the ketstate correlation coefficients.We may carry out exactly the same procedure for the bra-state in order to find Ẽ I , although the problem is linear with respect to {x I } in this case.These equations are thus rewritten conveniently for the ket state as and for the bra state as, Clearly, these equations may be solved for x I and xI by iterating them "directly" until convergence.Indeed, the local memory usage is vastly reduced because we do not need to store any Jacobian or other large matrix that scales in size with N 2 F .Furthermore, the computational problem posed by solving equations ( 22) and ( 23) via direct iteration may be solved using parallel processing.The different equations of equations ( 22) and ( 23) for different values of the index I are determined separately on different processors.The resulting data for these equations for the different values of I are then stored locally to each processor.At each iteration of the "direct iteration" method we find the right-hand sides of those relevant values of I allocated to each processor.We then collect the right-hand side into a single array and this forms our values for x I or xI for the next iteration.Again, we note that we must solve the ket-state equations of equation ( 22) first and then these values for the ket-state coefficients are used in the bra-state equations of equation ( 23).This approach is a simple "brute-force" method, although it has been found to be surprisingly successful at going to very high orders of approximation.Indeed, we may now treat of order 10 6 fundamental clusters using this approach and for approximately 10 2 − 10 3 processors used in parallel.Clearly, a similar approach may also be used to find the "generalized" expectation values via parallel processing.

F. The excited-state formalism
We now consider how the excited state may be treated using the CCM via a high-order approach.We begin by remarking that the excited-state wave function is given by The Schrödinger equation, E e |Ψ e = H|Ψ e and the equivalent equation for the ground state lead (after some simple algebra) to e X e |Φ = e −S [H, X e ]e S |Φ (≡ R|Φ ) , where e ≡ E e − E g is the excitation energy.We note that the excited-state correlation operator is written as, X e = I =0 Equation ( 26 We may now form the basic equations for the excited state, given by e X e I = Φ|C − I e −S [H, X e ]e S |Φ , ∀I = 0 , which is a generalized set of eigenvalue equations with eigenvalues e and corresponding eigenvectors X e I .We note that the choice of clusters for the excited-state may be different from those for the ground state.For example, the ground state for the model considered here is in the subspace s z T ≡ i s z i = 0, whereas the excited state has s z T ≡ i s z i = +1.The number of excited-state "fundamental" clusters that are distinct under the translational and point-group symmetries of the lattice and Hamiltonian is given by N fe .

G. High-order excited-state operators and commutations
In a similar manner as for the ground-state, we now define excited state operator via where the indices {i 1 , • • • , i l } run over all lattice sites.We assume explicitly again that there are (l!) orderings of the indices (even for s > 1 2 ).The index I corresponds to one of the choices of {i 1 , • • • , i l } for the fundamental set of configurations for the excited state, such that X e I ≡ X e ii,•••,i l .We now also define the further high-order operators for the excited state, given by The following commutation relations may also be proven:

H. Deriving and solving the excited state equations
We now wish to determine and solve the CCM excited-state equations given by equation ( 28) Specific terms in the Hamiltonian are now explicitly written in terms of the new excited-state high-order CCM operators as: (Note that s − |Φ = 0 is again implicitly assumed in equation (32) above.)Again, we now "patternmatch" the C − I operators (this time with respect to the fundamental set of the clusters in the excited state) to those relevant terms in the Hamiltonian from equation (28) above in order to form the CCM excited-state equations at a given level of approximation.By contrast to the case for the ground state, we see that the high-order operators of equation (30) are linear in those terms in equation (32).We choose the eigenvalue of the lowest value to be our result, and this method was found to provide good results in regions of the parameter space for which the model state was a good choice.Again we note that we have formed an eigenvalue problem, which is readily solved using a standard eigenvalue solver.However, the computational problem thus formed uses local memory that scales with the number of fundamental clusters used in the excited state, i. e., as N 2 fe .Again, this becomes prohibitive computationally for extremely large values of N fe and so we again use direct iteration methods.

Figure 1 .
Figure 1.CCM results at the LSUB14 level of approximation for the ground-state nearestneighbor ket-state correlation coefficients of the spin-half J1-J2 antiferromagnet on the linear chain.The nearest-neighbor coefficients S a 2 and S b 2 of the symmetry breaking dimerized solution are shown by the full lines.Results for the usual ('Néel-type') solution in the region J2/J1 > J2/J1|c 1 (where S a 2 = S b 2 ) are shown by the dotted line.Below the (bifurcation) CCM critical point at J2/J1|c 1 there is only the solution with S a 2 = S b 2 .A termination point J2/J1|t of the CCM equations for the dimerized solution, at which point the real solution to the CCM equations terminates, is indicated by the boxes.

Figure 2 .
Figure 2. CCM results at the LSUB14 level of approximation for the ground-state nearestneighbor bra-state correlation coefficients of the spin-half J1-J2 antiferromagnet on the linear chain.The nearest-neighbor coefficients Sa 2 and Sb 2 of the dimerized solution are shown by the full lines.Results for the usual ('Néel-type') solution for J2/J1 > J2/J1|c 1 (where Sa 2 = Sb 2 ) are shown by the dotted line.Below the critical point at J2/J1|c 1 both solutions coincide.Results for these bra-state correlation coefficients diverge at the critical point J2/J1|c 1 .Another termination point J2/J1|t is shown by the boxes on the right-hand side of the figure.

Figure 3 .
Figure 3. CCM results for the ground-state energy of the spin-half J1-J2 antiferromagnet with J1 = 1 on the linear chain for the LSUBm approximation with m = {6, 8, 10, 12, 14}.Results of exact diagonalizations for N = 28 and N = 32 number of sites are also shown.

Figure 4 .
Figure 4. CCM results for the ground-state energy of the spin-half J1-J2 antiferromagnet with J1 = 1 on the linear chain at the LSUB14 level of approximation.The dimerized and usual ('Néel-type') solutions are shown in this figure.Results of exact diagonalizations for N = 28 and N = 32 number of sites are also shown.The CCM bifurcation and termination points for the dimerized solution are shown by the boxes.

Figure 5 .
Figure 5. CCM results for the sublattice magnetization M of the spin-half J1-J2 antiferromagnet on the linear chain.Below the critical point at J2/J1|c 1 the results for the usual ('Néel-type') and the dimerized solution coincide.At J2/J1|c 1 the sublattice magnetization of the dimerized solution exhibits a jump, whereas we remark that M for the usual ('Néel-type' -not shown) solution is continuous.

Figure 7 .
Figure 7. CCM results for the excited-state energy gap, , of the spin-half J1-J2 antiferromagnet on the linear chain plotted with respect to J2/J1.

Table 1 .
CCM results for the positions of the range of the dimerized phase.