Solution (source code)

= Solution

Take the couplings in the <Hamiltonian> to include <inverse temperature>, so the <Boltzmann factor> is $e^{-H}$. If physical energies are used instead, first replace $J,L,D,K$ by their products with $1/(k_BT)$. Split the on-site term equally between its two bonds. For spin states $s,t\in\{1,0,-1\}$ the <spin-chain transfer matrix> has entries
$$
W_{st}=\exp\left[Jst+Ls^2t^2+\frac D2(s^2+t^2)+K\right].
$$
Using the specified coupling coordinates, in the order $(1,0,-1)$ this is
$$
\boxed{W=\frac{e^K}{z}\begin{pmatrix}1&x&y\\x&z&x\\y&x&1\end{pmatrix}.}
$$
For example, $e^{D/2}=x/z$ and $e^{-J+L+D}=y/z$. This is the nearest-neighbour <Blume–Emery–Griffiths model> with an additive constant and the paper's sign convention for $D$.

Summing the periodic spin chain gives the <partition function> $Z_N=\operatorname{tr}W^N=\sum_{a=1}^3\lambda_a^N$. For finite real couplings all entries of $W$ are strictly positive; the <Perron–Frobenius theorem> gives a unique positive <dominant eigenvalue> with $\lambda_1>|\lambda_a|$ for $a\ne1$. Since $W$ is a <symmetric matrix>, all its <eigenvalues> are real. Thus in the <thermodynamic limit>,
$$
\boxed{\frac{F_N}{Nk_BT}\longrightarrow-\log\lambda_1.}
$$
The positivity condition is stronger and more useful than mere ordering by signed value: subdominant <eigenvalues> can be negative. Also, the printed strict ordering between the other two is not guaranteed for all couplings. For instance, $J=L=D=0$ makes $W$ the positive constant <matrix> $e^K\mathbf1\mathbf1^T$ with two equal zero <eigenvalues>. Degeneracy there does not affect the largest-eigenvalue limit; no explicit generic <eigenvalues> are needed.

For <spin magnetization>, introduce a dimensionless field $h$ through $-h\sum_i\sigma_i$, or insert $S=\operatorname{diag}(1,0,-1)$ into the <trace>. Then
$$
\langle\sigma_i\rangle=\frac{\operatorname{tr}(SW^N)}{\operatorname{tr}W^N}\longrightarrow v_1^TSv_1=\left.\partial_h\log\lambda_1(h)\right|_{h=0},
$$
where $v_1$ is a normalized <eigenvector> of the largest <eigenvalue> and $W_{st}(h)=e^{h(s+t)/2}W_{st}(0)$. <Spin inversion symmetry> gives $v_{1,+}=v_{1,-}$, since the positive <eigenvector> of the largest <eigenvalue> is unique. Therefore
$$
\boxed{\langle\sigma_i\rangle=0}
$$
at zero field, both at finite $N$ and in the finite-coupling <thermodynamic limit>. The finite-$N$ result follows directly by pairing each configuration with its spin-reversed partner. <Eigenvalues> at a single fixed field do not determine a general observable: a field derivative of the largest <eigenvalue>, or its <eigenvector>, is required. Positivity and the <real analytic> dependence on the couplings exclude a finite-temperature <spontaneous symmetry breaking> transition in this one-dimensional finite-range chain; singular zero-temperature coupling limits require separate treatment.

For even $N$, <spin decimation> on alternate sites sums the middle spin of each two-bond segment, hence the coarse <spin-chain transfer matrix> is $W'=W^2$. Define
$$
A=1+x^2+y^2,\qquad B=x(1+y+z),\qquad C=x^2+2y,\qquad E=z^2+2x^2.
$$
Direct multiplication gives $W^2=e^{2K}z^{-2}\begin{pmatrix}A&B&C\\B&E&B\\C&B&A\end{pmatrix}$. Matching its entry ratios to the original parameterization yields the <spin-1 chain decimation recursion>
$$
\boxed{x'=\frac{x(1+y+z)}{1+x^2+y^2},\qquad y'=\frac{x^2+2y}{1+x^2+y^2},\qquad z'=\frac{z^2+2x^2}{1+x^2+y^2}.}
$$
The remaining overall positive factor is absorbed into $K'$. Keeping that factor preserves the <free energy> as well as normalized spin <probabilities>. Indeed $\operatorname{tr}(W^2)^{N/2}=\operatorname{tr}W^N$, exactly; on an odd ring an unmatched boundary segment needs separate handling rather than assuming a uniform two-site block decomposition.