{'fraction_non_alphanumeric': 0.04042983648966867, 'fraction_numerical': 0.010878577824947763, 'mean_word_length': 4.511241089380369, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multispinon excitations in the spin S=1/2 antiferromagnetic Heisenberg model

With the commutation relations of the spin operators, we first write out the equations of motion of the spin susceptibility and related correlation functions that have a hierarchical structure, then under the"soft cut-off"approximation, we give a set of equations of motion of spin susceptibilities for a spin S=1/2 antiferromagnetic Heisenberg model, that is independent of whether or not the system has a long range order in the low energy/temperature limit. Applying for a chain, a square lattice and a honeycomb lattice, respectively, we obtain the upper and the lowest boundaries of the low-lying excitations by solving this set of equations. For a chain, the upper and the lowest boundaries of the low-lying excitations are the same as that of the exact ones obtained by the Bethe ansatz, where the elementary excitations are the spinon pairs. For a square lattice, the spin wave excitation (magnons) resides in the region close to the lowest boundary of the low-lying excitations, and the multispinon excitations take place in the high energy region close to the upper boundary of the low-lying excitations. For a honeycomb lattice, we have one kind of"mode"of the low-lying excitation. The present results obey the Lieb-Schultz-Mattis theorem, and they are also consistent with recent neutron scattering observations and numerical simulations for a square lattice.


I. INTRODUCTION
Spins are neither bosons nor fermions, and their commutation relations make spin problems so difficult. A spin system is a subject in which there are few exactly solvable models which are nontrivial. Only a few of the models have solutions [1][2][3][4] which are well understood, in spite of the fact that many of them have been intensely studied [5][6][7]. The most challenge of a spin system is that there is absent of an analytical method directly applied for it beyond one dimension (1D) without the help of the slave particle representations.
Instead of directly studying a spin system, one usually maps it to a many-body problem which is a strongly interacting system. For example, a spin S = 1/2 chain may be mapped exactly into an interacting spinless fermion system, with the help of the Jordan-Wigner transformation [8,9] which is believed to be valid only for one dimension. While for a high-dimensional spin system, the spin operators are usually represented by slave bosons/fermions with some constraint conditions, then the system is mapped into an interacting boson/fermion system which has still not been very successful because of the strong interactions among bosons/fermions [10,11]. In contrast with a bosonic/fermionic system in which the basic low-lying excitation are quasiparticles obeying Bose-Einstein/Fermi-Dirac statistics, the low-lying excitations [12] of a spin S = 1/2 antiferromagnetic system are magnons with spin S = 1, where a magnon may be seen as a bound state (triplet) of two spinons that each spinon has a spin S = 1/2. However, in the 1D case, the spinons are nearly deconfined, and they become the elementary low-lying excitations [13,15] of the system. The calculations [18][19][20][21][22] based on the exact solution of the Bethe ansatz and the recent neutron scattering measurements on one-dimensional [16] spin S = 1/2 Heisenberg antiferromagnets strongly support the picture that the spinons are the elementary low-lying excitations. Recent neutron scattering experiment [17] shows that for two-dimensional (2D) spin S = 1/2 Heisenberg antiferromagnets, the spinons may be nearly deconfined in some short wave-length regions of the Brillouin zone (BZ).
Theoretically, the 2D Heisenberg antiferromagnet has been extensively studied by a variety of numerical approaches that try to completely understanding of these experimental observations with nearly deconfined spinons [17,23] or multi-magnon excitation [24][25][26], where there does not have a convincing unambiguous evidence to support which one of them as the low-lying excitations in short wave length region of the BZ. However, the numerical calculations based on the exact solution of the Bethe ansatz and the neutron scattering experimental observations unambiguously show that the low-lying excitations of the spin S = 1/2 systems are distributed a broad region in the frequency and momentum plane with the upper and the lowest boundaries, especially in the 1D case that the lowest boundary mainly represent the two spinon excitations, that has the same dispersion as that of usual spin wave, while the upper boundary represents the multispinon (pairs) excitation, that has a dispersion with a period twice that of the spin wave.
In contrast to usual slave particle methods, the equation of motion of Green's function approach can be a good candidate in studying of the low-lying excitations of a spin S = 1/2 antiferromagnetic system, in which one can easily write out the equations of motion of all high order correlation functions appearing in the equation of motion of Green's function.
The equations of motion of Green's function and related high order correlation functions are tightly coupled with each other, and they have a simply hierarchic structure that is an unclosed set of equations. In each level of this hierarchic structure there are many correlation functions where they construct a subset of equations.
In Refs. [17,18], the authors used the equation of motion of Green's function to have calculated the low-lying excitation of a 1D spin 1/2 Heisenberg model with usual cut-off approximations taken for high order correlation functions, and they had obtained the lowlying excitation spectrum that is consistent with the exact one [12] only for small momentum.
At larger momentum, however, their results heavily deviate from the exact ones. Moreover, it is hard to have the spin wave (the lowest boundary) of the low-lying excitations as that ones by the Bethe ansatz and the recent neutron scattering measurements, because in the previous calculations of Refs. [17,18] they cannot self-consistently calculated equations of motion of the high order correlation functions appearing in the same level of this hierarchic structure.
In this paper, we first write out of a complete hierarchic structure of the equations of motion of multiple-point correlation functions. Instead of taking usual cut-off approximations for high order correlation functions, we solve self-consistently the equations of motion of multiple-point correlation functions under "soft cut-off " approximations, then we can obtain the upper and lowest boundaries of the low-lying excitations of the magnons/(pair) spinons in the whole BZ for a spin S = 1/2 antiferromagnetic Heisenberg model in 1D and 2D. For 1D, the upper and lowest boundaries of the low-lying excitations have the same dispersion as that ones of the Bethe ansatz in the whole BZ, and for 2D, they are completely consistent with the recent neutron scattering experimental observations and numerical simulations. This paper is organized as follows. In section II, we give a detail explanation of our present method. Under the "soft cut-off" approximations [29], we write out the general expressions of equations of motion of the transverse and longitudinal spin susceptibilities in the N = 1 and N = 2 levels, respectively, in Section III. Then we apply these equations of motion of the transverse and longitudinal spin susceptibilities for the 1D and 2D cases, and calculate the low-lying excitation spectrums of the magnons/(pair) spinons in the whole BZ in Sections IV-VI. Finally we give our conclusions and discussions in Section VII. More technical calculations for the high order multiple-point correlation functions are put in the Appendixes.

II. BASIC IDEA OF THE ALGEBRAIC EQUATION OF MOTION APPROACH
For the spin operators of the spin S = 1/2, they have a SU(2) symmetry. As an unperturbable theory, we extend the hierarchical Green's function approach [29]  In the following sections, we use the bipartite sublattice representation to write out the EOMs of spin susceptibility and related multiple-point correlation functions, and all calculations about the EOMs of the multiple-point correlation functions are made on the lattice sites. The prominent advantage of the bipartite sublattice representation is that it can greatly simplify our calculating for the high order multiple-point correlation functions that are in fact the tensors whose indexes denoted by the lattice site coordinates. Finally, we only retain the results that are independent of the bipartite sublattice representation.

III. BASIC EQUATIONS OF MOTION OF THE SPIN SUSCEPTIBILITY
If we only consider the contributions of the related multiple-point correlation functions belonging to N = 1 level, under the "soft cut-off " approximation [29], we can obtain the following equations of motion (EOMs) of the transverse and longitudinal spin susceptibilities, These EOMs of the spin susceptibility are universal for a general spin S = 1/2 antiferromagnetic Heisenberg model, and they can be used to calculate its low-lying excitation spectrum on a variety of lattice sites, where the static constants can be self-consistently determined by a set of equations of equal-time spin susceptibility derived from the relation, s + i s − i = 1 2 + s z i , and sum rules. However, under this simple approximation, the above EOMs of the spin susceptibility can only give the reliable upper boundary of the low-lying excitation and they cannot give usual spin wave excitation that survive in the large momentum and low energy limit region.
In order to studying the lowest low-lying excitations residing in the large momentum region, we have to consider the contributions of the high order related multiple-point correlation functions belonging to N = 2 level.
Under the "soft cut-off" approximation, as including the contributions of the related multiple-point correlation functions belonging to the N = 2 level (see Appendix C), for example, we can obtain the following EOMs of the transverse spin susceptibility (for simplicity, The coefficients Γ lij (ω) and Π lij (ω) can be approximately written as that, Obviously, as ω > J, the coefficients Γ lij (ω) and Π lij (ω) are positive, and in the limit, ω/J → ∞, they go to zero, then the Eqs. (5,6) are reduced to the Eqs. (3,4), which means that the latter is the high energy limit of the former. To the contrary, in the limit, ω/J → 0, the coefficient . In this case, we can obtain the lowest boundary of the low-lying excitations of the spins.
The poles appearing in the coefficients Γ lij (ω) and Π lij (ω) are artificial, that originate from the approximations taken for the high order related multiple-point correlation functions belonging to the N = 2 level. The summation over the site variables of Γ lij (ω) and Π lij (ω) is very clear for a chain, while for the high dimensional lattice case, such as for a square lattice and a honeycomb lattice, this summation must be careful, due to the number of the next nearest neighbor sites becomes large. For a square lattice, the summation over the next nearest neighbor sites of the site x i is restricted as the sites x i±2ex and x i±2ey , where e x and e y are the x-axis and y-axis unit vectors, respectively, and we discard other sites, such as, x i±ex±ey and x i±ex∓ey , due to these sites can be reached from the site x i by two different ways.
In comparison with the EOMs of the transverse spin susceptibility in Eqs. (3,4), the ones in Eqs.(5,6) have a prominent character that there emerges the Γ lij (ω) term on the right hand side, which is contributed by the high order multiple-point correlation functions. In the bipartite sublattice representation, the Γ lij (ω) term represents the relation between the spin susceptibility on different sites of the same spin ingredient, which is survived in the low energy limit. While, the J U ij (ω) term describes the relation between the spin susceptibility on the nearest neighbor sites of the different spin ingredients, which is going to zero in the low energy limit. Under the condition of the locally short range antiferromagnetic correlation, the former one corresponds to the effect of twice spin-flipping process on different sites, which is a pair of kink and anti-kink in 1D, and the latter one is the effect of one spin-flipping process on the nearest neighbor sites, which is a kink in 1D. Based on these considerations, it is convincible to assume that the Γ lij (ω) term represents the low-lying excitations of magnons (pairs of spinons), and the J U ij (ω) term describes the low-lying excitations of nearly deconfined spinons. While, the magnons and nearly deconfined spinons are coexisting in the mid energy range where Γ lij (ω) and J U ij (ω) are finite. This picture is completely consistent with the exact one by the Bethe ansatz [13] in 1D.

IV. A SPIN CHAIN
For simplicity, we first consider a spin chain without the longitudinal coupling J z ij = 0, (XY model) to calculate the low-lying excitations of spins. According to the Eqs.(1,2), we have the low-lying excitation spectrum of the XY model, where choosing the lattice constant one. Based on the exact solution of the Bethe ansatz of the spin S = 1/2 chain XY model [14], the authors [15] had calculated the low-lying spectrum which is that, ε Bethe The spectrums ε XY 0k and ε Bethe k both have the same dispersion in the range of momentum, −π ≤ k ≤ π, and the difference between them is only their coefficients. The ε XY 0k is the upper boundary of spectrum of the XY model [12], and it represents the low-lying excitations of nearly deconfined spinons [13]. In order to have the lowest boundary of spectrum, like that in Eqs. (5,6), it needs to calculate the contributions coming from the high order multiple-point correlation functions to the spin susceptibility χ zz iq (ω) and χ zz iq (ω). After including the contributions of the multiplepoint correlation functions belonging to N = 2 level (see Appendix B), we obtain the lowest boundary of spectrum of the XY model, which is the low-lying excitations (spin wave) of magnons, and it possesses double periodicity of | sin (k) |, like the exact one [12,13,15].
In the case of the isotropic couplings J ⊥ = J z = J, with the Eqs.(1-4), we can obtain the following low-lying excitation spectrums, Here ε L/U 0k represent the lowest-and upper-boundary of the low-lying excitations only considering the contributions of the high order related multiple-point correlation functions belonging to N = 1 level. Notice that the low-lying excitation spectrum ε L 0k will disappear without using the bipartite sublattice representation, thus it may be an artificial result produced by the bipartite sublattice representation.
The spectrum ε U 0k has the same dispersion on the whole range of k as that of the S = 1 low-lying excitation of the Bethe ansatz, ε Bethe U k = πJ| sin k 2 |, while the spectrum ε L 0k is different from another one of the S = 1 low-lying excitation spectrum of the Bethe ansatz, After including the contributions of the multiple-point correlation functions belonging to N = 2 level, the momentum dependence in 0 ≤ k ≤ π of the upper boundary of the low-lying excitation spectrum is the same as that of ε U 0k , only their coefficients are modified, which describes the low-lying excitations of nearly deconfined spinons; However, the lowest boundary of the low-lying excitations is heavily modified, and it possesses double periodicity of | sin (k) |, like that for the XY model, and it has the same momentum dependence in 0 ≤ k ≤ π as the exact one ε Bethe Lk (see below).
With the help of the Eqs.(5,6) that included the contribution of the high order related multiple-point correlation functions belonging to N = 2 level, we obtain the following equations that can be used to determine the low-lying excitation spectrums of the system, where ∆ (ω) = J 2 + 4Γ (ω) + Here the Eq.(13) will disappear without using the bipartite sublattice representation. Thus, the physical low-lying excitations of the system is determined by the Eq.(12) which is independent of the bipartite sublattice representation. The momentum dependence of the low-lying excitations is controlled by the factors, 1 − cos (k) and 1 − cos (2k), that show different behavior around k ∼ 0 and k ∼ π. In the regime around k ∼ 0, the factors 1 − cos (k) and 1 − cos (2k) both go zero, and the possible excitation region of the spins is narrow. However, in the momentum regime around k ∼ π, there emerges a broad region of the low-lying excitation in the frequency and momentum plane, where its upper and lowest boundaries can be determined by the Eq.(12), since the factor 1 − cos (k) goes to 2, while the factors 1 − cos (2k) goes to zero.
In the high energy limit, ω/J → ∞, the coefficient ∆ (ω) takes the value, ∆ (ω) = J 2 , and the coefficient Γ (ω) goes to zero, Γ (ω) = 0, in which we can have the upper boundary of the low-lying excitation. In the low energy limit, ω/J → 0, the coefficient ∆ (ω) goes to zero, ∆ (0) = 0, and the coefficient Γ (ω) becomes a constant, Γ (0) = −J 2 /16, where the lowest boundary of the low-lying excitations can be determined. Under these two limits, we have the lowest and the upper boundaries of the low-lying excitations, Comparing with the Eq.(11), we find that the high order multiple-point correlation functions belonging to N = 2 level determines the lowest boundary ε L k of the low-lying excitation (spin wave) of the magnons, and they have little influence on the upper boundary ε U k of the lowlying excitations of the nearly deconfined spinons. In the whole region of the momentum, 0 < k < π, the lowest boundary ε L k and the upper boundary ε U k of the low-lying excitations both have the same dispersion as the exact ones ε Bethe Lk and ε Bethe U k of the Bethe ansatz. Moreover, there emerges a broad mixed region in the frequency ω and momentum k between ε L k and ε U k of the low-lying excitations of the magnons and nearly deconfined spinons determined by the Eq.(12).

V. THE SQUARE LATTICE
For the square lattice, according to the Eqs. (5,6), the low-lying excitation spectrums of the spins are determined by the following equations, where ζ k = cos k x +cos k y , and η k = cos 2k x +cos 2k y . These equations are similar to that ones for a chain, and the last equation (16) will disappear without using the bipartite sublattice representation. Thus the physical low-lying excitations of the system are determined by the Eq.(15) which is independent of the bipartite sublattice representation.
Like that for a chain, the low-lying excitation spectrums is determined by the Eq.(15), in which there are two kinds of the low-lying excitations represented by the factors 2 − η k , and 2 − ζ k , respectively. Around the momentum, k = (0, 0), the possible region of the low-lying excitations in the frequency ω axis is narrow, since both the factors 2 − ζ k and 2 − η k go to zero.
Thus we obtain the lowest boundary of the low-lying excitations, that represents the spin wave excitation of magnons (paired spinons). In the broad mixed region between ε L (k) and ε U (k), there are two kinds of modes of the low-lying excitations, that are represented by the factors, 2−ζ k and 2−η k , respectively. In practice, the coefficients ε L (k) and ε U (k) are modified due to both the coefficients ∆ (ω) and Γ (ω) are the function of the frequency ω.
According to the Eq.(15), the lowest boundary of the low-lying excitations around k = (π, 0) or k = (0, π) is generally different from that ones around k = π 2 , π 2 , even though they have the same upper boundary of the low-lying excitations represented by ε U (k). In the region around k = (π, 0) or k = (0, π), these two modes of the low-lying excitations have slowly varying momentum dependence of the forms cos ∆k x ± cos ∆k y and cos 2∆k x ± cos 2∆k y , where ∆k x and ∆k y are small quantities away from the point k = (π, 0) or k = (0, π). While in the region around k = π 2 , π 2 , these low-lying excitations have the momentum dependence of the forms sin ∆k x + sin ∆k y and cos 2∆k x + cos 2∆k y . The lowlying excitations in these two different regions have a distinct symmetry about ∆k x and ∆k y . On the other hand, in these two different regions, the difference between the upper and the lowest boundaries of the low-lying excitation is much less than that in the region around k = (π, π). These prominent characters of the low-lying excitations in the short wave-length regions of the Brillouin zone have been observed in recent neutron scattering observations [17] and numerical calculations [23] where this phenomenon is explained as nearly deconfined multispinon excitations.
The EOMs of spin susceptibility in the Eqs.(5,6) for a spin 1/2 Heisenberg model are independent of the dimensions of the system, and they valid for a chain and a square lattice.
According to the explanation of the elementary excitations for a 1D spin 1/2 Heisenberg model [13], it is convincible to believe that the lowest boundary of the low-lying excitations corresponds to the spin wave excitation of magnons (paired spinons), and the upper boundary of the low-lying excitations describes the excitations of nearly deconfied spinons, while in the broad mixed region between the lowest and upper boundaries the low-lying excitations there exists a cross-over from magnons (paired spinons) close to the lowest boundary to nearly deconfined spinons near the upper boundary. It is reasonable to conjecture that the coupling strength between two spinons is decreased as the frequency ω increasing from the lowest boundary to upper boundary of the low-lying excitations. This explanation of the low-lying excitations of the system is consistent with the experimental observations [17] and numerical calculations [23].
In the low energy limit, ω/J → 0, according to the Eqs. (5,6), we obtain the following equation that can be used to determine the lowest boundary of the low-lying excitations, The lowest boundary of excitations corresponds the usual spin wave which mainly takes place in the region around k = (π, π), while in other regions, such as around k = (π, 0) and k = π 2 , π 2 , it is very weak. Usually it is seen as the low-lying excitations of magnons (paired spinons). Another one mode resides mainly in the higher energy and/or short wave-length regions, such as around k = (π, 0) and k = π 2 , π 2 . It is called the low-lying excitations of nearly deconfined spinons. In the mixed region between the lowest and upper boundaries of the low-lying excitations there is a cross-over of the elementary excitations from magnons close to the lowest boundary to nearly deconfined spinons near the upper boundary. For a honeycomb lattice, the upper and the lowest boundaries of the low lying excitations can be described by one mode of the low-lying excitation which is represented by the factor ξ k , where the upper boundary of spectrum is gapful in whole BZ, and the lowest boundary of spectrum has the zero points at the high symmetry points k = (0, 0) and k = 2π 3 ±1, ± √ 3 of the honeycomb lattice.
All of these prominent characters of the low-lying excitation are independent of whether the system has a long range order that may modify the spectral weight of the low-lying excitations of the system, and they are completely determined by the local SU(2) symmetry of the spin 1/2 antiferromagnetic Heisenberg model and the structure of the lattice of the spins residing in, such as, a square lattice, a honeycomb lattice, or others. However, the parameter < s z i > appearing in the equation of motion of the spin susceptibility can be used to judge whether the system has a long range order. The recent neutron scattering observations [17] have clearly shown that in the regime around k = (π, π) the spectral function of the lowlying excitations has anomalously broad peaks, and this exotic behaviour of the low-lying excitations can be reasonably explained by the present calculations. Moreover, the present results obey the Lieb-Schultz-Mattis theorem and its generalizations, and for a chain, the upper and the lowest boundaries of the low-lying excitation have the same dispersion as that ones of the Bethe ansatz, only their coefficients are different, which can be modified by including the contributions of the high order multiple-point correlation functions.
Just as shown [13] for a chain, the elementary low-lying excitations are spinons, and they carry fractional spin (S = 1/2) which restricts them to being created in (multiple) pairs. The lowest boundary of the low-lying excitations is the spectrum of the two spinon excitation with the spin S = 1 (magnon). However, the upper boundary of the low-lying excitation has a different dispersion from that of this excitation spectrum of the two spinon excitation, where its period is twice that of the latter one. It is well known that the upper boundary of the low-lying excitations is the spectrum of nearly deconfined multispinon (pairs) excitations [18][19][20][21][22]. For a square lattice, it has the similar excitation spectrums like that for a chain, where the lowest boundary of the low-lying excitation around k = (π, π) represents the spin wave excitation represented by the factor 2 − η k , while as away from this region, the spin wave excitation is strongly suppressed, and it becomes very weak, since the coefficient of the factor 2 − η k rapidly decreasing as the frequency ω increasing. In a broad higher energy region between the upper and the lowest boundaries around the momentum k = (π, π) and other short wave-length regions, where the spin wave excitation nearly disappears, there emerges another low-lying excitation represented by the factor 2 − ζ k , where its coefficient rapidly increasing as the frequency ω increasing. that mainly contributes to the spectral weight of the low-lying excitations in these regions. These our results are completely consistent with the neutron scattering observations and numerical calculations for the square lattice. For the honeycomb lattice, the upper and the lowest boundaries of the low lying excitations of the spins can be described by one mode represented by the factor ξ k , that can be tested in the future neutron scattering experiments.
where the antiferromagnetic exchanges J ⊥ ij (here we take J x ij = J y ij = J ⊥ ij ) and J z ij are restricted to nearest neighbor spins < i, j >, J ⊥ ij = J ⊥ , J z ij = J z , other cases, they are zero, and for the isotropic case, J ⊥ = J z = J. The spin operators s i and τ i satisfy the commutation relations ( = 1), where δ ij is the Kronecker delta function, and ǫ µνλ is an antisymmetry tensor, ǫ xyz = 1.
With the Heisenberg representation, the time dependence of the spin operators is that, where the operator o i = s i , τ i .
According to the algebraic equation of motion approach, we need the following commutation relations, where X These commutation relations are the basic ingredients to writing equations of motion of high order correlation functions.
The transverse and longitudinal spin susceptibilities are defined as that, where µ, ν = ±, z. In order to tersely represent the equation of motion of the spin susceptibility, we define the following correlation functions that some of them appearing in the hierarchic series of EOM of the spin susceptibility, called related multiple-point correlation functions, ij , and the parameter N represents the level of the corresponding correlation function in the hierarchic structure of the EOM of the spin susceptibility. For the longitudinal spin susceptibility χ zz ij (t) and χ zz ij (t), we need to define another special related multiple-point correlation function, K ijq (t) = iθ(t) < X With the commutation relations in Eqs (24,25), we can write out the following the EOMs of the spin susceptibility after taking the Fourier transformation of time, ωχ −+ that are used to calculated the spin susceptibility of the spin 1/2 XY model.
With the help of the Eqs. (27,28), we can write out the EOMs of the multiple-point where we have neglected the static quantities appearing in these EOMs.
As writing out the summation over the lattice sites in the right side of the above EOMs, there may appear some multiple-point correlation functions belonging to the N = 3 level that have two same labels, such as, to simplify these multiplepoint correlation functions where there emerge some ones belonging to the N = 1 level, and finally discarding multiple-point correlation functions belonging to the N = 3 level [29] (called a "soft cut-off" approximation), we can rewrite out the Eqs.(42,43) as that, where there appear three spin susceptibilities in the above solutions defined on three neighbor sites.
where we have neglected the static quantities appearing in these EOMs.
Under the "soft cut-off" approximation that applying the relations ( s z for the correlation functions that having two same labels appearing in the summations of the right hand side of the Eqs.(47-52), and discarding the multiple-point correlation functions belonging to the N = 3 level, we can further simplify these equations as that, where A lij = (J il ) 2 2 (2 − δ jl ), B lij = (J ij ) 2 +(J jl) (1 − δ jl ).
The EOMs of the correlation functions F (X (−) ) mijq (ω) and F (X (−) ) imjq (ω) can be directly solved, and their solutions are written as that, where D X (ω) = (ω 2 − J 2 ) 2 − J 4 4 . The Eqs.(55-58) are a set of coupled equations that can be approximately solved. As a zeroth order approximation, we first decouple these equations by discarding the correlation functions F lijq (ω) appearing in the right hand side of these equations, in which they become independent with each other, and we can straight solve them. Then we use these approximation solutions to replace them that appearing in the right hand side of other equations, respectively. Under these approximations, we finally obtain the following solutions of the correlation functions F where Γ lij (ω) = These solutions of the multiple-point correlation functions F