
\subsection{Energy gap, critical point and scaling with system size}


\begin{figure}[h]
\includegraphics[trim=0cm 3cm 0cm 3cm,clip=true,width=0.5\textwidth]{supp1.pdf}
\caption{\label{suppfig1} (a) Energy gap $\Delta$ between ground and first excited state for a 10-spin chain. The energy gap decreases as $\Delta=1-g/g_c(p)+O(g^2)$, where $g_c(p)$ is the critical interaction strength separating the paramagnetic and ferromagnetic phases. In the finite-sized system, we define $g_c(p)$ to be the point where $\partial^2\Delta/\partial g^2$, shown in (b), is a maximum. Inset in (a) shows the dependence of $\Delta$ on $N$ for $g=0.4\omega$.}
\end{figure}
The energy gap $\Delta$ between the ground and first excited state is shown in Fig.~\ref{suppfig1}(a). The energy gap decreases approximately linearly with increasing $g/g_c(p)$. In an infinite system with nearest neighbour interactions, the exact energy gap is $\Delta=|\omega-g|$~\cite{pfeuty1970} and hence $\partial^2\Delta/\partial g^2=\delta(\omega-g)$. Finite-size effects regularize the divergence of $\partial^2\Delta/\partial g^2$, however we still observe a clear peak at a critical value of $g$, see Fig.~\ref{suppfig1}(b). We define $g_c(p)$ to be the value of $g$ corresponding to this peak. The critical value is close to $\omega$ for $p=\infty$ and decreases for decreasing $p$. The critical value will depend on $\omega$. In the main text, we fix $g_c(p)$ by defining this to be the critical value at $\omega=\omega_0$.

The dependence of work and efficiency on chain size $N$ is shown in Fig.~\ref{suppfig2}(a) for $g=0.4\omega_0$. For $p=\infty,3$ the performance scales extensively for $N\gtrsim 6$. Finite size effects are slightly more pronounced for $p=2$, however the scaling is close to extensive. For $p=1$ the performance scales super-extensively with $g$. This is due to a dependence of $\Delta$ on $N$, see inset to Fig.~\ref{suppfig1}(a). The critical point $g_c$ also depends on $N$ for $p=1$. In the cold, perturbative limit (see next section), $\Delta\sim \omega-g(\gamma+\ln N)$ with $\gamma$ the Euler-Mascheroni constant. The two-level approximation [Eq.~4 in the main text] then gives $W\sim N^{1+\beta_H g}$. This approximation captures the dependence of work on particle number for small $g/g_c$, see Fig.~\ref{suppfig2}(a,ii),(b,ii). Increasing either $g$ or $N$ increases $g/g_c$ and the perturbative approximation $W\sim N^{1+\beta_H g}$ breaks down. Scaling $g$ by $g_c$ effectively Kac renormalizes the interactions, rendering the long-range system extensive~\cite{kac1963,kac1969,cannas1996,gonzalezlazo2021,zhang2017}.

\begin{figure}[H]
\includegraphics[trim=0cm 8cm 0cm 8cm,clip=true,width=\textwidth]{supp2.pdf}
\caption{\label{suppfig2} The (a,i) work output and (a,ii) efficiency scale approximately extensively for $p=2,3,\infty$ whereas they scale non-extensively for $p=1$. The non-extensive scaling for $p=1$ agrees with the perturbative prediction $W\sim N^{1+\beta_H g}$ (black dashed lines). Results for $g=0.4\omega_0$. (b) The non-extensive scaling of $W$ and $\eta$ for $p=1$ and $g/\omega_0=0.2,0.3,0.4,0.5$ along with the perturbative prediction $W\sim N^{1+\beta_H g}$ (matching dashed black lines). The gray lines in (a,ii),(b,ii) are the Carnot efficiency $\eta_\mathrm{C}$. All results are for $\beta_H=10\omega_0^{-1}$, $\beta_C=2\beta_H$ and $r=r_\mathrm{NI}^\mathrm{max}$.}
\end{figure}


\subsection{Performance of the transverse Ising model}
The dimensionless free energy (free energy divided by $k_BT$) of the transverse Ising model confined to the first $N$ excited states is
\begin{equation}
\begin{split}
\ln Z&=\sum_{k=0}^{N-1} e^{-\beta\sqrt{\omega^2+g^2-2\omega g\cos(2\pi k/N)}},\\
&\approx\frac{N}{\pi}\int_0^\pi e^{-\beta\sqrt{\omega^2+g^2-2\omega g\cos\theta}}\,d\theta.
\end{split}
\end{equation}
Here and below we work in a frame where the ground state energy is zero. In the adiabatic limit, the system energy at the four points in the cycle in Fig.~1(a) of the main text are
\begin{equation}
\begin{split}
(1)&\hspace{0.5cm}\langle E\rangle=\frac{N}{\pi}\int_0^\pi \sqrt{r^2\omega_0^2+g^2-2r\omega_0 g\cos\theta}e^{-\beta_H\sqrt{r^2\omega_0^2+g^2-2r\omega_0 g\cos\theta}}\,d\theta,\\
(2)&\hspace{0.5cm}\langle E\rangle=\frac{N}{\pi}\int_0^\pi \sqrt{\omega_0^2+g^2-2\omega_0 g\cos\theta}e^{-\beta_H\sqrt{r^2\omega_0^2+g^2-2r\omega_0 g\cos\theta}}\,d\theta,\\
(3)&\hspace{0.5cm}\langle E\rangle=\frac{N}{\pi}\int_0^\pi \sqrt{\omega_0^2+g^2-2\omega_0 g\cos\theta}e^{-\beta_C\sqrt{\omega_0^2+g^2-2\omega_0 g\cos\theta}}\,d\theta,\\
(4)&\hspace{0.5cm}\langle E\rangle=\frac{N}{\pi}\int_0^\pi \sqrt{r^2\omega_0^2+g^2-2r\omega_0 g\cos\theta}e^{-\beta_C\sqrt{\omega_0^2+g^2-2\omega_0 g\cos\theta}}\,d\theta.
\end{split}
\end{equation}
The heat and work outputs can be computed from these, see Fig.~1(b) in the main text.


\subsection{Perturbative treatment}

We assume a large chain so that boundary effects can be ignored and hence we approximate the system as being translationally invariant (equivalently we can impose periodic boundary conditions). The spin operators can be converted to expressions in terms of bosonic operators via a Holstein-Primakoff transformation~\cite{holstein1940},
\begin{equation}
\hat{\sigma}_x^{(i)}\rightarrow \frac{\left(\sqrt{1-\hat{a}_i^\dagger\hat{a}_i}\right)\hat{a}_i+\hat{a}_i^\dagger\left(\sqrt{1-\hat{a}_i^\dagger\hat{a}_i}\right)}{2},\hspace{1cm}\hat{\sigma}_z^{(i)}\rightarrow\frac{1}{2}-\hat{a}_i^\dagger\hat{a}_i.
\end{equation}
In the low excitation regime we expand the Hamiltonian [Eq.~1 in the main text] to quadratic order in the bosonic operators $\hat{a}_i$, $\hat{a}_i^\dagger$. This gives~\cite{holstein1940,vogl2020,konig2021}
\begin{equation}
\hat{H}=-\frac{\omega N}{2}+\omega\sum_{i=1}^N\hat{a}_i^\dagger\hat{a}_i-\frac{g}{4}\sum_{\substack{i,j=1\\(j\ne i)}}^NJ_{ij}\left(\hat{a}_i+\hat{a}_i^\dagger\right)\left(\hat{a}_j+\hat{a}_j^\dagger\right).
\end{equation}
Fourier transforming and carrying out a Bogoliubov transformation gives~\cite{cevolani2016}
\begin{equation}
\hat{H}=\epsilon_0+\sum_{k=0}^{N-1} \left(\omega\operatorname{sign}\left(\omega-g\tilde{J}_k\right)\sqrt{1-\frac{2g}{\omega}\tilde{J}_k}\right)b_k^\dagger b_k,
\end{equation}
with $\epsilon_0$ the ground-state energy, $b_k$ bosonic operators for the quasiparticle modes and
\begin{equation}\label{Jk}
\tilde{J}_k=\sum_{m=1}^N \frac{\cos\left(2\pi k m/N\right)}{m^p}=C_p^N\left(\frac{2\pi k}{N}\right),
\end{equation}
where $C_p^N(\theta)=\sum_{m=1}^N\cos(m\theta)/m^p$ is the finite-$N$ generalized Clausen function~\cite{wu2010}. For nearest neighbour interactions, we have $\tilde{J}_k=\cos (2\pi k/N)$.

From hereon we assume $\operatorname{sign}(\omega-g\tilde{J}_k)>0$. The thermodynamic properties of the spin chain can be determined from the partition function,
\begin{equation}
Z=\prod_{k=0}^{N-1}\frac{1}{1-\exp\left(-\beta\omega\sqrt{1-\frac{2g}{\omega}\tilde{J}_k}\right)},
\end{equation}
with logarithm (dimensionless free energy),
\begin{equation}
\begin{split}
\ln Z&=\sum_{k=0}^{N-1} \ln \left[\frac{1}{1-\exp\left(-\beta\omega\sqrt{1-\frac{2g}{\omega}\tilde{J}_k}\right)}\right],\\
&\approx\sum_{k=0}^{N-1}\exp\left(-\beta\omega\sqrt{1-\frac{2g}{\omega}\tilde{J}_k}\right).
\end{split}
\end{equation}
The latter approximation assumes low temperature. We have ignored the ground state energy, which adds an inconsequential constant to $\ln Z$. To proceed analytically, we presume small $g\tilde{J}_k/\omega$ and approximate $\sqrt{1-2g\tilde{J}_k/\omega}\approx 1-g\tilde{J}_k\omega$. Hence
\begin{equation}\label{Zapprox}
\ln Z\approx e^{-\beta\omega}\sum_{k=0}^{N-1}\exp\left(m\beta g\tilde{J}_k\right).
\end{equation}

We reserve a discussion of $p=1$ for later and for now assume $p>1$. For large $N$, $C_p^N(\theta)\rightarrow C_p(\theta)=\sum_{m=1}^\infty \cos(m\theta)/m^p$ and we can replace the sum over $k$ in Eq.~\eqref{Zapprox} by an integral
\begin{equation}
\sum_{k=0}^{N-1}\exp\left(\beta g\tilde{J}_k\right)\rightarrow \frac{N}{\pi}\int_0^\pi \exp\left(\beta g C_p(\theta)\right)\,d\theta.
\end{equation}
For even integers $p$,
\begin{equation}
C_p(\theta)=-\frac{(-1)^{p/2}(2\pi)^p}{2p!}B_p\left(\frac{\theta}{2\pi}\right),
\end{equation}
with $B_p(x)$ the Bernoulli polynomials, which are polynomials of order $p$. We obtain the following results for $p=\infty$ and $p=2$,
\begin{equation}\label{p12}
\begin{split}
\frac{1}{\pi}\int_0^\pi\exp\left(\beta g C_\infty(\theta)\right)\,d\theta&=\frac{1}{\pi}\int_0^\pi\exp\left(\beta g \cos(\theta)\right)\,d\theta=I_0(\beta g)\overset{\beta\rightarrow\infty}\sim \frac{e^{\beta g}}{\sqrt{2\pi \beta g}},\\
\frac{1}{\pi}\int_0^1\exp\left(\beta g C_2(\theta)\right)\,d\theta&=\frac{2 e^{\beta g\pi^2/6}}{\sqrt{\beta g\pi^2}}D\left(\frac{\sqrt{\beta g\pi^2}}{2}\right)\overset{\beta\rightarrow\infty}\sim \frac{e^{\beta g\pi^2/6}}{3\beta g\pi^2/6},
\end{split}
\end{equation}
where $I_0(x)$ is the modified Bessel function and $D(x)$ is Dawson's function (expressible in terms of the error function $\operatorname{erf}(x)$ via $D(x)=-(\sqrt{\pi}/2)e^{-x^2}i\operatorname{erf}(ix)$).

In general, for $p>3$ the integral $\int_0^1 \exp((m\beta gC_p(2\pi x)))\,dx$ can be approximated at low temperatures using the method of steepest descent,
\begin{equation}
\frac{1}{\pi}\int_0^\pi \exp(\beta gC_p(\theta))\,d\theta=\sqrt{\frac{1}{2\pi \beta gC_{p-2}(0)}}\exp(\beta gC_p(0)),
\end{equation}
where $C_p(0)=\sum_{m=1}^\infty m^{-p}=\zeta(p)$. This does not work for $p=2$ since $C_2'(\theta)|_{\theta=0}\ne 0$. This reflects that the spectrum is linear rather than quadratic around the lowest energy state. Nor does it work for $p=3$, since $C_3^{''}(\theta)|_{\theta=0}$ diverges. We observe numerically that $\exp(\beta gC_3(\theta))$ is dominated by its small $\theta$ behaviour for large $\beta g$. Hence for $p=3$ we expand the Clausen function in a power series around the maximum $\theta=0$,
\begin{equation}\label{p3}
\begin{split}
\frac{1}{\pi}\int_0^\pi \exp(\beta gC_3(\theta))\,d\theta &\approx \frac{\exp(\beta g\zeta(3))}{\pi}\int_0^\pi\exp\left(-\frac{\beta g \theta^2}{4}\left(3-\ln\theta^2\right)\right)\,d\theta,\\
&=\frac{\exp(\beta g\zeta(3))}{\pi}\int_0^\pi\exp\left(-\frac{\beta g \theta^2}{4}\left(3+\ln(\beta g)-\ln(\beta g\theta^2)\right)\right)\,d\theta,\\
&=\frac{\exp(\beta g\zeta(3))}{\pi\sqrt{\beta g}}\int_0^{\pi\sqrt{\beta g}}u^{u^2/2}\exp\left(-\frac{1}{4}u^2\left(\ln(\beta g)+3\right)\right)\,du,\\
&\overset{\beta\rightarrow\infty}\sim \frac{\exp(\beta g\zeta(3))}{\pi\sqrt{\beta g}}\int_0^{\pi\sqrt{\beta g}}\exp\left(-\frac{1}{4}u^2(\ln(\beta g)+3)\right)\,du,\\
&=\frac{\exp(\beta g\zeta(3))}{\sqrt{\pi \beta g(3+\ln\beta g)}}\operatorname{erf}\left(\frac{\pi}{2}\sqrt{3\beta g+\beta g\ln\beta g}\right),\\
&\overset{\beta\rightarrow\infty}\sim \frac{\exp(\beta g\zeta(3))}{\sqrt{\pi \beta g(3+\ln\beta g})}.
\end{split}
\end{equation}
Hence we find a logarithmic correction $\ln \beta g$ to the partition function for $p=3$, which interpolates between $p=2$ and $p>3$. Note that $3>2\ln \pi\approx 2.29$ and so the expansion around $\theta=0$ creates no convergence problems for large $\theta>1$. The asymptotic behaviour in the fourth line follows by replacing $u^{u^2/2}$ by $\lim_{u\rightarrow 0^+}u^{u^2/2}=1$ due to $\exp(-\frac{1}{4}u^2(\ln(\beta g)+3))$ being sharply peaked around the origin for large $\beta g$. The validity of this approximation was confirmed numerically. Equations~\eqref{p12} and \eqref{p3} give $\ln Z=N\mathcal{G}_p(\beta g)e^{-\beta \Delta}$ with $\Delta=\omega-g\zeta(p)$ and $\mathcal{G}_p(\beta g)$ arising from thermal fluctuations and dependent on $p$.


For $p=1$, Eq.~\eqref{Jk} at $k=0$ diverges with $N$ as $\tilde{J}_0\approx \ln N+\gamma$, with $\gamma$ the Euler-Mascheroni constant. We separate out this term in Eq.~\eqref{Zapprox}. For $k\ne 0$, Eq.~\eqref{Jk} converges and we can use $\sum_{m=0}^\infty \cos(m\theta)/m=C_1(\theta)=-\ln |2\sin(\theta/2)|$. Hence
\begin{equation}
\ln Z\approx Ne^{-\beta\omega}\left(e^{\beta g (\ln N+\gamma)}+2^{1-\beta g}N^{-1}\sum_{k=1}^{(N-1)/2}\left(\sin \frac{\pi k}{N}\right)^{-\beta g}\right),
\end{equation}
where we have assumed $N$ is odd for simplicity (assuming even $N$ will give the same final result below). Noting that the sum is dominated by small $k$ terms, we can use the small-angle approximation $\sin(\theta)\approx\theta$,
\begin{equation}
\begin{split}
\ln Z&\approx Ne^{-\beta\omega}\left(e^{\beta g(\ln N+\gamma)}+\frac{1}{\pi}\left(\frac{N}{2\pi}\right)^{\beta g-1}\sum_{k=1}^{(N-1)/2}k^{-\beta g}\right),\\
&\approx Ne^{-\beta(\omega-g(\ln N+\gamma))}\left(1+\frac{2}{N}\left(\frac{1}{2\pi e^\gamma}\right)^{\beta g}\zeta\left(\beta g\right)\right),\\
&\approx Ne^{-\beta(\omega-g(\ln N+\gamma))},\\
&=e^{-\beta \omega}e^{\beta g\gamma}N^{1+\beta g}.
\end{split}
\end{equation}
Hence for $p=1$ fluctuations are suppressed for low temperatures and large $N$, and the dimensionless free energy is $\ln Z\approx Ne^{-\beta\Delta}$ with $\Delta\approx \omega-g(\ln N+\gamma)$. The low-temperature work output within this approximation is then (see Fig.~\ref{suppfig2}),
\begin{equation}\label{Wapprox}
W\approx N\omega_0(r-1)\left(e^{-\beta_H\Delta}-e^{-\beta_C\Delta}\right)\approx \omega_0(r-1)e^{-\beta_H\omega_0}e^{\beta_Hg\gamma}N^{1+\beta_Hg}.
\end{equation}


Summarising, to lowest order in $g/\omega$ and for large $N$ and $\beta$, we obtain the dimensionless free energies,
\begin{equation}
\ln Z\approx \left\{\begin{array}{ll}\sqrt{\frac{1}{2\pi\beta g\zeta(p-2)}}e^{-\beta\Delta},&\hspace{0.5cm}p>3,\\
\sqrt{\frac{1}{\pi\beta g(3+\ln\beta g)}}e^{-\beta\Delta},&\hspace{0.5cm}p=3,\\
\frac{1}{3\beta g\zeta(2)}e^{-\beta\Delta},&\hspace{0.5cm}p=2,\\
e^{-\beta\Delta},&\hspace{0.5cm}p=1,\end{array}\right.
\end{equation}
with $\zeta(s)=\sum_{m=1}^\infty 1/m^s$ the Riemann-zeta function and $\Delta(\omega)=\omega(1-\omega_0 g/\omega g_c)$.



\subsection{Calculation of entanglement}
Within the quadratic approximation in the previous section, the first excited state is a state with one excitation uniformly spread across all spins,
\begin{equation}
\ket{1}=\frac{1}{\sqrt{N}}\sum_{i=1}^N\hat{\sigma}_+^{(i)}\ket{0}.
\end{equation}
We consider a partition dividing the chain in half. The entanglement entropy with respect to this partition and state $\ket{1}$ is
\begin{equation}
\mathcal{S}=-\operatorname{Tr}_R(\rho_L\ln\rho_L),
\end{equation}
with $\rho_L=\operatorname{Tr}_L(\ket{1}\bra{1})$ and $\operatorname{Tr}_{R(L)}$ a partial trace over states in the right(left) half of the chain. It is straightforward to show that this gives $\mathcal{S}=\ln (N/2)$.

For low temperatures, we can qualitatively describe the system by the thermal state
\begin{equation}
\rho_T^{01}=\frac{1}{1+e^{-\beta\Delta}}\left(\ket{0}\bra{0}+e^{-\beta\Delta}\ket{1}\bra{1}\right),
\end{equation}
where $\Delta=\omega(1-g/g_c)$ within the quadratic approximation. The Peres-Horodecki criterion states that a separable density matrix has positive partial trace~\cite{peres1996,horodecki1996}. Hence, if the density matrix does not have positive partial trace, the state is entangled (``non-PPT entanglement''). Due to the symmetry of the partition and the state, a sufficient condition for $\rho_T^{01}$ to exhibit non-PPT entanglement is $\operatorname{Tr}(\hat{M}\otimes\hat{M}\rho_T^{01})<0$, with $\hat{M}$ any Hermitian operator acting on either the left ($\hat{M}\otimes$) or right ($\otimes \hat{M}$) partition of the chain~\cite{toth2009}. We choose $\hat{M}\otimes\hat{M}=\prod_{i=1}^N\hat{m}_i$ with
\begin{equation}
\hat{m}_i=\frac{1}{(1+\alpha)^{1/N}}\left[\frac{\hat{\mathcal{I}}}{2}-\hat{\sigma}_z^{(i)}-\alpha\left(\frac{\hat{\mathcal{I}}}{2}+\hat{\sigma}_z^{(i)}\right)\right],
\end{equation}
with $\hat{\mathcal{I}}$ the identity operator and $\alpha>0$ a parameter to be chosen. This gives
\begin{equation}
\operatorname{Tr}(\hat{M}\otimes\hat{M}\rho_T^{01})=\frac{1-\alpha e^{-\beta\Delta}}{(1+\alpha)(1+e^{-\beta\Delta})}.
\end{equation}
We can make $\operatorname{Tr}(\hat{M}\otimes\hat{M}\rho_T^{01})$ arbitrarily close to $-(1+e^{\beta\Delta})^{-1}$ by choosing $\alpha$ to be large ($\alpha\gg e^{\beta\Delta}$). This choice of $\hat{M}$ gives negative $\operatorname{Tr}(\hat{M}\otimes\hat{M}\rho_T^{01})$ and hence the state exhibits non-PPT entanglement.

\subsection{Calculation of counterdiabatic drive}
Here we derive Eq.~6 in the main text. We find the coefficients $C_{ij}$ in $\hat{H}_\mathrm{cd}=\sum_{ij(j\ne i)}C_{ij}\hat{\sigma}_x^{(i)}\hat{\sigma}_y^{(j)}$ by minimising the action
\begin{equation}
S=\operatorname{Tr}(G(\hat{H}_\mathrm{cd})^2),
\end{equation}
with $G=\frac{\partial \hat{H}}{\partial t}+i[\hat{H}_\mathrm{cd},\hat{H}]$~\cite{sels2017,kolodrubetz2017}. Hence we need to solve
\begin{equation}
\frac{\partial S}{\partial C_{ij}}=2\operatorname{Tr}\left(G\frac{\partial G}{\partial C_{ij}}\right)=0.
\end{equation}
For notational simplicity, we set $J_{ii}=C_{ii}=0$. We have
\begin{equation}
\begin{split}
G&=-\omega'(t)\sum_p\sigma_z^p-i\omega(t)\sum_{i,j,p}C_{ij}[\sigma_x^i\sigma_y^j,\sigma_z^p]-ig\sum_{i,j,p,q}C_{ij}J_{pq}\sigma_x^i[\sigma_y^j,\sigma_x^p\sigma_x^q],\\
&=-\omega'(t)\sum_p\sigma_z^p+\omega(t)\sum_{i,j,p}C_{ij}(\sigma_x^i\sigma_x^j\delta_{pj}-\sigma_y^i\sigma_y^j\delta_{pi})-g\sum_{i,j,p,q}C_{ij}J_{pq}\sigma_x^i(\sigma_z^p\sigma_x^q\delta_{jp}+\sigma_x^p\sigma_z^q\delta_{jq}),\\
&=-\omega'(t)\sum_p\sigma_z^p+\omega(t)\sum_{i,j}C_{ij}(\sigma_x^i\sigma_x^j-\sigma_y^i\sigma_y^j)-2g\sum_{i,p,q}J_{pq}C_{ip}\sigma_x^i\sigma_x^q\sigma_z^p.
\end{split}
\end{equation}
Hence
\begin{equation}
\frac{\partial G}{\partial C_{mn}}=\omega(t)(\sigma_x^m\sigma_x^n-\sigma_y^m\sigma_y^n)-2g\sum_{\ell\ne n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n.
\end{equation}
We now want to calculate the trace of $G\partial G/\partial C_{mn}$ ($m\ne n$). Considering the three terms in $G$ and the two terms in $\partial G/\partial C_{mn}$ gives a total of six terms. However three of these are zero, since the trace of terms with an odd number of spin-1/2 operators is zero. The three remaining terms are as follows.

Firstly,
\begin{equation}
\operatorname{Tr}\left[g\omega'(t)\sum_p\sigma_z^p\sum_{\ell}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right]=2^N\frac{g\omega'(t)J_{mn}}{16}.
\end{equation}


Next,
\begin{equation}
\operatorname{Tr}\left[\omega(t)^2\sum_{p,q}C_{pq}(\sigma_x^p\sigma_x^q-\sigma_y^p\sigma_y^q)(\sigma_x^m\sigma_x^n-\sigma_y^m\sigma_y^n)\right]=2^N\frac{\omega(t)^2C_{mn}}{4}.
\end{equation}


Finally,
\begin{equation}
\begin{split}
&\operatorname{Tr}\left[g^2\sum_{p,q,r(q,r\ne p)}J_{pq}C_{pr}\sigma_x^r\sigma_x^q\sigma_z^p\sum_{\ell\ne n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right],\\
&=\operatorname{Tr}\left[g^2\sum_{p,q,r(q,r\ne p)}J_{pq}C_{pr}\sigma_x^r\sigma_x^q\sigma_z^p\sum_{\ell\ne n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right],\\
&=\operatorname{Tr}\left[g^2\sum_{p,q}J_{pq}C_{pq}(\sigma_x^q)^2\sigma_z^p\sum_{\ell\ne n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right]+\operatorname{Tr}\left[g^2\sum_{p,q,r(\mathrm{distinct})}J_{pq}C_{pr}\sigma_x^r\sigma_x^q\sigma_z^p\sum_{\ell\ne n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right],\\
&=2^N \frac{g^2\sum_q J_{nq}C_{nq}J_{nm}}{64}+\operatorname{Tr}\left[g^2\sum_{p,q,r(\mathrm{distinct})}J_{pq}C_{pr}\sigma_x^r\sigma_x^q\sigma_z^p\sum_{\ell\ne m,n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right],\\
&=2^N \frac{g^2\sum_q J_{mn}J_{nq}C_{nq}}{64}+2^N\frac{g^2\sum_{\ell\ne m,n}\left(J_{n\ell}^2C_{mn}+J_{mn}J_{n\ell}C_{n\ell}\right)}{64},\\
&=2^N \frac{g^2(J_{mn}(JC)_{nn}+(J^2)_{nn}C_{mn}/2-J_{mn}^2C_{mn})}{32}.
\end{split}
\end{equation}


Combining terms gives the coupled linear equations that determine $C_{mn}$,
\begin{equation}
g\omega'(t)J_{mn}+2\omega(t)^2C_{mn}+g^2(J_{mn}(JC)_{nn}+(J^2)_{nn}C_{mn}/2-J_{mn}^2C_{mn})=0.
\end{equation}
This gives
\begin{equation}\label{Cmn}
C_{mn}=-\frac{g\omega'(t)J_{mn}+g^2J_{mn}(JC)_{nn}}{2\omega(t)^2+g^2(J^2)_{nn}/2-g^2J_{mn}^2}.
\end{equation}
Multiplying both sides by $g^2J_{nm}$ and summing over $m$ gives
\begin{equation}\label{Cnn}
g^2(JC)_{nn}=-(g\omega'(t)+g^2(JC)_{nn})f_n,
\end{equation}
with
\begin{equation}
f_n=\sum_m\frac{g^2J_{mn}^2}{2\omega(t)^2+g^2\left[(J^2)_{nn}/2-J_{mn}^2\right]}.
\end{equation}
Equation~\eqref{Cnn} can be rearranged to give
\begin{equation}
g^2(JC)_{nn}=-\frac{g\omega'(t)f_n}{1+f_n}.
\end{equation}
Substituting this into Eq.~\eqref{Cmn} gives
\begin{equation}
C_{mn}=-\frac{g\omega'(t)\left(\frac{1}{1+f_n}\right)J_{mn}}{2\omega(t)^2+g^2\left[(J^2)_{nn}/2-J_{mn}^2\right]}.
\end{equation}
For an infinite chain with $J_{mn}=1/|m-n|^p$, $(J^2)_{nn}/2=\zeta(2p)$ with $\zeta(s)$ the Riemann zeta function and $f_n=f$ is independent of $n$. Furthermore, we find that for integer values of $p$ we can evaluate $f$ analytically using Mathematica (for example, for $p=1$, $f=1-x\cot x$ with $x=g\pi/\sqrt{g^2\zeta(2)+64\omega(t)^2}$). For nearest neighbour interactions, $f=g^2/\omega(t)^2$, and
\begin{equation}
C_{mn}=-\frac{g\omega'(t)}{2(\omega(t)^2+g^2)}\delta_{n,m+1}.
\end{equation}





