Quantum kinetic perturbation theory for near-integrable spin chains with weak long-range interactions

TL;DR

This paper introduces a perturbative quantum kinetic approach to study the relaxation dynamics of weakly long-range interacting spin chains, enabling high-accuracy simulations of spin correlations and thermalization processes.

Contribution

It develops a novel perturbative scheme based on quantum kinetic equations for near-integrable spin chains with long-range interactions, including multiple truncation levels and benchmarking.

Findings

The method accurately captures prethermalization and thermalization dynamics.

It allows simulations of systems with around 100 sites for moderate times.

The approach can be extended to other weakly non-integrable spin models.

Abstract

For a transverse-field Ising chain with weak long-range interactions we develop a perturbative scheme, based on quantum kinetic equations, around the integrable nearest-neighbour model. We introduce, discuss, and benchmark several truncations of the time evolution equations up to eighth order in the Jordan-Wigner fermionic operators. The resulting set of differential equations can be solved for lattices with $O(10^2)$ sites and facilitates the computation of spin expectation values and correlation functions to high accuracy, at least for moderate timescales. We use this scheme to study the relaxation dynamics of the model, involving prethermalisation and thermalisation. The techniques developed here can be generalised to other spin models with weak integrability-breaking terms.

Equations

$$\mathcal{H}={\mathcal{H}}_{\text{int}}+{\mathcal{H}}_{\text{pert}},$$

$${\mathcal{H}}_{\text{int}}=J_x\sum _l{\mathcal{S}}_l^x{\mathcal{S}}_{l+1}^x+h\sum _l{\mathcal{S}}_l^z$$

$${\mathcal{H}}_{\text{pert}}=\frac{J_z}{2}\sum _{l,m}\frac{1}{d(m{)}^{\alpha }}{\mathcal{S}}_l^z{\mathcal{S}}_{l+m}^z$$

$${\mathcal{H}}_{\text{int}}=\sum _k{ϵ}_k\left({\eta }_k^{†}{\eta }_k-\frac{1}{2}\right)$$

$${ϵ}_k=\sqrt{h^2+h{J}_x\cos k+J_x^2/4},$$

$$\mathcal{H}=$$

$$+\sum _{\mathbf{k}}\left(B_{\text{mbox}I}(\mathbf{k}){\eta }_{-k_1}^{}{\eta }_{k_2}^{}{\eta }_{-k_3}^{}{\eta }_{k_4}^{}+B_{\text{mbox}II}(\mathbf{k}){\eta }_{k_1}^{†}{\eta }_{k_2}^{}{\eta }_{-k_3}^{}{\eta }_{k_4}^{}+...+B_{\text{mbox}V}(\mathbf{k}){\eta }_{k_1}^{†}{\eta }_{-k_2}^{†}{\eta }_{k_3}^{†}{\eta }_{-k_4}^{†}\right)$$

$$i{d}_t\mathcal{O}=\left[\mathcal{O},\mathcal{H}\right],$$

$$i{d}_t{\eta }_k^{†}{\eta }_k^{}=$$

$$+B_{\text{mbox}II}(\mathbf{k}){\Delta }_{\text{mbox}II}(\mathbf{k}){\eta }_{k_1}^{†}{\eta }_{k_2}^{}{\eta }_{-k_3}^{}{\eta }_{k_4}^{}+...+B_{\text{mbox}V}(\mathbf{k}){\Delta }_{\text{mbox}V}(\mathbf{k}){\eta }_{k_1}^{†}{\eta }_{-k_2}^{†}{\eta }_{k_3}^{†}{\eta }_{-k_4}^{†}),$$

$$\mathrm{C}_{2}^{1}=\bigl{\{}\eta^{{\dagger}}_{k}\eta_{k^{\prime}}^{\phantom{\dagger}}~{}|~{}\forall k,k^{\prime}\in\text{Br}\bigr{\}}\equiv\bigl{\{}\eta^{{\dagger}}\eta\bigr{\}},$$

$${\mathrm{C}}^p:=\bigcup _{\text{deg}\in [[p,2p]]}{\mathrm{C}}_{\text{deg}}^p,{\mathrm{C}}_{\text{deg}}:=\bigcup _{p\in [[\lfloor \frac{\text{deg}+1}{2}\rfloor ,\text{deg}]]}{\mathrm{C}}_{\text{deg}}^p.$$

$$\mathrm{T}=\mathrm{C}_{0}\cup\mathrm{C}_{1}\cup\mathrm{C}_{2}\equiv\bigl{\{}\mathbbm{1},~{}\eta,~{}\eta^{\dagger},~{}\eta^{{\dagger}}\eta,~{}\eta\eta,~{}\eta^{\dagger}\eta^{\dagger}\bigr{\}}$$

$$|\downarrow \cdots \downarrow ⟩={\mathcal{G}}_{\lfloor N/2\rfloor }|0⟩,$$

$${\mathcal{G}}_n:=\frac{1}{W_n}\left(1+\sum _{s=1}^n(-i{)}^s\sum _{0<k_1<\cdots <k_s<\pi }\frac{v_{k_1}}{u_{k_1}}\cdots \frac{v_{k_s}}{u_{k_s}}{\eta }_{-k_s}^{†}\cdots {\eta }_{-k_1}^{†}{\eta }_{k_1}^{†}\cdots {\eta }_{k_s}^{†}\right)$$

$$W_n^2=1+\sum _{s=1}^n\sum _{0<k_1<\cdots <k_s<\pi }{\left(\frac{v_{k_1}}{u_{k_1}}\cdots \frac{v_{k_s}}{u_{k_s}}\right)}^2.$$

$$|{\psi }^n⟩={\mathcal{G}}_n|0⟩$$

$$\bigl{\langle}\eta^{{\dagger}}_{-k}\eta^{{\dagger}}_{k}\eta_{-k^{\prime}}^{\phantom{{\dagger}}}\eta_{k^{\prime}}^{\phantom{{\dagger}}}\bigr{\rangle}-\bigl{\langle}\eta^{{\dagger}}_{-k}\eta^{{\dagger}}_{k}\bigr{\rangle}\bigl{\langle}\eta_{-k^{\prime}}^{\phantom{{\dagger}}}\eta_{k^{\prime}}^{\phantom{{\dagger}}}\bigr{\rangle}=\frac{v_{k}v_{k^{\prime}}}{u_{k}u_{k^{\prime}}}\left(W_{1}^{-4}-W_{1}^{-2}\right),$$

$$|\xi ⟩,|{\xi }^{\prime }⟩\in H^q:=\underset{i=0}{\overset{q}{⨁}}{H}^{(i)},$$

$${\tilde{\mathrm{T}}}_{\text{deg}}:=\bigcup _{j=0}^{\text{deg}}{\tilde{\mathrm{C}}}_j$$

$${\tilde{\mathrm{T}}}_4=\bigcup _{j=0}^4{\tilde{\mathrm{C}}}_j\equiv \left\{\text{mathbbm}1;{\eta }^{†}{\eta }^{†},{\eta }^{†}\eta ,\eta \eta ;{\eta }^{†}{\eta }^{†}{\eta }^{†}{\eta }^{†},{\eta }^{†}{\eta }^{†}{\eta }^{†}\eta ,{\eta }^{†}{\eta }^{†}\eta \eta ,{\eta }^{†}\eta \eta \eta ,\eta \eta \eta \eta \right\},$$

$$p=1+\max \left\{q\in \mathbb{N}|\forall |\xi ⟩,|{\xi }^{\prime }⟩\in H^q,⟨\xi |\mathcal{O}|{\xi }^{\prime }⟩=0\right\},$$

$${\tilde{\mathrm{T}}}^p:=\bigcup _{j=0}^p{\tilde{\mathrm{C}}}^j,$$

$$\begin{array}{rl}{\tilde{\mathrm{T}}}^4=\bigcup _{j=0}^4{\tilde{\mathrm{C}}}^j\equiv & \{\text{mathbbm}1;{\eta }^{†}\eta ;{\eta }^{†}{\eta }^{†},\eta \eta ,{\eta }^{†}{\eta }^{†}\eta \eta ;{\eta }^{†}{\eta }^{†}{\eta }^{†}\eta ,{\eta }^{†}\eta \eta \eta ,{\eta }^{†}{\eta }^{†}{\eta }^{†}\eta \eta \eta ;\\ & {\eta }^{†}{\eta }^{†}{\eta }^{†}{\eta }^{†},\eta \eta \eta \eta ,{\eta }^{†}{\eta }^{†}{\eta }^{†}{\eta }^{†}\eta \eta ,{\eta }^{†}{\eta }^{†}\eta \eta \eta \eta ,{\eta }^{†}{\eta }^{†}{\eta }^{†}{\eta }^{†}\eta \eta \eta \eta \}.\end{array}$$

$${\tilde{\mathrm{T}}}_{\text{deg}}^p={\tilde{\mathrm{T}}}_{\text{deg}-2}\cup {\tilde{\mathrm{C}}}_{\text{deg}}^p,$$

$${\tilde{\mathrm{T}}}_6^3\equiv \left\{\text{mathbbm}1;{\eta }^{†}{\eta }^{†},{\eta }^{†}\eta ,\eta \eta ;{\eta }^{†}{\eta }^{†}{\eta }^{†}{\eta }^{†},{\eta }^{†}{\eta }^{†}{\eta }^{†}\eta ,{\eta }^{†}{\eta }^{†}\eta \eta ,{\eta }^{†}\eta \eta \eta ,\eta \eta \eta \eta ;{\eta }^{†}{\eta }^{†}{\eta }^{†}\eta \eta \eta \right\}.$$

$$\min \left[\mathrm{Tr}(\rho \mathcal{D}),1-\mathrm{Tr}(\rho \mathcal{D})\right]\ll 1/2,$$

$$\mathcal{D}=\frac{1}{N}\sum _k{\eta }_k^{†}{\eta }_k^{}$$

$$\Delta ⟨{\mathcal{S}}_l^z{⟩}_{\tilde{\mathrm{T}}}(t)=\sqrt{\frac{{\int }_0^t{\left|⟨{\mathcal{S}}_l^z{⟩}_{\tilde{\mathrm{T}}}(u)-⟨{\mathcal{S}}_l^z{⟩}_{\text{ED}}(u)\right|}^{2}du}{1+{\int }_0^t{\left|⟨{\mathcal{S}}_l^z{⟩}_{\text{ED}}(u)\right|}^{2}du}}.$$

$${\tilde{\mathrm{T}}}_2⪯{\tilde{\mathrm{T}}}^2\prec {\tilde{\mathrm{T}}}_4⪯{\tilde{\mathrm{T}}}_6^3\prec {\tilde{\mathrm{T}}}_6^4⪯{\tilde{\mathrm{T}}}^4,$$

$${\mathcal{H}}_2={\mathcal{H}}_0+\sum _k\left(A_{\text{mbox}I}(k){\eta }_{-k}^{}{\eta }_k^{}+A_{\text{mbox}II}(k){\eta }_k^{†}{\eta }_k^{}+A_{\text{mbox}III}(k){\eta }_k^{†}{\eta }_{-k}^{†}\right)$$

Full text

1 Introduction

Equilibration and thermalisation are topics that link nonequilibrium physics to equilibrium physics, and they play a fundamental role for the validity and success of thermodynamics. These topics have a long history and have been studied in a variety of settings, including classical mechanics vs. quantum mechanics, closed systems vs. open systems, and others. Renewed interest in equilibration and thermalisation in isolated quantum systems was to a large extend triggered by experimental progress in preparing and manipulating assemblies of cold atoms that are extremely well isolated from their surroundings; see [1, 2] for reviews. Near-integrable systems, consisting of a dominant integrable part plus a small integrability-breaking perturbation, have been studied early on in some of these experiments, including the celebrated quantum Newton’s cradle by Kinoshita et al. [3]. The integrability-breaking perturbation ensures thermalisation to a microcanonical equilibrium, and the relaxation dynamics towards equilibrium in a near-integrable system takes place in two stages on widely separated timescales [4, 5, 6, 7]: A fast decay, termed prethermalisation, to a long-lasting nonequilibrium state that is characterised by a so-called generalised Gibbs ensemble (GGE) [8, 9]; and a second step, in which relaxation to thermal equilibrium, as described by the ordinary Gibbs ensemble, occurs on a much longer timescale, once the integrability-breaking perturbation becomes relevant.

Accurate and reliable calculations of these phenomena are challenging, at least when going beyond the small system sizes of $O(10)$ where exact diagonalisation is feasible. Perturbative techniques around the integrable limit suggest themselves for the problem at hand, and various types of such techniques have been employed in the context of prethermalisation, including a flow-equation methods [10], self-consistent mean-field techniques [11], and quantum kinetic theory [12, 13, 14, 15]. The notion of quantum kinetic theory subsumes a number of approximate methods based on identifying certain classes of operators (usually those of higher degree in the normal-ordered ladder operators; see Section 3 for more precise statements) as negligible, and deriving a reduced set of equations of motion for the remaining operators only [16]. In the abovementioned Refs. [12, 13, 14, 15] quantum kinetic theories are developed for studying bosons or fermions in one spatial dimension.

The quantum kinetic theory we develop in the present paper differs from those works in several important aspects. The integrable part of the Hamiltonian $\mathcal{H}_{\text{int}}$ we consider is an Ising spin chain with nearest-neighbour interaction and a transverse magnetic field, see Eq. (2). Our aim is to study the effect of weak long-range interactions, where we define long range as a power-law decay $|i-j|^{-\alpha}$ with the distance between lattice sites $i$ and $j$, where $\alpha$ is some nonnegative exponent.111The notion of long-range interactions is not unanimously defined. In some communities only exponents $\alpha$ smaller than the spatial dimension of the system are called long-range. Our terminology includes these cases, but is less restrictive. The specific long-range perturbation $\mathcal{H}_{\text{pert}}$ we consider is given in Eq. (3), but other types can be treated similarly. As is well known, the transverse-field Ising chain $\mathcal{H}_{\text{int}}$ can be mapped onto noninteracting fermions by a Jordan-Wigner transformation, followed by a Fourier and a Bogoliubov transformation [17, 18], which, one might think, should bring us back onto the familiar terrain of near-integrable fermionic models. However, applying the same sequence of transformations to $\mathcal{H}_{\text{pert}}$ generates complicated, non-number-conserving terms beyond those that are usually considered in fermionic models. As a consequence of these additional terms, quantum kinetic equations scale less favourably with the system size and the search for an optimised truncation scheme for those equations becomes a necessity. In Section 3 of this paper we introduce and discuss several such truncation schemes and benchmark them against exact results.

The long-range part $\mathcal{H}_{\text{pert}}$ of the Hamiltonian can be a small perturbation for one of two reasons: either because of a small prefactor $J_{z}$ in (3), or because of a large value of the long-range exponent $\alpha$. The quantum kinetic theory we develop in this paper applies to both cases, but the applications and results of Section 4 are for the latter case. To the best of our knowledge, this is the first example of a quantum kinetic perturbation theory that essentially uses $1/\alpha$ as a small parameter. The truncated set of quantum kinetic equations allows us to study the time evolution of spin expectation values and spin–spin correlation functions to high accuracy. Moreover, unlike some of the other kinetic equations techniques, our method does not require correlation functions to factorise as in the conditions of Wick’s first theorem. From our results we can distinguish different relaxation stages of the model, including prethermalisation due to the integrable part of the Hamiltonian, as well as the onset of thermalisation caused by the integrability-breaking terms.

2 Time evolution equations of a long-range spin chain

2.1 Near-integrable transverse-field Ising chain

We consider the Hamiltonian

[TABLE]

where

[TABLE]

describes an integrable transverse-field Ising chain, and

[TABLE]

is a long-range contribution. However, the methods developed in the following are expected to be applicable to a broader class of perturbations. Here, $l\in\llbracket 1,N\rrbracket$ labels the sites of a chain of length $N$. To each lattice site $l$ a spin-$1/2$ operator $\boldsymbol{\mathcal{S}}_{l}=(\mathcal{S}_{l}^{x},\mathcal{S}_{l}^{y},\mathcal{S}_{l}^{z})$ is associated, satisfying the commutation relations $\boldsymbol{\mathcal{S}}_{l}\times\boldsymbol{\mathcal{S}}_{q}=i\delta_{l,q}\boldsymbol{\mathcal{S}}_{l}$ (in units of $\hbar\equiv 1$). We assume periodic conditions $\boldsymbol{\mathcal{S}}_{N+1}\equiv\boldsymbol{\mathcal{S}}_{1}$, so that $\mathcal{H}$ is translationally invariant. Additionally, the Hamiltonian is invariant under the $\mathbb{Z}_{2}$ symmetry $x\to-x$. To account for the periodic boundary conditions, we define the distance between lattice sites $l$ and $l+m$ as the shortest connection around the circle, $d(m)=\min(m,N-m)$. The long-range interactions in (3) decay like a power law $d(m)^{-\alpha}$ with the distance, where $\alpha$ is some nonnegative exponent. To enforce that this term contains exclusively interactions beyond nearest neighbours, the sum over $m$ extends over $\llbracket 2,N-2\rrbracket$ only. The magnetic field strength is denoted by $h$, and $J_{x}$ and $J_{z}$ are pair coupling constants.

The integrable part $\mathcal{H}_{\text{int}}$ of the Hamiltonian is known to be exactly solvable by a Jordan-Wigner transformation, followed by a Fourier and a Bogoliubov transformation [17, 18] (see Section A.1). By means of this procedure, the integrable part (2) of the Hamiltonian can be brought into the quadratic form

[TABLE]

with dispersion relation222In case of a magnetic field reversal, such as the one we will use in Section C.2, this formula should be modified by replacing $\epsilon_{k}\to-\epsilon_{k}$ (c.f. Section A.1), which has an effect on the dynamics if $\mathcal{H}_{\mathrm{pert}}\neq 0$.

[TABLE]

where $k$ labels the momenta in the first Brillouin zone. $\eta_{k}^{\phantom{\dagger}}$ and $\eta_{k}^{\dagger}$ are fermionic operators satisfying the anticommutation relations \bigl{\{}\eta_{k}^{\phantom{\dagger}},\eta_{k^{\prime}}^{\dagger}\bigr{\}}=\delta_{k,k^{\prime}} and \bigl{\{}\eta_{k}^{\phantom{\dagger}},\eta_{k^{\prime}}^{\phantom{\dagger}}\bigr{\}}=0. It is crucial for what follows to also express the perturbation $\mathcal{H}_{\text{pert}}$ in this preferred fermionic basis in which $\mathcal{H}_{\text{int}}$ is quadratic and diagonal. This guarantees that, when making approximations by neglecting high-order terms, the error will be small. A proper definition of the notion of high-order operators is given in Section 2.2. The main steps of transforming $\mathcal{H}_{\text{pert}}$ into the fermionic quasi-particle basis $\eta^{{\dagger}}_{k},\eta_{k}^{\phantom{\dagger}}$ are reported in Section A.2, leading to the normal-ordered fermionic representation

[TABLE]

of the full Hamiltonian. Here we have employed the vector notation $\bm{k}=(k_{1},\dotsc,k_{4})$, and $\sum_{\bm{k}}$ indicates a summation over all momenta $k_{q}$ ($q=1,\dotsc,4$) in the Brillouin zone. The coefficients $A_{\mbox{\tiny I}},\dotsc,B_{\mbox{\tiny V}}$ are defined in Section A.2.2, and they contain contributions from the coupling constants in the original Hamiltonian, as well as combinatorial contributions that arise from normal-ordering. Normal-ordering leads to significantly more complicated expressions here, but it will be crucial for identifying negligible terms in the approximation scheme of Section 2.2. $\mathcal{H}_{0}$ in (2.1) denotes a term of degree zero in the fermionic operators, i.e., proportional to the identity. This term is irrelevant for the dynamics, but will be important for the definition of initial conditions. All terms in (2.1) are momentum conserving due to the translational invariance of the spin model, and of even degree in the fermionic operators because of the $\mathbb{Z}_{2}$ symmetry.

The Hamiltonian (2.1) is of quartic degree in the fermionic operators $\eta$, $\eta^{\dagger}$. This is different from the conventional long-range Ising model in a transverse field [19, 20, 21, 22] where, instead of $\mathcal{H}_{\text{pert}}$, a perturbation $J_{x}\sum_{l,m}\mathcal{S}_{l}^{x}\mathcal{S}_{l+m}^{x}d(m)^{-\alpha}/2$ with long-range couplings between the $x$-components of the spin operators is used, which leads to fermionic terms of arbitrarily high degree. The somewhat less conventional Hamiltonian (1)–(3) we chose is a convenient model for studying approximation methods for the dynamics of the spin chain: all deviations from the exact dynamics are expected to be genuine effects of the approximations made in the time-evolution equations, as no approximations have to be made on the level of the Hamiltonian.

2.2 Equations of motion

The time-evolution equation of an operator $\mathcal{O}$ in the Heisenberg picture is given by the von Neumann equation

[TABLE]

where $d_{t}\equiv\frac{d}{dt}$. To construct a quantum kinetic theory for the model (1)–(3), we require time-evolution equations of normal-ordered products of fermionic operators. For example, for $\mathcal{O}=\eta_{k}^{\dagger}\eta_{k}^{\phantom{\dagger}}$, a straightforward but tedious calculation yields

[TABLE]

where $\Delta_{\mbox{\tiny I}}(\bm{k})=\delta_{k}^{-k_{1}}+\delta_{k}^{k_{2}}+\delta_{k}^{-k_{3}}+\delta_{k}^{k_{4}}$, $\Delta_{\mbox{\tiny II}}(\bm{k})=\Delta_{\mbox{\tiny I}}(\bm{k})-\delta_{k}^{-k_{1}}-\delta_{k}^{k_{1}}$, $\Delta_{\mbox{\tiny III}}(\bm{k})=\Delta_{\mbox{\tiny II}}(\bm{k})-\delta_{k}^{k_{2}}-\delta_{k}^{-k_{2}}$, etc. The right-hand side of Eq. (7) generally involves time-evolved operators distinct from $\mathcal{O}$. Therefore, to solve the equation of motion (2.2), similar equations of motion have to be derived for the operators occurring on the right-hand side. In general, this will lead to a system of coupled differential equations whose number scales exponentially with the system size $N$. This is a problem of a complexity comparable to that of solving the von Neumann equation for the density operator in the Schrödinger picture, which is intractable already for moderate system sizes in most cases.

Our aim is to find a smaller differential system that is suitable for approximating the dynamics generated by the Hamiltonian (1), while being numerically tractable for larger system sizes. For this purpose, it will be convenient to classify operators according to their degree and their $p$-particle number.

Definition 1

Consider a product $\mathcal{O}=\mathcal{A}_{1}\cdots\mathcal{A}_{q}$ of fermionic operators $\mathcal{A}_{i}\in\{\eta_{k_{i}}^{{\dagger}},\eta_{k_{i}}\}$. Denote by $a\in\mathbb{N}$ the number of annihilation operators in $\mathcal{O}$, and by $c\in\mathbb{N}$ the number of creation operators. If all annihilation operators are to the right of all creation operators, the operator $\mathcal{O}$ is said to be normal-ordered. The degree of such a normal-ordered product is then defined as $\mathrm{deg}=a+c$, and we call the integer $p=\mathrm{max}(a,c)$ the $p$-particle number of $\mathcal{O}$.

We define the class $\mathrm{C}_{\text{deg}}^{p}$ as the unique set of normal-ordered products of fermionic operators with a given degree and a given $p$-particle number, e.g.

[TABLE]

where Br denotes the Brillouin zone, and the rightmost expression is a slightly abusive shorthand notation. We furthermore define superclasses of, respectively, fixed $p$-particle number and degree,

[TABLE]

The union $\mathrm{F}=\bigcup_{p=0}^{N}\mathrm{C}^{p}$, whose number of elements is exponentially large in the system size $N$, then spans the vector space of all fermionic operators acting on Fock space. To reduce the size of the system of coupled differential equations generated by (7), we introduce a truncation $\mathrm{T}\subsetneq\mathrm{F}$ as the union of, in general, several classes $\mathrm{C}_{\text{deg}}^{p}$. For example,

[TABLE]

corresponds to a truncation at the quadratic level, neglecting all terms of degree larger than two in the differential system. A first requirement on $\mathrm{T}$ to be a useful truncation is that it gives access to the observable(s) of interest. For instance, the spin component $\mathcal{S}_{l}^{z}$, expressed in terms of the fermionic operators $\eta_{k}$ and $\eta^{\dagger}_{k}$, is a linear combination of the identity and of some quadratic operators [see Eq. 55]. For simulating the dynamics of $\mathcal{S}_{l}^{z}$, the truncation must therefore contain at least the terms in (11), even though that selection may not be sufficient to obtain a good approximation. Other choices of observables may require larger truncations. Additionally, for the sake of numerical efficiency, we want $\mathrm{T}$ to contain only the most relevant terms to describe the dynamics, at least for the time-window and observable of interest, and for the level of precision required. In that sense, kinetic theory can be seen as a perturbation theory which aims at structuring the set of all operators into a hierarchy, and then truncates that hierarchy at a chosen level.

Such truncation schemes, as is evident from the definitions of the degree and the $p$-particle number, are based on the concept of normal ordering, and at least some kind of ordering is required for a consistent classification of operators and the establishment of a hierarchy. Even if both, the Hamiltonian and the observable are given in normal-ordered form, the commutator in (7) will usually create non-normally ordered terms in the system of coupled differential equations, which have to be normal-ordered before a truncation can be performed. For the applications considered in the present paper, the number of coupled differential equations typically scales like $N^{2}$ or $N^{3}$ with the system size $N$, and is therefore very large for system sizes of tens or even hundreds of spins. Hence, normal-ordering by hand is an arduous task. To avoid this, we have developed an algorithm, which we call the LKE (Linear Kinetic Equations) code, that takes care of the following tasks.

1. (i)

Symbolic calculation, for unspecified indices $k_{1},k_{2},\dotsc$, of the normal ordering of the commutators between all types of elements in $\mathrm{T}$. Technically this is equivalent to the derivation of Wick’s second theorem [23, 24]. 2. (ii)

Use the results of (i) to derive, from Eq. 7, the differential system $D$ for all operators $\mathcal{O}\in\mathrm{T}$. 3. (iii)

For a given initial density operator $\rho$, define $X_{0}$ as the vector composed of all $\langle\mathcal{O}\rangle=\mathrm{Tr}(\rho\mathcal{O})$. Numerically solve the coupled linear differential equations $\dot{X}=DX$ with initial condition $X(0)=X_{0}$. 4. (iv)

From $X(t)$, calculate the expectation value of the spin observable of interest, e.g. $\langle\mathcal{S}_{l}^{z}\rangle(t)$.

Our LKE code is different from other kinetic equations techniques, like the one developed in [14, 15] for Hubbard-type lattice models, in that it does not require the conditions of Wick’s first theorem to hold. Moreover, our approach gives direct access to correlation functions.

2.3 Initial states

In principle, the LKE code described above is not restricted to specific initial states, but specific choices may simplify the problem by reducing the size of the differential system $D$. In particular, spatially homogeneous initial states, which are invariant under discrete lattice translations, are a convenient choice, because they simplify the fermionic representation of observables like $\mathcal{S}_{l}^{z}$ (see Section A.3). This symmetry, as well as other ones, can be used to reduce the size of the differential system of kinetic equations, an issue that is discussed in detail in Appendix B. In principle, and for convenience, one could choose a homogeneous initial state that has a simple form in the fermionic basis. More relevant for physical applications, however, are initial states that have a simple form in the spin basis, as in this case it is more likely that such a state can be prepared experimentally.

A homogeneous initial state with a particularly simple form in the spin basis is a fully $z$-polarised state $\ket{\downarrow\cdots\downarrow}$, defined such that it satisfies $\mathcal{S}_{l}^{z}\ket{\downarrow\cdots\downarrow}=-1/2\ket{\downarrow\cdots\downarrow}$ for all $l$. Since the time evolution is calculated in the $\eta$-basis, the initial state needs to be transformed into that basis as well. Fortunately, fully polarised spin states have a convenient expression in the fermionic language. For instance, one can show that the fully down $z$-polarized spin state is transformed into the Bogoliubov basis according to

[TABLE]

where $\ket{0}$ denotes the Bogoliubov vacuum and

[TABLE]

for $n\in\llbracket 1,\lfloor N/2\rfloor\rrbracket$. $u_{k}$ and $v_{k}$, defined in Section A.1, are the coefficients of the Bogoliubov transformation that diagonalises the integrable part of the Hamiltonian, and $W_{n}$ is a normalization constant defined by

[TABLE]

We then define down truncated polarised states as

[TABLE]

with $n\in\llbracket 1,\lfloor N/2\rfloor-1\rrbracket$. For small sizes and/or large magnetic field amplitudes $|h|$, any of these states is a good approximation of the “proper” polarised state (12). However, for large systems or small magnetic fields, the LKE code is expected to perform well only for initial states (15) with small $n$, an effect that will become clearer in the context of the $p$-particle structure introduced in Section 3.1 and further discussed in Section D.1.

The symmetry properties of the truncated polarized states $\ket{\psi^{n}}$ can be used to further reduce the size of the differential system of kinetic equations. Firstly, these states belong to the even sector of the Fock space, and the time evolution under the Hamiltonian (2.1) preserves this evenness. Secondly, fermions created by the operator $\mathcal{G}_{n}$ in (13) always come in pairs with opposite momenta $k_{i}$ and $-k_{i}$, and one can show that a differential system restricted to products of operators that take into account this pair structure is sufficient to describe not only the initial state, but also the time evolution of a truncated polarized state. Similar to the classes of operators defined in Eqs. (9) and (10), we denote by $\tilde{\mathrm{C}}_{\text{deg}}^{p}$ the set of products of fermionic operators with a certain degree and $p$-particle number, with the additional constraints of satisfying momentum conservation, belonging to the even sector of the Fock space, and taking into account the pair structure of $\mathcal{G}_{n}$. A detailed account of these symmetries and a definition of the symmetry-reduced classes $\tilde{\mathrm{C}}_{\text{deg}}^{p}$ is given in Appendix B.

Lastly, it is worth noting that the truncated polarised states (15) do not satisfy the conditions of Wick’s first theorem. This is an interesting observation because of the fact that, different from other quantum kinetic equations that can be found in the literature, our LKE code does not rely on the validity of Wick’s theorem. For instance, for the state $\ket{\psi^{1}}$ one can show that

[TABLE]

and hence Wick’s first theorem does not apply. Similar conclusions can be drawn for all truncated polarized states.

3 Truncations and hierarchies

Our aim is to establish a hierarchy between operators according to their relevance for the time evolution, and then truncate that hierarchy at a certain level in order to reduce the size of the differential system of kinetic equations and render it numerically more manageable. For instance, in a weakly-interacting classical kinetic theory, one would first select the ballistic terms, for which the particles are non-interacting. If higher accuracy is needed one would include $2$-particle scattering terms, and so on. In this section we adapt this intuitive classical picture to the fermionic Hamiltonian (2.1) and comment on the role of initial conditions for selecting a suitable truncation scheme. In Section 3.3 we assess the quality of the approximations by benchmarking the results from different truncation schemes against exact results.

There is no rigorous theory that demands that hierarchies and truncation schemes be based on normal ordering, but on the more intuitive level one can reason as follows. Consider two fermionic states

[TABLE]

i.e. both states reside in the sector of Fock space that corresponds to at most $q$ fermions. Then, for any normal-ordered product of fermionic operators $\mathcal{O}\in\mathrm{C}_{\text{deg}}$, it follows that $\bra{\xi}\mathcal{O}\ket{\xi^{\prime}}=0$ if $\text{deg}>2q$, whereas such a matrix element can be nonzero if $\text{deg}\leq 2q$. The same is not true without normal-ordering, i.e. for a non-normal-ordered product $\mathcal{O}$ of fermionic operators, the matrix element $\bra{\xi}\mathcal{O}\ket{\xi^{\prime}}$ can be nonzero regardless of the degree of $\mathcal{O}$. This implies that, when disregarding non-normal-ordered operators of, say, $\text{deg}=4$, one is neglecting information not only about three and more fermions, but also about single fermions and pairs of fermions. This would contradict the intuitive, classical idea of a truncation scheme that we invoked at the beginning of this section. For a more detailed discussion of the reasoning behind normal ordering, see Chapter 4.1 of [24].

3.1 Definitions of truncations

Truncations based on the degree of operators

Based on the discussion in the preceding paragraph, it is natural to base a hierarchy of fermionic operators on their degree. Correspondingly, a truncation scheme

[TABLE]

is defined such that it contains only normal-ordered products of fermionic operators up to a certain degree (and which additionally meet the symmetry requirements discussed in Appendix B). The cardinality of $\tilde{\mathrm{T}}_{\text{deg}}$, and therefore the number of variables in $\dot{X}=DX$, scales like $N^{\text{deg}/2}$ with the number of sites $N$. For instance, at the quartic level, we have

[TABLE]

where the symbol ; is used to easily distinguish between the classes.

Truncations based on the $p$*-*particle number

Modifying the idea leading to the hierarchy (18), one can order the operators according to the integer

[TABLE]

which is precisely the $p$-particle number introduced in Section 2.2. The corresponding truncation is defined as

[TABLE]

which, for example, yields

[TABLE]

The two truncations $\tilde{\mathrm{T}}^{p}$ and $\tilde{\mathrm{T}}_{\text{deg}}$ are equivalent for Hamiltonians which obey fermion number conservation (assuming that $\text{deg}=2p$), but they differ in cases where, like in our fermionic Hamiltonian (2.1), terms like $\eta\eta$ or $\eta^{\dagger}\eta^{\dagger}$ create or destroy pairs of fermions. The cardinality of $\tilde{\mathrm{T}}^{p}$ scales like $N^{p}$ with the system size $N$.

Truncation based on degree and $p$*-*particle number

We can combine the ordering principles of the previous paragraphs in different ways. We introduce the truncation

[TABLE]

adding to the terms in $\tilde{\mathrm{T}}_{\text{deg}-2}$ only those of a specific degree and $p$-particle number. For example, for $\text{deg}=6$ and $p=3$ we have

[TABLE]

The number of normal-ordered products in $\tilde{\mathrm{T}}_{6}^{3}$ scales like $N^{3}$ with the number $N$ of spins.

3.2 Domain of validity of the truncations

The truncations introduced above are expected to yield good approximations of the dynamics for sufficiently short times. Longer times can be reached by tuning $\mathcal{H}$ and/or $\rho$ in a way such that expectation values of higher-degree fermionic operators are small. The main parameter for tuning the spin Hamiltonian (1)–(3) is the long-range variable $\alpha$. The larger $\alpha$, the closer the fermionic version (2.1) of the Hamiltonian is to the noninteracting integrable case. The smaller the interactions are, the longer it takes to build up correlations between fermions, and hence higher-degree fermionic operators remain close to their initial values for a longer time.

Another requirement for a truncation to yield a good approximation is that the initial state $\rho$ is uncorrelated or at most weakly correlated in the fermionic basis. Moreover, when using a truncation based on the $p$-particle hierarchy of Section 3.1, the approximation works particularly well for initial states with a small fermion density. By means of the particle–hole transformation of Section C.2 the validity can be extended to initial states having either a small fermion density or a small hole density,

[TABLE]

where

[TABLE]

is the (Bogoliubov) fermion density operator. Note that, as observed in Section D.1, for sufficiently small magnetic field amplitudes the fully-polarized state $\ket{\downarrow\cdots\downarrow}$ is expected to violate condition (25). However, for a given system size $N$, one can choose $n$ sufficiently small such that the truncated polarized state $\ket{\psi^{n}}$ falls into the range of validity of the truncation. Similarly, the validity of the condition (25) can be enforced by increasing, at fixed $n$, the system size $N$.

3.3 Benchmarking

In this section we assess the performance of the LKE code when using the truncations introduced in Section 3.1. We compare to exact diagonalization (ED) results [25] for system sizes up to $N=12$. As a measure for the accuracy, we use an indicator proportional to the time-integrated Euclidean distance between the LKE expectation value and the ED expectation value,

[TABLE]

Based on the results for various truncation schemes shown in Fig. 1, we make the following observations: For sufficiently short times, the accuracies of the truncations follow the hierarchy

[TABLE]

where the symbol $a\prec b$ means that the truncation $a$ is less accurate than $b$, and $\preceq$ means that the related truncations are equivalent in the large $\alpha$ limit. We note that the intuitive idea of a hierarchy based on the $p$-particle number and the degree is confirmed,333A more detailed benchmarking, which we do not show here, reveals that $\tilde{\mathrm{T}}^{p-1}\prec\tilde{\mathrm{T}}^{p}$ for $1\leq p\leq 4$ on the one hand, and $\tilde{\mathrm{T}}_{\text{deg}-2}\prec\tilde{\mathrm{T}}_{\text{deg}}$ for $2\leq\text{deg}\leq 6$ on the other hand, providing evidence of both, a degree hierarchy and a $p$-particle hierarchy. However, after the quartic level, such schemes are coarse, and it is the purpose of Fig. 1 to propose intermediate levels of approximation. but that “shortcuts” seem to exist, i.e. lower-order truncations that achieve more or less the same level of accuracy. For instance, $\tilde{\mathrm{T}}_{6}^{4}$ scales like $N^{3}$ with the system size, whereas $\tilde{\mathrm{T}}^{4}$ involves differential systems of size $O(N^{4})$, but for the model we study these two schemes become equivalent for large $\alpha$. Similarly, $\tilde{\mathrm{T}}_{4}$ and $\tilde{\mathrm{T}}_{6}^{3}$ give results of essentially the same accuracy, although the first truncation contains a significantly smaller number of operators, and is therefore numerically favourable. We found these observations to hold for all system sizes $N\leq 12$ for which we had ED results available for comparison, and we do not see any reason why the observed patterns should not remain valid for larger systems with otherwise similar parameter values.

As a rule of thumb, on a regular desktop computer we can deal with system sizes $\sim 10^{3}$ when using a truncation scheme for which the corresponding differential system in the LKE code scales linearly with $N$; system sizes of order $\sim 10^{2}$ when the scaling is quadratic in $N$; and sizes of order $\sim 40$ in the case of cubic scaling. Quadratic truncations, while scaling linearly with $N$, cannot capture effects beyond integrability, and hence are not suitable for our purposes. In the following we use the compromise $\tilde{\mathrm{T}}_{4}$, which scales quadratically in $N$,444Another promising quartic choice that scales quadratically with $N$ is $\tilde{\mathrm{T}}_{2}\cup\tilde{\mathrm{C}}_{4}^{2}\cup\tilde{\mathrm{C}}_{4}^{3}\preceq\tilde{\mathrm{T}}^{3}$. as our default truncation for the applications discussed in Section 4.

4 Prethermalisation and thermalisation in the long-range Ising chain

A nonintegrable isolated quantum system of large but finite size is expected to thermalise in a probabilistic sense, meaning that, at sufficiently late times, the expectation value $\langle\mathcal{O}\rangle(t)$ of a physically reasonable observable $\mathcal{O}$ is very close to its thermal equilibrium expectation value $\langle\mathcal{O}\rangle_{\text{th}}$ for most $t$ [26, 27]. Fluctuations around equilibrium are present, but their size is suppressed for large system sizes $N$; and while large deviations from equilibrium may occur, they are extremely rare.

In the transverse-field Ising chain with long-range interactions (1)–(3), the integrability of $\mathcal{H}_{\text{int}}$ is broken by the presence of $\mathcal{H}_{\text{pert}}$ which, for the large $\alpha$-values we are considering, is a weak perturbation. The relaxation to equilibrium of weakly nonintegrable systems, consisting of an integrable part plus a small nonintegrable perturbation, has been studied extensively in the literature (see [5] and references therein). For such systems, an out-of-equilibrium initial state typically approaches equilibrium in two stages [7]: On a rather short timescale, a long-lasting prethermalised nonequilibrium state is reached. This state is described by a so-called generalised Gibbs ensemble (GGE) [8, 9], which, in addition to conservation of energy, takes into account also all the other conserved local charges of the integrable part of the Hamiltonian. Proper thermal equilibrium, as described by the ordinary Gibbs ensemble, is expected to be approached only much later, once the integrability-breaking perturbation becomes relevant.

We expect a similar behaviour in the transverse-field Ising chain with long-range interactions, but with the difference that $\mathcal{H}_{\text{pert}}$ in (3) contains integrable as well as nonintegrable contributions, as is evident from the presence of quadratic as well as quartic terms in the fermionic Hamiltonian (2.1). Both types of contributions are of small magnitude, controlled by a combination of the parameters $J_{z}$ and $\alpha$. Notwithstanding the similar magnitudes of the integrable and nonintegrable contributions in $\mathcal{H}_{\text{pert}}$, the two terms will have different effects on the equilibration of the system. The integrable portion of $\mathcal{H}_{\text{pert}}$ will contribute a small shift to the GGE that is reached in the initial relaxation step due to $\mathcal{H}_{\text{int}}$. The nonintegrable portion of $\mathcal{H}_{\text{pert}}$ is generically expected to effect proper thermalisation to a Gibbs state on a timescale proportional to the squared inverse of the magnitude of the nonintegrable term [5].

In the following we make use of the LKE code with a suitable truncation scheme in order to probe the relaxation dynamics of the transverse-field Ising chain with long-range interactions. As our local observable of interest we choose $\mathcal{S}_{l}^{z}$, the $z$-component of the spin at site $l$. This observable has the advantage of being of a simple form not only in the spin framework, but also in the fermionic language, where it is a quadratic operator (55).

4.1 Quadratic fermionic Hamiltonian

To distinguish effects of nonintegrability from those of the integrable model, we define the Hamiltonian

[TABLE]

consisting of only the quadratic terms in the fermionic Hamiltonian (2.1). $\mathcal{H}_{2}$ differs from $\mathcal{H}_{\text{int}}$ for finite $\alpha$, owing to the fact that $\mathcal{H}_{\text{pert}}$ contains not only quartic, but also quadratic contributions. For any quadratic Hamiltonian, the truncation scheme $\tilde{\mathrm{T}}_{2}$ yields exact results,555This is a consequence of the fact that $id_{t}\mathcal{O}\in\operatorname{Span}\tilde{\mathrm{T}}_{2}$ for all $\mathcal{O}\in\tilde{\mathrm{T}}_{2}$; see Appendix B. and the use of such a low-order truncation scheme will allow us to deal with fairly large system sizes of $\mathcal{H}_{2}$ in the following.

Using that scheme in the LKE code, we show in Fig. 2 exact results for the time evolution of $\langle\mathcal{S}_{l}^{z}\rangle$ generated by the quadratic Hamiltonian $\mathcal{H}_{2}$ for various values of the long-range parameter $\alpha$. The most striking feature in this plot are the drastic changes, occurring periodically in the time evolution with a period of approximately 1200. These features have been termed traversals in [28], and they can be understood as a finite-size effect: The dynamics is controlled by pairs of quasiparticles travelling across the chain in opposite directions, and for $\mathcal{H}_{\text{int}}$ the maximum velocity of quasiparticle propagation is known [29]. Because of the periodic boundary conditions, the effect of returning quasiparticles that have travelled the full length of the circle will be felt after a time $\tau_{\text{trav}}\simeq N/|J_{x}|$ and multiples thereof. For $\alpha<\infty$ this timescale changes only slightly. As Fig. 2 illustrates, the traversals spoil the relaxation behaviour that is visible up to $t\simeq 1200$. In this way, finite system sizes limit the timescales that can be assessed. For the quest of observing equilibration, which occurs on the slowest relevant timescale of a system, this poses a challenge.

From now on we will focus exclusively on times $t$ up to $\tau_{\text{trav}}$. In that time window we observe in Fig. 2 a rapid rise of $\langle\mathcal{S}_{l}^{z}\rangle$ from approximately $-0.335$ [c.f. (101) and the discussion in Section D.2.2] to an $\alpha$-dependent value around $-0.3$. We estimate the corresponding timescale to be

[TABLE]

which yields a value $\tau\simeq 1$ for the parameters of Fig. 2, in agreement with the results shown in the plot. After that fast initial rise, a prethermalisation plateau is reached. We expect, but have not explicitly confirmed, that the attained long-time values agree with the GGE equilibrium values of $\mathcal{H}_{2}$ for the initial state used.

4.2 Beyond integrability

In this section we go beyond integrability by considering the full nonquadratic fermionic Hamiltonian (2.1), which is equivalent to the long-range spin model (1)–(3). In this case the quadratic truncation scheme $\tilde{\mathrm{T}}_{2}$ is not sufficient anymore, and we opt instead for using $\tilde{\mathrm{T}}_{4}$ as a compromise between accuracy and numerical efficiency (see Section 3.3). $\tilde{\mathrm{T}}_{4}$ scales quadratically with the system size, which restricts the system sizes we can deal with on a regular desktop computer to $N=120$. Because of the traversals discussed in Section 4.1, this will limit the timescales that can faithfully be observed to $\tau_{\text{trav}}\simeq 120$.

In Fig. 3 the time evolution of $\langle\mathcal{S}_{l}^{z}\rangle$ is shown for various values of the long-range parameter $\alpha$, for the full nonintegrable Hamiltonian $\mathcal{H}$ as well as for the quadratic Hamiltonian $\mathcal{H}_{2}$. The left panel of Fig. 3 shows that, for transverse magnetic field $h=-1$, the dynamics under $\mathcal{H}$ (solid lines) and $\mathcal{H}_{2}$ (dashed lines) are almost indistinguishable. Prethermal values are rapidly reached, and no subsequent drift towards the thermal values (indicated by straight solid lines) is evident on the accessible timescale. For $h=-0.51$ (right panel of Fig. 3), which is close to the quantum critical point of the model,666We have not studied the location of the critical point $|h_{\text{c}}|$ for finite $\alpha$. For a perturbation of the form $\sum_{l,m}\mathcal{S}_{l}^{x}\mathcal{S}_{l+m}^{x}d(m)^{-\alpha}$ this question has been addressed in [21]. Unlike in that case, our perturbation (3) couples spin components in the magnetic field direction, and for that reason we expect that the location of the critical point remains largely unaffected. For the parameter values we use, we hence expect that $h=-0.51$ is in the paramagnetic phase for all $\alpha$, which seems to be confirmed by numerical results. dashed and solid lines clearly differ, indicating that the quartic terms in the nonintegrable Hamiltonian $\mathcal{H}$ have a sizeable effect on the dynamics, at least for the smaller $\alpha$-values considered. Moreover, the presence of nonintegrable terms appears to promote thermalisation, shifting the time-evolving expectation values closer to their thermal equilibrium value.

The thermal equilibrium values shown in Fig. 3 are calculated according to

[TABLE]

where $Z(\beta)=\operatorname{Tr}(e^{-\beta\mathcal{H}})$ is the partition function. Since we are considering an isolated system where energy is conserved, the inverse temperature $\beta$ is fixed through the initial state $\ket{\psi^{1}}$ implicitly via

[TABLE]

at least under the idealisation of the thermodynamic limit. If $N$ is finite but sufficiently large, we expect (32) to still be valid. Based on this assumption we use exact diagonalisation (ED) for spin chains of up to 12 sites to determine $\beta$ and $\langle\mathcal{S}_{l}^{z}\rangle_{\text{th}}$. For the energy density $\bra{\psi^{1}}\mathcal{H}\ket{\psi^{1}}/N$, an exact expression is known, see (96). ED results for several small system sizes are then extrapolated to the system sizes of interest. For a magnetic field $h=-1$, we find an inverse temperature $\beta\simeq-5$, more or less independent of the value of $\alpha$. For $h=-0.51$, $\beta$ ranges from $\simeq-8.4$ to $\simeq-10.2$ for $\alpha$ between $4$ and $8$. Smaller $\alpha$ are not considered, as the validity of the approximations in our quantum kinetic theory become questionable in that case.

Negative inverse temperatures $\beta$ are known to occur in equilibrium systems with (upper and lower) bounded energy spectra if the entropy decreases as a function of energy in the high-energy region [30]. According to (65), in our model the transition from positive to negative $\beta$ takes place at the energy $\nu_{\text{th}}=0$. From Eq. (65) we furthermore find d_{\beta}\nu_{\text{th}}(\beta)\big{|}_{\begin{subarray}{c}\beta=0\end{subarray}}<0, and it is reasonable to assume that $\beta\mapsto\nu_{\text{th}}(\beta)$ is monotonous,777For instance, the condition $\mathrm{sp}(\mathcal{H})\subset\mathbb{R}_{-}^{*}$ is sufficient to obtain $d_{\beta}\nu_{\text{th}}(\beta)<0$. Moreover, from (45) we know that the eigenvalues of $\mathcal{H}_{\text{int}}$ are all strictly negative as long as $h<0$, which is always the case in this section. Perturbation theory therefore proves the strict positivity of the spectrum of the full Hamiltonian (2.1) for a sufficiently large (but finite) value of $\alpha$. which is consistent with the numerical results of Fig. 3. According to (98) the initial states we use correspond to positive energy densities, which, by virtue of the above monotonicity argument, imply negative temperatures. In Section C.1 we propose a method that allows us to tune the energy density of the initial state, and hence the effective temperature, while still using truncated polarised states and staying in the regime where our quantum kinetic theory remains valid.

4.3 Spreading of correlations

The linear quantum kinetic theory we developed in Section 2.2 does not make use of a Wick factorisation of correlations. This not only allows us to deal with correlated initial states, as discussed at the end of Section 2.3, but also provides access to a subset of fermionic correlation functions. At the level of the $\tilde{\mathrm{T}}_{4}$-truncation (19) that we use, we obtain all the fermionic correlations functions necessary for calculating the spin–spin correlations $\left\langle\mathcal{S}_{l}^{x}\mathcal{S}_{l+1}^{x}\right\rangle$, $\left\langle\mathcal{S}_{l}^{x}\mathcal{S}_{l+2}^{x}\right\rangle$, and $\left\langle\mathcal{S}_{l}^{z}\mathcal{S}_{l+m}^{z}\right\rangle$ for all $m\in\llbracket-\lfloor N/2\rfloor,\lfloor N/2\rfloor\rrbracket$ (see Section A.4 for explicit formulas). As an example, we show in Fig. 4 the time-evolution of the connected $zz$-correlation function

[TABLE]

for $h=-1$ and $h=-0.51$, starting from the initial state $\ket{\psi^{1}}$. This state has nonvanishing correlations in the spin basis, which, as is visible from the short-time behaviour in Fig. 4, are smaller for $h=-1$, and larger for $h=-0.51$. This $h$-dependence is a consequence of the fact that tuning $h$ not only modifies the Hamiltonian (2.1), but also changes the initial state (13), which affects the magnitude of the correlations (A.4). As discussed in Section D.1, $\ket{\psi^{1}}$ becomes a good approximation of $\ket{\downarrow\cdots\downarrow}$ in the limit of large negative magnetic fields, which is consistent with our observation of weaker initial correlations for a magnetic field of larger magnitude.

In Fig. 4 we show $\left\langle\mathcal{C}^{z}_{m}\right\rangle$ as a function of time $t$ and distance $m$ between lattice sites. We observe a rapid decay, on a timescale of the order one, of the nonlocal (i.e., $m$-independent) initial correlations, followed by a “lightcone”-like spreading of correlations in space and time [31, 32, 33, 34, 15]. In the presence of long-range interactions, a variety of analytical, numerical, as well as experimental results indicate that, at least for sufficiently small values of $\alpha$, the linear shape of the cone gets replaced by a curved shape [35, 36, 37, 38, 39, 40]. For $\alpha=4$ as used in Fig. 4 a curved shape is not visible, and it has in fact been conjectured that correlations spread strictly linearly for $\alpha$ larger than some critical value [41].

The spatial decay of the correlations, i.e. the $m$-dependence of $\left\langle\mathcal{C}_{m}^{z}\right\rangle$ at a given time $t$, appears significantly sharper in the left plot ($h=-1$) of Fig. 4 compared to the right plot ($h=-0.5$), a feature that is particularly striking at small $|m|$. This observation is consistent with the expectation that the correlation length diverges in the vicinity of the quantum critical point, which for $\alpha=4$ is expected to be close to the quantum critical point $h_{\text{c}}=-1/2$ of $\mathcal{H}_{\text{int}}$. Strictly speaking this argument is valid only for the groundstate and in equilibrium, but it is reasonable to expect signatures of quantum critically to persist at small Bogoliubov particle densities $\langle\mathcal{D}\rangle$. The initial state $\ket{\psi^{1}}$ has indeed a small $\langle\mathcal{D}\rangle$ (90) and, since the integrability breaking is weak, this remains true for fairly long times $t>0$. Furthermore, while global equilibrium has not yet been reached, regions that are some distance away from the edges of the lightcone seem to have equilibrated at least locally, with a correlation length that is presumably similar to that of the global equilibrium. Assuming all this heuristic reasoning to be valid, we interpret the qualitative differences between the two plots in Fig. 4 as consequences of the distance of the magnetic field values $h=-1$, respectively $h=-0.51$, from the quantum phase transition.

5 Conclusions

We have constructed quantum kinetic equations for describing the nonequilibrium dynamics of a transverse-field Ising chain with a weak integrability-breaking perturbation. The computational method we developed makes use of the Jordan-Wigner fermionic representation of the transverse-field Ising model and takes into account the integrability-breaking perturbation up to a certain degree in the time-evolution equations of operators. Which operators to include and which operators to neglect in the time-evolution equations is a crucial issue and strongly affects the accuracy of the approximation. In Section 3 we have introduced, discussed, and benchmarked several truncation schemes, all of which are based on the normal-ordering of products of fermionic operators. Based on the numerical benchmarking, we found the quartic truncation scheme $\tilde{\mathrm{T}}_{4}$ to be numerically efficient and at the same time adequate for studying effects beyond integrability. Truncation schemes involving sixth order terms can reduce errors in time-evolved expectation values by almost an order of magnitude, but become very costly in computation time. Using the truncation scheme $\tilde{\mathrm{T}}_{4}$, which scales quadratically in the system size $N$, we can reach sizes of up to $N=120$ on a desktop computer, but with more effort and/or high-performance computing facilities this value can certainly be pushed quite a bit further.

The model we have studied is the integrable transverse-field Ising chain with nearest-neighbour interactions (2), with an added integrability-breaking long-range perturbation (3). The perturbation can be made small by choosing either the coupling coefficient $J_{z}$ in (3) to be small, or the long-range exponent $\alpha$ to be large. The latter case, which we focus on in this paper, is, to the best of our knowledge, the first perturbative technique that uses $1/\alpha$ as a small parameter.

Research on systems with long-range interactions usually focusses on one of the following two cases: (i) Systems where the long-range exponent $\alpha$ is smaller than the spatial dimension of the system. In this case a number of unconventional thermodynamic and dynamic properties are known to occur, including thermal phase transitions in one-dimensional models [42], nonequivalence of statistical ensembles [43, 44, 45], and others. Our quantum kinetic theory applies to this regime when $J_{z}$ in (3) is sufficiently small, but we have not studied this case in detail in the present work. (ii) Values of the long-range exponent $\alpha$ that are relevant for recent experiments with ultracold atoms. For spin-$1/2$ models, these are in particular $\alpha=3$ (magnetic atoms, polar molecules, Rydberg atoms) and $\alpha=6$ (Rydberg atoms); see [46] for an overview. The latter case should certainly fall into the range of validity of our quantum kinetic perturbation theory when using $1/\alpha$ as a small parameter.

An important feature of our theory is that we do not assume the conditions of Wick’s first theorem to hold, i.e., unlike in some related work [14], we do not reduce expectation values of quartic fermionic terms into products of expectation values of quadratic terms. This comes with some advantages and some disadvantages. A disadvantage is that we need to solve a substantially larger set of coupled differential equations, which leads to restrictions on the accessible system sizes. This is attenuated to some extend by the fact that we deal with ordinary linear differential equations, whereas application of Wick’s theorem leads to nonlinearities. Moreover, since quartic terms are not broken up into quadratic ones, we have access to the corresponding quantum correlation functions. In the language of spin models, this gives us access not only to spin expectation values, but also to spin–spin correlation functions, as shown in Fig. 4.

As an application of our perturbative scheme we studied the influence of the small parameter $1/\alpha$ on prethermalisation and thermalisation in the weakly long-range transverse-field Ising chain. Finite-size effects restrict the accessible timescales to $t\simeq N$, which in turn implies a limitation on the relaxation phenomena one can observe. Relaxation due to the integrable part $\mathcal{H}_{2}$ occurs on a timescale of $O(1)$ and is easily observed. Thermalisation to a Gibbs state, induced by the nonintegrable part of $\mathcal{H}_{\text{pert}}$, takes place on a slower timescale. While we were not able to observe the full approach to thermal equilibrium in time, we do see that the presence of nonintegrable terms pushes the spin expectation values closer to their thermal values. This effect is more pronounced closer to the quantum critical point of the model, but we do not have a satisfactory explanation for this observation.

Acknowledgement

The authors benefited from helpful discussions with Fabian Essler, Johannes Kriel, and Stefan Kehrein. M. K. acknowledges financial support by the South African National Research Foundation through the Incentive Funding Programme and the Competitive Funding for Rated Researchers.

Appendix A Transformation into the diagonal basis of $\mathcal{H}_{\text{int}}$

As discussed in Section 2.1, we want to express the Hamiltonian as well as observables of interest as normal-ordered products of the fermionic operators that diagonalise the integrable part of the Hamiltonian. The reasoning behind this strategy is that high-order terms in those normal-ordered operator products are expected to be less relevant for the dynamics, and the kinetic equations derived in this paper are obtained by neglecting certain classes of normal-ordered fermionic operators. In Section A.1 we briefly recapitulate the standard result of diagonalising the transverse-field Ising chain with nearest-neighbour interactions by means of a Jordan-Wigner transformation, followed by a Fourier and a Bogoliubov transformation. In Section A.2 we express the long-range contribution $\mathcal{H}_{\text{pert}}$ in terms of the Bogoliubov fermions of Section A.1, and in the Appendices A.3 and A.4 we do the same for spin components $\mathcal{S}_{l}^{z}$ and spin–spin correlations, which will be our observables of interest.

A.1 Integrable part

The integrable part $\mathcal{H}_{\text{int}}$ of our Hamiltonian (1)–(3) describes a one-dimensional Ising chain in a transverse magnetic field with nearest-neighbour interactions. In this section we review the standard procedure of mapping this part of the Hamiltonian to noninteracting fermions by means of Jordan-Wigner, Fourier, and Bogoliubov transformations; see [18, 29, 47] for more detailed accounts.

A.1.1 Jordan-Wigner transformation

We consider the set of operators $\mu_{i}^{{\dagger}},\mu_{j}^{\phantom{{\dagger}}}$ satisfying the fermionic anticommutation relations

[TABLE]

Defining the Jordan-Wigner transformation

[TABLE]

with $S_{l}^{-}:=\left(S_{l}^{+}\right)^{\dagger}$ and $2S_{l}^{z}:=\left[S_{l}^{+},S_{l}^{-}\right]$, it is straightforward to verify that the fermionic anticommutation relations of $\mu_{i}^{{\dagger}},\mu_{j}^{\phantom{{\dagger}}}$ imply that $\mathcal{S}_{l}^{\pm},\mathcal{S}_{l}^{z}$ obey spin commutation relations, as required. Expressing the integrable part of the Hamiltonian (2) in terms of the fermionic operators $\mu_{i}^{{\dagger}},\mu_{j}^{\phantom{{\dagger}}}$ one obtains

[TABLE]

A.1.2 Fourier transformation

The spin Hamiltonian is invariant under discrete translations, which suggests to search for the eigenvectors of $\mathcal{H}_{\mathrm{int}}$ among the Fourier modes of a one-dimensional lattice with periodic boundary conditions. Our convention for the discrete Fourier transformation is

[TABLE]

The periodic boundaries impose conditions on the permissible momenta $k$ over which the sum in (37) extends, dependent on the total number of fermions on the chain. This number is obtained through the operator

[TABLE]

and we denote its eigenvalues by $M$. Then the permissible values of the momenta are given by $k=\frac{2\pi}{N}(q+1/2)$ with $q\in\llbracket 0,N-1\rrbracket$ if $M$ is even, and by $k=\frac{2\pi}{N}q$ if $M$ is odd.

The fermionic parity (i.e., the evenness or oddness of $M$) is conserved under the time evolution not only of the integrable part, but also of the full Hamiltonian, $\left[\mathcal{H},e^{i\pi\mathcal{M}}\right]=0$. We will restrict our attention to initial states from the even parity sector, and parity conservation will preserve that restriction for all later times. For convenience, within that sector we shift the Brillouin zone to be as symmetric as possible around zero by choosing the integers $q\in\llbracket-\lfloor N/2\rfloor,\lfloor N/2\rfloor-1\rrbracket$ in the above definition of the momenta.

To prepare for the Bogoliubov transformation to follow, we express $\mathcal{H}_{\text{int}}$ as a sum of matrix products,

[TABLE]

where

[TABLE]

with

[TABLE]

Since $R_{k}^{{\dagger}}=R_{k}^{\phantom{{\dagger}}}$, Eq. (39) can be diagonalized by means of a unitary transformation.

A.1.3 Bogoliubov transformation

We introduce the change of basis

[TABLE]

where $u_{k},v_{k}\in\mathbb{R}$. We denote by $\eta^{{\dagger}}_{k},\eta_{k}$ the image of the Fourier basis under this rotation, i.e.

[TABLE]

Since $\mathrm{det}(U_{k})=1$, there exists a real number $x_{k}$ such that $u_{k}=\cos(x_{k}),~{}v_{k}=\sin(x_{k})$. Requiring

[TABLE]

to be diagonal yields the Bogoliubov angle $x_{k}=\frac{1}{2}\tan^{-1}\left(b_{k}/a_{k}\right)$. Expressing $\mathcal{H}_{\mathrm{int}}$ in terms of the thus defined fermionic operators $\eta^{{\dagger}}_{k},\eta_{k}$ one obtains the diagonal Hamiltonian (4) with the dispersion relation

[TABLE]

This dispersion relation differs from the one in (5), and also from what is given in most papers and textbooks, by the factor of $\operatorname{sgn}(a_{k})$ [48]. This variant turns out to be useful in Section C.2, where we define a particle–hole mapping to reach the high-temperature regime, a transformation that, at least in the ferromagnetic phase, is equivalent (in the active viewpoint of symmetries) to a reversal of the magnetic field $h$.

A.2 Perturbation

In this section we express the long-range perturbation (3) of the Hamiltonian in terms of the Bogoliubov fermions $\eta^{{\dagger}}_{k},\eta_{k}$ defined in Section A.1.3.

A.2.1 Jordan-Wigner and Fourier transformations

Inserting the definitions of $\mu^{{\dagger}}_{l},\mu_{l}$ and $\tilde{\mu}^{{\dagger}}_{k},\tilde{\mu}_{k}$ into (3), we obtain

[TABLE]

where

[TABLE]

restricts the summation to momentum-conserving terms modulo $2\pi$. We have used the shorthand notation $\tilde{\mu}_{i}\coloneqq\tilde{\mu}_{k_{i}}$ and the truncated zeta-function

[TABLE]

A.2.2 Bogoliubov transformation

Performing the Bogoliubov transformation $\tilde{\mu}_{j}=u_{j}\eta_{j}-iv_{j}\eta^{{\dagger}}_{-j}$, the long-range part (46) of the Hamiltonian becomes

[TABLE]

where we have defined

[TABLE]

When normal-ordering the terms in the second sum of (A.2.2), further quadratic terms will emerge, which can be merged with the quadratic terms in $\mathcal{H}_{\mathrm{int}}$. The full Hamiltonian can finally be written as

[TABLE]

where $\mathcal{H}_{0}=H_{0}\,\mathbbm{1}$ with

[TABLE]

is proportional to the identity operator, and hence irrelevant for the dynamics. Here we have introduced the notations

[TABLE]

Note that $A_{\mbox{\tiny II}}(k)\in\mathbb{R}$ and $A_{\mbox{\tiny I}}(k),A_{\mbox{\tiny III}}(k)\in i\mathbb{R}^{2}$. Hermiticity imposes $\overline{A}_{\mbox{\tiny II}}=A_{\mbox{\tiny II}}$, which implies $A_{\mbox{\tiny II}}(k)\in\mathbb{R}$, but $A_{\mbox{\tiny I}}(k),A_{\mbox{\tiny III}}(k)\in i\mathbb{R}^{2}$ is not a necessary condition to fulfil $\overline{A}_{\mbox{\tiny I}}=A_{\mbox{\tiny III}}$. In addition, $A_{\mbox{\tiny I}}(-k)=-A_{\mbox{\tiny I}}(k)$, $A_{\mbox{\tiny III}}(-k)=-A_{\mbox{\tiny III}}(k)$ (consequence of the statistics), and $A_{\mbox{\tiny II}}(-k)=A_{\mbox{\tiny II}}(k)$. Similarly, $B_{\mbox{\tiny III}}(\bm{k})\in\mathbb{R}$, $B_{\mbox{\tiny II}}(\bm{k}),B_{\mbox{\tiny IV}}(\bm{k})\in i\mathbb{R}^{2}$, and $B_{\mbox{\tiny I}}(\bm{k}),B_{\mbox{\tiny V}}(\bm{k})\in~{}\mathbb{R}^{2}$. Hermiticity of the quartic part is guaranteed by the relations $B_{\mbox{\tiny I}}(1,2,3,4)=B_{\mbox{\tiny V}}(2,1,4,3)$, $B_{\mbox{\tiny II}}(1,2,3,4)=\overline{B}_{\mbox{\tiny IV}}(4,3,2,1)$, and $B_{\mbox{\tiny III}}(1,2,3,4)=B_{\mbox{\tiny III}}(4,3,2,1)$. Finally, $B_{\mbox{\tiny I}}(-\bm{k})=B_{\mbox{\tiny I}}(\bm{k})$, $B_{\mbox{\tiny II}}(-\bm{k})=-B_{\mbox{\tiny II}}(\bm{k})$, $B_{\mbox{\tiny III}}(-\bm{k})=B_{\mbox{\tiny III}}(\bm{k})$, with similar relations for $B_{\mbox{\tiny IV}},~{}B_{\mbox{\tiny V}}$.

A.3 Transformation of $\mathcal{S}_{l}^{z}$

Similarly, by performing Jordan-Wigner, Fourier, and Bogoliubov transformations, the $z$-component of the spin operator can be expressed in terms of the Bogoliubov fermions,

[TABLE]

In general, this expressions contains also terms of the form $\langle\eta_{1}^{{\dagger}}\eta_{2}\rangle$, $\langle\eta_{-1}\eta_{2}\rangle$, and $\langle\eta_{1}^{{\dagger}}\eta_{-2}^{{\dagger}}\rangle$, which, for $k_{1}\neq k_{2}$, are not momentum conserving. For (discrete) translationally invariant initial states, however, we show in Appendix B that such non-momentum-conserving terms have zero expectation values at all times. Therefore, the expectation value of (54) at time $t$ simplifies to

[TABLE]

This is an important simplification for the LKE code, as it allows us to restrict the set of operators considered in the code to momentum-conserving products of Bogoliubov fermions.

A.4 Transformation of correlation functions

In the spin picture we define the connected $xx$-correlation function as

[TABLE]

Because of the $\mathbb{Z}_{2}$ symmetry of Hamiltonian and initial state we are using, we have $\langle\mathcal{S}_{l}^{x}\rangle=0$ at all times and hence $\langle\mathcal{C}_{1}^{x}\rangle=\langle\mathcal{S}_{l}^{x}\mathcal{S}_{l+1}^{x}\rangle$. By performing Jordan-Wigner, Fourier, and Bogoliubov transformations and assuming a translationally invariant state, this correlation function turns out to be quadratic when expressed in terms of the Bogoliubov fermions,

[TABLE]

where

[TABLE]

Similarly, correlation functions $\langle\mathcal{C}^{x}_{m}\rangle=\langle\mathcal{S}_{l}^{x}\mathcal{S}_{l+m}^{x}\rangle$ can be expressed in terms of Bogoliubov fermions for any $m$, and they turn out to be of degree $2m$ in the fermionic basis. With the $\tilde{\mathrm{T}}_{4}$ truncation (19) we are using in Sections 4.2 and 4.3 we have access to all quartic fermionic operators, but not to all sixth-order terms, so that we can calculate $\langle\mathcal{C}^{x}_{2}\rangle$, but not the $xx$ correlations of spins that are further than two lattice sites apart.

For the $zz$-correlation function with respect to a translationally invariant state we find

[TABLE]

with

[TABLE]

Making use of $\langle\mathcal{S}_{l}^{z}\rangle$ (55) then allows us to express the connected $zz$-correlation function (33) in the fermionic basis. Since the resulting expression is quartic in the fermionic operators, it follows that the $\tilde{\mathrm{T}}_{4}$ truncation (19) is sufficient for calculating $\langle\mathcal{C}^{z}_{m}\rangle$ for arbitrary $m$.

Appendix B Symmetries

In Section 2.2 we have introduced truncations of the equations of motion based on the notions of the degree and the $p$-particle number of products of fermionic operators. The idea behind such a truncation is the knowledge (or belief, or hope) that the neglected classes of operators do not contribute significantly to the dynamics of the quantities of interest, at least on a certain time scale. Neglecting these operators drastically reduces the size of the differential system that needs to be considered. In addition to such approximately vanishing quantities, symmetries of the Hamiltonian and/or the initial state may lead to a strict decoupling of equations of motion, in the sense that the differential equations for a certain set of fermionic products is strictly independent of some other set of fermionic products (although the converse may in general not be true). Such a decoupling can, on top of the truncation that is essentially introduced by hand, further reduce the size of the differential system.

For instance, in the Schrödinger picture, if $\mathcal{O}$ is a normal-ordered product that does not commute with a given symmetry $\mathcal{U}$, i.e. $\mathcal{O}\in\mathrm{C}_{\cancel{\mathcal{U}}}$, and if the operators $\mathcal{U},\rho,\mathcal{H}$ verify the conditions

1. (a)

at time $t$, if $\left[\rho(t),\mathcal{U}\right]=0$, then $\operatorname{Tr}\left[\rho(t)\mathcal{O}\right]=0$ holds true, 2. (b)

$[\mathcal{H},\mathcal{U}]=0$, 3. (c)

$[\rho(0),\mathcal{U}]=0$,

then it follows888For infinitesimal $dt$, hypotheses (a) and (b) imply that $\rho(t+dt)=-i\left[\mathcal{H},\rho\right]dt+\rho(t)$ commutes with $\mathcal{U}$, so that $\mathrm{Tr}\left(\rho(t^{\prime})\mathcal{O}\right)=0$ for all $t^{\prime}>t$. that $\mathrm{Tr}\left(\rho(t)\mathcal{O}\right)=0$ for all $t$. As a result, operators in the complement $\mathrm{C}_{\cancel{\mathcal{U}}}$ can be safely ignored, without any approximation, when constructing the differential system of kinetic equations.

As pointed out in Section 2.3, the truncated polarised states (15) that we use as initial states have the following symmetries and resulting conservation laws:

1. (i)

Discrete translation invariance of the initial states $\ket{\psi^{n}}$ and the Hamiltonian $\mathcal{H}$, resulting in conservation of the total lattice momentum. For the example of the class of products of fermionic operators of degree 2 and $p$-particle number 1 in (9), the symmetry-reduced set is then given by \bigl{\{}\eta^{{\dagger}}_{k}\eta_{k}^{\phantom{\dagger}}~{}|~{}\forall k\in\text{Br}\bigr{\}}, which scales with system size as $N$ instead of $N^{2}$. Similarly, momentum conservation reduces the number of operators in any of the classes $\mathrm{C}_{4}^{p}$ to scale with system size like $N^{3}$ instead of $N^{4}$. 2. (ii)

Spin inversion symmetry in the $x$-direction of the spin Hamiltonian (1), leading to conservation of the fermionic parity \langle\exp\bigl{(}i\pi\sum_{k}\eta^{\dagger}_{k}\eta^{\phantom{{\dagger}}}_{k}\bigr{)}\rangle in the fermionic picture. Since the truncated polarized initial states $\ket{\psi^{n}}$ lie entirely in the even parity sector of the Fock space, all odd-degree normal-ordered products of fermionic operators strictly do not contribute and can be excluded from the differential system of kinetic equations.999 Many odd-degree normal-ordered products are also non-momentum-conserving, and have already been eliminated in (i), e.g. all $\eta_{k}$ for which $k\not\in\left\{0,\pi\right\}$. More generally, the set $S$ of operators of degree $2j+1\leq 2N$ that can be written as normal-ordered products of a momentum-conserving product of degree $2j$ times $\eta_{\pi}$ has also been eliminated before, and $S$ scales like $2j-1$. 3. (iii)

Fermions created by the operator $\mathcal{G}_{n}$ in (13) always come in pairs with opposite momenta $k_{i}$ and $-k_{i}$. Through the definition of the truncated polarized states $\ket{\psi^{n}}$ in Eq. (15), this implies that

[TABLE]

To exploit in the kinetic equations the symmetries (i)–(iii), we define, based on the classes $\mathrm{C}_{\text{deg}}^{p}$ of normal-ordered products of fermionic operators defined in Section 2.2, the symmetry-restricted classes

[TABLE]

This definition excludes, in agreement with (ii), all normal-ordered products of odd degree, and the construction via pair operators $\mathcal{P}_{i}$ in the second line of (62) guarantees momentum conservation (i) as well as the pair structure (iii). For example, the class $\tilde{\mathrm{C}}_{4}^{2}$ contains only elements of the form $\eta^{\dagger}_{-k}\eta^{\dagger}_{k}\eta^{\phantom{{\dagger}}}_{-k^{\prime}}\eta^{\phantom{{\dagger}}}_{k^{\prime}}$ and $\eta^{\dagger}_{k}\eta^{\dagger}_{k^{\prime}}\eta^{\phantom{{\dagger}}}_{k}\eta^{\phantom{{\dagger}}}_{k^{\prime}}$. We define $\tilde{\mathrm{C}}^{p}$ and $\tilde{\mathrm{C}}_{\text{deg}}$ analogous to (10), and the union

[TABLE]

contains all normal-ordered products of fermionic operators that satisfy the symmetry restrictions (i)–(iii). A proper subset $\tilde{\mathrm{T}}\subsetneq\tilde{\mathrm{F}}$ can then serve as a truncation of the differential system of kinetic equations.

For a truncated polarized initial state $\rho=\ket{\psi^{n}}\bra{\psi^{n}}$, we have $\mathrm{Tr}\left(\rho\mathcal{O}\right)=0$ for all $\mathcal{O}\in\mathrm{F}\setminus\tilde{\mathrm{F}}$. Moreover, for all $\tilde{\mathcal{O}}\in\tilde{\mathrm{F}}$, one finds $id_{t}\tilde{\mathcal{O}}\in\operatorname{Span}\tilde{\mathrm{F}}$.101010This follows from a similar property of the $(p\leq 2)$-subset $\tilde{\mathrm{T}}_{2}=\tilde{\mathrm{C}}_{0}\cup\tilde{\mathrm{C}}_{2}$, where it turns out that $\left[\mathcal{O},\mathcal{Q}\right]\in\operatorname{Span}\tilde{\mathrm{T}}_{2}$ for all $\mathcal{O},\mathcal{Q}\in\tilde{\mathrm{T}}_{2}$. Hence, a differential system that contains only products of operators from the set $\tilde{\mathrm{F}}$ is sufficient to describe not only the initial state, but also the time evolution of a truncated polarized state. The symmetry restrictions reduce the number of operators to be considered in the LKE code significantly, for example from $N^{4}$ to $N^{2}$ when going from $\mathrm{C}_{4}^{p}$ to $\tilde{\mathrm{C}}_{4}^{p}$.

Appendix C Dynamics and thermodynamics in the high-temperature limit

In Section 4.2 we employed exact diagonalisation of small, finite systems for calculating thermal expectation values, and extrapolated these values to large system sizes $N$. In this section we discuss how, without resorting to exact diagonalisation nor to conventional diagrammatic perturbation theory at equilibrium[15], thermal expectation values in the high-temperature limit can be obtained for large system sizes. We calculate, up to second order in the inverse temperature $\beta$, analytical expressions for the thermal values of both the energy density $\langle\mathcal{H}\rangle/N$ and the spin observable $\langle S_{l}^{z}\rangle$.

The truncated down-polarised states $\ket{\psi^{n}}$ that were introduced in (15), and for which our LKE code was tailored, do not usually fall into the regime where $\beta\ll 1$, and the same holds true for their up-polarized counterparts, which we denote by $\ket{\chi^{n}}$. However, by considering a suitable superposition of $\ket{\psi^{n}}$ and $\ket{\chi^{n}}$ the energy density can be tuned into the high-temperature regime. Under suitable conditions one can then argue that the dynamics of $\ket{\psi^{n}}$ and $\ket{\chi^{n}}$ decouples. Making use of a particle–hole transformation, the LKE code can be used for the calculation of the decoupled time evolution of $\ket{\chi^{n}}$, and the outcome can be compared to the thermal values obtained from a high-temperature expansion.

C.1 Thermal equilibrium results in the high-temperature limit

For sufficiently high temperatures, the Boltzmann factor can be expanded up to second order in $\beta$,

[TABLE]

Based on this expansion, and making use of $\operatorname{Tr}\left(\mathcal{H}\right)=0$, one can derive, up to the same order in $\beta$, the thermal expectation value of the energy density [49],

[TABLE]

where $K_{1}=2^{-N}\operatorname{Tr}\mathcal{H}^{2}/N>0$ and $K_{2}=-2^{-N-2}\operatorname{Tr}\mathcal{H}^{3}/N$. We find

[TABLE]

with the partial zeta-function $\zeta_{N}$ as defined in (48), which, for $\alpha>1/2$, converges in the large-$N$ limit. At next order, after a long calculation we find

[TABLE]

with $q\equiv m+n\bmod N$ and $L=\left\{N-1-m,N-m,N+1-m\right\}$. Equations (66) and (67) can be evaluated without too much effort for very large system sizes. For the purpose of benchmarking the LKE code, we are interested in $\alpha$ substantially larger than $1/2$. For such values we observe that $K_{1}$ and $K_{2}$ converge quickly with increasing $N$, being essentially indistinguishable from their infinite-system limit already for system sizes $\sim 10^{2}$. In that same regime of $\alpha$-values we also find that $K_{2}\ll K_{1}$. In fact, we expect more generally that odd orders in $\beta$ dominate the Taylor expansion (65), because $\nu_{\mathrm{th}}(\beta)$ is approximately an odd function. Along similar lines, we obtain

[TABLE]

for the observable of interest, and again the quadratic term in $\beta$ is negligible if $\alpha$ is appreciably larger than $1/2$. Based on these results, we can obtain accurate thermal expectation values $\upsilon_{\mathrm{th}}$ and $\langle S_{l}^{z}\rangle_{\text{th}}$ for the system sizes of order $\sim 10^{2}$ we want to benchmark against (or for much larger ones), in the regime of not-too-small $\alpha$ and small inverse temperatures $\beta$.

C.2 Particle–hole conjugation

The truncated polarized initial states $\ket{\psi^{n}}$, being highly ordered, are typically outside the high-temperature regime. Our strategy for obtaining a high-temperature initial state is to define, in the fermionic picture, a unitary particle–hole transformation, which maps fermions onto empty sites and vice versa. In the spin picture, this corresponds to transforming a fully down-polarized state $\ket{\downarrow\cdots\downarrow}$ into a fully up-polarized state $\ket{\uparrow\cdots\uparrow}$. Similarly, a truncated down-polarized state $\ket{\psi^{n}}$ is transformed into a truncated up-polarized state $\ket{\chi^{n}}$, and we find that a suitably chosen linear combination of the two states falls into the high-temperature regime, as demonstrated in Section D.2.1 for $n=1$ and $n=\lfloor N/2\rfloor$.

We define the particle–hole transformation as the unitary operator $\mathcal{U}$ that transforms Jordan–Wigner fermionic operators according to

[TABLE]

Under this transformation, all Jordan-Wigner fermionic states are mapped onto their particle–hole counterparts, e.g. $\mathcal{U}\ket{1011}=\ket{0100}$. From Eq. (69) it follows that $\mathcal{U}$ acts on Bogoliubov fermionic operators as

[TABLE]

By means of this transformation, we define truncated up-polarized states as

[TABLE]

To understand how the Hamiltonian is transformed under $\mathcal{U}$, we note that, in the spin framework, the spin operators are transformed like

[TABLE]

Transforming the Hamiltonian (1) under $\mathcal{U}$ corresponds to a sign reversal of constants,

[TABLE]

In the case of $2|h|>|J_{x}|$, to which we apply the LKE code in this paper, we have $\operatorname{sgn}(2h+J_{x}\cos k)=\operatorname{sgn}(h)$ for all momenta $k$, which implies that, according to (45), the magnetic field reversal effected by $\mathcal{U}$ amounts to replacing $\epsilon_{k}\to-\epsilon_{k}$ in $\mathcal{H}_{\mathrm{int}}$. For the perturbation $\mathcal{H}_{\mathrm{pert}}$, the dictionary (53c)–(53j) remains unchanged.

C.3 Decoupled dynamics

For truncated down-polarized states $\ket{\psi^{n}}$ with small $n$, we have by construction that the number of fermions is small, $\bra{\psi^{n}}\sum_{k}\eta^{{\dagger}}_{k}\eta^{\phantom{{\dagger}}}_{k}\ket{\psi^{n}}\ll N/2$. As discussed in Section 3.2, this is a requirement for our kinetic theory to provide a good approximation. For the truncated up-polarized states $\ket{\chi^{n}}$, in contrast, we have a large number of fermions, $\bra{\chi^{n}}\sum_{k}\eta^{{\dagger}}_{k}\eta^{\phantom{{\dagger}}}_{k}\ket{\chi^{n}}\gg N/2$, and the kinetic theory is expected to fail. However, such a state has only a small number of holes, and by applying the particle–hole transformation $\mathcal{U}$ introduced in Section C.2, it can be mapped to a state that satisfies the requirement of a small fermion number. Likewise by means of $\mathcal{U}$, the corresponding kinetic equations are obtained,

[TABLE]

where $\tilde{\mathcal{H}}\coloneqq\mathcal{H}(-h,-J_{x})$ is the transformed Hamiltonian and $\tilde{d}_{t}$ denotes the time-differential operator under $\tilde{\mathcal{H}}$. The time-evolution of a superposition

[TABLE]

of truncated up- and down-polarized states is then given by

[TABLE]

where we have defined $D_{t}\coloneqq|y_{1}|^{2}d_{t}-|y_{2}|^{2}\tilde{d}_{t}$. If $\ket{\psi^{n}}$ and $\ket{\chi^{n}}$ are few-fermion vectors, the second term on the right hand side of Eq. (76) is small. This is a consequence of the fact that the $p$-particle number of the terms in $\left[\mathcal{H},\mathcal{S}_{l}^{z}\right]$ is at most $p=4$, and hence does not couple few-fermion states to few-hole states in Fock space. Therefore, for $\alpha$ and $t$ sufficiently small, the dynamics of few-fermion and few-hole states approximately decouples,

[TABLE]

On the practical side, this implies that the time-evolution of a superposition $\ket{\phi^{n}}$ of truncated up- and down-polarized initial state can be computed with the LKE code by making use of the original Hamiltonian $\mathcal{H}$ as well as its field-inverted counterpart $\tilde{\mathcal{H}}$. We expect that the decoupling approximation (77) works well when $n$ is small, and, for fixed $n$, becomes better with increasing $N$. The excellent performance of the LKE code and the decoupling approximation is illustrated and benchmarked in Fig. 5.

Appendix D Initial conditions

D.1 Truncated polarized states and $p$-particle structure

In this subsection we sketch the main steps required in rewriting $\ket{\downarrow\cdots\downarrow}$ in the $\eta$-basis. We also comment on the link between this state and the $p$-particle structure, elaborating in particular on the role of the system size $N$ and the parameter $h/J_{x}$. These arguments will provide the main motivation for using small-$n$ truncated polarized states in Section 4.

Fully polarized state

Applying the Jordan-Wigner and Fourier transformations, the fully down-polarized spin state $\ket{\downarrow\cdots\downarrow}$ is mapped onto the Fourier vacuum $\ket{0^{\mathrm{F}}}$, defined as the only state satisfying

[TABLE]

The Bogoliubov transformation $\tilde{\mu}_{k}\mapsto u_{k}\eta_{k}-iv_{k}\eta_{-k}^{{\dagger}}$, however, mixes creation and annihilation operators in such a way that the Bogoliubov vacuum $\ket{0}$ for which

[TABLE]

is different from $\ket{0^{\mathrm{F}}}$. In the Fock space of Bogoliubov fermions, the Fourier vacuum can be expanded as

[TABLE]

where the $\lambda_{k_{1},\dotsc,k_{r}}$ are the coefficients we need to determine. Since the operator $u_{k}\eta_{k}-iv_{k}\eta^{{\dagger}}_{-k}$ is block anti-diagonal in the decomposition $H=H_{\mathrm{even}}\oplus H_{\mathrm{odd}}$ of the Bogoliubov-Fock space, condition (78) remains true for the restrictions $\tilde{\mu}_{k}|_{\mathrm{even}}$ and $\tilde{\mu}_{k}|_{\mathrm{odd}}$ of $\tilde{\mu}_{k}$ to the even and odd sectors, respectively. One can then prove that

[TABLE]

i.e. vectors in the odd sector do not contribute to $\ket{\downarrow\cdots\downarrow}$. To prove this result, we apply (78) to $\tilde{\mu}_{k}|_{\mathrm{odd}}$, which yields

[TABLE]

where the notation $\sum_{\bm{k}(q)}$ with $\bm{k}(q)=(k_{1},\dotsc,k_{2q+1})$ denotes a summation over all $k_{1}<\dotsb<k_{2q+1}$. Here we have defined $\mathcal{O}_{\bm{k}(q)}:=\eta^{{\dagger}}_{k_{1}}\cdots\eta^{{\dagger}}_{k_{2q+1}}$ and have labelled by $\underline{k},~{}\cancel{k}$ a coefficient that contains, respectively excludes, the momentum $k$ in its definition, e.g.

[TABLE]

with an analogous definition for $\Lambda_{\bm{k}(q)}^{\cancel{-k}}$. Projecting Eq. 82 onto $\ket{0}$ one finds that $\lambda_{k}=0$ for all $k$. For $q\geq 1$, we note that $\eta_{k}^{\phantom{{\dagger}}}\mathcal{O}_{\bm{k}(q)}$ and $\eta_{-k}^{\dagger}\mathcal{O}_{\bm{k}(q-1)}$ are of degree $2q$, so that the projection on any $2q$-excited state gives a recurrence relation between some of the $\Lambda_{\bm{k}(q)}^{\underline{k}}$ and some of the $\Lambda_{\bm{k}(q-1)}^{\cancel{k}}$. For instance, for $q=1$ one obtains

[TABLE]

It then follows that, for any $k_{2},k_{3}$ such that $-k\not\in\left\{k_{2},k_{3}\right\}$, we have $\Lambda^{\underline{k}}_{k,k_{2},k_{3}}=0$. Otherwise,

[TABLE]

and reasoning by induction leads to Eq. (81). Applying (78) to $\mu_{k}|_{\mathrm{even}}$ gives similar results for all the $\lambda_{k_{1},\dotsc,k_{2q}}$ coefficients, except that $\lambda_{0}\neq 0$. Finally one obtains $\ket{\downarrow\cdots\downarrow}=\mathcal{G}_{\lfloor N/2\rfloor}\ket{0}$, with

[TABLE]

where $W^{2}=1+\sum_{s=1}^{\lfloor N/2\rfloor}L_{s}$ with

[TABLE]

We prefer Eq. (86) over the exponential formulation $\mathcal{G}_{\lfloor N/2\rfloor}=W^{-1}\exp\left[i\sum_{0<k<\pi}\left(v_{k}/u_{k}\right)\eta_{-k}^{\dagger}\eta_{k}^{\dagger}\right]$ that has been used in the literature [29], as (86) naturally lends itself to the definition of truncated polarized states, as discussed in the next paragraph.

Truncated polarized states: properties and limitations

The $n$th truncated polarized state $\ket{\psi^{n}}$, as defined in Eq. (15), is obtained by truncating the first sum in (86) at the index $s=n$. Such a truncation preserves the reflection symmetry in momentum space, and therefore allows us to make use of the symmetry restrictions discussed in Appendix B. Also, the numerical computation of the vector $\ket{\psi^{n}}$ scales like $N^{n}$, whereas that of the fully polarised state scales like $\lfloor N/2\rfloor!$. Finally, truncated polarized states come with the additional advantage of allowing for a systematic tuning of the particle density such that, by making $n$ sufficiently small, the validity of the LKE code can be ensured. This can be useful for small magnetic fields $h$, where \bigl{\langle}\sum_{k}\eta_{k}^{\dagger}\eta_{k}^{\phantom{{\dagger}}}\bigr{\rangle}\ll N/2 in general does not hold for a (non-truncated) polarized state, as illustrated in Fig. 6.

A straightforward calculation shows that the expectation value of the fermionic particle density (26) with respect to any truncated polarized state is given by

[TABLE]

For the moment let us assume that, for sufficiently large $N$, $n$ being fixed,

[TABLE]

which will be justified towards the end of the section. Then it follows from (88) and (89) that

[TABLE]

for sufficiently large $N$. From the upper bound (90), the criterion (25) for the validity of the $p$-particle truncation for a truncated polarized initial state $\ket{\psi^{n}}$ immediately follows for $p\gtrsim n$. Moreover, since $\ket{\psi^{n}}=\ket{\downarrow\cdots\downarrow}$ for $n=\lfloor N/2\rfloor$, there must exist a smallest-possible $\nu$ for which, for a given system size $N$, $\langle\mathcal{D}\rangle_{\psi^{\nu}}\simeq\langle\mathcal{D}\rangle_{\downarrow}$ to a desired level of accuracy. It seems reasonable to assume that $\ket{\psi^{\nu}}$ is then a good approximation of the corresponding fully-polarized vector.

We chose to present in Section 4 LKE results only for $\ket{\psi^{1}}$, where (25) holds by construction, independently of the choice of $h/J_{x}$. $\ket{\psi^{1}}$ is not necessarily a good approximation of the fully polarized state, as for instance in Fig. 3 for parameter values $N=120$ and $h=-0.51$. While also in this case $\ket{\psi^{1}}$ is a perfectly legitimate choice as an initial state in the LKE code, it does not have a simple (approximate) representation in the spin picture, and the physical relevance of such an initial state is unclear.

In the remainder of this section we complete the above reasoning by providing a justification of the asymptotic property (89). We start by showing that $L_{1}$ grows linearly with $N$ asymptotically in the large-$N$ limit. From the definition (87) of $L_{s}$, together with the expressions of the Bogoliubov coefficients $u_{k}$ and $v_{k}$ derived in Section A.1.3, it follows that

[TABLE]

where

[TABLE]

and $\kappa=2h/J_{x}$ is the order parameter. For $\kappa>1$, each term in the sum of (91) is positive and smaller than one, which implies $L_{1}(N)\leq N$. To also prove a lower bound on $L_{1}$, we define, for a fixed $0<\epsilon\ll\pi/4$, the interval

[TABLE]

chosen such that, for a fixed $M>1$, for all $\kappa\in\left]1,M\right]$ and $\theta\in I$ we have

[TABLE]

Then it follows that

[TABLE]

where the first inequality follows from positivity of $g$, and the second from the fact that $g$ is an increasing function on its domain. Note that $g(A)>0$ because, by construction, $A>0$. The third inequality in (95) is then valid except for very small $N$, where the discreteness of the Brillouin zone may spoil it. Taking the upper and the lower bound together, it follows that $L_{1}(N)=\Theta(N)$. Along similar lines one can show that $L_{s}(N)=\Theta(N^{s})$, which implies (89).

D.2 Initial expectation values

In this section we gather expressions of the energy density $\nu=\langle\mathcal{H}\rangle/N$ and the spin observable $\langle\mathcal{S}_{l}^{z}\rangle$ for several specific states: $\ket{0}$, $\ket{\psi^{1}}$ and $\ket{\downarrow\cdots\downarrow}$, which are referred to in the main text, as well as their particle–hole counterparts $\ket{N}$, $\ket{\chi^{1}}$ and $\ket{\uparrow\cdots\uparrow}$, which are required in Appendix C.

D.2.1 Energy density

We obtain $\nu_{0}=\langle\mathcal{H}_{0}\rangle/N$ for the Bogoliubov vacuum $\ket{\psi^{0}}=\ket{0}$, and $\nu_{\lfloor N/2\rfloor}=-h/2+J_{z}\zeta/8$ for the fully $z$-polarised state $\ket{\psi^{\lfloor N/2\rfloor}}=\ket{\downarrow...\downarrow}$. For the truncated polarised state $\ket{\psi^{1}}$ one finds

[TABLE]

The particle–hole counterpart $\overline{\nu}_{1}:=\bra{\chi^{1}}\mathcal{H}\ket{\chi^{1}}/N$ is identical to (96) under magnetic field reversal $h\to-h$. For $N>8$ one finds $\bra{\chi^{1}}\mathcal{H}\ket{\psi^{1}}=0$, and hence the energy density of the superposition $\ket{\phi^{1}}=y_{1}\ket{\psi^{1}}+y_{2}\ket{\chi^{1}}$ simplifies to

[TABLE]

where $y_{1}$ and $y_{2}$ are complex coefficients normalized such that $|y_{1}|^{2}+|y_{2}|^{2}=1$. As a rule of thumb,

[TABLE]

so that for $h<0$ it follows with (45) that $\nu_{1}>0$ and $\overline{\nu}_{1}<0$. It is therefore possible to choose $y_{1},y_{2}$ such that $\nu_{\phi^{1}}\simeq 0$, which turns out to be useful in Section C.2 as a way of constructing initial states in the regime of small energy densities, which, according to Eq. 65, correspond to small inverse temperatures. To calculate the energy densities related to other truncated polarized states, we rely on the LKE code, which gives numerically exact initial energy densities for all truncations that contain $\tilde{\mathrm{T}}_{4}$ as a subset. Finally, for the state $\ket{\phi^{\lfloor N/2\rfloor}}$, which is used in Fig. 5, one can check that

[TABLE]

and thus it is once more possible to adjust $y_{1}$ and $y_{2}$ such that $\upsilon_{\phi^{\lfloor N/2\rfloor}}=0$.

D.2.2 Spin observable

We obtain $\langle\mathcal{S}_{l}^{z}\rangle=\pm 1/2$ for a fully $z$-polarized initial state, and

[TABLE]

for the Bogoliubov vacuum $\ket{0}$ and the anti-vacuum $\ket{N}$, with $\Gamma_{N}$ as defined in (53b). For $\ket{\psi^{1}}$ we find

[TABLE]

where the second term on the right-hand side is a correction to $\langle\mathcal{S}_{l}^{z}\rangle_{0}$ of order $O(1/N)$. Indeed,

[TABLE]

Moreover, (53b) [and hence (100)] can be regarded as a Riemann sum, and one can take the continuum limit $N\to\infty$. At finite $N$, the error made in the substitution by an integral is bounded by a term of order $O(1/N)$, and using (102) we obtain

[TABLE]

For $J_{x}=-1$ this implies that

[TABLE]

Since $h\mapsto\lim_{N\to+\infty}\langle\mathcal{S}_{l}^{z}\rangle_{1}(h)$ is strictly increasing on $]{-\infty},-1/2]$, Eq. (104) provides an analytical justification of the numerically observed initial values of $\langle\mathcal{S}_{l}^{z}\rangle$ in Figs. 2 and 3. In Fig. 3 for instance where $N=120$, the initial values predicted by (104) are $\langle\mathcal{S}_{l}^{z}\rangle_{1}\simeq-0.33$ for $h=-0.51$ and $\langle\mathcal{S}_{l}^{z}\rangle_{1}\simeq-0.46$ for $h=-1$, at order zero in $1/N$. One can then check that, with the exact expression (101), the error is indeed smaller than $10^{-2}$.