Confinement in the three-state Potts quantum spin chain in extreme ferromagnetic limit

TL;DR

This paper studies the three-state Potts quantum spin chain's dynamics in a ferromagnetic limit, revealing resonant excitations and post-quench evolution through a perturbative approach.

Contribution

It introduces a perturbative method to analyze the Potts model's dynamics, capturing phenomena beyond previous methods, especially in the oblique quench regime.

Findings

Analytic prediction of post-quench time evolution.

Identification of resonant excitations near stability thresholds.

Excellent agreement between analytical results and numerical simulations.

Abstract

We investigate the dynamics of the three-state Potts quantum spin chain in the extreme ferromagnetic limit using a perturbative expansion in the transverse magnetic field. We demonstrate that a perturbative approach provides access to important features beyond the reach of previous studies, most notably the description of resonant excitations and the analytic prediction of post-quench time evolution. A central focus is the oblique quench regime - a feature unique to the Potts model with no Ising counterpart - in which unconfined kink excitations hybridise with the two-kink bound states. We provide a detailed examination of the analytic structure of the two-kink scattering amplitude, tracing the transformation of stable excitations into resonances near stability thresholds. Our analytical results for the excitation spectra and magnetisation dynamics show excellent agreement with numerical simulations of non-equilibrium dynamics.

Tables

Table 1. Table 3.1 : Energy Spectra of Bound States and Resonances ( g = 0.2 , P = 0 g=0.2,P=0 ). This table lists the discrete energies ( E n E_{n} ) of stable collisionless bubbles ( h 2 > 0 h_{2}>0 ) and mesons ( h 2 < 0 h_{2}<0 ), alongside the complex energies of resonance states ( E res , l E_{\text{res},l} ). The real part of the resonance energy indicates the peak position in the Fourier spectrum, while the imaginary part represents the decay width resulting from hybridisation with the unconfined two-kink continuum.

$h_2$	$E_1$	$E_2$	$E_3$	$E_4$	$E_{res,1}$	$E_{res,2}$
0.1	1.70975	1.60194	1.50022	1.40001	1.8382 - 0.00898 i	2.0169 - 0.01813 i
0.2	1.60932	1.40079	1.20003	1.00000	1.8608 - 0.02246 i
-0.1	2.29025	2.39806	2.49978	2.59999	1.9831 - 0.01813 i	2.1618 - 0.00898 i
-0.2	2.39071	2.5992	2.79997	3.00000	2.1392 - 0.02246 i

Equations

$$H_{IC}(h_x,h_z)=-\frac{1}{2}\sum _{j=-\infty }^{\infty }({\sigma }_j^z{\sigma }_{j+1}^z+h_x{\sigma }_j^x+h_z{\sigma }_j^z).$$

$$[H_{IC}(h_x,0),U]=0,$$

$$U=\underset{j=-\infty }{\overset{j=+\infty }{⨂}}{\sigma }_j^x.$$

$${}^{(1)}⟨\text{vac}|{\sigma }_j^z|\text{vac}{⟩}^{(1)}=\bar{\sigma },{}^{(2)}⟨\text{vac}|{\sigma }_j^z|\text{vac}{⟩}^{(2)}=-\bar{\sigma },$$

$$H_{\text{IFT}}(\tau ,h)=H_{\ast }+\int dx\left(\tau ϵ(x)+h\sigma (x)\right)$$

$$\eta =2\pi \frac{\tau }{|h{|}^{8/15}}.$$

$$2m<M_1(\eta )<M_2(\eta )<\dots <M_{N(\eta )}(\eta )<2M_1(\eta ).$$

$$\underset{\eta \to {\eta }_N+0}{\lim }{M}_N(\eta )=2M_1({\eta }_N).$$

$$\underset{\eta \to -\infty }{\lim }{M}_1(\eta )=m.$$

$$H=-\sum _{j=-\infty }^{\infty }\sum _{\mu =1}^3(P_j^{\mu }{P}_{j+1}^{\mu }+h_{\mu }{P}_j^{\mu })-g\sum _{j=-\infty }^{\infty }{\tilde{P}}_j,$$

$$H(\mathbf{v},g)$$

$$H_0=-\sum _{j=-\infty }^{\infty }\sum _{\mu =1}^3\left(-\frac{2}{3}+P_j^{\mu }{P}_{j+1}^{\mu }\right),$$

$$V(\mathbf{v})=-\sum _{j=-\infty }^{\infty }\left({\tilde{P}}_j+\sum _{\mu =1}^3{v}_{\mu }\left(P_j^{\mu }-c^{\mu }\right)\right),$$

$$W_{\infty }=\underset{j=-\infty }{\overset{\infty }{⨂}}[{\mathbb{C}}^3{]}_j.$$

$$P_j^1=\left(\begin{array}{ccc}2/3& 0& 0\\ 0& -1/3& 0\\ 0& 0& -1/3\end{array}\right),P_j^2=\left(\begin{array}{ccc}-1/3& 0& 0\\ 0& 2/3& 0\\ 0& 0& -1/3\end{array}\right),$$

$$P_j^3=\left(\begin{array}{ccc}-1/3& 0& 0\\ 0& -1/3& 0\\ 0& 0& 2/3\end{array}\right),{\tilde{P}}_j=\left(\begin{array}{ccc}0& 1/3& 1/3\\ 1/3& 0& 1/3\\ 1/3& 1/3& 0\end{array}\right).$$

$$[H,{\hat{T}}_1]=0.$$

$${\hat{T}}_1|\mu {⟩}_j=|\mu {⟩}_{j-1},$$

$$W_N=\underset{j=-N/2+1}{\overset{N/2}{⨂}}[{\mathbb{C}}^3{]}_j.$$

$$H_0|0{⟩}^{(\mu )}=0,$$

$$|0{⟩}^{(\mu )}=\underset{j\in \mathbb{Z}}{⨂}|\mu {⟩}_j.$$

$$\mathcal{L}=\underset{\mu ,{\mu }^{\prime }=1}{\overset{3}{⨁}}{\mathcal{L}}_{\mu {\mu }^{\prime }}.$$

$${\mathcal{L}}_{11}={\mathcal{L}}_{11}^{(0)}\oplus {\mathcal{L}}_{11}^{(2)}\oplus {\mathcal{L}}_{11}^{(3)}\oplus \dots .$$

$$H_0|\Phi ⟩=n|\Phi ⟩,$$

$$|{\mathbf{K}}_{1,\nu }(j_1){\mathbf{K}}_{\nu ,1}(j_2)⟩=\left(\underset{j=-\infty }{\overset{j_1}{⨂}}|1{⟩}_j\right)\left(\underset{j=j_1+1}{\overset{j_2}{⨂}}|\nu {⟩}_j\right)\left(\underset{j=j_2+1}{\overset{\infty }{⨂}}|1{⟩}_j\right),$$

$$⟨{\mathbf{K}}_{1,\nu }(j_2){\mathbf{K}}_{\nu ,1}(j_1)|{\mathbf{K}}_{1,{\nu }^{\prime }}(j_1^{\prime }){\mathbf{K}}_{{\nu }^{\prime },1}(j_2^{\prime })⟩={\delta }_{\nu {\nu }^{\prime }}{\delta }_{j_1,j_1^{\prime }}{\delta }_{j_2,j_2^{\prime }},$$

$${\hat{T}}_1|{\mathbf{K}}_{1,\nu }(j_1){\mathbf{K}}_{\nu ,1}(j_2)⟩=|{\mathbf{K}}_{1,\nu }(j_1-1){\mathbf{K}}_{\nu ,1}(j_2-1)⟩.$$

$$|\mathbf{K}(j_1)\mathbf{K}(j_2){⟩}_{\iota }=\frac{1}{\sqrt{2}}(|{\mathbf{K}}_{1,2}(j_1){\mathbf{K}}_{2,1}(j_2)⟩+\iota |{\mathbf{K}}_{1,3}(j_1){\mathbf{K}}_{3,1}(j_2)⟩),$$

$$|K_{1,\nu }(p_1)K_{\nu ,1}(p_2)⟩=\sum _{j_1=-\infty }^{\infty }\sum _{j_2=j_1+1}^{\infty }\sum _{{\nu }^{\prime }=2}^3[e^{i(p_1{j}_1+p_2{j}_2)}{\delta }_{\nu ,{\nu }^{\prime }}+S_{\nu ,{\nu }^{\prime }}(p_1,p_2)e^{i(p_1{j}_2+p_2{j}_1)}]$$

$$\times |{\mathbf{K}}_{1,{\nu }^{\prime }}(j_1){\mathbf{K}}_{{\nu }^{\prime },1}(j_2)⟩,$$

$$|K(p_1)K(p_2){⟩}_{\iota }=\frac{1}{\sqrt{2}}(|K_{1,2}(p_1)K_{2,1}(p_2)⟩+\iota |K_{1,3}(p_1)K_{3,1}(p_2)⟩),\text{with }\iota =\pm$$

$$|j,P{⟩}_{\nu }=\sum _{j_1=-\infty }^{\infty }\exp \left[iP\left(j_1+\frac{j}{2}\right)\right]|{\mathbf{K}}_{1,\nu }(j_1){\mathbf{K}}_{\nu ,1}(j_1+j)⟩,$$

$${\hat{T}}_1|j,P{⟩}_{\nu }=e^{iP}|j,P{⟩}_{\nu },$$

Full text

1. 1 Introduction

2. 2 Spectrum of the three-state Potts quantum spin chain

3. 2.1 A brief review of the Ising case

4. 2.2 The three-state Potts model

5. 2.2.1 The Potts spin chain Hamiltonian

6. 2.2.2 The purely transverse case

7. 2.2.3 Confinement in the Potts model

8. 2.3 Structure of the Hilbert space

9. 2.4 The spectrum of two-kink states

10. 2.4.1 Pure transverse chain: kink excitations

11. 2.4.2 Two-kink states obtained from $\mathcal{L}_{11}^{(2)}$: the aligned case

12. 2.4.3 Two-kink states obtained from $\mathcal{L}_{11}^{(2)}$: the oblique case

13. 3 Resonance excitations in the oblique regime

14. 3.1 The bound state and resonance spectrum

15. 3.2 Trajectory of poles and zeroes of the scattering amplitude

16. 4 Evolution of magnetisation after a quantum quench

17. 4.1 Perturbation theory for Potts quench dynamics

18. 4.2 Comparison to numerical simulations

19. 4.2.1 Time evolution of the magnetisation

20. 4.2.2 Quench spectroscopy in aligned quenches

21. 4.2.3 Quench spectroscopy in oblique quenches

22. 5 Conclusions

23. A Details of the perturbative description of the excitation spectrum

24. A.1 The pure transverse chain: kink excitations

25. A.2 Two-kink states in the aligned case

26. A.3 Two-kink states in the oblique case

27. B True/false vacuum states for $g\ll 1$

28. C Perturbative expansion for Potts quench dynamics

29. D Quantum quench in the Ising spin chain at a weak transverse magnetic field

30. E Details of the numerical simulations

31. References

1 Introduction

Confinement occurs when particles cannot exist in isolation, but instead are observed only within bound states, a phenomenon most famously exemplified by quarks in hadrons [1, 2]. In systems with one spatial dimension, such as two-dimensional Quantum Field Theories (QFT) and quantum spin chains, the confinement of kink topological excitations is a common feature of models with discrete symmetry and order-disorder phase transitions [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14].

The most thoroughly investigated case is the Ising Field Theory (IFT) and its lattice counterpart, the mixed-field Ising spin chain. In the ferromagnetic phase, applying a longitudinal magnetic field induces a long-range attractive interaction between kinks, leading to their confinement into ”meson” (two-kink) bound states.

Beyond the spectrum, it is also interesting to consider non-equilibrium dynamics initialised from a fully polarised state. In the Ising model in the ferromagnetic phase, it was shown to lead to dynamical confinement when the longitudinal magnetic field is parallel to the initial magnetisation [15], and mesonic excitations dominate the dynamics. When the longitudinal magnetic field is anti-parallel to the initial magnetisation, it leads to Wannier-Stark localisation [16, 17], with the dynamics dominated by bubbles nucleating over a false vacuum. These bubbles undergo Bloch oscillations, leading to the suppression of false vacuum decay [18].

However, the $\mathbb{Z}_{2}$ symmetry of the Ising model significantly restricts its dynamics, limiting its ”hadron” spectrum exclusively to mesons. In contrast, the three-state Potts model possesses an $S_{3}$ permutation symmetry that allows for both meson and ”baryon” (three-kink) excitations in the ferromagnetic phase [19, 9, 20, 21, 22]. While recent studies of the Potts quantum spin chain have successfully used semiclassical quantisation to explain two-kink spectra [23], these methods face inherent limitations. Specifically, semiclassical approaches cannot account for the time evolution when the system is driven out of equilibrium. A particularly challenging non-equilibrium regime is the oblique quench, which has no Ising counterpart. In these cases, localisation is only partial, leading to unconfined kink excitations [23]. The hybridisation between the unconfined two-particle continuum and the bound states leads to resonances in the spectrum that also require an analytical treatment beyond semiclassical approximations.

In this work, we address these challenges by examining the Potts model in the extreme ferromagnetic regime using a perturbative expansion in the transverse magnetic field. This approach provides several key advantages. It allows the calculation of the excitation spectrum and the post-quench time evolution of magnetisation up to second order, respectively. Specifically, the calculation successfully determines the spectrum of the resonance excitations via the analytic structure of scattering amplitudes, thus overcoming the limitations of the semiclassical and exact diagonalisation methods. Furthermore, it provides a predictive framework for post-quench dynamics that captures features such as high-frequency oscillations, well beyond the reach of previous methods.

The paper is organised as follows. Section 2 introduces the Potts Hamiltonian and discusses its spectrum. Section 3 explores the analytic structure of scattering amplitudes and the transformation of stable bound states into resonances. Section 4 presents our results on quench dynamics, comparing perturbative predictions with iTEBD numerical simulations. We conclude in Section 5. Appendices A-C contain the details regarding the perturbative calculations of the excitation spectrum, vacuum structure and quench time evolution for the Potts spin chain, Appendix D demonstrates the perturbative approach to quench time evolution using the Ising case, while Appendix E gives the technical details of the iTEBD simulations.

2 Spectrum of the three-state Potts quantum spin chain

2.1 A brief review of the Ising case

As a warm-up, we recall the essential features of confinement in the Ising spin chain. The mixed-field Ising quantum spin chain is defined by the Hamiltonian:

[TABLE]

Here $\sigma^{x,z}$ are Pauli matrices, $j$ enumerates the chain sites, $h_{x}$ is the transverse and $h_{z}$ is the longitudinal magnetic field.

At $h_{z}=0$, the Hamiltonian (2.1) has the $\mathbb{Z}_{2}$-symmetry

[TABLE]

induced by the spin-flip operator

[TABLE]

The model (2.1) is integrable at $h_{z}=0$, and reduces to a free-fermion model by means of the Jordan-Wigner transformation. The quantum phase transition takes place on the axis $h_{z}=0$ at $h_{x}=1$, which separates the paramagnetic ($h_{x}>1$) and ferromagnetic ($0<h_{x}<1$) phases.

In the ferromagnetic phase $0<h_{x}<1$ at $h_{z}=0$, the Hamiltonian $\mathbb{Z}_{2}$-symmetry is spontaneously broken, and model (2.1) has two degenerate ferromagnetic vacua $|\text{vac}\rangle^{(\mu)}$, $\mu=1,2$, which are distinguished by the sign of the spontaneous magnetisation:

[TABLE]

with $\bar{\sigma}=(1-h_{x}^{2})^{1/8}>0$. Elementary excitations in this case are the kinks $K_{12}$ and $K_{21}$, which interpolate between the two vacua.

Applying a longitudinal magnetic field $h_{z}>0$ breaks the model’s integrability and explicitly breaks the $\mathbb{Z}_{2}$ symmetry of its Hamiltonian. As a result, the vacuum $|\text{vac}\rangle^{(1)}$ decreases in energy and becomes the true ground state, while the vacuum $|\text{vac}\rangle^{(2)}$ increases in energy and transforms into a metastable false vacuum. The energy difference between the true and false vacua induces a long-range attractive interaction between kinks, which, in turn, leads to their confinement: isolated kinks do not exist anymore in the system, and the kinks are bound into compound particles called ‘mesons’, the spectrum of which has been studied extensively [8, 15, 18, 24].

The scaling limit of the spin chain (2.1) to the Ising field theory (IFT) is defined in the vicinity of the second-order phase transition point $h_{x}=1$, i.e. in the limit $|h_{x}-1|\to 0$ with scaling $h_{z}\to 0$ appropriately [25]:

[TABLE]

where ${H}_{*}$ is the fixed-point Hamiltonian (corresponding to a massless Majorana fermions), $\epsilon$ is the energy operator (fermion mass term), and $\sigma$ is the continuum field corresponding to the order parameter $\sigma^{z}_{i}$. The physics of the IFT is determined by a single dimensionless scaling parameter $\eta$, defined as [7]:

[TABLE]

Confinement in the IFT was first examined by McCoy and Wu [3], and the mass spectrum of mesons was determined by combining analytic and numerical methods [3, 6, 7, 26, 8, 27, 21]. An important feature of the spectrum is the existence of resonant excitations, first noted in [3].

In the limits $\eta\to\pm\infty$, corresponding to zero magnetic field $h=0$, the IFT reduces to a relativistic QFT describing free spinless Majorana fermions with mass $m=2\pi|\tau|$. These fermions represent ordinary particles in the paramagnetic phase ($\eta\to-\infty$), whereas in the ferromagnetic phase ($\eta\to+\infty$) they become topological excitations - kinks interpolating between two degenerate ferromagnetic vacua. To describe the evolution of the particle content and physical properties of the IFT as the scaling parameter $\eta$ varies between the limiting values $\pm\infty$, McCoy and Wu proposed the following scenario [3].

The case $0<\eta<\infty$ corresponds to the ferromagnetic IFT in the presence of a magnetic field $h$, which lifts the degeneracy between the two ferromagnetic ground states and leads to confinement of kinks into the ‘meson’ bound states. The mass $M_{1}(\eta)$ of the lightest meson approaches from above the two-kink mass value $2m$ as $\eta\to+\infty$, and increases as $\eta$ decreases. For a fixed positive value of $\eta$, there exists a finite number $N(\eta)$ of stable meson modes with masses $M_{n}(\eta)$ satisfying $M_{n}(\eta)<M_{n+1}(\eta)$ for $n=1,\ldots,N(\eta)-1$. The number $N(\eta)$ diverges as $\eta\to+\infty$ and decreases as $\eta$ decreases.

The mass $M_{N(\eta)}(\eta)$ of the heaviest stable meson does not exceed the threshold value $2M_{1}(\eta)$. Thus, for generic $\eta>0$, the masses of the $N(\eta)$ stable mesons satisfy

[TABLE]

When the parameter $\eta$ approaches a certain value $\eta_{N}$ from above, the mass of the $N$-th meson reaches the stability threshold $2M_{1}(\eta_{N})$:

[TABLE]

Upon decreasing $\eta$ below $\eta_{N}$, the $N$-th meson becomes unstable with respect to decay into two lightest ($n=1$) mesons and transforms into a resonance state.

Further decreasing $\eta$ to negative values leads to a continued reduction in the number of stable particles $N(\eta)$. Eventually, for $\eta<\eta_{2}$ (with some negative $\eta_{2}$), only a single stable particle remains in the theory. In the limit $\eta\to-\infty$, this particle becomes a free Majorana fermion, and

[TABLE]

The McCoy-Wu scenario outlined above was subsequently confirmed and elaborated in numerous analytical and numerical investigations; see [7, 27] and references therein. Two particularly striking aspects of this scenario are the intrinsic presence of resonances in the IFT and the transformation of heavy stable mesons into resonances as the scaling parameter $\eta$ is tuned and the meson mass crosses the stability threshold $2M_{1}(\eta)$. The detailed mechanism of this transformation, which was not specified in [3], was later proposed by Fonseca and Zamolodchikov [7]. Unfortunately, a consistent analytical treatment remains elusive as it requires knowledge of the elastic two-particle scattering amplitude of the lightest ($n=1$) meson. Since these mesons are themselves composite particles, an analytic calculation of their mutual scattering amplitude remains a highly non-trivial and unresolved problem.

There are strong reasons to expect that the meson-resonance transformation scenario suggested in [7] is rather general and also occurs in other QFTs and spin-chain models exhibiting confinement. It turns out that resonance excitations also play an important role in the non-equilibrium dynamics of the Potts spin chain, and we present an analysis of their spectrum in Section 3. In the present section, we focus on the meson spectrum.

2.2 The three-state Potts model

2.2.1 The Potts spin chain Hamiltonian

The three-state Potts quantum spin chain is defined by the Hamiltonian:

[TABLE]

where the index $j\in\mathbb{Z}$ enumerates the sites of the infinite chain. It is convenient to modify the Hamiltonian (2.10) to the following form by adding a suitable constant:

[TABLE]

where $\mathbf{v}=(v_{1},v_{2},v_{3})$ is composed of the rescaled longitudinal magnetic fields $v_{\mu}=h_{\mu}/g$ and the constant shift is specified by $c^{1}=2/3$, and $c^{2}=c^{3}=-1/3$. The Hamiltonian acts on the infinite tensor product Hilbert space

[TABLE]

The three-dimensional local vector space $[\mathbb{C}^{3}]_{j}$, associated with the site $j$, has the basis $|\mu\rangle_{j}$, with $\mu=1,2,3$. In this basis, the action of operators $P^{\mu}_{j}$, $\mu=1,2,3$, and $\tilde{P}_{j}$ in the space $[\mathbb{C}^{3}]_{j}$ is determined by the $3\times 3$ - matrices:

[TABLE]

The parameter $g>0$ is the transverse magnetic field, and parameters $h_{\mu}$, with $\mu=1,2,3$, constitute the three components of the longitudinal magnetic field.

The Hamiltonian (2.11) commutes with the unit-step translation operator $\hat{T}_{1}$:

[TABLE]

The latter maps the three-dimension space $[\mathbb{C}^{3}]_{j}$ onto $[\mathbb{C}^{3}]_{j-1}$, acting on the basis vectors of the former as follows:

[TABLE]

with $j\in\mathbb{Z}$, and $\mu=1,2,3$.

The Hamiltonian (2.11) corresponds to the thermodynamic limit of the Potts spin chain, which is the main subject of our interest. However, in certain intermediate calculations, we also address the finite Potts spin chain with an even number of sites $N$. The Hamiltonian $H_{N}$ of the latter is obtained from the right-hand side of (2.11), in which the index $j$ runs from $-N/2+1$ to $N/2$. The Hamiltonian $H_{N}$ acts on the $3^{N}$-dimensional vector space

[TABLE]

The spectrum of the model (2.11) at $h_{\mu}=0$ can be constructed perturbatively in the extreme ferromagnetic limit $g\ll 1$ [28]. It turns out that the spectrum of this model can also be computed in the extreme ferromagnetic limit $g\ll 1$ at any fixed values of the components $v_{\mu}={h_{\mu}}/{g}\neq 0$ of the rescaled longitudinal magnetic field. A similar technique has been used previously for studying the kink confinement in the extreme anisotropic limit of the mixed-field Ising spin chain [8, 29], and of the antiferromagnetic XXZ spin chain [30, 31, 10, 32, 33]. Perturbative expansion in the transverse field was also found to be useful for the study of transport in the confining Ising spin chain [16].

2.2.2 The purely transverse case

The Hamiltonian $H(\mathbf{v},g)$ at $\mathbf{v}=0$ is known as the transverse field Potts quantum spin chain. It is invariant under the group $S_{3}$ of permutations of three ‘colours’, and displays a continuous order-disorder phase transition. In the ferromagnetic phase $0<g<1$, it has three energetically degenerate polarised vacua $|\text{vac}\rangle^{(\mu)}$, which for $g=0$ take the form:

[TABLE]

enumerated by the index $\mu=1,2,3$. The quasi-particle spectrum is made up from six types of kinks $K_{\mu\nu}$ (’quarks’) interpolating between $|\text{vac}\rangle^{(\mu)}$ and $|\text{vac}\rangle^{(\nu)}$, with $\mu,\nu=1,2,3$.

The scaling limit of the purely transverse spin chain, called the Potts field theory (PFT), is integrable, and the kink scattering matrix [34] and form factors of physically relevant operators [35] in this model are exactly known. Note, however, that the Potts quantum spin chain with a finite $g$ is not integrable, even in the purely transverse case.

2.2.3 Confinement in the Potts model

Similarly to the IFT, application of a longitudinal magnetic field breaks the integrability of the three-state Potts field theory (PFT) and leads to confinement of kinks, as shown by Delfino and Grinza [19]. These authors performed the symmetry analysis of the kink bound states in the $q$-state PFT and showed that in contrast to the $q=2$ Ising model, both mesons and baryons are allowed at $q=3$. The meson masses in the $q$-state PFT in the weak confinement regime were analytically calculated to leading order in $h$ in [9], while the masses of several lightest baryons in the three-state PFT in the leading order in $h$ have been calculated in [20]. These predictions were later confirmed by direct numerical calculations [21].

Recently, confinement of kinks was studied away from the scaling limit, i.e. in the three-state Potts spin chain defined by the Hamiltonian (2.11) [23]. Similarly to the field theory version, switching on a (weak) longitudinal field in the ferromagnetic phase leads to confinement, with the spectrum containing meson and baryon excitations.***We note that the realisation of baryon excitations in condensed matter systems has attracted recent interest [36, 37, 38]. It was found that the spectra of two-kink excitations can be explained well using semiclassical quantisation [23].

In the Potts model, the non-equilibrium dynamics following quantum quenches is largely analogous to the Ising case when the longitudinal field is either parallel or antiparallel to the initial magnetisation. The main difference from the Ising case is the presence of baryonic oscillation signatures. However, for the Potts model, it is also possible to realise oblique quenches in which the longitudinal field is in a direction different from the initial magnetisation, and such quenches have no Ising counterparts [23]. Unlike the Ising model, where the longitudinal field confines all excitations, the $S_{3}$ symmetry of the Potts model allows specific field orientations to retain a residual degeneracy. In these cases, confinement is only partially effective, which leads to the presence of unconfined kink excitations. Additionally, the hybridisation of the corresponding unconfined two-kink continuum with the bound states leads to the appearance of resonances in the spectrum, which cannot be treated using semiclassical methods.

Due to the appearance of false vacuum states, the excitation spectrum also contains bubbles which correspond to the nucleation of the true vacuum from the false vacuum according to the Coleman scenario [39]. Contrary to field theory, these excitations can be localised in a lattice setting (Wannier-Stark localisation) [16, 17, 40, 24]. Similarly to mesons, the bubble excitations can also hybridise with an unconfined two-kink continuum, leading to the appearance of resonances [23]. Once again, while the bubble spectrum admits a semiclassical treatment, the resonances cannot be incorporated into this description.

In the present work, we use a perturbative approach to overcome the aforementioned limitations of the methods used in [23] by determining the excitation spectrum of the resonances and describing the post-quench dynamics of the model.

2.3 Structure of the Hilbert space

The structure of the vector space of the low-energy excitations of the unperturbed Hamiltonian $H_{0}$ is, to a great extent, similar to the infinite antiferromagnetic XXZ spin-1/2 chain in the anisotropic limit [41].

Denote by $\mathcal{L}$ the space of states for the Hamiltonian $H_{0}$, i.e. the Hilbert space spanned by the eigenvectors $\otimes_{j=-\infty}^{\infty}|\mu_{j}\rangle_{j}$ of the Hamiltonian $H_{0}$ that have finite energies (eigenvalues). The Hilbert space $\mathcal{L}$ splits into the direct sum of nine subspaces (sectors):

[TABLE]

The subspace $\mathcal{L}_{\mu\mu^{\prime}}$ is spanned by the basis vectors $\otimes_{j=-\infty}^{\infty}|\mu_{j}\rangle_{j}$, which satisfy the requirement (i) that there exist two integers $j_{1},j_{2}\in\mathbb{Z}$, such that $j_{1}<j_{2}$; and $\mu_{j}=\mu$ for all $j<j_{1}$, while $\mu_{j}=\mu^{\prime}$, for $j>j_{2}$. This means that the relevant configurations interpolate between the spin polarisations $\mu$ and $\mu^{\prime}$. The subspaces $\mathcal{L}_{\mu\mu^{\prime}}$ with $\mu=\mu^{\prime}$, and with $\mu\neq\mu^{\prime}$ are called the topologically neutral, and topologically charged sectors, respectively. Since the local spin operators $P_{j}^{\mu}$ map each sector $\mathcal{L}_{\mu\mu^{\prime}}$ onto itself, different sectors are separate from each other.

In our studies of non-equilibrium dynamics in Section 4, we start from a state polarised in direction $1$. As a result, we are interested in the topologically neutral sector $\mathcal{L}_{11}$, whose structure we describe below in detail.

The Hilbert space $\mathcal{L}_{11}$ is the direct sum of the subspaces $\mathcal{L}_{11}^{(n)}$, with $n=0,2,3,4,\ldots$:

[TABLE]

The basis vector $|\Phi\rangle=\otimes_{j=-\infty}^{\infty}|\mu_{j}\rangle_{j}$ belongs to the subspace $\mathcal{L}_{11}^{(n)}$, if it satisfies the requirement (i) with $\mu=\mu^{\prime}=1$, and

[TABLE]

i.e., it contains $n$ domain walls. The space $\mathcal{L}_{11}^{(0)}$ is one-dimensional spanned by the fully polarised state $|0\rangle^{(1)}\in\mathcal{L}_{11}^{(0)}$, while all subspaces $\mathcal{L}_{11}^{(n)}$ with $n>0$ have infinite dimensions.

The subject of our particular interest is the two-kink subspace $\mathcal{L}_{11}^{(2)}$. The basis in this subspace is formed by the localised two-kink states $|\mathbf{K}_{1,\nu}(j_{1})\mathbf{K}_{\nu,1}(j_{2})\rangle$, with $\nu=2,3$ and $-\infty<j_{1}<j_{2}<\infty$. These states correspond to fully polarised domains joined at single sites and are defined as follows:

[TABLE]

which satisfy the normalisation condition:

[TABLE]

where $\langle\mathbf{K}_{1,\nu}(j_{2})\mathbf{K}_{\nu,1}(j_{1})|$ is the ’bra’-vector corresponding to the ‘ket’-vector $|\mathbf{K}_{1,\nu}(j_{1})\mathbf{K}_{\nu,1}(j_{2})\rangle$. The unit-step translation operator $\hat{T}_{1}$ defined by equation (2.15) acts on the states (2.22) in the following way:

[TABLE]

We also use the parity eigenbasis in $\mathcal{L}_{11}^{(2)}$ formed by the states

[TABLE]

with $\iota=\pm$, denoting the even and odd states respectively.

For the subsequent analysis, it is also essential to introduce another basis in $\mathcal{L}_{11}^{(2)}$ given by the momentum eigenstates (‘Bethe states’)

[TABLE]

where $-\pi<p_{1},p_{2}\leq\pi$ are the quasimomenta of the two kinks, $\nu=2,3$ and $S_{\nu,\nu^{\prime}}(p_{1},p_{2})$ are the two-kink scattering amplitudes defined in Eq. (A.2). We also use the parity-even/odd two-kink Bethe states

[TABLE]

which diagonalise the two-kink scattering matrix (A.4). Yet another basis is obtained by transforming to the centre of mass frame:

[TABLE]

with $\nu=2,3$, $j=1,2,\ldots$, which satisfy

[TABLE]

where $-\pi<P\leq\pi$ is the overall quasimomentum. We can also introduce even/odd combinations of these states by :

[TABLE]

Due to equation (2.21), all states in $|\Phi\rangle\in\mathcal{L}_{11}^{(n)}$ are the eigenstates of the unperturbed Hamiltonian $H_{0}$ with the same eigenvalue $n$. Therefore, to the first order in the transverse magnetic field $g$, it is only necessary to diagonalise the perturbed Hamiltonian (2.11) inside the subspace $\mathcal{L}_{11}^{(2)}$. We note that while higher-order corrections generally mix the two-kink subspace with states of higher kink numbers, previous results show that this mixing is small and a two-kink approximation [8, 25, 40, 24] provides a very accurate description of the spectrum as shown in [23], similarly to the Ising case [8, 25, 40, 24].

Denoting the projection operator onto the two-kink subspace $\mathcal{L}_{11}^{(2)}$ as $\mathcal{P}_{11}^{(2)}$, the restriction of the Hamiltonian $H(\mathbf{v},g)$ to the two-kink subspace can be defined as

[TABLE]

Now we consider the eigenvalue problem

[TABLE]

with $|\Psi(\mathbf{v},g=0)\rangle\in\mathcal{L}_{11}$. Since the Hamiltonian $H(\mathbf{v},g)$ commutes with the translation operator $\hat{T}_{1}$, we can require that

[TABLE]

with quasimomentum $-\pi<P\leq\pi$. The formal expansion of both sides of equation (2.32) in integer powers of $g$ yields:

[TABLE]

Subsequently equating the coefficients in the resulting power series leads to the Rayleigh-Schrödinger perturbation theory [42, 43]. For the deformed vacuum state $|\text{vac}(\mathbf{v},g)\rangle^{(1)}$, the zero-order terms in the perturbation expansions (2.34) read as follows:

[TABLE]

The perturbative calculation of the deformed vacuum states (2.35) for $g\neq 0$ is described in Appendix B.

Turning to the two-kink excitations, one finds that any two-kink state $|\Psi\rangle\in\mathcal{L}_{11}^{(2)}$ is the eigenstate of the unperturbed Hamiltonian $H_{0}$ corresponding to the same energy $E_{0}=2$:

[TABLE]

Switching on $g\neq 0$ in the Hamiltonian (2.11) lifts this degeneracy at the first order in $g$, and the ‘good’ zero-order eigenstates [42, 43] $|\Psi_{0}(\mathbf{v})\rangle$ can be obtained by diagonalising the Hamiltonian $H^{(2)}(\mathbf{v},g)$ defined by (2.31). In the following sections, we diagonalise the Hamiltonian $H^{(2)}(\mathbf{v},g)$ at three different choices of the vector $\mathbf{v}$.

• •

In Subsection 2.4.1, we put $\mathbf{v}=\mathbf{0}$, corresponding to the purely transverse case. In this case, the excitations are two-kink states of the form (2.26).

• •

The case $\mathbf{v}=(v_{1},0,0)$ is studied in Section 2.4.2. In this case, the longitudinal magnetic field $\mathbf{h}=g\,\mathbf{v}$ is either parallel (at $v_{1}>0$) or antiparallel (at $v_{1}<0$) to the orientation of the magnetisation of the first vacuum $|0\rangle^{(1)}$. Following [23], we refer to these two cases as positively and negatively aligned, and the two-kink eigenstates are called mesons/bubbles, respectively. For $v_{1}>0$, the interaction between the two-kinks is attractive, and the excitations are true bound states (a.k.a. mesons), while for $v_{1}<0$, the interaction is repulsive and the excitations correspond to bubbles of the true vacuum forming inside a bulk false vacuum, which are stabilised by Wannier-Stark localisation.

• •

In Subsection 2.4.3, we address the oblique regime $\mathbf{v}=(0,v_{2},0)$, in which the field $\mathbf{v}$ is neither parallel nor anti-parallel to the magnetisation of the state $|0\rangle^{(1)}$. As in [23], we distinguish positive and negative oblique regimes, according to the sign of $v_{2}$. Besides the presence of mesonic excitations for $v_{2}<0$ and bubbles for $v_{2}>0$, we also find unconfined kink excitations in this case. Hybridisation of these states with the meson/bubble states leads to the appearance of resonances in the spectrum. Resonance excitations were observed in [23] using quench spectroscopy; however, they cannot be described by the semiclassical methods used in that paper. We determine the two-kink scattering amplitudes, which we use in Section 3 to calculate the resonance spectrum exploiting their analytic structure following the Fonseca-Zamolodchikov approach [7].

2.4 The spectrum of two-kink states

2.4.1 Pure transverse chain: kink excitations

The Hamiltonian $H(\mathbf{v},g)$ (defined in Eq. (2.11)) with the parameters $\mathbf{v}=\mathbf{0}$ corresponds to the infinite three-state Potts spin chain with purely transverse magnetic field. Here, we discuss the results of the perturbative analysis of this model in the extreme ferromagnetic limit $g\ll 1$ in the two-kink sector. Additional details of the calculation are provided in Appendix A.1.

We start from the eigenvalue problem (2.32) with $\mathbf{v}=\mathbf{0}$. To find the ‘good’ zero-order wave function $|\Psi_{0}(\mathbf{0})\rangle\in\mathcal{L}_{11}^{(2)}$ that stands in the expansion (2.34), one has to diagonalise the Hamiltonian $H^{(2)}(\mathbf{0},g)$. The solution to the latter problem is given by the two-kink ‘Bethe states’ $|{K}_{1,\nu}(p_{1}){K}_{\nu,1}(p_{2})\rangle$ defined by (2.26) or $|{K}(p_{1}){K}(p_{2})\rangle_{\iota}$ defined by (2.27). Both of these two-kink states diagonalise simultaneously the operators $H^{(2)}(\mathbf{0},g)$, and $\hat{T}_{1}$; for example

[TABLE]

where

[TABLE]

is the kink dispersion law to the linear order in the small parameter $g$.

We note that the dynamical properties of the $q$-state Potts spin chain in the presence of the pure transverse magnetic field were previously studied perturbatively by Rapp and Zaránd [28]. Here we go beyond their results by presenting the explicit formulas (A.2) for the two-kink scattering matrix in the case $q=3$, which is essential for the analysis of the resonance spectrum performed in Section 3.

2.4.2 Two-kink states obtained from $\mathcal{L}_{11}^{(2)}$: the aligned case

Application of the longitudinal magnetic field $\mathbf{v}=(v_{1},0,0)$, with $v_{1}=h_{1}/g$, explicitly breaks the $S_{3}$ symmetry of the Hamiltonian (2.11), with the effect of this field on the energies of the vacua of the model shown schematically in Fig. 2.1. For $h_{1}>0$, the energy of the first vacuum $|\text{vac}(\mathbf{v},g)\rangle^{(1)}$ decreases with respect to the energies of the second and third ones $|\text{vac}(\mathbf{v},g)\rangle^{(\nu)}$, $\nu=2,3$. As a result, the state $|\text{vac}(\mathbf{v},g)\rangle^{(1)}$ becomes the true ground state of the spin chain, while for $h_{1}<0$ it becomes the metastable (false) vacuum, with an energy larger than those of the two degenerate ground states $|\text{vac}(\mathbf{v},g)\rangle^{(2,3)}$. In both cases, the application of the weak longitudinal magnetic field leads to the formation of localised two-kink states a.k.a. ’mesons’ and ’bubbles’. In [23], these excitations were described in the semiclassical approximation. For $h_{1}>0$, the mesons correspond to the confinement of the kinks, while the bubbles appearing for $h_{1}<0$ correspond to the nucleation of the true vacuum. These states are fully analogous to the mesons and bubbles obtained in the Ising spin chain [8, 40, 24].

Here we briefly discuss the spectrum of these two-kink bound states to the leading order in the weak transverse magnetic field $g$, with the details of the calculation relegated to Appendix A.2. We use the short-cut notation $H^{(2)}(v_{1},g)$ for the Hamiltonian (2.31). The eigenstates can be expressed in terms of the basis $|j,P\rangle_{\iota}$ defined by (2.30) as follows:

[TABLE]

Together with the boundary condition that the wave-function vanishes as $j\rightarrow\infty$, the above equation determines the energy eigenvalues. The resulting energy spectrum is discrete, and we obtain in Appendix A.2 the dispersion laws $E_{\iota,n}(P,h_{1})$ of mesons at $h_{1}>0$, and bubbles at $h_{1}<0$ to the first order in the small parameter $g$.

The energy spectra $E_{\pm,n}(P,h_{1})$ of the three lightest mesons at $g=0.2$, $h_{1}=0.05$ are shown in Fig. 2.2. The boundaries of the two-kink continuum at $h_{1}=0$ are given by

[TABLE]

using Eq. (2.39).

A critical feature of the two-kink spectrum is the distinction between collisional and collisionless states [24, 23]. This distinction is most intuitively understood through a semiclassical lens.

• •

Mesonic excitations ($h_{1}>0$): For the lightest (lowest-energy) mesons, the constituent kinks follow trajectories that cross, resulting in periodic mutual collisions. As one moves to higher-energy mesonic states, the attractive longitudinal field induces Wannier-Stark localisation [16, 17]. This effect confines the constituent kinks to separate regions, preventing them from colliding. In this collisionless regime, the kinks undergo independent Bloch oscillations [18]. Consequently, the energy levels become uniformly spaced, with the gap determined by the Bloch oscillation frequency.

• •

Bubble excitations ($h_{1}<0$): The above scenario is mirrored for bubbles, with the primary difference being that the bubble spectrum is bounded from above [40]. This results in an inverted energy hierarchy in the bubble sector, where the highest-energy states are collisional, whereas the lower-energy states correspond to collisionless configurations stabilised by Wannier-Stark localisation.

The collisional particles lie inside the interval of the two-kink continuum at $h_{1}=0$ [23], i.e., between $E_{\text{min}}(P)$ and $E_{\text{max}}(P)$.

2.4.3 Two-kink states obtained from $\mathcal{L}_{11}^{(2)}$: the oblique case

This section addresses the oblique regime, defined by the rescaled longitudinal field $\mathbf{v}=(0,v_{2},0)$. This configuration is uniquely characteristic of the three-state Potts model and has no counterpart in the Ising spin chain. The presence of a non-zero $v_{2}$ deforms the three ferromagnetic vacua and shifts their energies, as illustrated in Fig. 2.3. A critical feature of this regime is that it allows for the coexistence of localised two-kink states (mesons for $v_{2}<0$ and bubbles for $v_{2}>0$) with unconfined kink excitations [23], which interpolate between the degenerate vacuum states. The hybridisation of these bound states with the two-particle continuum leads to the formation of resonances, which are considered in detail in Section 3.

Throughout this Section, the notation $H^{(2)}(v_{2},g)$ for the Hamiltonian (2.31) is used instead of $H^{(2)}(\mathbf{v},g)$. The details of the perturbative calculation are relegated to Appendix A.3.

To determine the spectrum, we solve the eigenvalue problem restricted to the two-kink subspace $\mathcal{L}_{11}^{(2)}$. At total quasimomentum $P$, we expand the eigenstate in terms of the basis $|j,P\rangle_{\nu}$ determined by (2.28):

[TABLE]

As shown in Appendix A.3, the above equation determines that the wavefunction coefficients have the form

[TABLE]

and

[TABLE]

where

[TABLE]

The excitations below the two-kink continuum ($E<E_{\text{min}}(P)$) correspond to collisionless bubbles, while those above the two-kink continuum ($E>E_{\text{max}}(P)$) are collisionless mesons.

Eq. (2.42) furthermore implies the boundary conditions

[TABLE]

where for simplicity we suppressed $E$ and $P$ from the arguments of the function $\psi_{\nu}(j;E,P)$. At $j\to\infty$, the wave-function $\psi_{2}(j)$ must vanish, while the absolute value of $\psi_{3}(j)$ must remain bounded.

In the interval $E_{\text{min}}(P)<E<E_{\text{max}}(P)$, substitution of (2.4.3), (2.44) into the equations of the boundary conditions (2.46) leads to two linear equations, determining the coefficients $A$ and $B$. For $P=0$, the wavefunction $\psi_{3}(j;E,P)$, given by Eq. (2.44), can be represented in the form:

[TABLE]

Then Eq. (2.46) yields

[TABLE]

where $z=e^{ip}$ with $p$ taking real values, $v_{2}=\frac{h_{2}}{g}$, and $c_{0}(z,v_{2})=\frac{2(z+z^{-1})}{3v_{2}}$. We remark that the symmetry (A.33) implies

[TABLE]

using that $v_{2}$ is real. Finally, we note that the wave-function $\psi_{2}(j;E,P=0)$ can be simply written as

[TABLE]

Note that the states in this range are two-kink scattering states and therefore their energy spectrum is continuous. While standard semiclassical methods fail in this regime, our perturbative approach allows us to identify resonance peaks appearing in the quench spectroscopy analysis performed in Section 4, by analysing the analytic structure of the two-kink scattering amplitude, as detailed in Section 3.

When the energy $E$ is outside the two-kink continuum, substituting Eqs. (2.4.3,2.44) into the boundary condition equation (2.46) yields two linear secular equations that determine the coefficient $C$ and the discrete energy spectrum $\{E_{n}(P,h_{2})\}_{n=1}^{\infty}$. In the negative oblique regime $h_{2}<0$, this energy spectrum corresponds to collisionless mesons,

[TABLE]

with $E_{1}(P,h_{2})>E_{\text{max}}(P)$, while in the positive oblique regime $h_{2}>0$, it corresponds to collisionless bubbles,

[TABLE]

with $E_{1}(P,h_{2})<E_{\text{min}}(P)$.

3 Resonance excitations in the oblique regime

As discussed above, a remarkable feature of the three-state Potts quantum spin chain is the coexistence of localised two-kink bound states—mesons and bubbles—with a continuum of stable, unconfined kinks in the oblique regime.

In this section, we provide a detailed analysis of the two-kink scattering amplitude in the small-$g$ limit to describe the interactions within this mixed particle content. We calculate the energies of stable bound states and the positions of resonances by analysing the poles of the scattering amplitude in the complex plane. We then demonstrate how resonances arise from the hybridisation of the unconfined two-particle continuum with bound states that would otherwise be stable in a purely confining regime. By tracking the trajectories of these poles as the longitudinal field parameter $v_{2}$ is varied, we describe the physical transformation of stable mesons and bubbles into unstable resonances. This analysis confirms that the three-state Potts model follows the general scenario proposed by Fonseca and Zamolodchikov for the Ising field theory, where heavy stable mesons cross a stability threshold and decay into the lightest available particles [7]. These results provide the theoretical basis for the quench spectroscopy features observed in Section 4, which remain beyond the reach of semiclassical approximations.

3.1 The bound state and resonance spectrum

In the oblique regime — defined by $h_{2}\neq 0$ and $h_{1}=h_{3}=0$ — the system’s dynamics is governed by localised two-kink states and a continuum of unconfined kinks. At zero total quasimomentum ($P=0$), we consider the wave-function $\psi_{3}(j;E,P=0)$ as defined in (2.47). The scattering amplitude $S(z,v_{2})=\frac{B_{in}(z^{-1},v_{2})}{B_{in}(z,v_{2})}$ of two unconfined kinks is given by:

[TABLE]

The scattering states correspond to $z=e^{ip}$ with the relative momentum variable $p$ taking a real value, i.e. $|z|=1$.

To find the spectrum of bound states, we analytically continue the scattering amplitude to the complex $z$-plane, where it is a meromorphic function for a fixed real value of $v_{2}\neq 0$. It has two classes of poles:

• •

Stable bound states: poles located on the real interval $z_{n}\in(-1,1)$ correspond to stable two-kink bound states, with energies:

[TABLE]

• –

In the positive oblique regime ($v_{2}>0$), these poles are positive ($z_{n}\in(0,1)$), representing collisionless bubbles.

• –

In the negative oblique regime ($v_{2}<0$), these poles are negative ($z_{n}\in(-1,0)$), representing collisionless mesons.

• •

Resonance states: poles $z_{l}$ that lie in the upper half-plane ($\text{Im }z_{l}>0$) outside the unit circle ($|z_{l}|>1$) determine the positions and widths of resonances

[TABLE]

where the real/imaginary parts give the energy/width, respectively.

Numerical values for these energies at $g=0.2$ and $P=0$ are detailed in Table 3.1. The energy dependences of the wave-function amplitude $\sqrt{A^{2}+B^{2}}$ and the scattering phase $\alpha=-i\log S$ are illustrated in Figure 3.1. As shown in Fig. 3.1, the stable collisionless particles appear as discrete energy levels (red gridlines) outside the continuous two-kink spectrum (shaded blue region). Conversely, resonance positions (blue gridlines) coincide with broader peaks within the continuum, where hybridisation between bound states and unconfined kinks occurs.

The complex energies $E_{\text{res},l}$ presented in Table 3.1 provide a quantitative basis for the quench spectroscopy features observed in the iTEBD simulations as we discuss later in Section 4. While the stable bound states ($E_{n}$) appear as sharp, delta-like peaks in quench spectroscopy, the resonance states are characterised by a finite decay width, determined by the imaginary part of the pole position in the complex energy plane.

3.2 Trajectory of poles and zeroes of the scattering amplitude

To understand the physical transformation of stable excitations into resonances, it is instructive to trace the evolution of the poles $p_{n}(v_{2})$ and zeroes $k_{n}(v_{2})$ of the two-kink scattering amplitude $\mathcal{S}(p,v_{2})=S(z,v_{2})|_{z=\exp(ip)}$ as the oblique parameter $v_{2}$ is tuned. Since the locations of poles $p_{n}(v_{2})$ of the scattering amplitude determine the energies $E_{n}(v_{2},g)$ of the two-kink bound states,

[TABLE]

this gives us insight into the behaviour of these energies near the stability threshold, and elucidates the transformation of the stable mesons (at $v_{2}<0$) and bubbles (at $v_{2}>0$) into the resonances.

We represent the scattering amplitude in terms of the Jost function $F(p,v_{2})$:

[TABLE]

where the Jost function is defined as:

[TABLE]

This function is $2\pi$-periodic in $p$, and also analytic in the complex $p$-plane at any fixed real $v_{2}\neq 0$. The poles of the scattering amplitude (3.5) are located at zeroes of the Jost function (3.6):

[TABLE]

In the following, we focus on the positive oblique regime $v_{2}>0$. The analysis of the negative oblique regime $v_{2}<0$ is entirely parallel, since the location of poles $p_{n}(v_{2})$ in the negative oblique regime $v_{2}<0$ can be found from the equality

[TABLE]

that follows from the symmetry relation $F(p+\pi,-v_{2})=-F(p,v_{2})$.

At sufficiently large values of $v_{2}$, the zeroes of the Jost function $F(p,v_{2})$ in the strip $-\pi<\mathrm{Re}\,p\leq\pi$ are exclusively on the imaginary $p$-axis, as illustrated in Figures 2(a), 2(b) for $v_{2}=2$. The left subfigure 2(a) shows the momentum dependences of the Jost functions $F(p,v_{2})$, and $F(-p,v_{2})$ at purely imaginary momenta $p=iq$. According to (3.5), the zeroes of these functions determine the locations of the poles and the zeroes of the scattering amplitude $\mathcal{S}(p,v_{2})$, which are schematically displayed in Fig. 2(b). These figures demonstrate that the poles and zeroes of the function $\mathcal{S}(p,v_{2})$ alternate on the imaginary $p$-axis forming pairs. In the ‘physical’ upper half-plane $\mathrm{Im}\,p>0$, the poles $p_{n}(v_{2})$ and zeroes $k_{n}(v_{2})$ are ordered as:

[TABLE]

In the lower half-plane $\mathrm{Im}\,p<0$, poles and zeroes of the scattering amplitude $\mathcal{S}(p,v_{2})$ are located at $\{-k_{n}(v_{2})\}_{n=1}^{\infty}$, and at $\{-p_{n}(v_{2})\}_{n=1}^{\infty}$, respectively, as required by the unitarity relation:

[TABLE]

Upon decreasing the parameter $v_{2}$, all poles and zeroes of the scattering amplitude move along the imaginary $p$-axis towards the origin $p=0$. When the parameter $v_{2}$ reaches its first critical value $a_{1}=1.47112\ldots$, the pole $p_{1}(v_{2})$ arrives at the origin $p_{1}(a_{1})=0$, where it merges with the zero of the function $\mathcal{S}(p,v_{2})$ arriving from the lower half-plane. This configuration of poles and zeroes of the function $\mathcal{S}(p,a_{1})$ is shown in Fig. 2(d). Fig. 2(c) illustrates the $q$-dependencies of the Jost functions $F(iq,v_{2})$ and $F(-iq,v_{2})$ at $v_{2}=a_{1}$.

Note, that at $v_{2}=a_{1}$, the energy of the first bubble $E_{1}(v_{2},g)$ reaches its stability threshold, which coincides with the lower bound of the two-kink continuum spectrum:

[TABLE]

Since the function $p_{1}(v_{2})$ remains analytic in $v_{2}$ in the vicinity of the critical point $a_{1}$, it varies at $v_{2}\to a_{1}$ as:

[TABLE]

with some $c_{1}>0$. Due to (3.4), the difference $E_{1}(v_{2},g)-E_{\text{min}}(g)$ vanishes quadratically in $(v_{2}-a_{1})$, when $v_{2}$ approaches the critical value $a_{1}$ from above:

[TABLE]

Fig. 2(e) displays the profile of the Jost functions at $v_{2}=1.468$, which is slightly below the first critical value $a_{1}$, with the corresponding configuration of poles and zeroes shown in Fig. 2(f). The pole $p_{1}(v_{2})$ drops into the lower half-plane $\mathrm{Im}\,p<0$ remaining purely imaginary. Formulas (3.12), and (3.13) still can be used in this case, but no stable excitation exists at $v_{2}<a_{1}$, with the energy given by Eq. (3.13) corresponding to the so-called ‘virtual level’, see e.g [44].

As the parameter $v_{2}$ is decreased further to the second critical value $v_{2}=b_{1}=1.4646\ldots$, two poles $p_{1}(v_{2})$ and $-k_{1}(v_{2})$ of the scattering amplitude $\mathcal{S}(p,v_{2})$ merge at the negative $p$-half-axes, as it is shown in Fig. 2(h), with the functions $p_{1}(v_{2})$, $k_{1}(v_{2})$, and $E_{1}(v_{2},g)$ having a square-root branching point at $v_{2}=b_{1}$.

Decreasing $v_{2}$ further leads to the situation shown in Figure 2(j), where the two poles shown by red crosses, which are located in the lower half-plane, are shifted away from the imaginary $p$-axis by gaining non-zero real parts $\sim\pm\sqrt{b_{1}-v_{2}}$. The corresponding energy $E_{1}(v_{2},g)$, in turn, gains a negative imaginary part $\sim-\sqrt{b_{1}-v_{2}}$, indicating that the virtual state transformed into a resonance.

As the parameter $v_{2}$ decreases further, the above scenario repeats and the bubbles enumerated by the indices $n=2,3,\ldots$ transfer one by one into resonances. Fig. 3.3 shows the resulting configuration of poles and zeroes of the scattering amplitude $\mathcal{S}(p,v_{2})$ at $v_{2}\simeq 0.5$ schematically.

4 Evolution of magnetisation after a quantum quench

In this section, we apply a perturbative analysis to the post-quench dynamics of the three-state Potts spin chain. Our goal is to calculate the time-dependent local magnetisation, $M_{\mu}(t)$, up to the second order in the transverse field $g$.

Before delving into the calculation itself, we note that perturbative approaches to quantum quenches have already been explored in the context of massive quantum field theories [45, 46, 47, 48, 49, 50]. This approach works for small quenches, where the quench magnitude can be defined as the post-quench energy density relative to the gap scale. One of the main drawbacks of this method is that it is based on the pre-quench Hilbert space, whereas the time evolution of the quench is determined by the post-quench spectrum.

An alternative approach is to expand in the post-quench Hilbert space using a form-factor approach [51, 52], which has two main advantages. First, the time evolution has the correct frequency spectrum of the post-quench energy system [47]. Second, resummation of secular contributions to the time evolution reproduces the long-time asymptotic behaviour of observables, most notably their exponential relaxation, and also gives a systematic perturbative expansion for the relaxation time. The relaxation time obtained this way for the Ising field theory was found to agree well with numerical simulations [25]. The drawback of the post-quench form factor expansion is that it requires knowledge of the energies of eigenstates of the post-quench Hamiltonian, their overlaps with the initial states, as well as the matrix elements of observables in the eigenstate basis. This information is available when both the post-quench Hamiltonian is integrable, and the initial state corresponds to a so-called integrable quench [53, 54], which is the case considered in [51, 52].

In order to demonstrate the efficiency of our perturbation technique, we apply it in Appendix D to describe the evolution of the magnetisation in the Ising spin chain in the limit of the weak transverse magnetic field in the post-quench Hamiltonian. The comparison of the perturbative results with the well-known form-factor result clarifies the way in which our perturbative expansion works. Essentially, it can only capture the time evolution on sufficiently short time scales, provided the quench itself is suitably small.

Since the Potts model is non-integrable, we can only resort to a perturbative approach to construct its time evolution. Despite its limitations, we demonstrate that it agrees reasonably well with the numerically computed time evolution in Ref. [23]. In addition, we also demonstrate that the spectrum in the two-kink sector, constructed above in Sections 2 and 3, shows an excellent agreement with the Fourier spectrum obtained by quench spectroscopy.

4.1 Perturbation theory for Potts quench dynamics

We explore essentially the same quench protocol as that studied in [23]. We start from the finite Potts spin chain having an even number of sites $N$, whose Hamiltonian is given by a simple modification of (2.11):

[TABLE]

where the constant shift is specified by $c^{1}=2/3$, and $c^{2}=c^{3}=-1/3$. The initial state is purely ferromagnetic along the direction $1$:

[TABLE]

For further calculations, let us denote by $|\text{vac}(\mathbf{v},g|N)\rangle$ the deformation of the vacuum $|0\rangle^{(1)}$ of the Hamiltonian $H_{N}(\mathbf{0},0)$:

[TABLE]

It is convenient to represent the deformed vacuum $|\text{vac}(\mathbf{v},g|N)\rangle$ in the form:

[TABLE]

where the vector $|\Omega(\mathbf{v},g|N)\rangle$ admits the following asymptotic expansion at $g\to 0$:

[TABLE]

such that

[TABLE]

The details of the perturbative analysis of the true/false vacuum states are relegated to the Appendix B.

At $t>0$, the average value $A_{\alpha}(t|N)$ of the observable corresponding to the operator $\hat{A}_{\alpha}$ reads:

[TABLE]

We consider the following three operators, which are basically the three magnetisations shifted by a constant:

[TABLE]

where $\hat{M}_{\alpha}$ are the three magnetisation operators with $\alpha=1,2,3$. Since

[TABLE]

only two of the operators are independent. The operators $\hat{A}_{\alpha}$ are preferred over the magnetisations $\hat{M}_{\alpha}$ because they satisfy the relation

[TABLE]

which leads to a substantial simplification of the calculation of the expectation value (4.6). The perturbative calculation of the aforementioned expectation value is discussed in detail in Appendix C. Here, we only summarise the key ideas and the results of the computation.

To calculate the $N\rightarrow\infty$ limit, necessary for the subsequent comparison with the numerical results, we introduce the following notations as well

[TABLE]

After taking the thermodynamic limit, we apply a further key simplification, provided by the two-kink approximation, which essentially corresponds to replacing the Hamiltonian $\Delta H(\mathbf{v},g)$ by its restriction to the subspace $\mathcal{L}_{11}^{(2)}$:

[TABLE]

Using the results of Appendix C, this approximation implies that

[TABLE]

The first term on the right-hand side of (4.12) is time independent and can be expressed as shown in Eq. (C.16). The second term describes the resonant oscillations of the magnetisations $A_{\alpha}^{(2)}(t)$ with high frequencies $\omega\simeq 2+O(g)$, while the last one represents the low-frequency contribution to the quantities $A_{\alpha}^{(2)}(t)$.

We expect that at $g\ll 1$, the difference between $A_{\alpha}(t)$ and $A_{\alpha}^{(2)}(t)$ remains small during the first stage of the evolution, i.e. for not too large $t$. However, at sufficiently large $t$, the multi-kink contributions to the state $|\Psi(t)\rangle=e^{-it\,\Delta H(\mathbf{v},g)}|\Omega_{1}(\mathbf{v})\rangle$ should become considerable, as the full Hamiltonian does not conserve the number of kinks. This leads to a substantial deviation of $A_{\alpha}(t)$ from its two-kink approximation $A_{\alpha}^{(2)}(t)$. We return to this issue later.

In the rest of this Section, we concentrate on the oblique regimes and put $\mathbf{v}=(0,h_{2}/g,0)$. We use the notation $H^{(2)}(h_{2},g)$ for the Hamiltonian, instead of $H^{(2)}(\mathbf{v},g)$. The calculation of the second and third terms of the right-hand side of the formula (4.12) in this regime is explained in detail in Appendix C. The resulting post-quench evolution of the magnetisations $M_{2}(t)$, and $M_{3}(t)$ is displayed in Fig. 4.1 at $g=0.2$ in the positive (a), and negative (b) oblique regimes with $h_{2}=0.1$, and $h_{2}=-0.1$, respectively. As one can see from Fig. 4.1, the magnetisation $M_{3}(t)$ linearly increases at large $t$ in the initial stage of the post-quench evolution, which is described by our perturbative approach:

[TABLE]

The coefficient $c$ that stands in this formula depends, of course, on the Hamiltonian parameters $h_{2}$ and $g$. By the means of Appendix C, the explicit form of $c$ reads as

[TABLE]

where the wave-functions $\psi_{\nu}(j;E,P=0)$ are determined by the equations (2.47) and (2.51), $E(p)=2-\frac{4g}{3}\cos p$ and $B(p)$ is defined as

[TABLE]

where the function $B_{in}(z,v_{2})$ is given by (2.48). For the values of the parameters of the Hamiltonian chosen in Fig. 4.1, one finds:

[TABLE]

For later use, we also include the Fourier transform of the expression (4.12). Here, we only focus on the second term on the right-hand side of (4.12), which is the only contributing term in the relevant frequency domain, as explained in the next section.

Using the results of Appendix C, the real part of the Fourier transform

[TABLE]

for the relevant frequency domain $\omega\in(2-\frac{4g}{3},2+\frac{4g}{3})$, takes the form:

[TABLE]

where

[TABLE]

4.2 Comparison to numerical simulations

4.2.1 Time evolution of the magnetisation

To validate the perturbative results, we realised the quench protocol described at the beginning of Section 4 using the infinite Time Evolving Block Decimation (iTEBD) method following Ref. [23]. The details of the simulations are provided in Appendix E.

The initial state is chosen as $|0\rangle^{(1)}$ given by the $N\rightarrow\infty$ limit of (4.2), which is evolved in time by the Hamiltonian (2.11) with exactly one non-zero finite longitudinal field in each case, as discussed previously. We examine four distinct scenarios: $h_{1}>0$ (positively aligned), $h_{1}<0$ (negatively aligned), $h_{2}>0$ (positive oblique) and $h_{2}<0$ (negative oblique). During the time evolution, we calculated the local magnetisations

[TABLE]

In aligned quenches, symmetry reduces the system to have a single independent magnetisation (chosen as $M_{1}$), whereas oblique quenches require two (typically chosen as $M_{2}$ and $M_{3}$). The underlying reason is that the three magnetisations sum up to zero, and in the aligned cases, the transformation swapping the spin directions 2 and 3 (leaving 1 intact) is a symmetry of the quench.

In Fig. 4.2, we show the independent magnetisations $M_{2}(t)$ and $M_{3}(t)$ in the case of oblique quenches calculated numerically with the iTEBD method for a specific parameter setting in each case, together with the analytical prediction for the magnetisations in the two-kink approximation that was determined by (4.12). The perturbative analytical results show clear convergence toward the full iTEBD numerical data as the transverse field $g$ is reduced from $0.2$ to $0.05$. However, this convergence is notably slower in the Potts model than in the Ising spin chain. (see D.1). The accuracy of the perturbative results for the time evolution of $M_{\mu}(t)$ is highly sensitive to the magnitude of $g$ (as illustrated in Fig. 4.2), while the spectral features (quench spectroscopy) remain significantly more robust, as demonstrated later in this section.

4.2.2 Quench spectroscopy in aligned quenches

The quasi-particle content relevant for the quench dynamics can be obtained from the Fourier transform of the magnetisations [15, 23]:

[TABLE]

where $N_{s}$ is the total number of the discrete time steps of the simulation, $\delta t$ is the time step, and $k$ is an integer such that $k\in[0,\ldots N_{s}-1]$. The specific values used in the simulations are provided in Appendix E.

In the following, we plot the numerically calculated $|M_{\mu}(\omega)|^{2}$ power spectra. Due to translational invariance, only eigenstates with total momentum $P=0$ contribute to the Fourier transform. In the following plots, the spectrum of the $P=0$ eigenstates obtained by perturbation theory is shown by vertical lines. Besides comparing the peak locations with the analytical results, we also explain the selection of the eigenstates expected to contribute in the different cases. A more detailed discussion, together with exact diagonalisation results for the spectrum, can be found in Ref.[23].

Positive aligned quenches -

Fig. 4.3 shows the typical Fourier spectrum of $M_{1}(\omega)$ for the case of a positively aligned quench.

Turning on a positive $h_{1}$ results in a doubly degenerate false and a single true vacuum state. As the initial state $|0\rangle^{(1)}$ is favoured by the longitudinal field, we expect excitations built upon the true vacuum to contribute to the quench with the largest overlaps. These are predominantly meson states, described in Section 2.4.2, with their energy determined by the solutions of the perturbative secular equation (A.24). Their overlap with the initial state is expected to decrease with their energy [40, 24]. Furthermore, because of the $\mathbb{Z}_{2}$ symmetry of the quench, only the $\iota=+$ (even) meson states contribute, and therefore we show the first few $\iota=+$ meson masses ($P=0$ energy eigenvalues), calculated from (A.24) using black vertical lines in Fig. 4.3. The locations of the peaks in the Fourier spectrum match the meson masses with reasonable precision. The grey region indicates the range of the continuous spectrum of two-kink states at zero longitudinal field. The boundaries of the region are determined by relation (2.41) at $P=0$:

[TABLE]

As already mentioned before, the meson excitations falling in this interval are collisional, i.e., the constituent kinks scatter on each other, while the ones outside are non-collisional, which means they consist of a pair of kinks localised by Bloch oscillations [24, 23]. We note that the exact diagonalisation (ED) results for the meson masses (with $L=10$ and periodic boundary conditions) match better with the peaks in $M_{1}(\omega)$, and they even allow for the identification of the additional small, high-energy peaks in the Fourier transform as baryons (three-kink bound states) [23]. However, the exact diagonalisation method suffers from other limitations: it is not sufficient for a complete classification of the eigenstates; moreover, it cannot be applied to locate resonances in oblique cases, as discussed later.

Negative aligned quenches -

In the negatively aligned case, the typical Fourier spectrum of $M_{1}(\omega)$ resulting from the simulations is illustrated in Fig. 4.4.

In this case, there are two degenerate true vacua and a single false vacuum. The initial state $|0\rangle^{(1)}$ is disfavoured by the longitudinal field; therefore, the excitations built upon the false vacuum are expected to dominate the quench. These are bubble states that correspond to the nucleation of a true vacuum domain, with their energies determined by the solutions of Eq. (A.24). Again, due to the aforementioned symmetry of the initial state, only the $\iota=+$ (even) bubbles contribute to the time evolution. The energy of the bubbles decreases with their size, i.e., with their species label; therefore, their overlap with the initial state increases with their energy [40, 24]. The first few (highest energy) $\iota=+$ bubble masses, calculated from (A.24), are plotted in Fig. 4.4 by black vertical lines which again coincide well with the locations of the peaks in $M_{1}(\omega)$. The grey-coloured region in the figure again denotes the two-kink continuous spectrum at $h_{1}=h_{2}=0$ with its boundaries determined by (4.22). Just as before, the excitations inside this interval are collisional, while the ones outside are non-collisional bubbles. Again, we remark that the bubble masses calculated via ED (with $L=10$ and periodic boundary conditions) match better with the peaks in the Fourier transform. Nevertheless, we repeatedly emphasise that the bubble energy eigenstates cannot be identified from the ED states without the help of the perturbative results or semiclassical quantisation, as done in Ref.[23]. The ED energy spectrum additionally enables us to identify the additional small peaks in the Fourier spectrum as baryonic bubbles (3-kink localised states built upon the false vacuum) [23].

4.2.3 Quench spectroscopy in oblique quenches

In Ref. [23], the positively and negatively aligned quenches considered above were also referred to as (standard) confining and anticonfining quenches, respectively, based on the main mechanism of localisation. In contrast, the oblique regimes are characterised as ”partial confining” or ”partial anticonfining” due to the coexistence of localised bound states and unconfined kink excitations. This mixed particle content leads to a complex spectral structure that cannot be captured by semiclassical methods but is precisely described by our perturbative framework. As mentioned, aligned quenches can also be realised in the Ising case; however, the oblique quenches have no Ising analogues.

Positive oblique quenches -

In the positive oblique case, we chose the two independent magnetisations as $M_{1}(t)$ and $M_{2}(t)$. The typical Fourier spectrum for $M_{1}(t)$ is visible in Fig. 4.5. $M_{2}(\omega)$ is qualitatively very similar and exhibits peaks at the same locations.

The presence of $h_{2}>0$ leads to a single true vacuum and doubly degenerate false vacua. The initial state $|0\rangle^{(1)}$ is disfavoured by the longitudinal field, and therefore, the elementary excitations built upon the false vacua have the largest overlaps with the initial state. The calculation of the energy spectrum of such excitations is discussed in Sections 2.4.3 and 3. As discussed previously in detail, outside the interval $(E_{\text{min}}(0),E_{\text{max}}(0))$ (where $E_{\text{min}}(0)$ and $E_{\text{max}}(0)$ are given by (4.22)) we get a discrete energy spectrum (2.53) describing collisionless bubbles. These excitations appear in the Fourier spectrum as sharp peaks. In the interval $(E_{\text{min}}(0),E_{\text{max}}(0))$, the resulting spectrum is continuous and admits resonances (unstable collisional bubbles) as explained in Section 3. Nevertheless, at the resonance locations, we still expect wide resonant peaks. Note that the resonance peaks are notably wider and less symmetric than their collisionless counterparts. As the resonance energy moves deeper into the two-kink continuum, the increasing imaginary part causes the peak to blend into the background spectral density, eventually becoming a broad hump or a cusp at the stability threshold. For instance, at $h_{2}=0.1$, the results in Table 3.1 show that the first resonance $E_{\text{res},1}$ has a smaller imaginary part compared to the second resonance $E_{\text{res},2}$. This accurately predicts why the first resonance peak in Fig. 5(a) appears sharper than the second one.

In Fig. 4.5, we show the analytically calculated collisionless bubble masses and the locations of resonances by black gridlines. The grey region again represents the two-kink continuous spectrum at $h_{1}=h_{2}=0$, lying within the interval ($E_{\text{min}}(0),E_{\text{max}}(0)$). The black lines within this region denote the resonances (unstable collisional bubbles), while the lines outside this interval stand for the collisionless bubbles. The positions of the lines match with the locations of the peaks in $M_{1}(\omega)$ (and $M_{2}(\omega)$) with good precision. The boundary (threshold) of the continuous energy spectrum corresponds to the appearance of a cusp in the Fourier spectrum, which is observed numerically, and its location approximately agrees with the analytical prediction. We remark that the collisionless bubble masses and the location of the cusp can be determined more precisely by ED; however, the ED method cannot locate the resonances.

Negative oblique quenches -

In the negative oblique case, the dynamics and its explanation is similar to the positive oblique case, with the bubbles replaced by mesons. The typical Fourier spectrum of $M_{1}(t)$ is shown in Fig. 4.6, while the other independent magnetisation, chosen as $M_{2}(t)$, has the same qualitative behaviour [23].

Setting $h_{2}<0$ results in a single false vacuum and doubly degenerate true vacua, meaning that $|0\rangle^{(1)}$ is favoured by the longitudinal field, thus the quasi-particles built upon the true vacua contribute dominantly to the quench. The corresponding energy spectrum was calculated in Sections 2.4.3 and 3. Outside the interval $(E_{\text{min}}(0),E_{\text{max}}(0))$ (marked by a grey region in the figure), the energy spectrum is discrete (2.52) and describes collisionless mesons which are present in the quench spectroscopy as sharp peaks. In the interval $(E_{\text{min}}(0),E_{\text{max}}(0))$, the spectrum is continuous and exhibits resonances (unstable collisional mesons). Similarly to our analysis in the positive oblique case, for $h_{2}=-0.1$, the results in Table 3.1 show that the first resonance $E_{\text{res},1}$ has a larger imaginary part compared to the second resonance $E_{\text{res},2}$ which explains why the second resonance peak in Fig. 6(a) appears sharper than the first one.

The edge of the continuous spectrum again corresponds to a cusp in the Fourier spectrum. Once again, the ED meson masses (and the location of the cusp computed by ED) are more precise, but the method fails to determine the positions of the resonances for reasons analogous to those encountered in the positively aligned case.

In the following, we compare the Fourier transform of the magnetisation $M_{1}(\omega)$ calculated numerically in the simulations for the oblique cases with the analytical results in Eq. (4.18). We expect the Fourier transform of $A_{\alpha}(t)$ given by (4.12) to match the corresponding numerically calculated magnetisation $M_{\alpha}(\omega)$ in the relevant frequency domain of the low-energy excitation spectrum up to a constant shift $\text{d}M_{\alpha}$ and a constant multiplicative factor $a$. The constant shift was set manually, while the multiplicative factor originates from the different definitions of the Fourier transforms as discussed in Appendix E. The difference between the definition of $A_{\alpha}(t)$ and $M_{\alpha}(t)$ (4.7) becomes irrelevant as they differ by a constant, which results in an additional peak at zero frequency in the Fourier spectrum, which is outside of the interval of interest. Furthermore, the first term on the right-hand side of (4.12) is time-independent, which again leads to a peak at zero frequency. Meanwhile, the last term represents the low-frequency contribution to $A_{\alpha}(\omega)$ and also gives rise to peaks outside the investigated interval. Therefore, in the quench spectroscopy, we neglect these terms and only focus on the second term, whose Fourier transform is given by (4.18). In Figs. 4.7 and 4.8 we show typical Fourier spectra for $a\,\text{Re}F^{(1)}_{1}(\omega)+\text{d}M_{1}$ together with the numerically calculated $\text{Re}M_{1}(\omega)$ in the oblique cases; further examples are included in Appendix E. The grey vertical lines mark the quasi-particle masses and resonance frequencies determined from perturbation theory. The plots demonstrate an excellent match between the analytic and numerical results, despite the presence of numerical noise in the iTEBD results.

5 Conclusions

In this work, we have developed a comprehensive perturbative description of the dynamics of the three-state Potts quantum spin chain in the extreme ferromagnetic regime, utilising the transverse magnetic field as the expansion parameter.

The primary achievement of this study is to provide analytical access to spectral and dynamical features beyond the reach of standard semiclassical and exact diagonalisation methods. While we recovered the known meson and bubble excitation spectra in the aligned cases, our approach proved particularly powerful in the oblique regime. In this case, switching on the longitudinal field in the direction of one vacuum leaves the remaining two degenerate, allowing for unconfined kink excitations [23].

Beyond the bound state spectrum, we explored the analytic structure of the two-kink scattering amplitude. By tracking the trajectory of bound-state poles near stability thresholds, we demonstrated that the three-state Potts model follows the general scenario proposed by Fonseca and Zamolodchikov for the Ising field theory, in which stable excitations transform into resonances [7]. We successfully reproduced the resonances arising from the hybridisation of bound states with the unconfined continuum, matching quench-spectroscopy results with high precision.

Furthermore, we provided an explicit evaluation of post-quench magnetisation evolution. Although this perturbative method is inherently time-limited and converges more slowly in the Potts case than in the Ising model, it provides a remarkably accurate description of short- to medium-term dynamics compared with numerical iTEBD simulations.

To conclude, this work provides a robust analytical foundation for studying non-equilibrium dynamics in non-integrable spin chains. An essential next step for this research is to extend the perturbative framework to include three-kink bound states (baryons). As baryons are a distinguishing feature of the Potts model over the simpler Ising chain, capturing their contribution to the post-quench spectrum remains an intriguing challenge for future investigations.

Acknowledgments

This work was partially supported by the Ministry of Culture and Innovation and the National Research, Development and Innovation Office (NKFIH) through the OTKA Grant K 138606. AK was also partially supported by the HUN-REN Hungarian Research Network through the HUN-REN-BME-BCE Quantum Technology Research Group, and by the Doctoral Excellence Fellowship Programme (DCEP), funded by the National Research Development and Innovation Fund of the Ministry of Culture and Innovation and the Budapest University of Technology and Economics, under a grant agreement with the National Research, Development and Innovation Office. GT was also supported by the Quantum Information National Laboratory of Hungary (Grant No. 2022-2.1.1-NL-2022-00004). We also acknowledge KIFÜ (Governmental Agency for IT Development, Hungary) for awarding us access to the Komondor HPC facility based in Hungary.

Appendix A Details of the perturbative description of the excitation spectrum

A.1 The pure transverse chain: kink excitations

Here we discuss the details of the perturbative analysis of the model described by the Hamiltonian $H(\mathbf{v},g)$ with $\mathbf{v}=\mathbf{0}$, which corresponds to the infinite three-state Potts spin chain with purely transverse magnetic field. The perturbative calculations describe the extreme ferromagnetic limit $g\ll 1$ in the two-kink sector.

As stated in the main text in Section 2.4.1, the first step in solving the eigenvalue problem (2.32) perturbatively is the diagonalisation of the Hamiltonian $H^{(2)}(\mathbf{0},g)$, which is given by the ’Bethe states’ (2.26):

[TABLE]

where $p_{1},p_{2}\in\mathbb{R}/2\pi\mathbb{Z}$ are the quasimomenta of the two kinks, $\nu=2,3$ and $S_{\nu,\nu^{\prime}}(p_{1},p_{2})$ are the entries of the two-kink scattering matrix:

[TABLE]

The Bethe states $|{K}_{1,\nu}(p_{1}){K}_{\nu,1}(p_{2})\rangle$ satisfy the Faddeev-Zamolodchikov commutation relation:

[TABLE]

Such states with $-\pi<p_{2}<p_{1}\leq\pi$ form the basis in the subspace $\mathcal{L}_{11}^{(2)}$, and satisfy the normalisation condition:

[TABLE]

for $-\pi<p_{2}<p_{1}\leq\pi$, $-\pi<p_{2}^{\prime}<p_{1}^{\prime}\leq\pi$, and $\nu,\nu^{\prime}=2,3$. We also use the basis (2.27) of the two-kink Bethe states:

[TABLE]

with $\iota=\pm$. From the second line of the above equation, it follows that these states diagonalise the two-kink scattering matrix:

[TABLE]

Most importantly, the two-kink Bethe states defined by Eqs. (A.1) and (A.3), diagonalise the operators $H^{(2)}(\mathbf{0},g)$ and $\hat{T}_{1}$ simultanoeusly:

[TABLE]

where

[TABLE]

is the kink dispersion law to the linear order in the small parameter $g$.

The scattering matrix defined by (A.2) satisfies the unitarity relations:

[TABLE]

where $p_{1},p_{2}\in\mathbb{R}$, $\nu,\nu^{\prime}=2,3$, and $\iota=\pm$.

A.2 Two-kink states in the aligned case

Here we consider the perturbative analysis of the low-energy spectrum of the Hamiltonian $H(\mathbf{v},g)$ with $\mathbf{v}=(v_{1},0,0)$ in the subspace $\mathcal{L}_{11}^{(2)}$. For the Hamiltonian, we use the short-cut notation $H^{(2)}(v_{1},g)$ and calculate its spectrum to the leading order in $g\ll 1$. We use the basis in the two-kink subspace $\mathcal{L}_{11}^{(2)}$:

[TABLE]

where $\iota=\pm$, $j=1,2\ldots$, and $-\pi<P\leq\pi$ is the quasimomentum. It follows from (2.24), (2.25), that

[TABLE]

Denote by $\mathcal{L}^{(2)}(P)=\mathcal{L}^{(2)}(P,+)\oplus\mathcal{L}^{(2)}(P,-)$ the subspace of $\mathcal{L}_{11}^{(2)}$ spanned by the vectors $|j,P\rangle_{\iota}$, with $j=1,2\ldots$, and $\iota=\pm$. For all $|\psi\rangle\in\mathcal{L}^{(2)}(P)$ we have

[TABLE]

It is convenient to define a modified scalar product $\langle\ldots|\ldots\rangle$ in the subspace $\mathcal{L}^{(2)}(P)$ which satisfies:

[TABLE]

with $\iota,\iota^{\prime}=\pm$, and $j,j^{\prime}=1,2,\ldots$.

Consider now the eigenvalue problem

[TABLE]

for $|\pi_{\iota}\rangle\in\mathcal{L}^{(2)}(P,\iota)$, and expand its eigenvector $|\pi_{\iota}\rangle$ in the basis $|j,P\rangle_{\iota}$:

[TABLE]

We simplify notations by omitting the labels $E_{\iota}$ and $P$ in $\psi_{\iota}(j;E_{\iota},P)$, and use $\psi_{\iota}(j)$. Eq. (A.16) implies the following second-order linear difference equation:

[TABLE]

where $j=2,3,\ldots$. Furthermore, the wave-function coefficients $\psi_{\iota}(j)$ must vanish for $j\to+\infty$, and also satisfy the boundary condition:

[TABLE]

Note, that the difference equation (A.18) and the boundary condition (A.19) are invariant under the transformation

[TABLE]

Using the following identity satisfied by the Bessel function of the first kind $J_{c}(Z)$:

[TABLE]

the solution for the wave-function reads

[TABLE]

where

[TABLE]

Substituting of (A.22) into (A.19) leads to the secular equation

[TABLE]

where $Z$ is given by (A.23). For $h_{1}>0$, solutions of this equation determine the dispersion laws $E_{\iota,n}(P,h_{1})$ of mesons, while for $h_{1}<0$ they yield the spectrum of bubbles, to the first order in the small parameter $g$. The resulting dispersion laws $E_{\iota,n}(P)$ of the few lightest mesons are depicted in Fig. 2.2 for $g=0.2$ and $h_{1}=0.05$.

The energies of mesons and bubbles furthermore satisfy the following equality:

[TABLE]

which follows directly from the symmetry property (A.20), valid to the first order in the parameter $g\ll 1$.

A.3 Two-kink states in the oblique case

We now turn to the calculation of the energy spectrum of the Hamiltonian $H(\mathbf{v},g)$ in the subspace $\mathcal{L}_{11}^{(2)}$ in the oblique regime $\mathbf{v}=(0,v_{2},0)$, to leading order in $g\ll 1$. Once again, we introduce the shorthand notation $H^{(2)}(v_{2},g)$ for the corresponding Hamiltonian.

We choose the basis

[TABLE]

in the subspace $\mathcal{L}_{11}^{(2)}$, with $\nu=2,3$, $j=1,2,\ldots$, and $-\pi<P\leq\pi$. These basis vectors are simply related to the vectors $|j,P\rangle_{\pm}$ given in (A.12):

[TABLE]

At a fixed $P$, vectors $|j,P\rangle_{\nu}$, with $\nu=2,3$, and $j=1,2\ldots$, form the basis of the subspace $\mathcal{L}^{(2)}(P)$ introduced in the previous Section, and satisfy the normalisation condition:

[TABLE]

To solve the eigenvalue problem

[TABLE]

we expand the eigenvector $|\Psi\rangle$ in the basis $|j,P\rangle_{\nu}$:

[TABLE]

For the sake of simplicity, we again suppress the variables $E$ and $P$ in the wave-functions $\psi_{\nu}(j;E,P)$ and use $\psi_{\nu}(j)$ instead. The eigenvalue problem (A.29) implies the following system of second-order linear difference equations:

[TABLE]

with $j=2,3,\ldots$. In addition, the wave function coefficients satisfy the boundary conditions

[TABLE]

Furthermore, the coefficients $\psi_{2}(j)$ must vanish for $j\to+\infty$, while the coefficients $\psi_{3}(j)$ must be at least bounded in absolute value.

Equations (A.31) and boundary conditions (A.32) do not change upon the transformation:

[TABLE]

The solution of equation (A.31a) reads

[TABLE]

where $j=1,2\ldots$,

[TABLE]

while the solution of equation (A.31b) can be written as

[TABLE]

where

[TABLE]

The discussion of the resulting spectra is included in the main text. Here, we only note that due to the symmetry property (A.33), the resulting meson and bubble energies in the oblique regimes satisfy the equality

[TABLE]

Appendix B True/false vacuum states for $g\ll 1$

Here we apply the standard Rayleigh-Schrödinger perturbation theory to describe the properties of the vacuum states of the Potts spin chain (2.11) for $g\ll 1$. A similar analysis for the antiferromagnetic XXZ spin chain in the anisotropic limit was described by Jimbo and Miwa on Pages 14-18 of their monograph [41].

Consider the Hamiltonian $H_{N}(\mathbf{v},g)$ given by (4.1) for some even number of sites $N$ with periodic boundary conditions. The deformation $|\text{vac}(\mathbf{v},g|N)\rangle$ of the vacuum $|0\rangle^{(1)}$ satisfies the following properties:

[TABLE]

The deformed vacuum $|\text{vac}(\mathbf{v},g|N)\rangle$ can be represented as

[TABLE]

where the vector $|\Omega(\mathbf{v},g|N)\rangle$ admits the following asymptotic expansion:

[TABLE]

such that

[TABLE]

The coefficients $|\Omega_{n}(\mathbf{v}|N)\rangle$ in this expansion can be perturbatively determined term by term. In particular:

[TABLE]

where we used the notation (2.28). The norm of the vector (B.3) can be expanded as

[TABLE]

Denote by ${\mathfrak{N}}$ the linear operator acting on the space

[TABLE]

in the following way:

[TABLE]

Then we can compute

[TABLE]

This equality shows that in the state $|\text{vac}(\mathbf{v},g|N)\rangle$, the spin chain of length $N$ contains about $\frac{Ng^{2}}{18}$ sites where the spin is oriented not in the first, but in the second or third direction. Therefore, such flipped spins are well-separated from one another for $g\ll 1$. E.g., at $g=0.2$, the distance between the neighbouring flipped spins is about $18/0.2^{2}=450$ spin-chain sites.

Similar reasoning leads to the relations

[TABLE]

for $j=-N/2+1,\ldots,N/2$.

Appendix C Perturbative expansion for Potts quench dynamics

Using the results of Appendix B, here we determine the post-quench time evolution of the expectation value $A_{\alpha}(t|N)$ of the operator $\hat{A}_{\alpha}$, defined by (4.6,4.7), in the Potts spin chain. Note that

[TABLE]

for $\alpha=1,2,3$. Rewriting equation (B.1c) as

[TABLE]

expanding $|\delta\Phi(\mathbf{v},g|N)\rangle$ to the second order in $g$ yields

[TABLE]

After substitution of the right-hand side of the above equality into (4.6), straightforward manipulations yield

[TABLE]

Taking the thermodynamic limit $N\to\infty$ results in

[TABLE]

where following the main text we introduced the notations

[TABLE]

The first term on the right-hand side of (C.5) does not depend on time, and gives the expectation value of the operator $\hat{A}_{\alpha}$ in the deformed vacuum state $|\text{vac}(\mathbf{v},g)\rangle$. It follows from (B.10), (4.7) that

[TABLE]

The vector $|\Omega_{1}(\mathbf{v})\rangle$ in the second and third terms on the right-hand side of (C.5) is given by (C.6c), and corresponds to a translation invariant two-kink state:

[TABLE]

Operators $\hat{A}_{\alpha}$ act as homomorphisms in the subspace $\mathcal{L}^{(2)}(0)$, whose basis is given by vectors $\{|j,0\rangle_{2,3}\}_{j=1}^{\infty}$ determined by (2.28). The matrix elements of the operators $\hat{A}_{\alpha}$ in this basis are

[TABLE]

where no summation over repeated indices is implied on the right-hand sides.

Note that the (modified) scalar product in the subspace $\mathcal{L}^{(2)}(0)$ is defined according to (A.28):

[TABLE]

with $j=1,2,\ldots$, and $\nu=2,3$. Accordingly, the completeness relation for the subspace $\mathcal{L}^{(2)}(0)$ reads

[TABLE]

where $\mathbf{1}$ is the identity operator acting on $\mathcal{L}^{(2)}(0)$.

In general, the Hamiltonian $\Delta H(\mathbf{v},g)$ determined by (C.6d) does not conserve the number of kinks. Therefore at $t>0$ the state $|\Psi(t)\rangle=e^{-it\,\Delta H(\mathbf{v},g)}|\Omega_{1}(\mathbf{v})\rangle$ gains contributions of the $n$-kink states with $n=3,4,\ldots$:

[TABLE]

where $|\Psi^{(n)}(t)\rangle\in\mathcal{L}_{11}^{(n)}$. This prevents the exact analytical calculation of the second and third terms on the right-hand side of (C.5) at $t>0$. The key simplification is provided by the two-kink approximation, which corresponds to replacing the Hamiltonian $\Delta H(\mathbf{v},g)$ by the restriction of the Hamiltonian $H(\mathbf{v},g)$ defined by (2.11) to the subspace $\mathcal{L}_{11}^{(2)}$:

[TABLE]

Note, that the Hamiltonians $H(\mathbf{v},g)$ and $\Delta H(\mathbf{v},g)$ differ only by the constant term $\sim g^{2}$, which arises from the vacuum energy (B.9). Substitution of (C.14) into (C.4), (C.6a) leads to the result

[TABLE]

where the scalar product is understood in accordance with (C.11). In the above expression, only the second and third terms exhibit time dependence, while the first (time-independent) term is determined by (B.2,B.3) and (C.7,C.8) as

[TABLE]

For the remainder of this section, we focus on the oblique regime $\mathbf{v}=(0,h_{2}/g,0)$, and use the short-hand notation $H^{(2)}(h_{2},g)$ for $H^{(2)}(\mathbf{v},g)$. We start by considering the matrix element

[TABLE]

that stands on the right-hand side of (C.15).

The identity operator (C.12) acting on the subspace $\mathcal{L}^{(2)}(0)$ admits the alternative representation:

[TABLE]

Here $|\psi_{n}\rangle$ and $|\psi(p)\rangle$ are the eigenstates of the Hamiltonian $H^{(2)}(h_{2},g)$ corresponding to the discrete and continuous energy spectra, respectively

[TABLE]

where $|\psi_{n}\rangle,|\psi(p)\rangle\in\mathcal{L}^{(2)}(0)$, and $E(p)=2-\frac{4g}{3}\cos p$. These are essentially the same states as given by Eq. (A.30) but reduced to the subspace $P=0$. Thus, with the notation $\psi_{\nu}(j;E,P=0)$ introduced in the main text for the wave-function, these states can be expressed as

[TABLE]

for the continuous spectra and similarly as

[TABLE]

for the discrete case. The Hamiltonian eigenstates corresponding to the continuous spectrum are normalised by the condition:

[TABLE]

where $p,p^{\prime}\in(0,\pi)$, and $\mathrm{B}(p)=B_{in}(e^{ip},h_{2}/g)$, with the function $B_{in}(z,v)$ defined by (2.48).

Combining (C.17) with (C.18), one finds:

[TABLE]

All scalar products and matrix elements on the right-hand side are real. Therefore:

[TABLE]

Applying the expansions (C.21,C.22), the expression of $|\Omega_{1}(h_{2})\rangle$ given by (B.4) and the relations (C.10), the relevant terms can be determined as

[TABLE]

with similar results for the continuous part of the spectrum. The Fourier transformation of the second term on the right-hand side of (C.15) reads:

[TABLE]

For the real part of $F_{\alpha}^{(1)}(\omega)$, this yields:

[TABLE]

For $\omega\in(2-\frac{4g}{3},2+\frac{4g}{3})$, we get, in particular:

[TABLE]

where

[TABLE]

For $\alpha=1,2$, the explicit form of formula (C.29) reads:

[TABLE]

Here we recall that

[TABLE]

where the function $B_{in}(z,v_{2})$ is given by (2.48).

Let us turn now to the matrix element

[TABLE]

which appears in the second line on the right-hand side of (C.15). Exploiting the completeness identity in the form (C.12), together with formulas (C.10), the right-hand side of (C.34) can be evaluated as:

[TABLE]

where

[TABLE]

Using the completeness identity in the form (C.18), together with formula (C.6c), the latter function can be represented in the explicit form:

[TABLE]

where

[TABLE]

and of course, once again $E(p)=2-\frac{4g}{3}\cos p$.

The post-quench time evolution of the magnetisations can be then calculated by the means of equations (4.7) and (C.15,C.16,C.25,C.26,C.36,C.38). The resulting time evolution of the magnetisations $M_{2}(t)$, and $M_{3}(t)$ are displayed in Fig. 4.1 at $g=0.2$ in the positive (a), and negative (b) oblique regimes with $h_{2}=0.1$, and $h_{2}=-0.1$, respectively.

Our perturbative approach predicts

[TABLE]

that is also reflected in Fig. 4.1 for the first stage of the time evolution. The explicit form of $c$ can be found by means of the asymptotical analysis of equations (C.36), (C.38) at $t\to\infty$. The final result reads:

[TABLE]

Appendix D Quantum quench in the Ising spin chain at a

weak transverse magnetic field

Here, we test our perturbative approach using quantum quenches in the transverse-field Ising chain (TFIC), defined by the Hamiltonian

[TABLE]

where $N$ is even, and the periodic boundary condition $\sigma_{N/2+1}^{z}=\sigma_{-N/2+1}^{z}$ is imposed. This Hamiltonian acts in the $2^{N}$-dimensional Hilbert space

[TABLE]

We denote the basis of the two-dimensional local vector space $[\mathbb{C}^{2}]_{j}$ associated with the site $j$ as $|\mu\rangle_{j}$, with $\mu=1,2$.

As mentioned in the Introduction, this model is integrable for any value of the transverse magnetic field $h_{x}$ and admits a representation in terms of non-interacting fermions via the Jordan-Wigner transformation. The free-fermionic representation allows one to describe the dynamics in the TFIC in a very simple way [55, 56, 57], making it a paradigmatic model for quantum quenches.

The small parameter of our perturbation theory is $h_{x}\ll 1$. We consider the quench protocol where the initial state is the pure ferromagnetic state along direction 1:

[TABLE]

This state represents one of two degenerate ground states of the Hamiltonian (D.1) at $h_{x}=0$:

[TABLE]

Let us tune the transverse magnetic field $h_{x}$ to some small positive value, and consider the eigenvalue problem:

[TABLE]

in which $E_{\text{vac}}(h_{x}|N)$ and $|\Omega(h_{x}|N)\rangle$ are understood as formal power series in $h_{x}$. In particular:

[TABLE]

We require $|\Omega_{0}(N)\rangle=|0\rangle^{(1)}$, and $\phantom{.}{}^{(1)}\langle 0|\Omega_{n}(N)\rangle=0$, for all $n=1,2,\ldots$. The normalised deformed vacuum state associated with (D.6) reads:

[TABLE]

Let us consider the operator $\hat{A}$, which differs by the constant shift $\frac{1}{2}$ from the magnetisation operator $\hat{M}$:

[TABLE]

The time evolution of its expectation value is given by

[TABLE]

Following the procedure described in Section 4, one obtains in the thermodynamic limit $N\to\infty$ the small-$h_{x}$ asymptotic formula for the magnetisation $M(t)$:

[TABLE]

which is analogous to (C.5), where we introduced the notations

[TABLE]

The next crucial step of the calculations is provided by the two-kink approximation, which replaces the Hamiltonian $\Delta H(h_{x})$ in the asymptotic formula (D.11) by its restriction $\Delta H^{(2)}(h_{x})$ to the subspace $\mathcal{L}_{11}^{(2)}$ spanned by the translation invariant two-kink states:

[TABLE]

which yields

[TABLE]

where

[TABLE]

and $\mathcal{P}_{11}^{(2)}$ is the projection operator onto the subspace $\mathcal{L}_{11}^{(2)}$. Taking into account the evident equality

[TABLE]

the result (D.14) can be rewritten as

[TABLE]

where

[TABLE]

The function $A_{1}^{(2)}(t)$ determined by (D.18) describes high-frequency oscillations of the magnetisation $M^{(2)}(t)$. The explicit representation for this function reads:

[TABLE]

where $J_{j}(\tau)$ is the Bessel function, and $\tau=2h_{x}t$.

In turn, the function $A_{2}^{(2)}(t)$ determined by (D.19) contains the slowly varying contribution to the magnetisation, and its explicit form reads:

[TABLE]

where

[TABLE]

The asymptotic behaviour of $A_{2}^{(2)}(t)|_{t=\tau/(2h_{x})}$ at large $\tau$ is given by:

[TABLE]

This result can be matched against the analytic results of Ref. [56], where the large-time asymptotic of the magnetisation is obtained as†††Note that our convention for the Hamiltonian and the magnetisation both differ by a factor of $1/2$, which accounts for the difference of Eq. (D.24) from Eq. (15) of Ref. [56].

[TABLE]

Expanding this expression to order $h_{x}^{3}$ gives

[TABLE]

which correctly reproduces the linear term in (D.23). The full exponential relaxation of the magnetisation, as manifest in (D.24), can be obtained by a form-factor resummation [56]; the perturbative method we use only gives its lowest-order expansion, which is linear in time. However, as noted in [56], the form factor resummation cannot reproduce the oscillating contributions which are clearly absent from (D.24, D.25). These oscillating terms were computed [51] and numerically verified [25] in the scaling field theory limit, which cannot be compared to the perturbative result valid for $h_{x}\ll 1$, since the scaling limit is obtained in the vicinity of the critical point $h_{x}=1$.

Fig. D.1 illustrates the post-quench time evolution of the magnetisation determined by equations (D.17), (D.20), and (D.21) at $h_{x}=0.2$, compared to the results of the numerical simulation of the time evolution. The simulation was performed using the infinite volume Time-Evolving Block Decimation (iTEBD) method with second-order Trotterisation, time step $\delta t=0.005$ and maximal bond dimension $\chi_{\text{max}}=300$. The Fourier images $F_{\alpha}(\omega)$ of the functions $A_{\alpha}^{(2)}(t)$ are defined according to:

[TABLE]

Here we present only the remarkably simple formula for the real part of the function $F_{1}(\omega)$:

[TABLE]

which holds at $\omega>0$.

Appendix E Details of the numerical simulations

The numerical simulations of the time evolution after quantum quenches were performed using the infinite volume Time-Evolving Block Decimation (iTEBD) method. In the simulations, we used second-order Trotterisation with time step $\delta t=0.005$, evaluating the magnetisation(s) after every 20th step. The maximal bond dimension was set at $\chi_{\text{max}}=300$ for simulations running to longer times ($t=1000$), and at $\chi_{\text{max}}=800$ for simulations running up to shorter times (approximately $t\approx 200$). Larger bond dimensions ensure higher precision; however, they consume much more computer resources and run much more slowly. Nevertheless, to achieve sufficient resolution for quench spectroscopy, long-time runs are required. The resulting Fourier spectra of the long-time simulations with $\chi_{\text{max}}=300$ and the short-time ones with $\chi_{\text{max}}=800$ essentially agree, as shown in Ref. [23].

After the quenches, we observe persistent oscillations in the magnetisation(s), as illustrated in the main text in Fig. 4.2. In the following, we supplement these results by showing the entanglement entropy $S$ for two qualitatively distinct cases in Fig. E.1 for both settings of the maximum bond dimension.

In Fig. E.1 (a), it is demonstrated that for the positively aligned quench, the entropy exhibits an apparent saturation up to short times but manifests a slow drift in the long-time simulations. In the negatively aligned case, the time evolution of the entropy shows qualitatively the same features. This behaviour results from the confining/anticonfining nature of the quenches. However, in the oblique cases, demonstrated by the positive oblique case in Fig. E.1 (b), the entropy grows fast due to the contribution from the unconfined kink degrees of freedom [23].

To compare the numerical and the analytical Fourier spectra, it is necessary to take into account the different definitions of the Fourier transform (C.27) and (4.21). As stated in the main text, these can differ by a multiplicative factor and a constant shift. By comparing the definitions, the multiplicative factor $a$ turns out to be given by

[TABLE]

In the long-time simulations, we time-evolved the initial state until $t=1000$ with $\delta t=0.005$, as already stated, yielding $N_{s}=200001$ steps, including the $t=0$ point. However, as we noted previously, we only calculated the magnetisation(s) at every 20th time step, which means that the number of steps relevant for the Fourier transformation is $N_{FT}=10001$, while the relevant time step is $\delta t_{\text{FT}}=0.1$, which can be substituted into Eq. (E.1) yielding $a\approx 0.10$.

To supplement the main text, we also show typical Fourier spectra of the analytically calculated $a\,\text{Re}\,F^{(1)}_{2}(\omega)+\text{d}M_{2}$ together with $\text{Re}\,M_{2}(\omega)$ determined in the iTEBD simulations for oblique quenches in Figs. E.2 and E.3. The numerical and analytical data show excellent agreement, with the shift $\text{d}M_{2}$ vanishing.