text:
# A New Method of the High Temperature Series Expansion

## I Introduction

The high-temperature expansion can be used in any dimension and has provided significant information on variety of models.Domb3  It is based on the Taylor expansion of the Boltzmann factor $`e^\beta `$ in $`\beta `$, around the high-temperature limit. Although the concept of the high-temperature expansion is quite simple, its high-order calculation requires sophisticated strategies. A standard method of the high-temperature expansion is called the finite cluster methodDomb3 ; oitmaa ; Gelfand  (FCM). It has been used for a long time,Domb3  and progress of the high-temperature series has been due mainly to improvement of computers. In the FCM, the series expansion is reduced to calculation in connected finite-size clusters. However, the FCM is inefficient for calculating a correlation function between rather distant sites: (i) The clusters have to include those two sites and calculation in such large clusters consumes a lot of time. (ii) The FCM does not use a valuable piece of information. That is, the high-temperature series of a long-range correlation starts at a certain finite order, namely, low-order coefficients are equal to zero —– physically, it means that a long-range correlation develops at low temperature.

In this paper, we formulate a new method oriented to the long-range correlation, namely, it overcome (i) and (ii) above. The new method is valid for the spin-1/2 Heisenberg model, or, more generally, for SU($`n`$) Heisenberg model with arbitrary $`n`$. High-temperature series coefficients are written in terms of cumulants. We start from the lowest-order nonzero cumulant in the correlation-function series. Next, we consider a deviation from it. Then, the deviation can be regarded as a sort of moment, which we call a quasi-moment. The corresponding quasi-cumulant enables us to calculate a number of terms higher than the lowest-order nonzero cumulant. In Ref. FK1d, , the Fourier transform of the correlation function is calculated using the new method and the FCM complementarily. That is, we have used the new method for the long-range correlations and the FCM for the short-range correlations.

In addition, also for the specific heat, we have adopted a similar strategy in Ref. FK1d, . That is, contributions from the required largest cluster are calculated by a new technique, and the FCM is used only for the smaller clusters. We utilize cumulants relevant to the ordering of quantum variables. By choosing nonzero cumulants from them, some of the series coefficients are calculated simply. We explain also this new technique in this paper.

The structure of this paper is as follows. In Sec. II, we present the model. We show a representation of SU($`n`$) Heisenberg model in terms of exchange operators. The relation of a correlation function and an expectation value of the exchange operator is mentioned. In Sec. III, we review the high-temperature expansion. Although we use this model for illustration, what is written here can generally be applied to other models. To clarify the purpose of our new method, we review the FCM in detail using a one-dimensional system as an explicit example. In Sec. IV, we describe the new method for the correlation function. Here, we derive some properties of the cumulants. Owing to those properties, we can define quasi-moments/cumulants. How to calculate series coefficients using those quasi-cumulants is explained. In Sec. V, we show the new technique for the specific heat. Finally, we give a summary in Sec. VI.

## II model

Let us consider the spin-1/2 Heisenberg model,

$$_{\mathrm{SU}(2)}^{(p)}:=2\underset{(i,j)}{}J_{i,j}𝒔_i𝒔_j,$$
(1)
where $`𝒔_i`$ is the spin operator at site $`i`$, and $`J_{i,j}`$ is the coupling constant. Here, $`(i,j)`$ represents that the summation is performed only once for each $`(i,j)`$ pair. The spin operators are rewritten as

$$2𝒔_i𝒔_j=P_{i,j}\frac{1}{2},$$
(2)
where $`P_{i,j}`$ is an exchange operator defined by

$`P_{i,j}|\mathrm{}\stackrel{i}{\alpha }\mathrm{}\stackrel{j}{\beta }\mathrm{}=|\mathrm{}\stackrel{i}{\beta }\mathrm{}\stackrel{j}{\alpha }\mathrm{},`$ (3)
for arbitrary $`\alpha =,`$, and $`\beta =,`$. Then, we can rewrite Eq. (1) as

$$_{\mathrm{SU}(2)}^{(p)}=\underset{(i,j)}{}J_{i,j}\left(P_{i,j}\frac{1}{2}\right).$$
(4)
This Hamiltonian has SU(2) symmetry. Since our method is valid more generally for an SU($`n`$) symmetric case with arbitrary positive integer $`n`$, we define the SU($`n`$) Heisenberg model below.

Let each site take one of the $`n`$ colors, and denote them as $`|\alpha `$ with $`\alpha =1,2,\mathrm{},n`$. An exchange operator is defined using $`X^{\alpha \beta }:=|\alpha \beta |`$ by

$$P_{i,j}:=\underset{\alpha =1}{\overset{n}{}}\underset{\beta =1}{\overset{n}{}}X_i^{\alpha \beta }X_j^{\beta \alpha },$$
(5)
for $`ij`$. Colors at sites $`i`$ and $`j`$ are exchanged by $`P_{i,j}`$, namely, Eq. (3) is satisfied for arbitrary colors $`\alpha `$ and $`\beta `$. Furthermore, for a later convenience, let us define

$$p_{i,j}:=P_{i,j}\frac{1}{n}.$$
(6)
In the case of SU(2),

$$p_{i,j}=2𝒔_i𝒔_j$$
(7)
is satisfied. Although we mostly use expressions for general $`n`$ in this paper, we expect that Eq. (7) helps instant understanding of explanations. In some cases, we explicitly use the SU(2) notation for a simple explanation.

The Hamiltonian of the SU($`n`$) Heisenberg model is given by

$$^{(p)}:=\underset{(i,j)}{}J_{i,j}\left(P_{i,j}\frac{1}{n}\right)=\underset{(i,j)}{}J_{i,j}p_{i,j}.$$
(8)
The constant term just yields an energy shift. Hence, for calculating thermal averages, we can use

$$:=\underset{(i,j)}{}J_{i,j}P_{i,j}.$$
(9)
We define the partition function as

$$Z:=\mathrm{Tr}e^\beta ,$$
(10)
and the free energy as

$$F:=\beta ^1\mathrm{ln}Z.$$
(11)
The average in this system is denoted by

$$\widehat{O}:=\mathrm{Tr}\left(\widehat{O}e^\beta \right)/Z,$$
(12)
where $`\widehat{O}`$ is an arbitrary operator.

In the case of spin-1/2 system, there are three independent interacting components, $`s^x,s^y,s^z`$. In general, when there are $`n`$ states per site, the maximum number of independent interacting components is $`n^21`$. The model, isotropic with respect to these components, is the SU($`n`$) Heisenberg model as explicitly shown in Appendix A.

What we present in this paper is the high temperature expansion of $`P_{i,j}`$. A correlation function $`X_i^{\alpha \beta }X_j^{\beta \alpha }`$ is calculated using a relation,

$`X_i^{\alpha \beta }X_j^{\beta \alpha }`$ $`=`$ $`{\displaystyle \frac{1}{n^21}}\left(P_{i,j}{\displaystyle \frac{1}{n}}\right),`$ (13)
for arbitrary $`\alpha \beta `$, $`ij`$. In addition, when $`i=j`$,

$$X_i^{\alpha \beta }X_i^{\beta \alpha }=\frac{1}{n},$$
(14)
for arbitrary $`\alpha \beta `$. These relations are shown also in Appendix A. In the case of SU(2), $`X_i^{\alpha \beta }X_j^{\beta \alpha }=2s_i^zs_j^z`$. The uniform susceptibility is equal to inverse temperature times the Fourier transform of the correlation function with wave-number zero. Furthermore, the specific heat is calculated from the internal energy $`=_{(i,j)}J_{i,j}P_{i,j}`$.

## III high-temperature expansion

In this section, we review high-temperature expansion. Although the Heisenberg model is used for explanation, this review is more general. That is, one can easily replace the Heisenberg model with a different model.

### III.1 Moments and cumulants for classical variables

Fundamental properties of moments and cumulants are frequently used in this paper. In order to conveniently refer to those properties, let us review moments and cumulants of classical variables $`x_i(i=1,2,\mathrm{})`$. We put most of the fundamental details of the review in Appendix B, and let us here just show one important formula used frequently in this paper. Here, a moment is denoted by $`x_jx_k\mathrm{}x_l_x`$, and a cumulant is denoted by $`[x_jx_k\mathrm{}x_l]_x`$.

Moments can be expanded using cumulants, and vice versa. However, in those relations, the number of terms in the expansion drastically increases as the order of moments or cumulants becomes higher. To avoid it, another relation,

$`[x_i\mathrm{}x_{\mathrm{}}]_x`$ $`=`$ $`x_i\mathrm{}x_{\mathrm{}}_x`$ (15)
$`{\displaystyle \underset{𝒫(\xi ,2)}{}}[x_i\mathrm{}x_j]_xx_k\mathrm{}x_m_x,`$
can be used.Domb3  Here, $`\xi `$ is the number of $`x`$-variables in the bracket in the l.h.s. Then, the summation denoted by $`𝒫(\xi ,2)`$ is taken over every partition of $`\xi `$ elements into two groups on the condition that one of the variables, for example $`x_i`$, must always be included in $`[\mathrm{}]_x`$. In addition, each bracket includes at least one $`x`$-variable. In this paper, we call this formula “the mixed expansion pivoting on $`x_i`$”. For example, the mixed expansion of $`[x_1x_2x_3]_x`$ pivoting on $`x_1`$ is written as

$`[x_1x_2x_3]_x`$ $`=`$ $`x_1x_2x_3_x[x_1x_2]_xx_3_x[x_1x_3]_xx_2_x`$
$`[x_1]_xx_2x_3_x.`$
### III.2 Moments and cumulants for quantum variables

As a reference system, we take a non-interacting system such that $`J_{i,j}=0`$ for every $`(i,j)`$ pair. Its partition function is written as $`Z_0:=\mathrm{Tr}\mathrm{\hspace{0.17em}1}`$, and the average in the non-interacting system is denoted by

$$\widehat{O}_0:=\frac{\mathrm{Tr}(\widehat{O})}{\mathrm{Tr}\mathrm{\hspace{0.17em}1}},$$
(16)
where $`\widehat{O}`$ is an arbitrary operator. This average $`\mathrm{}_0`$ plays a role of $`\mathrm{}_x`$ in the previous subsection.

There are a number of ways to define moments and cumulants of quantum variables because of their noncommutativity.Kubo  One of the ways is defining a moment by a symmetrized product,

$`P_{i_1,j_1}P_{i_2,j_2}\mathrm{}P_{i_k,j_k}_\mathrm{s}`$
$`:={\displaystyle \frac{1}{k!}}{\displaystyle \underset{\sigma }{}}P_{i_{\sigma (1)},j_{\sigma (1)}}P_{i_{\sigma (2)},j_{\sigma (2)}}\mathrm{}P_{i_{\sigma (k)},j_{\sigma (k)}}_0,`$ (17)
where $`\sigma `$ represents a permutation of the indeces. Note that the ordering of exchange operators $`P_{i,j}`$ is unimportant in $`\mathrm{}_\mathrm{s}`$. Hereafter, we use $`\mathrm{}_0`$ rather than $`\mathrm{}_\mathrm{s}`$ if they are obviously equivalent to each other, e.g., $`P_{i,j}_\mathrm{s}=P_{i,j}_0`$; further examples are commented in Appendix C.

The partition function, Eq. (10), is rewritten as

$$Z=\mathrm{Tr}\mathrm{\hspace{0.17em}1}\frac{\mathrm{Tr}e^\beta }{\mathrm{Tr}\mathrm{\hspace{0.17em}1}}=Z_0e^\beta _0.$$
(18)
This averaged Boltzmann factor works as the generating function of the symmetrized moments mentioned above, namely, with $`\lambda _{i,j}=\beta J_{i,j}`$,

$$P_{i,j}\mathrm{}P_{k,l}_\mathrm{s}=\frac{}{\lambda _{i,j}}\mathrm{}\frac{}{\lambda _{k,l}}e^\beta _0|_{\lambda =0}.$$
(19)
The corresponding cumulants are denoted by $`[\mathrm{}]_\mathrm{s}`$, and given by

$$[P_{i,j}\mathrm{}P_{k,l}]_\mathrm{s}=\frac{}{\lambda _{i,j}}\mathrm{}\frac{}{\lambda _{k,l}}\mathrm{ln}e^\beta _0|_{\lambda =0}.$$
(20)
These moments and cumulants have the properties explained in the previous subsection and Appendix B.

The free energy is rewritten as

$$F=\frac{1}{\beta }\mathrm{ln}Z_0\frac{1}{\beta }\underset{m=1}{\overset{\mathrm{}}{}}\frac{(\beta )^m}{m!}[^m]_\mathrm{s},$$
(21)
$`[^m]_\mathrm{s}`$ $`=`$ $`{\displaystyle \underset{(i_1,j_1)}{}}\mathrm{}{\displaystyle \underset{(i_m,j_m)}{}}J_{i_1,j_1}\mathrm{}J_{i_m,j_m}`$ (22)
$`\times \left[P_{i_1,j_1}\mathrm{}P_{i_m,j_m}\right]_\mathrm{s}.`$
Many terms in the summation of Eq. (22) are equivalent because the ordering of operators in $`[\mathrm{}]_\mathrm{s}`$ is unimportant. That is, in the case of the cumulant $`[P_{1,2}^{k_{12}}P_{2,3}^{k_{23}}\mathrm{}]_\mathrm{s}`$, then, $`m!/(k_{12}!k_{23}!\mathrm{})`$ terms are equivalent.

From Eqs. (21) and (22), we can derive

$`P_{i,j}={\displaystyle \frac{F}{J_{i,j}}}`$
$`={\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\beta )^{m1}}{(m1)!}}{\displaystyle \underset{(i_1,j_1)}{}}\mathrm{}{\displaystyle \underset{(i_{m1},j_{m1})}{}}`$
$`J_{i_1,j_1}\mathrm{}J_{i_{m1},j_{m1}}\left[P_{i,j}P_{i_1,j_1}\mathrm{}P_{i_{m1},j_{m1}}\right]_\mathrm{s}`$
$`={\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\beta )^m}{m!}}\left[P_{i,j}^m\right]_\mathrm{s}.`$ (23)
Even if $`J_{i,j}=0`$ in the original Hamiltonian (e.g., for non-nearest neighbors in the case of nearest-neighbor interaction), one can start from the general formulation written above. After deriving all the formulae, one can put $`J_{i,j}=J`$ or 0 according to the original Hamiltonian.

Let us consider cumulants for $`p_{i,j}`$ instead of $`P_{i,j}`$. See the definition of cumulants, Eq. (20). If one replaces $`P_{i,j}`$ by $`p_{i,j}`$, the generating function of cumulants changes only by $`1/n_{(i,j)}\lambda _{i,j}`$, and thus only the first-order cumulants are affected by this replacement. In other words,

$$[P_{i,j}]_\mathrm{s}=P_{i,j}_0=1/n$$
(24)
differs from

$$[p_{i,j}]_\mathrm{s}=p_{i,j}_0=0.$$
(25)
On the other hand, the higher order cumulants are equal. Namely,

$$[p_{i_1,j_1}p_{i_2,j_2}\mathrm{}p_{i_{\mathrm{}},j_{\mathrm{}}}]_\mathrm{s}=[P_{i_1,j_1}P_{i_2,j_2}\mathrm{}P_{i_{\mathrm{}},j_{\mathrm{}}}]_\mathrm{s}$$
(26)
when $`[\mathrm{}]_\mathrm{s}`$ has two or more $`p_{i,j}`$ operators. For practical calculation, we mainly use $`P_{i,j}`$ for convenience. In some cases, however, the property $`p_{i,j}_0=0`$ helps simplification. Hereafter, we use $`P_{i,j}`$ and $`p_{i,j}`$ interchangeably with keeping in mind Eqs. (24)-(26).

### III.3 Properties of the cumulants

Let us introduce a diagrammatic representation of the moments and the cumulants. Although we mainly use those diagrams for later calculation, let us define them here in advance because they are convenient also for simple explanation in this review.

Figure 1 shows examples of the diagrams for $`P_{i,j}_\mathrm{s}`$ and $`[P_{i,j}]_\mathrm{s}`$. We call them bond-diagrams. Each dot represents a site. A segment connecting sites $`i`$ and $`j`$ represents $`P_{i,j}`$, and we call it a bond. Note that a product of two equivalent operators appearing in $`\mathrm{}_\mathrm{s}`$ \[e.g., $`P_{3,4}P_{3,4}`$ in the r.h.s. of Fig. 1(a)\] cannot be reduced to the identity because these are not always next to each other in the symmetrization. Figure 1(a) is a moment bond-diagram, and (b) is a cumulant bond-diagram.

A cumulant is equal to zero if the variables in $`[\mathrm{}]_\mathrm{s}`$ can be partitioned into two groups which are independent of each other on averaging $`\mathrm{}_0`$. For example, $`[P_{1,2}P_{3,4}]_\mathrm{s}=0`$ because the trace of sites 1,2 is taken independently of that of sites 3,4 in the non-interacting system. In other words, a cumulant is equal to zero if the bond-diagram is not connected as shown in Fig. 2(a). Here, a solid rectangle in Fig. 2 represents arbitrary bonds. Therefore, as for cumulants, we consider only connected diagrams hereafter.

Figure 2(b) is a first-order cumulant. In a higher-order cumulant bond-diagram, more than one bonds must meet at each dot. If not, it is equal to zero as shown in Fig. 2(c) and (d). Furthermore, a cumulant bond-diagram is equal to zero if two parts in it are connected by only one bond as shown in Fig. 2(e). We prove these properties in Appendix D.

### III.4 An example of the finite cluster method

The finite cluster method (FCM) is a standard method for the high-temperature expansion.Domb3 ; oitmaa  It is also used for perturbative expansions starting from exactly solvable models.Gelfand  In order to explain it in detail, we give an explicit example using the one-dimensional nearest-neighbor interacting system with a periodic boundary condition, namely, $`J_{i,j}=J\delta _{j,i+1}`$, $`(1iN)`$, where the site $`N+1`$ is equivalent to site 1. In this case, the procedure of the FCM is reduced to a very simple form.

In the FCM, series coefficients in the thermodynamic limit are exactly obtained by summing contributions from finite-size clusters. The clusters are subsets of the infinite lattice, and we need only connected clusters because cumulant bond-diagrams must be “connected”. Namely, the one dimensional system requires only open chains. The cluster of size $`\mathrm{}`$ is defined as

$$_{\mathrm{}}^{(p)}:=J\underset{i=1}{\overset{\mathrm{}1}{}}\left(P_{i,i+1}\frac{1}{n}\right)=J\underset{i=1}{\overset{\mathrm{}1}{}}p_{i,i+1},$$
(27)
where we have chosen to use $`p_{i,j}`$ rather than $`P_{i,j}`$ in order to make use of $`[p_{i,j}]_\mathrm{s}=0`$. The expectation value in this cluster is given by

$`p_{i,j}_{\mathrm{}}`$ $`:=`$ $`{\displaystyle \frac{\mathrm{Tr}(p_{i,j}e^{\beta _{\mathrm{}}^{(p)}})}{\mathrm{Tr}(e^{\beta _{\mathrm{}}^{(p)}})}}`$ (28)
$`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\beta )^m}{m!}}\left[p_{i,j}\left(_{\mathrm{}}^{(p)}\right)^m\right]_\mathrm{s}.`$
We assume $`i<j`$ hereafter. Note that $`p_{i,j}_{\mathrm{}}=0`$ when $`i<1`$ or $`j>\mathrm{}`$ because of a property shown in Fig. 2(c) and (d). Since a finite cluster is used instead of the original system, $`p_{i,j}_{\mathrm{}}`$ contains only a subset of terms appearing in $`p_{i,j}`$. The important point is that $`p_{i,j}_{\mathrm{}}`$ includes contributions from clusters smaller than $`\mathrm{}`$, and thus simple summation of $`p_{i,j}_{\mathrm{}}`$ over $`\mathrm{}`$ yields multiple counting of those terms. To avoid this, the FCM requires subtraction of contributions from the smaller clusters. We define a net contribution from the $`\mathrm{}`$-site cluster,

$`p_{i,i+x}_{\mathrm{}}^{}`$ $`:=`$ $`p_{i,i+x}_{\mathrm{}}`$ (29)
$`{\displaystyle \underset{\mathrm{}_1=0}{\overset{i1}{}}}{\displaystyle \stackrel{}{}_{\mathrm{}_2=0}^\mathrm{}j}p_{(i\mathrm{}_1),(i\mathrm{}_1)+x}_{\mathrm{}\mathrm{}_1\mathrm{}_2}^{},`$
where the summation $`^{}`$ excludes $`\mathrm{}_1=\mathrm{}_2=0`$. Here, $`\mathrm{}_1`$ ($`\mathrm{}_2`$) represents the number of removed sites from the side of 1 ($`\mathrm{}`$). When $`i=1`$ and $`i+x=\mathrm{}`$, this expression is reduced to $`p_{1,\mathrm{}}_{\mathrm{}}^{}=p_{1,\mathrm{}}_{\mathrm{}}`$. Let $`p_{i,i+x}_{\mathrm{}}^{}`$ denote the total net contribution to $`p_{i,i+x}`$ from clusters smaller than or equal to $`\mathrm{}`$,

$`p_{i,i+x}_{\mathrm{}}^{}:={\displaystyle \underset{l=x+1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{j=1}{\overset{lx}{}}}p_{j,j+x}_l^{}.`$ (30)
Here, $`p_{i,i+x}_{\mathrm{}}^{}`$ does not depend on $`i`$ because of the translational symmetry. Note that $`p_{i,i+x}=p_{i,i+x}_{\mathrm{}}^{}`$. When $`\mathrm{}`$ is finite, the series for $`p_{i,i+x}_{\mathrm{}}^{}`$ is correct up to a finite order. To think about the “finite order”, let us take a cumulant bond-diagram appearing in the $`m`$-th order. The diagram has a bond $`P_{i,i+x}`$ and $`m`$ bonds from the Hamiltonian. In order to give a nonzero contribution, there must be at least one bond per nearest-neighbor pair between $`i`$ and $`j`$. Furthermore, nearest-neighbor pairs not between $`i`$ and $`i+x`$ —– less than $`i`$ or greater than $`i+x`$ —– must have at least two bonds per pair. Therefore, there is a contribution from $`(\mathrm{}+1)`$-site cluster when $`mx+2(\mathrm{}x)`$. That is, $`p_{i,i+x}_{\mathrm{}}^{}`$ is correct up to $`O[(\beta J)^{2\mathrm{}x1}]`$.

The above is the usual method of the FCM. However, in one dimension, these equations can be reduced to simpler forms as shown in Appendix E.

## IV A new method for the correlation function

Let us think about calculating the Fourier transform of the correlation function,

$$S^z(q):=\frac{1}{n^21}\underset{x}{}p_{j+x,j}e^{iqx}.$$
(31)
Series of $`S^z(q)`$ up to $`O[(\beta J)^M]`$ requires series of each $`p_{j+x,j}`$ up to $`O[(\beta J)^M]`$. Here, the important point is that the series coefficients of $`p_{j+x,j}`$ up to $`O[(\beta J)^{|x|1}]`$ are equal to zero as explained in Fig. 2(c),(d) and in Sec. III.4. Hence, we need $`p_{j+x,j}`$ only for $`1xM`$. In addition, that property provides us a weak point of the FCM: When $`M`$ is fixed, large $`x`$ requires a larger cluster than small $`x`$ does. Thus, the most time-consuming part is calculation in the largest cluster, $`x=M`$, namely, calculation of Tr$`{}_{}{}^{(1+M)}\{p_{1,1+M}(_{1+M}^{(p)})^M\}`$. However, what we need for $`x=M`$ is only the $`M`$-th order coefficient because we already know that the lower orders are equal to zero. Therefore, if only the FCM is used, the most time-consuming calculation yields very little information.

Our goal in this paper to formulate a new method is as follows. If the contributions from large clusters can be calculated by another method, the FCM can be used only for smaller clusters.Singhcomment  Accordingly, one can calculate up to high orders. In Ref. FK1d, , we have calculated $`S^z(q)`$ up to $`O[(\beta J)^{19}]`$. However, we have used the FCM only for $`\mathrm{}13`$. The corrections to the results have been calculated by the method explained in this section.

### IV.1 Traces using combinatorics

We calculate traces of products of exchange operators by decomposing every permutation into a product of independent cyclic permutationshand ; Chen  as explained in the following, (an explicit example is given in Appendix F). Let us consider a trace of $`P:=P_{i_1,j_1}P_{i_2,j_2}\mathrm{}P_{i_m,j_m}`$, with $`\mathrm{Tr}^{(\mathrm{})}`$ denoting the trace in the $`\mathrm{}`$-site system,

$`\mathrm{Tr}^{(\mathrm{})}P`$ (32)
$`:=`$ $`{\displaystyle \underset{\alpha _1=1}{\overset{n}{}}}\mathrm{}{\displaystyle \underset{\alpha _{\mathrm{}}=1}{\overset{n}{}}}\stackrel{1}{\alpha _1}\stackrel{2}{\alpha _2}\mathrm{}\stackrel{\mathrm{}}{\alpha _{\mathrm{}}}|P|\stackrel{1}{\alpha _1}\stackrel{2}{\alpha _2}\mathrm{}\stackrel{\mathrm{}}{\alpha _{\mathrm{}}}`$
$`=`$ $`{\displaystyle \underset{\alpha _1=1}{\overset{n}{}}}\mathrm{}{\displaystyle \underset{\alpha _{\mathrm{}}=1}{\overset{n}{}}}\stackrel{1}{\alpha _1}\stackrel{2}{\alpha _2}\mathrm{}\stackrel{\mathrm{}}{\alpha _{\mathrm{}}}|\stackrel{1}{\alpha _{P1}}\stackrel{2}{\alpha _{P2}}\mathrm{}\stackrel{\mathrm{}}{\alpha _P\mathrm{}}`$
$`=`$ $`{\displaystyle \underset{\alpha _1=1}{\overset{n}{}}}\mathrm{}{\displaystyle \underset{\alpha _{\mathrm{}}=1}{\overset{n}{}}}\delta _{\alpha _1,\alpha _{P1}}\delta _{\alpha _2,\alpha _{P2}}\mathrm{}\delta _{\alpha _{\mathrm{}},\alpha _P\mathrm{}},`$
where $`\alpha _{Pi}`$ refers to the color at position $`i`$ after the permutation $`P`$. The summation of $`\alpha _i`$ makes a contribution only when $`\alpha _i=\alpha _{Pi}`$ for every $`i`$. Consider using this relation successively starting from $`i`$. That is, $`\alpha _i`$ is equal to $`\alpha _{Pi}`$, and then $`\alpha _{Pi}`$ is equal to $`\alpha _{P^2i}`$, …, one can repeat this procedure until coming back to $`\alpha _i`$ at a certain power of $`P`$, namely, $`\alpha _i=\alpha _{Pi}=\alpha _{P^2i}=\mathrm{}=\alpha _i`$. In other words, all the variables whose subscript belong to one cyclic permutation in $`P`$ have to be equal. Since any permutation can be decomposed into a product of independent cyclic permutations, the number of independent variables of the summation is the number $`Y(P)`$ of cyclic permutations of $`P`$. Therefore the trace is given by, $`\mathrm{Tr}^{(\mathrm{})}P=n^{Y(P)}`$, and accordingly,

$$P_{\mathrm{}}=\frac{\mathrm{Tr}^{(\mathrm{})}P}{\mathrm{Tr}^{(\mathrm{})}1}=n^{Y(P)\mathrm{}}.$$
(33)
In fact, this method allows us to calculate Tr$`{}_{}{}^{(\mathrm{})}(P_{i,j}P)`$ at the same time as Tr$`{}_{}{}^{(\mathrm{})}P`$ for an arbitrary $`(i,j)`$ pair.Chen  Let us remember the permutation $`P`$ is decomposed into cycles in calculating $`\mathrm{Tr}^{(\mathrm{})}P`$. When the cycle which $`i`$ belongs to is different from what $`j`$ belongs to, the next operation $`P_{i,j}`$ unites the two cycles into one. Then, the number of cycles decreases by one, namely, $`\mathrm{Tr}^{(\mathrm{})}(P_{i,j}P)=n^{Y(P)1}`$. On the contrary, if $`i`$ and $`j`$ belong to one cycle of $`P`$, then $`P_{i,j}`$ breaks the cycle into two. Then, the number of cycles increases by one, namely, $`\mathrm{Tr}^{(\mathrm{})}(P_{i,j}P)=n^{Y(P)+1}`$.

### IV.2 Explicit calculation of cumulants

To represent $`P_{i,j}_0`$, we introduce an unsymmetrized version of a bond-diagram. Figure 3 shows its examples. Since this diagram is equivalent to an “amida lottery”, used in Japan to peacefully decide how to distribute a fixed number of objects among an equal number of people, here we call it an amida-diagram. Each vertical line represents a site, and a horizontal line between site $`i`$ and $`j`$ represents $`P_{i,j}`$. The order of horizontal lines from bottom to top should coincide with the order of $`P_{i,j}`$ in $`\mathrm{}_0`$ from right to left. We bracket an amida-diagram by $`\mathrm{}`$ in order to distinguish it from a cumulant version of the amida-diagram introduced later in Sec. V.

Since expectation values are not changed by transpositions of neighboring commutative operators, the ordering of $`P_{1,2}`$ and $`P_{3,4}`$ is irrelevant to the results. On the other hand, two exchange operators sharing one of the site-indeces, such as $`P_{1,2}`$ and $`P_{2,3}`$, are not commutative, and in general the ordering of those bonds affects the result. For example, the value of Fig. 3(a) is different from that of (b). However, in some cases, even the order of such “site-sharing” bonds does not affect the result. Let us think about the value of an amida-diagram in Fig. 4(a) by counting cycles. The “initial” state is the identity permutation, which has $`\mathrm{}`$ cycles. Then, every time $`P_{i,j}`$ is applied, it unites two cycles into one. After all the $`P_{i,j}`$ operators are applied, there is only one cycle, which includes all the $`\mathrm{}`$ sites. Hence the expectation value is $`n^1\mathrm{}`$, and it does not depend on the ordering of those $`P_{i,j}`$ operators. Therefore, after the symmetrization, the moment bond-diagram Fig. 4(b) also gives $`n^1\mathrm{}`$.

Let us think about adding $`P_{1,\mathrm{}}`$ to Fig. 4(a) on the top. Since sites 1 and $`\mathrm{}`$ belong to the same cycle, the operation of $`P_{1,\mathrm{}}`$ cuts the cycle. Namely, the permutation has two cycles and the diagram gives $`n^2\mathrm{}`$. Again, after the symmetrization, the corresponding moment bond-diagram also gives $`n^2\mathrm{}`$.

In Fig. 5(a), we calculate a cumulant using the values derived above. This is the lowest-order nonzero cumulant in the $`p_{1,\mathrm{}}`$ series. We use the mixed expansion, Eq. (15), pivoting on $`P_{1,\mathrm{}}`$. The expression of the expansion is very simple because $`[P_{1,\mathrm{}}]_\mathrm{s}`$ is the one and only nonzero cumulant which includes $`P_{1,\mathrm{}}`$.

Next, let us consider adding one more bond, $`P_{1,2}`$. The mixed expansion pivoting on $`P_{1,\mathrm{}}`$ is shown in Fig. 5(b). The coefficient 2 of the last term is due to the two ways of choosing which one of the two $`P_{1,2}`$ is put in $`[\mathrm{}]_\mathrm{s}`$. The first term of the r.h.s. is calculated as shown in Fig. 6. The symmetrization is equivalent to averaging over all the possible configurations of bonds. In fact, however, the value of each amida-diagram is determined only by relative configurations of two $`P_{1,2}`$ and one $`P_{2,3}`$. There are only two patterns, namely, Fig. 6(b) and (c). Their r.h.s. are obtained by simplification $`P_{1,2}P_{1,2}=1`$ and $`P_{1,2}P_{2,3}P_{1,2}=P_{1,3}`$. Here, an “$`\mathrm{\Omega }`$-shape” bond in Fig. 6(c) represents that the bond is not connected to the site at “$`\mathrm{\Omega }`$”. We already know how to calculate the r.h.s. of Fig. 6(b) and (c). We also know that configuration of $`P_{i,i+1}`$ ($`i3`$) does not affect the value of the amida-diagram. That is, the value of the amida-diagram is determined only by the position of $`P_{2,3}`$.

In order to obtain the l.h.s. of Fig. 6(a), we have to count how often each of the configurations Fig. 6(b) and (c) appears in the symmetrization. For example, the prefactor 1/3 in the last term of Fig. 6(a) is obtained by

$$\frac{\{\text{the number of configurations of Fig. }\text{6}\text{(c)}\}}{\{\text{the total number of possible configurations}\}}.$$
(34)
It can be considered like this: The vertical line of site 2 is partitioned into three regions by the two $`P_{1,2}`$. The probability of finding $`P_{2,3}`$ in the intermediate region is 1/3. More explicitly, it can be given also by integrals,

$$\frac{_0^1d\tau _{12}_{\tau _{12}}^1d\tau _{12}^{}_{\tau _{12}}^{\tau _{12}^{}}d\tau _{23}}{_0^1d\tau _{12}_{\tau _{12}}^1d\tau _{12}^{}_0^1d\tau _{23}}=\frac{1}{3},$$
(35)
where $`\tau _{12}`$, $`\tau _{12}^{}`$, $`\tau _{23}`$ represent positions of lower $`P_{1,2}`$, upper $`P_{1,2}`$, $`P_{2,3}`$, respectively. Here, $`\tau =0`$ is the bottom, $`\tau =1`$ is the top. On the other hand, the “probability” of the other configuration Fig. 6(b) is 2/3. Consequently, the symmetrized value is calculated as shown in Fig. 6(a). This formula is true only when $`\mathrm{}3`$. If $`\mathrm{}=2`$, we cannot take an average over position of $`P_{2,3}`$, and Fig. 6(c) does not appear.

Examples of higher order results are shown in Fig. 7. We can find simple properties in these cumulants.

1. A cumulant is reduced by a factor $`1/n`$ as $`\mathrm{}`$ increases by one.
2. A double bond in Fig. 7(a) yields a factor of $`(\frac{n}{3}\frac{2}{n})`$, two distant double bonds in Fig. 7(d) yield square of this factor. In fact, we can make a more general statement: When multiple bonds are distant, we can simply make a product of contributions from each of multiple bonds.

We more rigorously derive property (i) in the next subsection and property (ii) in the subsection after the next. After that, we formulate the new method using those properties.

### IV.3 The “contractible” property

First, we consider a contraction of an amida-diagram. In an amida-diagram, let us focus on a part that includes two sites as shown in Fig. 8. Because of a property of a trace, the ordering of operators can be cyclically rotated. In other words, amida-diagrams have a periodic boundary condition. Hence we can put a bond $`P_{j,k}`$ at the bottom of the amida-diagram without loss of generality. Let us compare the l.h.s. and the r.h.s. at each position of bonds. In the l.h.s., the bottommost bond $`P_{j,k}`$ unites two cycles at sites $`j`$ and $`k`$ of the initial state, and makes a cycle of the two sites. At this stage, nothing happens to the r.h.s., and there is a cycle at $`j`$. Then, the rest of the operations of the l.h.s. and the r.h.s. are the same. That is, $`P_{i,j}`$ unites this cycle and another cycle; when $`P_{k,l}`$ in the l.h.s. unites cycles or breaks a cycle, $`P_{j,l}`$ in the r.h.s. does the same thing. Therefore, the number of cycles of the l.h.s. is equal to that of the r.h.s. As for the denominator of the trace, the l.h.s. has one more site than the r.h.s. Hence, in total, the r.h.s. needs a factor $`1/n`$. The logic above is true even if site $`i`$ is the same as site $`l`$.

The symmetrization of all the operators except $`P_{j,k}`$ is equivalent to the symmetrization of all the operators as explained in Appendix C. Hence, in all the terms in the symmetrization, one can put $`P_{j,k}`$ at the bottom of amida-diagrams and apply the contraction rule above. Therefore, we can obtain a contraction property also for a moment bond-diagram as shown in Fig. 9.

Finally, we prove that a cumulant bond-diagram has the same contraction property as the corresponding moment bond-diagram does. We use a mathematical induction for the proof.

1. The simplest nonzero cumulant, the l.h.s. of Fig. 5(a), has this contraction property.
2. Let us think about a certain cumulant bond-diagram. We assume that “the lower-order cumulant bond-diagrams have this contraction property.” Then, we use the mixed expansion, Eq. (15), pivoting on one of the bonds as shown Fig. 10. We use this assumption for the third term. The property of moments can be applied to the first and second terms.

Because of (i) and (ii), the contraction property is true for any cumulant bond diagram that satisfies the condition. We call this relation the Theorem I.

### IV.4 The “detachable” property

Some of the cumulants can be decomposed into a “background” and local quantities. We give an example in Fig. 11.

As already appeared several times, an intersite that has only one bond plays a special role in this paper. We call it a single bond here. On the other hand, an intersite that has more than one bond is called a multiple bond.

In this subsection, we prove Theorem II shown in Fig. 12. Here, a solid rectangle represents arbitrary bonds. A dotted line stands for a sequence of single bonds. The bond-diagram in the r.h.s. of (a) has the same length as that of the l.h.s.; the bond-diagram in a denominator of (b) has the same length as that of the numerator. This decomposition is possible even if the bond-diagram has more than two blocks of multiple bonds, on condition that the number of single bonds should be greater than that of blocks.

In other words, the Theorem II is described as follows: a bond-diagram can be rewritten by the product of a background and local quantities. The background is obtained by replacing each multiple bond with a single bond. One of the local quantities is obtained by replacing each multiple bond with a single bond except those in that block, and dividing it by the background.

#### IV.4.1 Proof of a lemma

For the proof, we introduce some notation. It is schematically shown in Fig. 13. A diagram is denoted by $`G`$. Here, solid rectangle $`A`$ ($`B`$) represents “rectangle I (II) and a sequence of single bonds” in Fig. 12 unitedly. We define $`G_A`$ ($`G_B`$) as a diagram such that each multiple bond in $`B`$ ($`A`$) of $`G`$ is replaced with a single bond. We define the “background” $`\mathrm{\Delta }_G`$ as a cumulant bond-diagram such that each multiple bond in $`A`$ and $`B`$ of $`G`$ is replaced with a single bond. Hence, there is a relation, $`\mathrm{\Delta }_G=\mathrm{\Delta }_{G_A}=\mathrm{\Delta }_{G_B}`$.

Let us consider the cumulant of $`G=P_{i,j}_{k,l}P_{k,l}`$. As is also shown in Fig. 14, the mixed expansion pivoting on $`P_{i,j}`$ is written as,

$$[G]_\mathrm{s}=C_1(G)+C_2(G),$$
(36)
$$C_1(G):=G_\mathrm{s}[P_{i,j}]_\mathrm{s}\underset{k,l}{}P_{k,l}_\mathrm{s},$$
(37)
$$C_2(G):=\underset{\mathrm{else}}{}\left[P_{i,j}\underset{k^{},l^{}}{}P_{k^{},l^{}}\right]_\mathrm{s}\underset{k^{\prime \prime },l^{\prime \prime }}{}P_{k^{\prime \prime },l^{\prime \prime }}_\mathrm{s}.$$
(38)
That is, the second term of $`C_1`$ is decomposition into $`P_{i,j}`$ and the rest. The other decompositions are included in $`C_2`$. In $`C_2`$, each cumulant includes more than one exchange operators, and thus, to give a nonzero cumulant, each site index has to appear at least twice as explained in Fig. 2(c) and (d). In other words, each cumulant bond-diagram in $`C_2`$ has to have a loop.

In order to prove the Theorem II, we first prove a lemma,

$$\frac{C_1(G)}{\mathrm{\Delta }_G}=\frac{C_1(G_A)}{\mathrm{\Delta }_G}\frac{C_1(G_B)}{\mathrm{\Delta }_G},$$
(39)
in the following.

First, we give an explicit example of $`C_1`$ using the cumulant in Fig. 5(b). In Sec. IV.2, we have symmetrized each of the terms in the mixed expansion individually. Here, however, we postpone the symmetrization procedure. That is, as shown in Fig. 15, we first sum amida-diagrams of the terms in $`C_1`$ with a certain configuration of bonds. After that, we symmetrize by averaging over all possible configuration of bonds. We carry out the symmetrization of all the operators except $`P_{1,\mathrm{}}`$ —– this is equivalent to the symmetrization of all the operators as shown in Appendix C. Let us consider Fig. 15(a). In the amida-diagram in the second term, sites 1 and $`\mathrm{}`$ belong to different cycles. Hence, when one more bond $`P_{1,\mathrm{}}`$ is added (to make the amida-diagram of the first term), the two cycles are united into one. Thus, the ratio of these amida-diagrams is $`1/n`$, which is equal to $`[P_{1,\mathrm{}}]_\mathrm{s}`$ in the second term. Accordingly, these terms cancel out each other. Therefore, there is a nonzero contribution to $`C_1`$ only when two sites 1 and $`\mathrm{}`$ belong to one cycle in the “amida-diagram without $`P_{1,\mathrm{}}`$” in the second term of $`C_1`$. Hereafter, such a configuration is called a $`C_1`$-finite configuration. Figure 15(b) is a $`C_1`$-finite configuration. Its contribution to $`C_1(G_{\mathrm{Fig}.\text{5}(\mathrm{b})})`$ is determined by “deviation from the background”. There is an “extra” bond compared to the background, and here it increases the number of cycles by one to make a factor $`n`$. On the other hand, the “probability” for this configuration to occur is 1/3 as shown in Eqs. (34) and (35). Therefore we obtain $`C_1(G_{\mathrm{Fig}.\text{5}(\mathrm{b})})=\frac{1}{3}n\mathrm{\Delta }_{G_{\mathrm{Fig}.\text{5}(\mathrm{b})}}`$.

The example above is a simple case. In general, “the probability of the configuration” times “the deviation from the background” is summed over all the $`C_1`$-finite configurations, namely,

$$C_1(G)=\mathrm{\Delta }_G\underset{\stackrel{C_1\text{-finite}}{\mathrm{config}.}}{}(\mathrm{probability})\times (\mathrm{deviation}).$$
(40)
Let us consider a diagram $`G`$ shown in Fig. 13. Then, in order to put sites 1 and $`\mathrm{}`$ in one cycle in the “amida-diagram without $`P_{1,\mathrm{}}`$” of $`G`$, each of $`G_A`$ and $`G_B`$ must also have a $`C_1`$-finite configuration. Then, the deviation from the background as a whole is given by (deviation in $`A`$)$`\times `$(deviation in $`B`$).

In fact, “probability” is also a product of those of $`A`$ and $`B`$ as shown below. For example, the probability of configuration in Fig. 16(a) is written as

$`2^2{\displaystyle _0^1}d\tau _{12}{\displaystyle _{\tau _{12}}^1}d\tau _{12}^{}{\displaystyle _{\tau _{12}}^{\tau _{12}^{}}}d\tau _{23}`$
$`{\displaystyle _0^1}d\tau _{45}{\displaystyle _{\tau _{45}}^1}d\tau _{45}^{}{\displaystyle _{\tau _{45}}^{\tau _{45}^{}}}d\tau _{34}.`$ (41)
The first three integrals are independent of the last three integrals, and thus the probability is written by the product of two parts. Next, let us think about Fig. 16(b). If we naively use $`\tau `$, the expression of the integrals is complicated. However, if we utilize the periodic boundary condition of $`\tau `$ and shift the integration variables as $`\rho :=\tau \tau _{23}`$ for $`\tau =\tau _{34},\tau _{34}^{},\tau _{45}`$, we obtain,

$`2^2{\displaystyle _0^1}d\tau _{12}{\displaystyle _{\tau _{12}}^1}d\tau _{12}^{}{\displaystyle _{\tau _{12}}^{\tau _{12}^{}}}d\tau _{23}`$
$`{\displaystyle _0^1}d\rho _{34}{\displaystyle _{\rho _{34}}^1}d\rho _{34}^{}{\displaystyle _{\rho _{34}}^{\rho _{34}^{}}}d\rho _{45}.`$ (42)
This expression is equivalent to Eq. (41). Hence the probability is written by the product of two parts again. The important point in deriving Eq. (42) is existence of $`P_{4,5}`$. Suppose the diagram does not have $`P_{4,5}`$, and has $`P_{1,4}`$ instead of $`P_{1,5}`$. The uppermost bond $`P_{1,4}`$ is always at $`\tau =1`$. If the integration variables are shifted, $`1\tau _{23}`$ appears in lower or upper bounds of the integrals. That is, the integrals are not written by a product of two parts any more in that case.

The shift of the integration variables mentioned above can generally be applied to any blocks $`A`$ and $`B`$. Hence, the probability in Eq. (40) is also a product of those of $`A`$ and $`B`$. The summation over $`C_1`$-finite configuration is double summation over $`C_1`$-finite configuration of $`A`$ and over $`C_1`$-finite configuration of $`B`$. Therefore,

$`{\displaystyle \frac{C_1(G)}{\mathrm{\Delta }_G}}`$ $`=`$ $`{\displaystyle \underset{\stackrel{C_1\text{-finite}}{\mathrm{config}.A}}{}}(\mathrm{probability}A)(\mathrm{deviation}A)`$ (43)
$`\times {\displaystyle \underset{\stackrel{C_1\text{-finite}}{\mathrm{config}.B}}{}}(\mathrm{probability}B)(\mathrm{deviation}B)`$
$`=`$ $`{\displaystyle \frac{C_1(G_A)}{\mathrm{\Delta }_{G_A}}}{\displaystyle \frac{C_1(G_B)}{\mathrm{\Delta }_{G_B}}}.`$
By using $`\mathrm{\Delta }_G=\mathrm{\Delta }_{G_A}=\mathrm{\Delta }_{G_B}`$, we obtain the lemma, Eq. (39).

Note that the shift of the integration variables mentioned above can be applied even if a diagram has more than two blocks of multiple bonds. A block and a single bond is regarded as one group whose members have a common shift of the integration variables. This is the reason for the restriction of the Theorem II; the number of single bonds should be greater than that of blocks.

#### IV.4.2 Proof of the Theorem II

In order to focus on the “multiple bond” region, we introduce another notation for moments and cumulants. Let us think about $`G_A`$. A set of bonds, $`A`$, is partitioned into two subsets in calculating $`C_2(G_A)`$. Let $`a`$ denote a subset of $`A`$. Then, $`Q(a)`$ denotes the cumulant in which $`A`$ is replaced by $`a`$. In choosing $`a`$, at least one bond per intersite must be chosen to make $`Q(a)`$ nonzero. A moment for bonds $`Aa`$ is denoted by $`(Aa)`$. Then, Eq. (36) for $`G_A`$ can be rewritten as

$`Q(A)=C_1(G_A){\displaystyle \underset{aA}{}}Q(a)(Aa).`$ (44)
The l.h.s. is equal to the excluded term $`a=A`$ in the summation in the r.h.s. because $`(\mathrm{})=1_0=1`$. Therefore, we obtain

$`C_1(G_A)={\displaystyle \underset{aA}{}}Q(a)(Aa).`$ (45)
The corresponding equation for graph $`G`$ is

$`{\displaystyle \frac{C_1(G)}{\mathrm{\Delta }}}={\displaystyle \underset{aA}{}}{\displaystyle \underset{bB}{}}{\displaystyle \frac{Q(a+b)}{\mathrm{\Delta }}}(Aa+Bb).`$ (46)
Since $`A`$ is not connected to $`B`$,

$$(Aa+Bb)=(Aa)(Bb).$$
(47)
On the other hand, according to the lemma and Eq. (45),

$`{\displaystyle \frac{C_1(G)}{\mathrm{\Delta }}}`$ $`=`$ $`{\displaystyle \frac{C_1(G_A)}{\mathrm{\Delta }}}{\displaystyle \frac{C_1(G_B)}{\mathrm{\Delta }}}`$ (48)
$`=`$ $`{\displaystyle \underset{aA}{}}{\displaystyle \underset{bB}{}}{\displaystyle \frac{Q(a)}{\mathrm{\Delta }}}{\displaystyle \frac{Q(b)}{\mathrm{\Delta }}}(Aa)(Bb).`$
Here, we use a mathematical induction again for the proof.

1. The Theorem II is true when there are two extra bonds on the background as explicitly calculated in Fig. 7(d).
2. We assume that the Theorem II is true when the number of extra bonds on the background is less than that of $`G`$. According to the assumption,

$$\frac{Q(a+b)}{\mathrm{\Delta }}=\frac{Q(a)}{\mathrm{\Delta }}\frac{Q(b)}{\mathrm{\Delta }},(a+bA+B).$$
(49)
Combine the r.h.s. of both Eqs. (46) and (48), use Eqs. (47) and (49). Then, all the terms except for $`(a,b)=(A,B)`$ cancel out and we obtain

$$\frac{Q(A+B)}{\mathrm{\Delta }}=\frac{Q(A)}{\mathrm{\Delta }}\frac{Q(B)}{\mathrm{\Delta }}.$$
(50)

According to (i) and (ii), the Theorem II is true for any $`G`$ that satisfies the conditions.

The proof mentioned above can be applied even if block $`A`$ or $`B`$ includes single bonds. Then, the Theorem II for a diagram with two blocks is repeatedly used to prove the Theorem II for a diagram with more than two blocks. Therefore, also a diagram with more than two blocks has the “detachable” property as mentioned first.

#### IV.4.3 Theorem I’

Now that the Theorem II is proved, the contraction can be written in a more general form. Although the Theorem I is for a sequence of three single bonds, sequences of two single bonds also can be contracted as shown in Fig. 17 as Theorem I’, on condition that each block must be separated by at least one single bond and two of the blocks must be separated by a sequence of at least two single bonds. The Theorem I’ is easily proved using Fig. 12(a); contract $`\mathrm{\Delta }_G`$ in the r.h.s. and reversely use the same relation to the contracted diagram.

### IV.5 The quasi-cumulant method

Apart from the restriction, the Theorem II is equivalent to a property of moments —– a disconnected moment is equal to a product of connected moments. In order to effectively use the Theorem II, we define quasi-moments and corresponding quasi-cumulants below. Our definition of a quasi-moment is a factor times a “local quantity” in the Theorem II. Although another definition can be the local quantity itself, it is less convenient for our purpose as discussed in Appendix G.

Let us think about the one-dimensional system again. We define the generating function of quasi-moments as

$`𝒢_{\mathrm{qm}}^{(i,j)}`$ $`:=`$ $`{\displaystyle \frac{1}{\lambda _{i,i+1}\lambda _{i+1,i+2}\mathrm{}\lambda _{j1,j}}}`$ (51)
$`\times {\displaystyle \frac{p_{i,j}}{[P_{i,j}P_{i,i+1}P_{i+1,i+2}\mathrm{}P_{j1,j}]_\mathrm{s}}}.`$
Only up to a finite order, the Theorem II can be applied and $`𝒢_{\mathrm{qm}}`$ works as a generating function. The larger $`|ij|`$ is, the higher order the Theorem II is valid up to. Here, we assume that $`|ij|`$ is large enough, and forget about this restriction for the time being. It will be commented on later. Then, the quasi-moments are defined as

$$P_{k,l}\mathrm{}P_{y,z}_\mathrm{q}:=\frac{}{\lambda _{k,l}}\mathrm{}\frac{}{\lambda _{y,z}}𝒢_{\mathrm{qm}}^{(i,j)}|_{\lambda =0}.$$
(52)
The zeroth-order quasi-moment satisfies $`1_\mathrm{q}=1`$. Some examples of quasi-moments are shown in Fig.18. More explicitly, the quasi-moments are obtained by

$`P_{1,2}^{m_1}P_{2,3}^{m_2}\mathrm{}P_{k1,k}^{m_{k1}}_\mathrm{q}={\displaystyle \frac{\left[P_{1,2}^{m_1+1}P_{2,3}^{m_2+1}\mathrm{}P_{k1,k}^{m_{k1}+1}P_{k,k+1}P_{1,k+1}\right]_\mathrm{s}}{(m_1+1)(m_2+1)\mathrm{}(m_{k1}+1)\left[P_{1,2}P_{2,3}\mathrm{}P_{k1,k}P_{k,k+1}P_{1,k+1}\right]_\mathrm{s}}},`$ (53)
where we have assumed that all the bonds are on the “background”. We need a special care when $`P_{k,l}`$ is outside the background, namely, $`k<i`$ or $`j<l`$ in the case of $`𝒢_{\mathrm{qm}}^{(i,j)}`$. Such a operator is explicitly denoted by $`P_{k,l}^{(\mathrm{o})}`$. Then, we can more generally write as,

$$(P_{1,2}^{(\mathrm{o})})^{m_1}\mathrm{}(P_{i1,i}^{(\mathrm{o})})^{m_{i1}}P_{i,i+1}^{m_i}\mathrm{}P_{k1,k}^{m_{k1}}_\mathrm{q}=\frac{\left[(P_{1,2}^{(\mathrm{o})})^{m_1}\mathrm{}(P_{i1,i}^{(\mathrm{o})})^{m_{i1}}P_{i,i+1}^{m_i+1}\mathrm{}P_{k1,k}^{m_{k1}+1}P_{k,k+1}P_{i,k+1}\right]_\mathrm{s}}{(m_i+1)\mathrm{}(m_{k1}+1)\left[P_{i,i+1}\mathrm{}P_{k1,k}P_{k,k+1}P_{i,k+1}\right]_\mathrm{s}}.$$
(54)
Next, using quasi-moment defined above, the generating function of quasi-cumulants is defined as,

$$𝒢_{\mathrm{qc}}^{(i,j)}:=\mathrm{ln}𝒢_{\mathrm{qm}}^{(i,j)}.$$
(55)
Then, the quasi-cumulants are defined by,

$$[P_{k,l}\mathrm{}P_{y,z}]_\mathrm{q}:=\frac{}{\lambda _{k,l}}\mathrm{}\frac{}{\lambda _{y,z}}𝒢_{\mathrm{qc}}^{(i,j)}|_{\lambda =0}.$$
(56)
According to the Theorem II, when none of $`k`$ or $`l`$ is equal to $`y`$ or $`z`$,

$`\mathrm{exp}\left[{\displaystyle \underset{k,l}{}}\lambda _{k,l}P_{k,l}+{\displaystyle \underset{y,z}{}}\lambda _{y,z}P_{y,z}\right]_\mathrm{q}=`$
$`\mathrm{exp}\left[{\displaystyle \underset{k,l}{}}\lambda _{k,l}P_{k,l}\right]_\mathrm{q}\mathrm{exp}\left[{\displaystyle \underset{y,z}{}}\lambda _{y,z}P_{y,z}\right]_\mathrm{q},`$ (57)
and thus, as discussed in Appendix B,

$$\left[\underset{k,l}{}P_{k,l}\underset{y,z}{}P_{y,z}\right]_\mathrm{q}=0.$$
(58)
Therefore, only “connected” quasi-cumulants remain nonzero. If only the quasi-moments are used, we have to classify connected and disconnected quasi-moments. However, if the quasi-cumulants are used, we need only connected quasi-cumulants and can avoid the complicated classification.

Suppose we set $`\lambda _{k,l}=\lambda `$ for every $`(k,l)`$ pair. Then, $`𝒢_{\mathrm{qc}}^{(i,i+x)}`$ $`(x>0)`$ is a function of $`x`$ and $`\lambda `$. The quasi-cumulants have a translational symmetry in the background. For example, $`[P_{i,i+1}^m]_\mathrm{q}=[P_{i+1,i+2}^m]_\mathrm{q}=\mathrm{}=[P_{i+x1,i+x}^m]_\mathrm{q}`$; and thus, these “length-1” quasi-cumulants can be summed up to $`x[P_{i,i+1}^m]_\mathrm{q}`$ in $`𝒢_{\mathrm{qc}}^{(i,i+x)}`$. If we in the same way simplify $`𝒢_{\mathrm{qc}}^{(i,i+x)}`$, length-2 quasi-cumulants are multiplied by $`x1`$, length-3 quasi-cumulants are multiplied by $`x2`$, and so on. The quasi-cumulants including $`P_{k,l}^{(\mathrm{o})}`$ are located at both ends, and thus, they are multiplied by two. Namely,

$$𝒢_{\mathrm{qc}}^{(i,i+x)}=𝒢_m^{(x)}+𝒢_m^{\mathrm{end}}+O(\lambda ^{m+1}),$$
(59)
$`𝒢_m^{(x)}`$ $`:=`$ $`x{\displaystyle \underset{1m1m}{}}{\displaystyle \frac{\lambda ^{m_1}}{m_1!}}[P_{i,i+1}^{m_1}]_\mathrm{q}`$ (60)
$`+(x1){\displaystyle \underset{\stackrel{1m_1,m_2}{m_1+m_2m}}{}}{\displaystyle \frac{\lambda ^{m_1+m_2}}{m_1!m_2!}}[P_{i,i+1}^{m_1}P_{i+1,i+2}^{m_2}]_\mathrm{q}`$
$`+\mathrm{}\mathrm{}\mathrm{}`$
$`+(xm+1)\lambda ^m[P_{i,i+1}\mathrm{}P_{i+m1,i+m}]_\mathrm{q},`$
$`𝒢_m^{\mathrm{end}}`$ $`:=`$ $`2{\displaystyle \underset{2m_1m}{}}{\displaystyle \frac{\lambda ^{m_1}}{m_1!}}[(P_{i1,i}^{(\mathrm{o})})^{m_1}]_\mathrm{q}`$ (61)
$`+2{\displaystyle \underset{\stackrel{2m_1,1m_2}{m_1+m_2m}}{}}{\displaystyle \frac{\lambda ^{m_1+m_2}}{m_1!m_2!}}[(P_{i1,i}^{(\mathrm{o})})^{m_1}P_{i,i+1}^{m_2}]_\mathrm{q}`$
$`+2{\displaystyle \underset{\stackrel{2m_1,m_2}{m_1+m_2m}}{}}{\displaystyle \frac{\lambda ^{m_1+m_2}}{m_1!m_2!}}[(P_{i2,i1}^{(\mathrm{o})})^{m_1}(P_{i1,i}^{(\mathrm{o})})^{m_2}]_\mathrm{q}`$
$`+\mathrm{}\mathrm{}\mathrm{}.`$
Hence, this expression for arbitrary $`x`$ can be derived if we obtain all the quasi-cumulants up to $`m`$-th order, which are calculable in small systems. Therefore, we can obtain $`p_{i,i+x}`$ up to $`(x+m)`$-th order for arbitrary $`x`$ by tracing back Eqs. (59)-(61), (55) and (51). We call this method the quasi-cumulant method (QCM).

The QCM has been done by a brute-force program of Mathematica. The quasi-cumulants can be calculated by Eq. (56), that is, differentiating the series of their generating function with respect to corresponding $`J_{i,i+1}`$, where $`J_{i,i+1}`$ are coupling constants of the Hamiltonian of a rather general form, $`_iJ_{i,i+1}P_{i,i+1}`$. We take the system as small as possible on condition that it be large enough to correctly calculate those quasi-cumulants. That is, quasi-cumulants of length $`l`$ can be obtained in the system of length $`l+1`$. First, we calculate the series of $`p_{i,j}`$ to obtain $`𝒢_{\mathrm{qm}}^{(i,j)}`$ up to $`m`$-th order. At each order of calculating $`^k`$, we exclude terms irrelevant to the derivative in order to save computational time and memory. Let us rewrite $`𝒢_{\mathrm{qm}}^{(i,j)}`$ as

$$𝒢_{\mathrm{qm}}^{(i,j)}=\underset{k=0}{\overset{m}{}}\frac{(\beta )^k}{k!}^k_\mathrm{q}+O(\beta ^{m+1}),$$
(62)
then, Eq. (15), the mixed expansion pivoting on one of $``$ can be used,

$$[^k]_\mathrm{q}=^k_\mathrm{q}\underset{l=1}{\overset{k1}{}}\frac{(k1)!}{(l1)!(kl)!}[^l]_\mathrm{q}^{kl}_\mathrm{q},$$
(63)
$$𝒢_{\mathrm{qc}}^{(i,j)}=\underset{k=1}{\overset{m}{}}\frac{(\beta )^k}{k!}\left[^k\right]_\mathrm{q}+O(\beta ^{m+1}).$$
(64)
Now we are ready for the derivatives. The point is that we need at most $`m+2`$ sites here, while we can reconstruct the generating function of larger systems from the quasi-cumulants and thus we can obtain $`p_{i,i+x}`$ up to $`(x+m)`$-th order arbitrary $`x`$.

Finally, we comment on the restriction of the QCM. The QCM is based on the Theorem II. Let us consider diagrams shown in Fig. 19. We cannot apply the Theorem II to these diagrams, and thus Eq. (57) is not valid. However, we can apply the Theorem II if the order is lower than these. Hence, we can use Eqs. (59)-(61), only for $`mx/2`$ ($`x`$: even) or $`m(x1)/2`$ ($`x`$: odd).

Note that the QCM is useful for contributions from large clusters rather than those from small clusters. Hence, by combining the QCM with the FCM, we can calculate the high-temperature series coefficients up to a high order. Figure 20 shows which of the FCM and the QCM was used to obtain the series coefficients in Ref. FK1d, . We calculated $`p_{i,i+x}`$ up to $`(x+6)`$-th order for $`x13`$ by the QCM. Note that the QCM results are valid for arbitrary $`x13`$, including the $`x\mathrm{}`$ limit as indicated in Fig. 20. The reason why we stopped the QCM at this order was mainly due to the FCM. If the FCM can treat larger systems, the restriction of the Theorem II is less important, and the QCM can go further. Our FCM program was for general SU($`n`$); we expect that more terms can be obtained if the interest is restricted to SU(2) in making the FCM program.

In Ref. FK1d, , the largest cluster we used for the FCM was composed of $`\mathrm{}=13`$ sites. As noted in the context of Eq. (30), $`p_{i,i+x}_{\mathrm{}}^{}`$ is correct only up to $`O[(\beta J)^{2\mathrm{}x1}]`$, and higher orders include contributions from the larger clusters. Hence, for high-order coefficients of $`x12`$, corrections to the FCM results must be calculated. The QCM can be used also for this purpose because these corrections are a subset of the terms obtained by the QCM. We put a fictitious coefficient on each quasi-cumulant located at both ends in Eq. (61). These fictitious coefficients carry information about the cluster-size of a cumulant, and we can distinguish terms needed for the corrections to the FCM result. Coefficients obtained by this procedure in Ref. FK1d,  are indicated by “FCM+QCM” in Fig. 20.

Since these series coefficients are obtained as explicit functions of $`n`$, the data of the series coefficients are too numerous to be fully listed in this paper. Hence, we list here only the results of the QCM, $`p_{i,i+x}`$ up to $`(x+6)`$-th order, and the other data will be provided on request. The coefficients of the series

$$p_{i,i+x}=n^{1x}\left(1+n^2\right)\underset{m}{}\pi _{x,m}(\beta J)^m$$
(65)
are given in the following:

$$\pi _{x,x}=1(x1),$$
(66)
$$\pi _{x,x+1}=\frac{\left(6+n^2\right)x}{6n}(x2),$$
(67)
$`\pi _{x,x+2}={\displaystyle \frac{n^2}{360}}\left[5x^2\left(6+n^2\right)^26x\left(30+15n^2+2n^4\right)+n^2\left(150+7n^2\right)\right]`$
$`(x4),`$ (68)
$`\pi _{x,x+3}`$ $`=`$ $`{\displaystyle \frac{n^3}{45360}}[35x^3(6+n^2)^3126x^2(180120n^2+3n^4+2n^6)`$ (69)
$`+2x\left(7560+17010n^22961n^4+302n^6\right)`$
$`+24n^2(945+21n^220n^4)](x6),`$
$`\pi _{x,x+4}`$ $`=`$ $`{\displaystyle \frac{n^4}{5443200}}[175x^4(6+n^2)^41260x^3(6+n^2)^2(30+15n^2+2n^4)`$ (70)
$`+2x^2\left(12474001965600n^2+519750n^448210n^6+6817n^8\right)`$
$`6x\left(226800+1020600n^2297990n^436795n^6+5468n^8\right)`$
$`+9n^2(352800+212520n^214850n^4+3287n^6)](x8),`$
$`\pi _{x,x+5}`$ $`=`$ $`{\displaystyle \frac{n^5}{359251200}}[385x^5(6+n^2)^54620x^4(6+n^2)^3(30+15n^2+2n^4)`$ (71)
$`+44x^3\left(2381400+3572100n^21260630n^4+147915n^616461n^8+2021n^{10}\right)`$
$`66x^2\left(22680006350400n^2+2891700n^4101940n^653013n^8+6496n^{10}\right)`$
$`+2x(35925120+264448800n^2182037240n^4`$
$`+17223030n^63830673n^8+517472n^{10})`$
$`36n^2(6652800+7983360n^260720n^4165099n^6+27776n^8)]`$
$`(x10),`$
$`\pi _{x,x+6}`$ $`=`$ $`{\displaystyle \frac{n^6}{5884534656000}}[175175x^6(6+n^2)^6`$ (72)
$`3153150x^5\left(6+n^2\right)^4\left(30+15n^2+2n^4\right)`$
$`+5005x^4\left(6+n^2\right)^2\left(38556004725000n^2+982800n^4+18870n^6+18947n^8\right)`$
$`6006x^3(306180000+731430000n^2413229600n^4+58660200n^6`$
$`+2763810n^81239825n^{10}+126908n^{12})`$
$`+13x^2(172260950400745940210400n^2+577975305600n^4101282227200n^6`$
$`+6937664580n^82448983250n^{10}+265586033n^{12})`$
$`6x(163459296000+1777619844000n^22058679022400n^4`$
$`+396577981800n^6+31224022830n^811642353905n^{10}+1392291386n^{12})`$
$`+135n^2(32691859200+61509127680n^215679984320n^4+695398704n^6`$
$`446293302n^8+62451523n^{10})](x12).`$
Even for $`\mathrm{}12`$, the QCM is valid for a number of series coefficients, and should coincide with the results of the FCM or the “FCM+QCM”. We have used this property for checking whether the computer codes are correct.

## V A technique for the specific heat

Our strategy for the specific heat in Ref. FK1d,  was similar to that for the correlation functions. That is, the most time-consuming part of the FCM is calculated by another method instead. If only the FCM is used, the series for $`P_{i,i+1}`$ up to $`O[(\beta J)^M]`$ requires systems with $`\mathrm{}M/2+1`$ ($`M`$: even) or $`\mathrm{}(M+1)/2+1`$ ($`M`$: odd). In Ref. FK1d, , we calculated $`P_{i,i+1}`$ up to $`O[(\beta J)^{22}]`$. However, the FCM was used only for $`\mathrm{}11`$; contributions from the required largest cluster, namely, the lowest- and the second-lowest-order nonzero contribution from a cluster, was directly calculated by a new technique explained below.

In fact, we can use $`\mathrm{}_0`$ as a moment and define a cumulantFulde  without the symmetrization. The generating functions of the moments and the cumulants are defined as

$$𝒢_{0\mathrm{m}}^{i_1,j_1,\mathrm{},i_k,j_k}:=e^{\lambda _1P_{i_1,j_1}}e^{\lambda _2P_{i_2,j_2}}\mathrm{}e^{\lambda _kP_{i_k,j_k}}_0,$$
(73)
$$𝒢_{0\mathrm{c}}^{i_1,j_1,i_2,j_2,\mathrm{},i_k,j_k}:=\mathrm{ln}𝒢_{0\mathrm{m}}^{i_1,j_1,i_2,j_2,\mathrm{},i_k,j_k},$$
(74)
respectively. Then, $`P_{i_1,j_1}\mathrm{}P_{i_k,j_k}_0`$ works as a moment, and the corresponding cumulant is define by

$`[P_{i_1,j_1}\mathrm{}P_{i_k,j_k}]_0:={\displaystyle \frac{}{\lambda _1}}\mathrm{}{\displaystyle \frac{}{\lambda _k}}𝒢_{0\mathrm{c}}^{i_1,j_1,\mathrm{},i_k,j_k}|_{\lambda =0}.`$ (75)
Here, the ordering of exchange operators in $`[\mathrm{}]_0`$ is important. In other words, the result can be changed by exchanging positions of uncommutative operators. Note that simplification such as $`P_{1,2}P_{1,2}=1`$ is not possible in $`[\mathrm{}]_0`$ while it is possible in $`\mathrm{}_0`$. In expanding the cumulants using the moments, the ordering of operators in $`\mathrm{}_0`$ must be the same as that in $`[\mathrm{}]_0`$, for example, $`[P_{1,2}P_{2,3}]_0=P_{1,2}P_{2,3}_0P_{1,2}_0P_{2,3}_0`$. This cumulant is represented by the cumulant version of an amida-diagram. Hereafter, we distinguish between a moment amida-diagram and a cumulant amida-diagram by bracketing the diagram by $`\mathrm{}`$ and $`[\mathrm{}]`$, respectively.

The free energy is rewritten as

$$F=\frac{1}{\beta }\mathrm{ln}Z_0\frac{1}{\beta }\underset{m=1}{\overset{\mathrm{}}{}}\frac{(\beta )^m}{m!}[^m]_0,$$
(76)
$`[^m]_0`$ $`=`$ $`{\displaystyle \underset{(i_1,j_1)}{}}\mathrm{}{\displaystyle \underset{(i_m,j_m)}{}}J_{i_1,j_1}\mathrm{}J_{i_m,j_m}`$ (77)
$`\times \left[P_{i_1,j_1}\mathrm{}P_{i_m,j_m}\right]_0.`$
These are equivalent to Eqs. (21) and (22) except for the subscript “0” replacing “s”. Here, however, the cumulants are generally dependent on the ordering of the operators while “s” makes many terms equivalent.

We can use general property of moments and cumulants again. If,

$`{\displaystyle \underset{k}{}}\mathrm{exp}\left[\lambda _kP_{i_k,j_k}\right]{\displaystyle \underset{l}{}}\mathrm{exp}\left[\lambda _lP_{i_l,j_l}\right]_0=`$
$`{\displaystyle \underset{k}{}}\mathrm{exp}\left[\lambda _kP_{i_k,j_k}\right]_0{\displaystyle \underset{l}{}}\mathrm{exp}\left[\lambda _lP_{i_l,j_l}\right]_0,`$ (78)
then, the derivative of its logarithm gives,

$$\left[\underset{k}{}P_{i_k,j_k}\underset{l}{}P_{i_l,j_l}\right]_0=0.$$
(79)
When none of $`i_k`$ or $`j_k`$ is equal to $`i_l`$ or $`j_l`$, this relation is obviously satisfied. Furthermore, the cumulant is equal to zero also when an amida-diagram is separated into two diagrams by cutting only one point on a vertical line as shown in Fig. 21. For example, the cumulant versions of the amida-diagrams in Fig. 3(a) and (b) are equal to zero. Let us think about the moment amida-diagram, Fig. 21(a), where bonds are located only in the painted regions. Bonds located left of $`\xi `$ are denoted by $`P_{i_l,j_l}`$, and those located right of $`\xi `$ are denoted by $`P_{i_r,j_r}`$. The diagram can be written as $`_lP_{i_l,j_l}_rP_{i_r,j_r}_0`$ by exchanging positions of commutative $`P_{i_l,j_l}`$ and $`P_{i_r,j_r}`$. Let us imagine calculating each of $`_lP_{i_l,j_l}`$ and $`_rP_{i_r,j_r}`$ by counting cycles, and make the product of these two at the end. This final operation unites the two cycles relevant to $`\xi `$ into one cycle because $`_lP_{i_l,j_l}`$ and $`_rP_{i_r,j_r}`$ are sharing only one site $`\xi `$. Hence, the number of cycles in the l.h.s. is less than that in the r.h.s. “by one”. However, the number of sites in the l.h.s. is also less than that in the r.h.s. “by one”. As a result, the l.h.s. is equal to the r.h.s. We can make the same argument to $`_le^{\lambda _lP_{i_l,j_l}}_re^{\lambda _rP_{i_r,j_r}}_0`$. Therefore, $`[_lP_{i_l,j_l}_rP_{i_r,j_r}]_0=0`$.

Here, we must be careful about the periodic boundary condition of amida-diagrams. If a diagram satisfies the condition above by using the periodic boundary condition, the amida-diagram is also equal to zero. We give an example in Fig. 22. At a sight, $`[P_{1,2}P_{2,3}P_{2,3}P_{1,2}]_0`$ does not satisfy the condition above. However, it is equivalent to $`[P_{1,2}P_{1,2}P_{2,3}P_{2,3}]_0=0`$.

Suppose there are two bonds per nearest-neighbor pair as shown in Fig. 24. It is a contribution from cluster-size $`\mathrm{}`$ for $`O[(\beta J)^{2(\mathrm{}1)1}]`$ of $`P_{i,i+1}`$, or equivalently, $`O[(\beta J)^{2(\mathrm{}1)}]`$ of the free energy. Because of the property mentioned above, many configurations of bonds make cumulant amida-diagrams equal to zero. In order to give a nonzero contribution, at each vertical line bonds at left and right have to appear alternately. We can make such diagrams by repeating a enlarging procedure shown in Fig. 23. That is, there are two possibilities for next bonds. In fact, the value of such a cumulant amida-diagram does not depend on those configurations as explained in Appendix H. Hence, we can just multiply the value of the cumulants by the probability for such nonzero configurations. We have calculated the probability by Mathematica using integrals like Eqs. (35) and (41).

The values of the cumulants are shown in Fig. 24. Namely, the cumulant is equal to $`(1n^2)`$ when $`\mathrm{}`$ is an even number, and $`(1n^2)`$ when $`\mathrm{}`$ is an odd number. We can prove this by using the mixed expansion pivoting on one of $`P_{1,2}`$ shown in Fig. 25. The first term is the moment. The second term is $`[P_{1,2}]`$ times the moment of the others. The rest of the terms are obtained by separating the diagram at sites from 2 to $`\mathrm{}1`$. In fact, only the last term survives and every term else cancels out one another. When $`\mathrm{}`$ is an even number, the sum of the first, second and third terms, the sum of the fourth and fifth terms, …, are equal to zero. When $`\mathrm{}`$ is an odd number, the sum of the first and second terms, the sum of the third and fourth terms, …, are equal to zero. The details are explained in Appendix H.

We can calculate also the next order by the same technique as above. Now, one of the nearest-neighbor pairs has three bonds. The value of the cumulant amida-diagram does not depend on the position of this triple bond as shown in Fig. 26. That is, the cumulant is equal to $`\frac{2}{n}(1\frac{1}{n^2})`$ when $`\mathrm{}`$ is an even number, and equal to $`\frac{2}{n}(1\frac{1}{n^2})`$ when $`\mathrm{}`$ is an odd number.

## VI summary

We have formulated a new method of high-temperature series expansion using the SU($`n`$) Heisenberg model. It is designed for efficiently calculating contributions from large clusters, which is actually the most time-consuming part of a standard method. The net contributions from a cluster to high-temperature series start from a certain order. We have focused on this property, and considered deviation from the lowest-order nonzero contribution.

Although we have mainly used one-dimensional systems to explain the new method, many of the techniques can be used also in higher dimensions. In contrast to one-dimensional systems, however, the calculation requires various clusters, and some of the clusters are unsuitable for our method. Figure 27 shows an example. The main advantage to use our method is that diagrams can be reduced to smaller diagrams. However, we cannot apply it to a cluster in Fig. 27 because it has many sites at which more than two bonds meet, and chains connecting those sites are too short. Nevertheless, the number of embeddings of such unsuitable clusters in the infinite lattice is small, and thus contributions from them are usually small. In fact, the number of embeddings of open chains, the most preferable clusters for our method, is the largest. Hence, the dominant contributions can be calculated by our method.

We expect that similar approaches are possible also for other models, namely, for Hamiltonians of lower symmetries.

###### Acknowledgements.

The author is grateful to M. Arikawa for indicating him the equivalence of an amida-lottery and a permutation. He also thanks N. Shannon, K. Penc and Y. Kuramoto for useful comments on the manuscript. This work was financially supported by Visitors’ Program of the MPI–PKS and the Japan Society for the Promotion of Science.

## Appendix A relations derived from the isotropy

For example, using $`n^21`$ generators such as

$$𝒳^{\alpha \beta }:=\sqrt{\frac{n}{2}}\left(X^{\alpha \beta }+X^{\beta \alpha }\right)(1\alpha <\beta n),$$
(80)
$$𝒴^{\alpha \beta }:=i\sqrt{\frac{n}{2}}\left(X^{\alpha \beta }X^{\beta \alpha }\right)(1\alpha <\beta n),$$
(81)
$`D^{(\alpha )}`$ $`:=`$ $`\sqrt{{\displaystyle \frac{n}{\alpha (\alpha 1)}}}\left({\displaystyle \underset{\gamma =1}{\overset{\alpha 1}{}}}X^{\gamma \gamma }(\alpha 1)X^{\alpha \alpha }\right)`$ (82)
$`(2\alpha n),`$
the exchange operator is rewritten as

$`nP_{i,j}`$ $`=`$ $`{\displaystyle \underset{1\alpha <\beta n}{}}\left(𝒳_i^{\alpha \beta }𝒳_j^{\alpha \beta }+𝒴_i^{\alpha \beta }𝒴_j^{\alpha \beta }\right)`$ (83)
$`+{\displaystyle \underset{2\alpha n}{}}D_i^{(\alpha )}D_j^{(\alpha )}+1.`$
Let us prove that each component contributes equally, namely,

$`𝒳_i^{\alpha \beta }𝒳_j^{\alpha \beta }=𝒴_i^{\alpha \beta }𝒴_j^{\alpha \beta }=D_i^{(\alpha )}D_j^{(\alpha )}`$
$`={\displaystyle \frac{nP_{i,j}1}{n^21}}.`$ (84)
1. The Hamiltonian is invariant under re-labeling of colors. Therefore, $`𝒳_i^{\alpha \beta }𝒳_j^{\alpha \beta }`$ does not depend on $`\alpha `$ or $`\beta `$.
2. $`𝒳_i^{\alpha \beta }𝒳_j^{\alpha \beta }=𝒴_i^{\alpha \beta }𝒴_j^{\alpha \beta }`$ for arbitrary $`\alpha `$, $`\beta `$, because these $`𝒳`$, $`𝒴`$ operators are $`x`$, $`y`$ components of Pauli matrix $`(\times \sqrt{n/2})`$ regarding only $`\alpha `$ and $`\beta `$.
3. We can apply the same argument also for $`z`$-component, and thus $`𝒳_i^{\alpha \beta }𝒳_j^{\alpha \beta }=D_i^{(2)}D_j^{(2)}`$.
4. By a direct calculation using

$$X_i^{\alpha \alpha }X_j^{\alpha \alpha }=X_i^{11}X_j^{11}$$

for arbitrary $`\alpha `$, and

$$X_i^{\alpha \alpha }X_j^{\beta \beta }=X_i^{11}X_j^{22}$$

for arbitrary $`\alpha \beta `$, we can derive

$`D_i^{(\alpha )}D_j^{(\alpha )}`$ $`=`$ $`n\left(X_i^{11}X_j^{11}X_i^{11}X_j^{22}\right)`$
$`=`$ $`D_i^{(2)}D_j^{(2)}`$
for arbitrary $`\alpha `$.

From (i)-(iv), we obtain Eq. (84).

Using Eq. (84), we can derive

$`X_i^{\alpha \beta }X_j^{\beta \alpha }`$ $`=`$ $`{\displaystyle \frac{1}{n}}𝒳_i^{\alpha \beta }𝒳_j^{\alpha \beta }`$ (85)
$`=`$ $`{\displaystyle \frac{1}{n^21}}\left(P_{i,j}{\displaystyle \frac{1}{n}}\right),`$
for $`\alpha \beta ,ij`$.

In addition, for $`i=j`$, since the average of $`|\alpha \alpha |`$ does not depend on $`\alpha `$,

$$X_i^{\alpha \beta }X_i^{\beta \alpha }=|\stackrel{i}{\alpha }\stackrel{i}{\alpha }|=\frac{1}{n}\underset{\gamma =1}{\overset{n}{}}|\stackrel{i}{\gamma }\stackrel{i}{\gamma }|.$$
(86)
Then, using

$$\underset{\gamma =1}{\overset{n}{}}|\stackrel{i}{\gamma }\stackrel{i}{\gamma }|=1,$$
(87)
we obtain

$$X_i^{\alpha \beta }X_i^{\beta \alpha }=\frac{1}{n}1=\frac{1}{n}.$$
(88)
## Appendix B Moments and cumulants for classical variables

In this section, we review fundamental properties of moments and cumulants of classical variables $`x_i(i=1,2,\mathrm{})`$. Note that some other properties are written in Sec. III.1. Here, the expectation value with respect to the distribution of $`\{x_i\}`$ is denoted by $`\mathrm{}_x`$. The generating function of moments is defined as

$$e^{_i\lambda _ix_i}_x.$$
(89)
That is, the moments are derived by

$$x_jx_k\mathrm{}x_l_x=\frac{}{\lambda _j}\frac{}{\lambda _k}\mathrm{}\frac{}{\lambda _l}e^{_i\lambda _ix_i}_x|_{\lambda =0},$$
(90)
where the subscript $`\lambda =0`$ represents that all the $`\lambda `$-variables are set to zero. On the other hand, cumulants are denoted by $`[\mathrm{}]_x`$ and defined as

$$[x_jx_k\mathrm{}x_l]_x=\frac{}{\lambda _j}\frac{}{\lambda _k}\mathrm{}\frac{}{\lambda _l}\mathrm{ln}e^{_i\lambda _ix_i}_x|_{\lambda =0}.$$
(91)
Namely, the generating function of the cumulants is written as

$$\mathrm{ln}e^{_i\lambda _ix_i}_x=[e^{_i\lambda _ix_i}1]_x.$$
(92)
A cumulant is equal to zero if the variables in $`[\mathrm{}]_x`$ can be partitioned into two groups which are independent of each other in averaging $`\mathrm{}_x`$. Suppose $`\{x_i\}`$ is independent of $`\{x_j\}`$. Let us in advance set $`\lambda _k`$ to zero if $`x_k`$ is not included in either $`\{x_i\}`$ or $`\{x_j\}`$. Then,

$$e^{_i\lambda _ix_i+_j\lambda _jx_j}_x=e^{_i\lambda _ix_i}_xe^{_j\lambda _jx_j}_x.$$

By taking the logarithm, we obtain

$$\mathrm{ln}e^{_i\lambda _ix_i}_x+\mathrm{ln}e^{_j\lambda _jx_j}_x.$$

That is, the generating function of the cumulants does not have any term which has a product of $`\lambda _i`$ and $`\lambda _j`$. Thus, the corresponding cumulant is equal to zero.

There are relations between moments and cumulants. Moments can be expanded using cumulants, for example,

$`x_1_x`$ $`=`$ $`[x_1]_x,`$
$`x_1x_2_x`$ $`=`$ $`[x_1x_2]_x+[x_1]_x[x_2]_x,`$
$`x_1x_2x_3_x`$ $`=`$ $`[x_1x_2x_3]_x+[x_1x_2]_x[x_3]_x+[x_1x_3]_x[x_2]_x`$
$`+[x_1]_x[x_2x_3]_x+[x_1]_x[x_2]_x[x_3]_x.`$
These relations are generally written as

$$x_i\mathrm{}x_l_x=\underset{\zeta }{}\underset{𝒫(\xi ,\zeta )}{}[x_i\mathrm{}x_j]_x\mathrm{}[x_k\mathrm{}x_m]_x.$$
(93)
These summations include every partition of $`x`$-variables. Here, $`\xi `$ is the number of $`x`$-variables in the bracket in the l.h.s. Then, the summation denoted by $`𝒫(\xi ,\zeta )`$ is taken over every partition of $`\xi `$ elements into $`\zeta `$ groups ($`1\zeta \xi `$). In other words, $`\zeta `$ is the number of brackets $`[\mathrm{}]_x`$ in the r.h.s. Here, each $`[\mathrm{}]_x`$ includes at least one $`x`$.

Conversely, cumulants can be expanded using moments, for example,

$`[x_1]_x`$ $`=`$ $`x_1_x,`$
$`[x_1x_2]_x`$ $`=`$ $`x_1x_2_xx_1_xx_2_x,`$
$`[x_1x_2x_3]_x`$ $`=`$ $`x_1x_2x_3_xx_1x_2_xx_3_xx_1x_3_xx_2_x`$
$`x_1_xx_2x_3_x+2x_1_xx_2_xx_3_x.`$
In general, Eq. (93) is inverted to

$`[x_i\mathrm{}x_{\mathrm{}}]_x`$ $`=`$ $`{\displaystyle \underset{\zeta }{}}{\displaystyle \underset{𝒫(\xi ,\zeta )}{}}(1)^{\zeta 1}(\zeta 1)!`$ (94)
$`\times x_i\mathrm{}x_j_x\mathrm{}x_k\mathrm{}x_m_x.`$
## Appendix C symmetrization

Let us consider expanding $`^m`$. It includes all the orderings of $`P_{i,j}`$ operators. In other words, it is already symmetrized. Therefore,

$$^m_\mathrm{s}=^m_0.$$
(95)
Accordingly,

$$e^\beta _\mathrm{s}=e^\beta _0.$$
(96)
Because of a a property of trace, one can cyclically rotate variables in the average. In each term of the symmetrization, let us put a certain operator, for example $`P_{i_k,j_k}`$, at the leftmost. Then, $`k`$ terms give the same ordering. Namely,

$`{\displaystyle \underset{\sigma }{}}P_{i_{\sigma (1)},j_{\sigma (1)}}P_{i_{\sigma (2)},j_{\sigma (2)}}\mathrm{}P_{i_{\sigma (k)},j_{\sigma (k)}}_0=`$
$`kP_{i_k,j_k}{\displaystyle \underset{\stackrel{~}{\sigma }}{}}P_{i_{\stackrel{~}{\sigma }(1)},j_{\stackrel{~}{\sigma }(1)}}P_{i_{\stackrel{~}{\sigma }(2)},j_{\stackrel{~}{\sigma }(2)}}\mathrm{}P_{i_{\stackrel{~}{\sigma }(k1)},j_{\stackrel{~}{\sigma }(k1)}}_0,`$
where $`\stackrel{~}{\sigma }`$ is a permutation of $`k1`$ elements. The factor $`k`$ cancels out $`k`$ in $`k!`$ in Eq. (17). That is, the symmetrization of all the operators except one is equivalent to the symmetrization of all the operators. Therefore, $`P_{i,j}^m_0`$ is also already symmetrized, and thus

$$P_{i,j}e^\beta _\mathrm{s}=P_{i,j}e^\beta _0.$$
(98)
## Appendix D proofs for cumulants

Here, we prove properties of cumulants in Fig. 2. For a simplicity, we use the SU(2) notation here. However, for general SU($`n`$), one can easily replace spin operators below by $`𝒳_i^{\alpha \beta }`$, $`𝒴_i^{\alpha \beta }`$, $`D_i^{(\alpha )}`$ in Appendix A.

In order to consider Fig. 2(c) and (d). Let us expand the l.h.s. of Eq. (26) using moments, namely, in terms of Eq. (94). Then, let us rewrite the $`p_{i,j}`$ operators using spin operators in terms of Eq. (7). If a site index $`i`$ appears only once in $`[\mathrm{}]`$, every term in the expansion is equal to zero because the average at site $`i`$ is evaluated as $`s_i^x_0=s_i^y_0=s_i^z_0=0.`$ Consequently, the cumulant is equal to zero.

Note that, in contrast to the Ising model, a moment bond-diagram in which an odd number of bonds meet at a site can give nonzero value because a product of three operators can give nonzero average, for example, $`s_i^xs_j^xs_i^ys_j^ys_i^zs_j^z_0`$.

Next, we prove Fig. 2(e). After Sec. IV.4.1 is explained, this property can be easily proved. Let us think about the mixed expansion pivoting on the bond in the center. Then, $`C_1`$ is equal to zero because this middle bond always unites two cycles. Furthermore, $`C_2`$ in Eq. (36) is equal to zero because every cumulant in it is equal to zero.

## Appendix E simple FCM formulae in one dimension

In one dimension, the equations in Sec. III.4 can be reduced to simpler forms as shown below. From Eq. (29), we can derive

$`p_{i,j}_{\mathrm{}}^{}`$ $`=`$ $`p_{i,j}_{\mathrm{}}p_{i,j}_\mathrm{}1p_{i1,j1}_\mathrm{}1`$ (99)
$`+p_{i1,j1}_\mathrm{}2.`$
In other words, Eq. (29) is inverted into Eq. (99) using full perimeter lattice constants explained in Ref. DombDomb, . Combining Eqs. (30) and (99), contributions from small clusters are canceled out and we obtain a simple formula,

$$p_{i,i+x}_{\mathrm{}}^{}=\underset{j=1}{\overset{\mathrm{}x}{}}p_{j,j+x}_{\mathrm{}}\underset{j=1}{\overset{(\mathrm{}1)x}{}}p_{j,j+x}_\mathrm{}1.$$
(100)
There are formulae also for the free energy and we have used it for calculating the specific heat of the SU($`\mathrm{}`$) limit in Ref. FK1d, . The free energy from the interaction part for $`_{\mathrm{}}^{(p)}`$ is $`f_{\mathrm{}}:=\beta ^1\mathrm{ln}\mathrm{exp}(\beta _{\mathrm{}}^{(p)})_0`$. The net contribution from the $`\mathrm{}`$-site cluster is defined as,

$$f_{\mathrm{}}^{}:=f_{\mathrm{}}\underset{\mathrm{}_1=1}{\overset{\mathrm{}1}{}}(\mathrm{}_1+1)f_\mathrm{}\mathrm{}_1^{}.$$
(101)
The total net contribution from clusters smaller than or equal to $`\mathrm{}`$ is defined as

$$f_{\mathrm{}}^{}:=\underset{l=1}{\overset{\mathrm{}}{}}f_l^{}.$$
(102)
The simplified version of these relations also exists. Equation (101) is rewritten asDombDomb

$$f_{\mathrm{}}^{}=f_{\mathrm{}}2f_\mathrm{}1+f_\mathrm{}2.$$
(103)
Finally, Eq. (102) is reduced to

$$f_{\mathrm{}}^{}=f_{\mathrm{}}f_\mathrm{}1.$$
(104)
Here, $`f_{\mathrm{}}^{}`$ is correct up to $`O[(\beta J)^{2\mathrm{}1}]`$. Note that $`f_{\mathrm{}}^{}=F_{\mathrm{int}}/N`$, where $`F_{\mathrm{int}}:=T\mathrm{ln}\mathrm{exp}(\beta ^{(p)})_0`$. In fact, Eq. (104) is reduced to $`f_{\mathrm{}}^{}=\mathrm{d}F_{\mathrm{int}}/\mathrm{d}N`$ in the $`\mathrm{}\mathrm{}`$ limit.

## Appendix F  an example of the trace

Let us consider a five-site system. Then,

$$\mathrm{Tr}^{(5)}\mathrm{\hspace{0.33em}1}=\underset{\alpha _1=1}{\overset{n}{}}\mathrm{}\underset{\alpha _5=1}{\overset{n}{}}\alpha _1\alpha _2\alpha _3\alpha _4\alpha _5|\alpha _1\alpha _2\alpha _3\alpha _4\alpha _5=n^5.$$
(105)
We show the calculation for $`P=P_{1,4}P_{3,5}P_{1,2}`$ below. Here, $`P_{i,j}`$ exchanges $`i`$-th state and $`j`$-th state.

$`\mathrm{Tr}^{(5)}P`$ (106)
$`=`$ $`{\displaystyle \underset{\alpha _1=1}{\overset{n}{}}}\mathrm{}{\displaystyle \underset{\alpha _5=1}{\overset{n}{}}}\alpha _1\alpha _2\alpha _3\alpha _4\alpha _5|P|\alpha _1\alpha _2\alpha _3\alpha _4\alpha _5`$
$`=`$ $`{\displaystyle \underset{\alpha _1=1}{\overset{n}{}}}\mathrm{}{\displaystyle \underset{\alpha _5=1}{\overset{n}{}}}\alpha _1\alpha _2\alpha _3\alpha _4\alpha _5|\alpha _4\alpha _1\alpha _5\alpha _2\alpha _3`$
$`=`$ $`{\displaystyle \underset{\alpha _1=1}{\overset{n}{}}}\mathrm{}{\displaystyle \underset{\alpha _5=1}{\overset{n}{}}}\delta _{\alpha _1,\alpha _4}\delta _{\alpha _2,\alpha _1}\delta _{\alpha _3,\alpha _5}\delta _{\alpha _4,\alpha _2}\delta _{\alpha _5,\alpha _3}`$
Let us consider a condition to give a nonzero contribution. Here, $`\delta _{\alpha _1,\alpha _4}`$ requires $`\alpha _1=\alpha _4`$, then $`\delta _{\alpha _4,\alpha _2}`$ requires $`\alpha _4=\alpha _2`$, then $`\delta _{\alpha _2,\alpha _1}`$ requires $`\alpha _2=\alpha _1`$. These are reduced to $`\alpha _1=\alpha _4=\alpha _2`$. In fact, this procedure is equivalent to finding a cyclic permutation $`1421`$. By the same argument, variables in another cyclic permutation $`353`$ also should be equal to each other, namely, $`\alpha _3=\alpha _5`$. Hence, there are only two independent variables in the multiple summation above. Accordingly, we obtain

$$\mathrm{Tr}^{(5)}P=n^2.$$
(107)
In other words, permutation $`(12345)(41523)`$ is composed of cycles $`(142)(421),(35)(53)`$, and then, the number of independent variables is given by the number of cycles, namely, two.

Let us operate one more exchange operator to $`P`$. When $`P_{1,2}`$ is operated, both 1 and 2 belong to a cycle $`(142)`$, and this operation breaks this cycle. That is, $`P_{1,2}P`$ has three cycles (1)(24)(35). In the case of $`P_{1,3}P`$, however, 1 belongs to (142), and 3 belongs to (35), and these cycles are united, namely, $`P_{1,3}P`$ has only one cycle (15342).

## Appendix G another definition of quasi-moments/cumulants

In fact, a quasi-moment can be defined without a prefactor in Eqs. (53) and (54). In this case, Fig. 12(b) is equal to a quasi-moment. The generating function for this quasi-moments is defined by a derivative as,

$`\stackrel{~}{𝒢}_{\mathrm{qm}}`$ $`:=`$ $`{\displaystyle \frac{}{\lambda _{i,i+1}}}{\displaystyle \frac{}{\lambda _{i+1,i+2}}}\mathrm{}{\displaystyle \frac{}{\lambda _{j1,j}}}`$ (108)
$`\times {\displaystyle \frac{p_{i,j}}{[P_{i,j}P_{i,i+1}P_{i+1,i+2}\mathrm{}P_{j1,j}]_\mathrm{s}}}.`$
Consequently, the generating function for the corresponding quasi-cumulants is written as,

$$\stackrel{~}{𝒢}_{\mathrm{qc}}:=\mathrm{ln}\stackrel{~}{𝒢}_{\mathrm{qm}}.$$
(109)
However, we have to integrate the quasi-moment generating function $`ji1`$ times to come back to the correlation function. Hence, it is difficult to calculate general dependence on $`ji`$, and thus this definition is not suitable for our purpose. This is the reason why we have adopted Eq. (51).

## Appendix H amida-diagrams

First, let us think about moment amida-diagrams. In order to give a nonzero contribution, a configuration at a nearest-neighbor pair is topologically equivalent to one of the two types shown in Fig. 28. A pair of $`P_{x,y}`$ can be simplified by exchanging subscripts $`x`$ and $`y`$ of operators between them. We can simplify left and right pairs of two bonds to give the r.h.s. of Fig. 28. In the parentheses, we schematically represent which sites are included in which cycle; sites belonging to a cycle are connected by a line.

In the moment versions of amida-diagrams in Fig. 24, we simplify bonds at every other intersite, namely, $`P_{1,2}`$, $`P_{3,4}`$, …. The simplified diagrams are made of the pieces in Fig. 28 and special pieces for the ends as shown in Fig. 29. When the end pieces are absent, any combination of the pieces in Fig. 28 gives two cycles. When $`\mathrm{}`$ is an even number, the end pieces are simple vertical lines, and do not do anything to cycles. Hence, the diagram has two cycles and gives $`n^2/n^{\mathrm{}}`$. When $`\mathrm{}`$ is an odd number, one of the end pieces is not a simple vertical line as shown in Fig. 29(b), and it unites two cycles into one. Therefore, the diagram has only one cycle and gives $`n/n^{\mathrm{}}`$. The important point is that the amida-diagrams of odd $`\mathrm{}`$ are equal to those of $`\mathrm{}+1`$.

Then, we can prove Fig. 24 using a mathematical induction. In Fig. 30, we explicitly calculates the lowest order. Next, let us consider the mixed expansion in Fig. 25 for the higher orders. Suppose $`\mathrm{}`$ is an odd number. The first and second terms cancel out each other by the same argument as noted in the context of Eq. (39) and Fig. 15. If Fig. 24 is true for the lower orders, the sum of the third and fourth terms, the sum of the fifth and sixth terms, …, are equal to zero. As a result, only the last term survives. Since the moment amida-diagram in it is equal to unity, we obtain $`(1n^2)`$. Next, suppose $`\mathrm{}`$ is an even number. The sum of the first and the second terms is $`\frac{n^2}{n^{\mathrm{}}}(1\frac{1}{n^2})`$. The third term is $`(1\frac{1}{n^2})`$ times $`\frac{n}{n^\mathrm{}1}`$. Hence, the first, second and third terms cancel out. Furthermore, if Fig. 24 is true for the lower order, the sum of the fourth and fifth terms, the sum of the sixth and seventh terms, …, are equal to zero. Namely, again, only the last term survives. Since the moment amida-diagram in it is equal to unity, we obtain $`(1n^2)`$.