# Exploring Two-dimensional Coherent Spectroscopy with Exact Diagonalization:
Spinons and Confinement in 1D Quantum Magnets

Yoshito Watanabe  Institute for Theoretical Physics, University of Cologne,
50937 Cologne, Germany  Simon Trebst  Institute for Theoretical Physics,
University of Cologne, 50937 Cologne, Germany  Ciarán Hickey  School of
Physics, University College Dublin, Belfield, Dublin 4, Ireland  Centre for
Quantum Engineering, Science, and Technology, University College Dublin,
Dublin 4, Ireland

###### Abstract

Two-dimensional coherent spectroscopy (2DCS) with terahertz radiation offers a
promising new avenue for the exploration of many-body phenomena in quantum
magnets. This includes the potential diagnosis of fractionalized excitations,
for which linear response often struggles due to the indistinguishability of a
continuum of fractional excitations from that caused by disorders or
impurities. However, the interpretation of the complex results produced by
2DCS remains a challenge, and a general prediction of the spectral
characteristics of different types of excitations has not yet been
established. In this paper, we develop a numerical approach based on exact
diagonalization (ED) to push our understanding of 2DCS towards different
scenarios. We first validate our approach by comparing numerical ED and exact
analytical results for the spectroscopic signatures of spinons in one-
dimensional transverse field Ising model and develop how to deal with the
inherently small system sizes in ED calculations. Augmenting the model by a
longitudinal field, we demonstrate significant changes to the 2DCS spectrum
upon the field-induced spinon pair confinement, which can be rationalized in
our ED calculations and from a “two-kink” model (in the absence of
integrability). One advantage of our ED approach is its possible extension to
finite temperatures, which we explore using thermally pure quantum states and
demonstrate to change the intensity and spectroscopical patterns of 2DCS when
going beyond the integrable model. Our numerically exact results provide a
benchmark for future experiments and theoretical studies relying on
approximation methods, and pave the way for the exploration of fractionalized
excitations in quantum magnets.

##  I Introduction

Terahertz (THz) spectroscopy, in which a THz pulse incident on a system is
used to excite and probe its dynamics, is widely employed in the study of
quantum materials [1, 2, 3, 4]. In particular, the THz range aligns well with
the typical excitation energies of quantum magnets. The technique can be
understood within the framework of linear response, providing information akin
to the dynamical structure factor, but restricted to zero momentum. In
conventional magnets, this form of linear response typically reveals sharp
excitation modes, such as magnons or triplons. However, in systems with
fractionalized excitations, linear response reveals a broad continuum of
excitations, reflecting the underlying fact that local operators create
multiple (deconfined) excitations [5, 6, 7]. Consequently, such systems
present a profound challenge in distinguishing between a continuum arising
from fractionalized excitations and one resulting from, for example, thermal
or impurity-induced disorder [8, 9, 10, 11, 12]. Furthermore, even in the
absence of disorder, extracting precise information from the continuum is
notoriously difficult. This complexity highlights the need for complementary
methods or approaches to accurately characterize dynamics in quantum
materials, particularly in distinguishing the nuances of its excitation
spectra.

A recently developed method to gain deeper insights into the dynamics of
unconventional magnets is to expand THz spectroscopy to include _multiple_
excitation pulses [13, 14, 15]. If one employs two THz pulses to probe the
system and, in addition, varies the timing between the two pulses, one can
perform two-dimensional coherent spectroscopy (2DCS) – a technique that, in
the optical range, has already been extensively used in fields such as quantum
chemistry, to probe the structure of complex molecules [16], and semiconductor
physics, to probe the dynamics of excitons [17]. At its heart, such a multi-
pulse approach allows one to extract the _non-linear_ response of the system
by subtracting the single pulse responses from the multi-pulse response. Then,
by performing a two-dimensional Fourier transform over the time arguments
$\tau$, the time delay between the two pulses, and $t$, the measurement time
after the second pulse, the frequency spectrum of these non-linear responses
can be revealed. With THz sources, it has been argued that 2DCS has the
potential to become one of the more potent tools to explore many-body
phenomena in quantum magnets, such as in studying spin waves [18] or, more
enticingly, fractionalized quasiparticles such as spinons [19, 20, 21] or
Majorana fermions [22, 23, 24, 25, 26].

2DCS, being a probe of higher-order dynamical correlation functions,
inherently contains more information than in linear response. One of the key
capabilities of non-linear response is the detection of
interactions/statistics between quasiparticles within the same species [27,
18], as well as between different modes [23, 14, 15]. Furthermore, recent
theoretical work [19], has shown that 2DCS can clearly identify one of the
most elementary forms of fractionalized excitations – spinon excitations in
the (exactly solvable) one-dimensional transverse field Ising model (1D-TFIM).
The “spinon echo” signal, represented by a sharp anti-diagonal line along
$\omega_{t}=-\omega_{\tau}$, arises from the interference of phases
accumulated by spinon pairs created by consecutive THz pulses. Weak
perturbations have been argued to introduce a finite lifetime to the spinon
pairs, which in turn manifests as an energy-dependent broadening of the spinon
echo signal, thereby revealing the lifetime of the individual spinon
excitations [28]. The appearance of a similar (sharp) anti-diagonal feature as
a signature of fractionalized excitations has also been shown in the Kitaev
honeycomb model [22], distinguishing it from broad continua due to more
trivial effects.

Figure 1: Schematic of the two-pulse measurement setup. Two THz pulses (in
red) are applied at times $t^{\prime}=0$ and $t^{\prime}=\tau$, and the
magnetization $m_{0\tau}$ (in blue) is measured at $t^{\prime}=t+\tau$.
Subtracting the single pulse responses, $m_{0}$ and $m_{\tau}$, from
$m_{0\tau}$ isolates the non-linear magnetization $m_{\text{NL}}$.

Despite the usefulness of 2DCS, theoretically, the calculation of the non-
linear response functions and the interpretation of the 2DCS spectrum can be
quite challenging. Here, we explore the utility of using exact diagonalization
(ED) to study 2DCS in 1D quantum magnets. ED is well-suited to achieving long-
time evolution, and correspondingly high-frequency resolution, but it is
obviously rather limited in terms of system sizes. As such, it offers a
different numerical trade-off than recent exploratory studies [29, 30] using a
time-evolution of matrix product states (MPS), which can deal with much larger
system sizes than ED but can be limited to shorter time scales due to the
growth of entanglement under time evolution [31]. To showcase the
applicability of our ED approach, we consider the 1D-TFIM and critically
discuss finite-size effects in its exactly solvable limit and then extend our
study to the case where the integrability is broken by a longitudinal field.
We demonstrate that the fine resolution in the frequency domain achieved by ED
reveals a detailed structure of the 2DCS spectrum characterized by a sequence
of spinon-pair bound states. Finally, we will extend our approach to finite
temperatures, a crucial developement for the comparison with experiments.

Our discussion in the following is structured as follows: We start, in Sec.
II, with a general overview of 2DCS, details on how we numerically compute the
associated non-linear response functions with ED, and the model Hamiltonian
that we focus on in this work. Next, in Sec. III, we investigate the origin of
finite-size effects in the 1D-TFIM and discuss how to interpret the resulting
2DCS spectra. With the understanding of the finite-size effects in hand, we
then turn to the case of a longitudinal field in Sec. IV. Combining our ED
results with a “two-kink” approximation, we show the confinment of spinon
pairs results in a significant change of the “spinon-echo” signal. Finally, in
Sec. V, we discuss the ED approach to finite temperatures and demonstrate how
the finite-temperature effects can be captured in 2DCS. In Sec. VI, we
summarize our results and discuss the experimental relevance of our findings.

##  II Methods

Before discussing any results, we first provide a brief overview of magnetic
2DCS and its theoretical underpinnings, details on how we numerically compute
the associated non-linear response functions with ED, and the model
Hamiltonian that we focus on in this work.

#### Formalism of 2DCS

We consider a simple time-dependent magnetic field consisting of two delta-
function pulses at times $t^{\prime}=0$ and $t^{\prime}=\tau$. The
magnetization is measured after the second pulse at a time
$t^{\prime}=t+\tau$. Thus, $\tau$ is the delay time between the two pulses and
$t$ is the measurement time after the second pulse. The $\alpha$-component of
the resulting time-dependent applied field can be written as
$B^{\alpha}(t^{\prime})=B^{\alpha}_{0}\delta(t^{\prime})+B^{\alpha}_{\tau}\delta(t^{\prime}-\tau)$.

We define $m^{\alpha}_{0\tau}(t,\tau)$ as the $\alpha$-component of the
induced magnetization at time $t^{\prime}=t+\tau$ after two successive pulses
at $t^{\prime}=0$ with strength $B_{0}$ and at $t^{\prime}=\tau$ with strength
$B_{\tau}$. Analogously, $m^{\alpha}_{0}(t,\tau)$ and
$m^{\alpha}_{\tau}(t,\tau)$ represent the induced magnetization after single
pulses at $t^{\prime}=0$ with strength $B_{0}$ and at $t^{\prime}=\tau$ with
strength $B_{\tau}$ respectively. Subtracting these single-pulse magnetization
responses ($m^{\alpha}_{0}$ and $m^{\alpha}_{\tau}$) from the two-pulse
magnetization $m^{\alpha}_{0\tau}$ isolates the non-linear magnetization as

$\begin{split}m_{\text{NL}}^{\alpha}(t,\tau)=&\,\,m^{\alpha}_{0\tau}(t,\tau)-[m^{\alpha}_{0}(t,\tau)+m^{\alpha}_{\tau}(t,\tau)]\\\
=&\,\,\chi_{\alpha\beta\gamma}^{(2)}(t,\tau)B^{\beta}_{\tau}B_{0}^{\gamma}+\\\
&\,\,\chi_{\alpha\beta\gamma\gamma}^{(3)}(t,\tau,0)B_{\tau}^{\beta}B_{0}^{\gamma}B_{0}^{\gamma}+\\\
&\,\,\chi_{\alpha\beta\beta\gamma}^{(3)}(t,0,\tau)B_{\tau}^{\beta}B_{\tau}^{\beta}B_{0}^{\gamma}+\dots,\end{split}$
(1)

wherein only mixed cross-terms survive (for example, second-order terms
proportional to $B_{0}^{2}$ and $B_{\tau}^{2}$ cancel out). This setup is
schematically illustrated in Fig. 1. Note that, in this two-pulse setup, there
is a single second-order susceptibility but two distinct third-order
contributions proportional to $B_{\tau}B_{0}^{2}$ and $B_{\tau}^{2}B_{0}$
respectively.

Concrete expressions for the non-linear susceptibilities can be derived using
time-dependent perturbation theory [32]. In general, they are written as the
equilibrium expectation value of nested commutators of the magnetization
operators. As an example, the relevant third-order susceptibility can be
neatly expressed in its most general form as

$\chi^{(3)}_{\alpha\beta\gamma\delta}(t_{3},t_{2},t_{1})=\frac{2}{N}\textrm{Im}[R_{1}+R_{2}+R_{3}+R_{4}],$
(2)

where $t_{1}$, $t_{2}$ and $t_{3}$ are the time delays for a general multi-
pulse setup, and the $R_{a}$ encode the contributions from the different
possible orderings of the operators (we have suppressed their arguments and
labels here for simplicity). As an example, $R_{1}$ is given by

$R_{1}=\left<M^{\gamma}(t_{1})M^{\beta}(t_{2}+t_{1})M^{\alpha}(t_{3}+t_{2}+t_{1})M^{\delta}(0)\right>,$
(3)

where $M^{\alpha}(t)=\sum_{i}S_{i}^{\alpha}(t)$. By inserting resolutions of
the identity, $R_{1}$ can be rewritten in terms of the energy eiegenstates
$\ket{P}$, $\ket{Q}$ and $\ket{R}$ as

$\begin{split}R_{1}=\sum_{PQR}&m^{\gamma}_{0R}m^{\beta}_{RQ}m^{\alpha}_{QP}m^{\delta}_{P0}\\\
&\times
e^{-\text{i}E_{P}t_{1}}e^{-\text{i}(E_{P}-E_{R})t_{2}}e^{-\text{i}(E_{P}-E_{Q})t_{3}},\end{split}$
(4)

where we have defined the magnetization matrix elements
$m^{\alpha}_{fi}=\bra{f}M^{\alpha}\ket{i}$. See Appendix A for the full
expressions for $\chi^{(2)}$ and $\chi^{(3)}$ and their associated $R_{a}$.

Setting $t_{1}=0,t_{2}=\tau,t_{3}=t$ and $t_{1}=\tau,t_{2}=0,t_{3}=t$ yields
$\chi_{\alpha\beta\gamma\gamma}^{(3)}(t,\tau,0)$ and
$\chi_{\alpha\beta\beta\gamma}^{(3)}(t,0,\tau)$ respectively. The two-
dimensional Fourier transform of $\chi^{(3)}$ over positive $t$ and $\tau$
generates the 2DCS frequency spectrum. We define
$\chi^{(3;1)}(\omega_{t},\omega_{\tau})\equiv\text{FT}[\theta(t)\theta(\tau)\chi^{(3)}(t,\tau,0)]$
and
$\chi^{(3;2)}(\omega_{t},\omega_{\tau})\equiv\text{FT}[\theta(t)\theta(\tau)\chi^{(3)}(t,0,\tau)]$.
It’s important to note here that the positive time constraint, $t>0,\tau>0$,
imposed by the form $\theta(t)\theta({\tau})\chi(t,\tau)$, results in a
convolution of
$[\delta(\omega_{t})+1/(\text{i}\pi\omega_{t})][\delta(\omega_{\tau})+1/(\text{i}\pi\omega_{\tau})]$
with $\tilde{\chi}(\omega_{t},\omega_{\tau})$, where
$\tilde{\chi}(\omega_{t},\omega_{\tau})$ is the unconstrained Fourier
transform of $\chi(t,\tau)$. Unfortunately, this convolution means that both
the real and imaginary parts of the Fourier transform contain artificial
broadening, e.g. an additional $1/(\omega_{t}\omega_{\tau})$ term for the real
part. This distortion of the spectrum of the pure
$\tilde{\chi}(\omega_{t},\omega_{\tau})$ is commonly referred to as “phase
twisting” [33, 28], and is a known impediment to the clean interpretation of
2DCS spectra [16, 34]. Nevertheless, we focus here only on
$\text{Re}[\chi^{(2)}(\omega_{t},\omega_{\tau})]$ and
$\text{Im}[\chi^{(3;1,2)}(\omega_{t},\omega_{\tau})]$, the parts which, in the
absence of the positive time constraint, would be non-vanishing.

#### Evaluation of non-linear susceptibilities with exact diagonalization

We individually compute the non-linear susceptibilities $\chi^{(2)}$ and
$\chi^{(3)}$ using the approach developed in Ref. [30]. Similar to the
experimental setup, this involves applying two Dirac-$\delta$ pulses to the
system and analyzing the resulting time evolution to determine the non-linear
response. In the remainder, we focus on the case in which the field pulses and
the measured magnetization are all aligned, so $\gamma=\beta=\alpha$. The
action of a single pulse at time $t^{\prime}$ on a state $\ket{\psi}$ is
described by

$\ket{\psi^{\prime}}=\exp\left(\text{i}B_{t^{\prime}}M^{\alpha}\right)\ket{\psi},$
(5)

where $B_{t^{\prime}}$ signifies the pulse magnitude. The wavefunction after
two successive pulses, applied at times $t^{\prime}=0$ and $t^{\prime}=\tau$
is

$\ket{\psi(t,\tau)}=e^{-\text{i}\mathcal{H}t}e^{\text{i}B_{\tau}M^{\alpha}}e^{-\text{i}\mathcal{H}\tau}e^{\text{i}B_{0}M^{\alpha}}\ket{\psi}.$
(6)

The magnetization per site after two successive pulses is thus given by

$m^{\alpha}(t,\tau,B_{\tau},B_{0})=\frac{1}{N}\bra{\psi(t,\tau)}M^{\alpha}\ket{\psi(t,\tau)}.$
(7)

We can then calculate the individual non-linear susceptibilities as follows

$\begin{split}\chi^{(2)}_{\alpha\alpha\alpha}(t,\tau)&=\left.\frac{\partial^{2}m^{\alpha}(t,\tau,B_{\tau},B_{0})}{\partial
B_{\tau}\partial B_{0}}\right|_{B_{0}=B_{\tau}=0},\\\
\chi^{(3)}_{\alpha\alpha\alpha\alpha}(t,\tau,0)&=\left.\frac{\partial^{3}m^{\alpha}(t,\tau,B_{\tau},B_{0})}{\partial
B_{\tau}\partial B_{0}^{2}}\right|_{B_{0}=B_{\tau}=0},\\\
\chi^{(3)}_{\alpha\alpha\alpha\alpha}(t,0,\tau)&=\left.\frac{\partial^{3}m^{\alpha}(t,\tau,B_{\tau},B_{0})}{\partial
B_{\tau}^{2}\partial B_{0}}\right|_{B_{0}=B_{\tau}=0}.\end{split}$ (8)

In practice, we numerically compute the derivatives using the central
difference method. For example, $\chi^{(2)}_{\alpha\alpha\alpha}(t,\tau)$ can
be computed using four different $(B_{0},B_{\tau})$ combinations. Setting
$B_{0}=B_{\tau}=B$, it can be written as

$\begin{split}\chi^{(2)}_{\alpha\alpha\alpha}(t,\tau)=\frac{1}{4B^{2}}&\left[m^{\alpha}(t,\tau,B,B)-m^{\alpha}(t,\tau,-B,B)\right.\\\
&\left.\\!\\!\\!\\!-m^{\alpha}(t,\tau,B,-B)+m^{\alpha}(t,\tau,-B,-B)\right].\end{split}$
(9)

This method is similarly applied to the calculation of the third-order
derivatives for $\chi^{(3)}_{\alpha\alpha\alpha\alpha}(t,\tau,0)$ and
$\chi^{(3)}_{\alpha\alpha\alpha\alpha}(t,0,\tau)$, using six combinations of
$(B_{0},B_{\tau})$. Throughout, we use $B_{0}=B_{\tau}=0.001$.

To make the Fourier transform well-behaved with finite time windows, a filter
function $e^{-\eta(t^{2}+\tau^{2})}$ is applied, which broadens the signal in
the frequency domain. We set a rather small broadening of $\eta=0.001$,
resulting in rather sharp peaks in the frequency domain. With this filter, we
time evolve to a maximum time of $t_{\text{max}}=\tau_{\text{max}}=150$. The
sampling intervals $\delta t$ and $\delta\tau$ determine the energy range as
$\omega^{\text{max}}_{t}=\pi/\delta t$ and
$\omega^{\text{max}}_{\tau}=\pi/\delta\tau$. We set $\delta
t=\delta\tau=0.25$. For the time-evolution and application of the pulses, both
of which are achieved by exponential operator multiplications, we use the
package “Expokit.jl”, a Julia implementation of EXPOKIT [35] that efficiently
executes these calculations using Lanczos routines.

Figure 2: Revival of the signal in the small system size in one-dimensional
transverse field Ising (1D-TFIM). Two-dimensional coherent spectroscopy (2DCS)
of 1D-TFIM with the model parameters $J=0.7$, $h_{x}=0.3$, i.e., ferromagnetic
phase. (a-e) Analytically calculated lnear-response
$\text{Im}[\chi^{(1)}_{xx}(\omega_{t})]$, third-order response
$\chi^{(3)}_{xxxx}(t,t+\tau,t+\tau)$ and $\chi^{(3)}_{xxxx}(t,t,t+\tau)$ for
$L=500$, and 2d-FT of them. Spinon pairs manifest as continuous pump-probe
(PP) signal at $\omega_{\tau}=0$ in $\chi^{(3;1)}_{xxxx}$ and spinon-
echo/rephasing (R) signal at $\omega_{t}=-\omega_{\tau}$ in
$\chi^{(3;2)}_{xxxx}$. (f-j) The same quantities obtained by exact
diagonalization (ED) for $L=24$. Though the small $t$ and $\tau$ data are
similar to the that of large system size, the signal starts to revive as $t$
and $\tau$ increase. The discrete character of the PP/R signal is evident. The
non-rephasing (NR) signal at $\omega_{t}=\omega_{\tau}$ is also visible in
$\chi^{(3;2)}_{xxxx}$, which is absent in the thermodynamic limit.

For the finite temperature simulations, we use the canonical thermal pure
quantum state $\left|\phi\right>=\exp(-\beta H/2)\left|\phi_{0}\right>$ [36,
37], where $\left|\phi_{0}\right>$ is a random vector whose norm is
initialized to one, as the initial state. We average the results over $N=10$
random initial states. See Appendix B for the details of the finite
temperature simulation.

Figure 3: Schematics of the rephasing (R) and non-rephasing (NR) signal. In
the experiment, signals are measured at $t,\tau>0$ represented by the shaded
area. (a) In the thermodynamic limit, the signal can be described by a
superposition of continuous poles; therefore, it takes the form of the
integral. R signal decays as $t+\tau$ increases, and NR signal decays as
$t-\tau$ increases. FT of the signal for $t,\tau>0$ only reveals R signal. (b)
In the finite size system, the signal takes the form of the sum due to the
discrete nature of possible $k$. R and NR signals are visible at $t,\tau>0$ as
a revival of the signal. FT of the signal yields discrete R and NR signals.

#### 1D Transverse Field Ising Model (1DTFIM)

The model we consider is the 1D-TFIM with the Hamiltonian

$\mathcal{H}=-J\sum_{i=1}^{L}\sigma_{i}^{z}\sigma_{i+1}^{z}-h_{x}\sum_{i=1}^{L}\sigma_{i}^{x}-h_{z}\sum_{i=1}^{L}\sigma_{i}^{z},$
(10)

where $\sigma_{i}^{\alpha}$ $(\alpha=x,y,z)$ are Pauli matrices at site $i$,
$J>0$ is the nearest neighbor ferromagnetic (FM) coupling, $h_{x}$ is the
transverse field, and $h_{z}$ the longitudinal field. We assume periodic
boundary conditions throughout. When $h_{z}=0$, the model is exactly solvable
via a Jordan-Wigner transformation, and, in the low-field ferromagnetically
ordered phase, the elementary excitations are the kinks, or “spinons”, of the
FM order (domain walls between blocks of aligned spins). These form a gapped
dispersive band of excitations with energies $\lambda_{k}$. However, a single
spin flip $\sigma_{i}^{x}$ generates a pair of spinons, with equal and
opposite momenta. This form of fractionalization results in a continuum
response at zero momentum when probed within the linear response regime.

A finite longitudinal field $h_{z}$ breaks the integrability of the model, and
generates a confining potential for the spinon pairs. Though the model is no
longer exactly solvable, the resulting dynamics can be well described by a
“two-kink approximation”, in which an effective Hamiltonian acting within the
subspace of two spinon states is constructed [38, 39]. At large $h_{z}$, the
spinons are tightly bound, and connect to the single spin-flip excitations of
the high-field limit.

##  III Finite-size effects

Since the available system sizes are limited in ED, we first need to have a
solid understanding of potential finite-size effects. Here, we investigate
finite size effects in two distinct scenarios, namely, (i) a magnetic field
pulse creates a pair of excitations, meaning that linear response exhibits a
continuum of excitations, and (ii) a magnetic field pulse creates only single
excitations, meaning that linear response exhibits a discrete set of
excitation modes. We will see that finite-size effects are more significant in
the former scenario.

###  III.1 Revival of signals

If a system of interest has a continuous spectrum at zero momentum, then
simulating 2DCS with a finite system size can lead to a spurious “revival of
signals” in the time domain. As a result, the Fourier-transformed signal in
the frequency domain deviates from the true behavior expected in the
thermodynamic limit. Here, we illustrate that the origin of this deviation is
a combination of the discrete nature of the spectrum in finite-sized systems
and the positive time constraint, $t>0,\tau>0$, enforced in the experimental
setup.

In Ref. [19], the full 2DCS spectrum for the exactly solvable 1D-TFIM was
obtained (see Appendix C for the derivation). The third-order susceptibility
$\chi^{(3)}_{xxxx}$ reads:

$\begin{split}&\chi^{(3)}_{xxxx}(t_{3},t_{2},t_{1})=\frac{1}{L}\sum_{0<k<\pi}\left[A_{k}^{(1)}+A_{k}^{(2)}+A_{k}^{(3)}+A_{k}^{(4)}\right],\end{split}$
(11)

where

$\begin{split}A_{k}^{(1)}&=-8\sin^{2}\theta_{k}\cos^{2}\theta_{k}\sin(2\lambda_{k}(t_{3}+t_{2}+t_{1})),\\\
A_{k}^{(2)}&=8\sin^{2}\theta_{k}\cos^{2}\theta_{k}\sin(2\lambda_{k}(t_{2}+t_{1})),\\\
A_{k}^{(3)}&=-4\sin^{4}\theta_{k}\sin(2\lambda_{k}(t_{3}+t_{1})),\\\
A_{k}^{(4)}&=-4\sin^{4}\theta_{k}\sin(2\lambda_{k}(t_{3}-t_{1})),\end{split}$
(12)

with $\lambda_{k}$ the single-spinon energies, and the angle $\theta_{k}$
defined by $\tan\theta_{k}=(2J\sin k)/(2J\cos k-2h_{x})$.

Figure 2 compares the analytical and numerical $\chi^{(1)}_{xx}$ and
$\chi^{(3;1,2)}_{xxxx}$ of the 1D-TFIM in the FM phase. The analytical data is
obtained for a large system size of $L=500$ using Eq. (12), representing the
expected behavior in the thermodynamic limit, and the numerical data is
obtained for $L=24$ with ED. For the ED data, we checked that the error
between the analytical result and the ED result for the same system size is
less than $10^{-5}$ in the time domain data, and no apparent increasing trend
is evident during the chosen time evolution window.

There are unsurprisingly a number of similarities between the two system
sizes. In $\chi^{(3;1)}_{xxxx}$, there is a pump-probe (PP) signal along the
$\omega_{t}$ axis at $(\omega_{t},\omega_{\tau})=(2\lambda_{k},0)$. In
$\chi^{(3;2)}_{xxxx}$, there is an anti-diagonal signal at
$\omega_{t}=-\omega_{\tau}=2\lambda_{k}$ corresponding to the “spinon-echo” or
rephasing (R) signal, originating from $A_{k}^{(4)}$. For the smaller system
size, the discrete nature of these signals simply results from the
discreteness of the momentum $k$. This is most evident in the linear response
$\chi^{(1)}_{xx}$, where the $L=24$ system contains a visibly discrete set of
$L/2$ peaks, corresponding to the $L/2$ allowed pairs of spinons with momenta
$k$ and $-k$, while the larger system size effectively has a continuum of
excitations due to its much larger number of $k$ points.

There are two obvious differences though between the two system sizes. In the
time domain, there are periodic revivals of the signal in the small system
size used for the ED that do not appear in the larger system size, and, in the
frequency domain, there are additional diagonal non-rephasing (NR) signals in
the small system size which are absent in the larger system size. To
understand this, let’s first consider a one-dimensional time domain, and for
simplicity assume that the energies of the spinons are equally spaced with a
uniform frequency spacing of $\Delta\omega\propto 1/L$. If one superposes
multiple sine waves, $\sum_{k>0}\sin(\omega_{k}t^{\prime})$, with a uniform
frequency spacing of $\Delta\omega$, then the wavelength of the resulting
beating pattern will be proportional to $1/\Delta\omega$. In the thermodynamic
limit, $L\rightarrow\infty$, $\Delta\omega\rightarrow 0$, and the spinon
spectrum becomes continuous, which means that the wavelength diverges and the
periodic signal disappears from the positive time axis, leaving a maxima at
$t^{\prime}=0$.

In the two-dimensional time domain of relevance here, a similar scenario plays
out. As illustrated in Fig. 3, for $L\rightarrow\infty$, we need to superpose
a continuum of sine waves, $\int dk\,g(k)\sin[\omega_{k}(t\pm\tau)]$, with
$g(k)$ simply representing matrix element factors. For the $t-\tau$ case,
there is a maxima along the line $t=\tau$, and for the $t+\tau$ case, the
maxima is along $t=-\tau$. Focusing on the relevant quadrant with
$t>0,\tau>0$, only the $t=\tau$ line can be observed which, after Fourier
transforming, is exactly the rephasing R signal. The $t=-\tau$ line is absent
from the purely positive time window, and thus the absence of the
corresponding non-rephasing NR signal is generally expected when the system
has a continuous spectrum in the thermodynamic limit. On the other hand, for
$L$ finite, we are adding a finite number of frequencies,
$\sum_{k}\,g(k)\sin[\omega_{k}(t\pm\tau)]$, resulting in periodic signals with
wavelengths roughly proportional to $L$. This means that, in the positive time
quadrant, $t>0,\tau>0$, both rephasing and non-rephasing lines are present,
and thus the Fourier transformed spectra have both signals present. The
appearance of the NR signal is thus a consequence of the discreteness of the
spectrum. In the ED, our time window for the FT in both $t$ and $\tau$ extends
to several times $L$, so we indeed observe both R and NR signals and,
similarly, we are able to resolve the discrete nature of the signals.

It should be noted that there is a straightforward solution in this particular
case to avoid such revivals of signals in the time domain. One can simply
choose the time window small enough such that no revivals occur within it
(alternatively, one can choose a large damping factor $\eta$ in the filter
function $e^{-\eta(t^{2}+\tau^{2})}$). Thus, the NR signal will be absent and
the R signal will be continuous, both as expected for the thermodynamic limit.
For large system sizes, this is relatively straightforward as the periodic
revivals are anyway spread far apart in time, but for the small system sizes
accessible with ED this heavily constrains the available time window and hence
leads to an extremely poor resolution in frequency space. More importantly,
it’s, in general, not possible to know a priori whether the revivals are an
artifact, to be avoided, or a feature, to be included, i.e. whether the true
physics is really a continuum or simply a dense set of discrete excitations
(as, for example, observed in the confining scenario of the next section).

In Fig. 3(a), it is the positive time constraint, $t>0,\tau>0$, that ensures
that it is only the R signal that is present in the thermodynamic limit. If we
were to include both positive and negative $\tau$, then both the NR and R
lines would be visible in the time domain, and hence one would observe both
diagonal NR and anti-diagonal R signals in the 2DCS spectrum. In addition, the
artificial broadening induced by the positive time constraint in the Fourier
transform would be absent (in the case of the imaginary part). Taken together,
this means that, for the relaxed constraint, $t>0$ only, the small system size
ED result is simply a discretized version of the results in the thermodynamic
limit, with all of the same qualitative features and signals (see Appendix D
for an illustration of this).

Figure 4: Twisted Kitaev Model. (a) Analytical result for $L=500$. In addition
to the R signal, there are streak features in both first and fourth quadrant,
originating from the interplay between different spinon modes. (b) ED result
for $L=24$. Additional diagonal NR signal at $\omega_{t}=\omega_{\tau}$ is
observed.

#### Twisted Kitaev model

As another example in which identical finite size effects are present, we show
results for another exactly solvable model, the twisted Kitaev model (TKM)
[40, 41, 42, 43]. The Hamiltonian is

$\mathcal{H}=-J\sum_{i}^{L^{\prime}}\left[\tilde{\sigma}_{2i-1}(\theta)\tilde{\sigma}_{2i}(\theta)+\tilde{\sigma}_{2i}(-\theta)\tilde{\sigma}_{2i+1}(-\theta)\right],$
(13)

where $L^{\prime}=L/2$,
$\tilde{\sigma}_{i}(\theta)=\sigma_{i}^{z}\cos(\theta/2)+\sigma_{i}^{y}\sin(\theta/2)$,
and $\theta$ is the “twist” angle. For $0\leq\theta<\pi/4$, the ground state
is a doubly degenerate FM state, polarized along the $z$-direction. The
elementary excitations are again spinons (kinks, or domain walls, of the FM
state), but now there are two kinds of spinons, with dispersions $l_{k}$ and
$\lambda_{k}$, where $k=2\pi n/L^{\prime}$ and $n=1,3,\cdots,L^{\prime}$ [41].

Analytical expressions for the second and third-order susceptibilities for the
TKM have been derived in Ref. [43], with the result for $\chi^{(3;2)}_{xxxx}$
shown in Fig. 4(a) for a large system size of $L=500$, again representing the
behavior of the thermodynamic limit. Similar to the 1D-TFIM, a magnetic pulse
along the $x$-direction only excites spinon pairs, with equal and opposite
momenta; therefore there appears again a spinon-echo R signal. In addition,
due to the interplay between the $l_{k}$ spinons and $\lambda_{k}$ spinons,
there are additional streak features for both positive and negative
$\omega_{\tau}$ [41]

In ED, on an $L=24$ chain, we observe all of the features of the analytical
result in the thermodynamic limit, but with discretized streaks due to the
discrete number of $k$ points. Again, as in the TFIM, the main difference, in
the frequency domain, is the additional diagonal NR signal, which is absent in
the thermodynamic limit. The small system sizes available in ED are again able
to capture the key qualitative features of the 2DCS spectrum, with the caveat
that an additional NR signal presents itself.

Figure 5: Single spin flip excitations in paramagnetic phase (PM). (a and b)
The result for $J=0$, where the system size independent exact result is
available. R and NR are observed at
$(\omega_{t},\omega_{\tau})=(2h_{x},-2h_{x})$ and $(2h_{x},2h_{x})$,
respectively. (c and d) ED result for $L=24$ with $J=0.1$. (e) Schematics of
the processes responsible for the NR signal. The intermediate state Q is
either the ground states (red) or two-spin flip state (blue). In the presence
of the interaction $J$, the two-spin flip states with zero total momentum
constitute a continuum, which ends up with discrete poles in the finite size
system. (f) Actual position of the poles for $L=24$. The data is obtained by
Fourier transforming both positive and negative time data, which eliminates
the phase twisting effect. Slightly different frequencies of the signals with
opposite signs are superposed, which manifests diverging behavior in the time
domain data (c). Figure 6: Spinon-confinement via an applied longitudinal
field with parameters $J=0.7,h_{x}=0.06$ in (a, b) the weakly confined limit,
$h_{z}=0.06$, and (c, d) the strongly confined limit, $h_{z}=0.4$. All
relevant susceptibilities, $\chi^{(1)}_{xx},\chi^{(2)}_{xxx}$, and
$\chi^{(3;1,2)}_{xxxx}$ are shown. The results in (a), (c) (upper row) are
obtained analytically using the two-kink approximation in the thermodynamic
limit, while the results in (b), (d) (lower row) are obtained with ED for
$L=24$. Figure 7: Evolution of the spinon-echo rephasing signal. A comparison
of $\chi^{(3;2)}_{xxxx}$, which contains the spinon-echo signal, for different
longitudinal fields computed with the two-kink approximation (top row) and ED
(bottom row). In both cases, the sharp spinon-echo signal at $h_{z}=0$ is
broadened as $h_{z}$ increases, ending as a set of discrete peaks in the
strongly confined regime. Note that since the two-kink approximation does not
include four-kink states, it overestimates the broadening of the spinon-echo
signal.

###  III.2 Discrete spin-flip excitations in the paramagnetic phase

Next, we discuss finite size effects in the paramagnetic (PM) phase of the
1D-TFIM, i.e. $h_{x}>J$. In this phase, a magnetic field pulse along the
$z$-direction creates a single spin flip excitation with momentum $k=0$, and
the spectrum is expected to be discrete even in the thermodynamic limit. We
thus expect that ED should be able to produce qualitatively similar results as
those expected in the thermodynamic limit. Also, unlike the case of the FM
ground state, we expect to observe an NR signal corresponding to single spin-
flip excitations.

Figure 5(a) and (b) show $\chi^{(3;2)}$ for the trivial limit of $J=0,h_{x}=1$
($h_{x}/J\rightarrow\infty$), in which the result can be obtained
analytically. In this transverse field-only limit, the ground state can be
written as $\ket{0}=\ket{\rightarrow\rightarrow\cdots\rightarrow}$.
Excitations involving either a single or double spin-flip at sites $i$ or
$(i,j)$ are represented by $\ket{\leftarrow}_{i}=\sigma_{i}^{z}\ket{0}$ for a
single site, and
$\ket{\leftarrow}_{i}\ket{\leftarrow}_{j}=\sigma_{i}^{z}\sigma_{j}^{z}\ket{0}\,(i\neq
j)$ for two sites, respectively. These are energy eigenstates with energies
$E=2h_{x}$ and $E=4h_{x}$ respectively.

Keeping in mind that the matrix elements involved for $\chi^{(3)}$ can be
written as $m^{z}_{0R}m^{z}_{RQ}m^{z}_{QP}m^{z}_{P0}$, we can think of two
types of processes: one with $\ket{Q}=\ket{0}$, simplifying the matrix element
to $|m^{z}_{0R}|^{2}|m^{z}_{0P}|^{2}$, and another with $\ket{Q}$ involving
double spin-flip states. In both cases, $\ket{P}$ and $\ket{R}$ are single
spin-flip states. Noting that $(\sigma_{z}^{i})^{2}$ is the identity operator,
the number of contributions scales as $L^{2}$ for $\ket{Q}=\ket{0}$ processes,
and as $2L(L-1)$ for processes with $\ket{Q}$ involving double spin-flip
states. Combining these together, it can be shown that the summation of the
different $R_{a}$ in Eq. (2) for both types of processes results in a
cancellation of the terms proportional to $L^{2}$, leaving only the terms
proportional to $L$. Dividing by $L$, we thus obtain system size-independent
peaks at $(\omega_{t},\omega_{\tau})=(2h_{x},2h_{x})$ and $(2h_{x},-2h_{x})$,
corresponding to NR and R signals respectively. Unlike in the low-field FM
phase, the NR signal does not vanish here as the pole is now discrete.

Next, we consider the effects of a finite Ising interaction $J$, by which the
spin-flip excitation obtains a dispersion. Figure 5(c) and (d) shows the ED
result with $L=24$ for $J=0.1,h_{x}=0.9$. Similar to the non-interacting case,
we observe NR and R signals at
$(\omega_{t},\omega_{\tau})\approx(2h_{x}-2J,2h_{x}-2J)$ and
$(2h_{x}-2J,-2h_{x}+2J)$, respectively, whose energies corresponds to that of
the single spin-flip excitation at $k=0$. In the time domain data, the
amplitude of the signal increases as $t$ increases, see Fig. 5(c). Similar to
the non-interacting case, the contribution from the processes with
$\ket{Q}=\ket{0}$ should be partially canceled out by the processes with
$\ket{Q}$ involving two-spin flip states. However, since the double spin-flip
states with $k=0$ now form a continuum in the thermodynamic limit, this not
only leads to the cancellation of the $\ket{Q}=\ket{0}$ processes, but also
results in a tail of signals in the frequency domain [33], as schematically
shown in Fig. 5(e). Indeed, in the zoomed in data of Fig. 5(f), both the
numerical NR and R signals exhibit a tail in the frequency domain. The
superposition of slightly different signal frequencies with opposite signs
generates the diverging behavior in the time domain data, in line with the
“long-time divergences” discussed in Ref. [18].

##  IV Breaking integrability and Confinement

Having understood the impact of finite size effects, we now consider the
addition of a longitudinal field, which introduces a confining potential for
the spinon pairs. This explicitly breaks the integrability of the model,
meaning it is no longer exactly solvable via a Jordan-Wigner transformation.
However, it is possible to obtain analytical results using the two-kink
approximation, derived by projecting the full Hamiltonian into the subspace of
two-spinon states [39, 44, 38] (see Appendix E for details of the analytical
calculation). In Ref. [29], they explored the impact of confinement on the
second-order response $\chi^{(2)}$. Here, we focus on $\chi^{(3)}$, and in
particular the evolution of the spinon-echo rephasing signal, predicted as a
key fingerprint of fractionalization [19], with increasing longitudinal field.

Figure 6 compares the analytical and ED results for $\chi_{xx}^{(1)}$,
$\chi_{xxx}^{(2)}$, and $\chi_{xxxx}^{(3;1,2)}$ for both the weakly confined,
with $h_{z}=0.06$, and strongly confined, with $h_{z}=0.4$, cases. The
confining potential leads to a set of spinon bound states with a discrete
energy spectrum $E_{n}$, clearly visible in the linear response susceptibility
$\chi^{(1)}_{xx}$. The energy spacing between these states increases as
$h_{z}$ increases, with for example only the lowest three visible within the
frequency range $0<\omega_{t}<6$ in the strongly confined limit, $h_{z}=0.4$.
The second-order susceptibilities are consistent with Ref. [29] and are
discussed in depth there.

In the weakly confined limit, $h_{z}=0.06$, we see that even a small
longitudinal field produces significant deviations from the integrable (zero
confinement) results. First, in $\chi_{xxxx}^{(3;1)}$, there are a sequence of
additional peaks surrounding the on-axis pump-probe signal. More
interestingly, in $\chi_{xxxx}^{(3;2)}$, the spinon-echo R signal now consists
of a grid of discrete peaks spread over a wide frequency range
$(2\leq\omega_{t},-\omega_{\tau}\leq 4)$, in contrast to the sharp anti-
diagonal line in the integrable case. The nature of these peaks can be
discerned by examining in more detail the structure of the results from the
two-kink approximation. The matrix elements involved can be written as
$m^{x}_{0R}m^{x}_{RQ}m^{x}_{QP}m^{x}_{P0}$, with the grid of rephasing peaks
corresponding to the case where the ground state is the intermediate state
$\ket{Q}=\ket{0}$, and $\ket{P}$ and $\ket{R}$ are spinon bound states (see
Fig. 13 in Appendix E for a plot of the contribution of just the
$\ket{Q}=\ket{0}$ case). Thus, we observe peaks at energies
$(\omega_{t},\omega_{\tau})=(E_{n},-E_{l})$, with intensities proportional to
$\left|m^{x}_{n0}\right|^{2}\left|m^{x}_{l0}\right|^{2}$ (peaks on the anti-
diagonal have $n=l$, and on the off-diagonal $n\neq l$). These peaks contain
similar information to linear response, as it is matrix elements of the form
$\left|m_{n0}\right|^{2}$ that appear there as well.

In the strongly confined limit, $h_{z}=0.4$, we see that the number of peaks
visible has sharply reduced, primarily due to the increased spacing in energy
between bound states. In $\chi_{xxxx}^{(3;2)}$, the peaks along the anti-
diagonal, as well as their associated cross-peaks, again originate from the
case where the ground state is the intermediate state $\ket{Q}=\ket{0}$, and
$\ket{P}$ and $\ket{R}$ are spinon bound states. On the other hand, the
terahertz rectification (TR) peaks on, and just off, the $\omega_{\tau}$ axis
originate from the case where the intermediate state is a two-kink state (see
Fig. 13 in Appendix E for a plot of just the intermediate two-kink
contributions).

It’s important to note that here the discreteness of the spectrum is not a
finite-size effect, as it was in the integrable FM phase, but rather it has a
clear physical origin. Therefore, the NR signals observed in
$\chi_{xxxx}^{(3;1,2)}$ are expected to remain even in the thermodynamic
limit. Indeed, the analytical results, qualitatively consistent with the ED
data, are computed in the thermodynamic limit (except for system size-
dependent intensity factors). Thus, in this particular case, the appearance of
a finite NR signal is a direct consequence of confinement converting the
spinon continuum into a discrete set of spinon bound states.

In Fig. 7, we show the full evolution of $\chi_{xxxx}^{(3;2)}$ with increasing
longitudinal field, computed both within the two-kink approximation and ED.
The sharp line features of the zero field case, corresponding to free,
deconfined spinon pairs, quickly decompose into a discrete set of peaks due to
the confining potential. The rephasing signal decomposes into a grid pattern
of discrete peaks, a non-rephasing signal appears with a similar grid of
diagonal and cross peaks, and the terahertz rectification signal along the
$\omega_{\tau}$ axis also decomposes into discrete peaks which expand along
both directions in frequency space. Also note that the two-kink approximation
naturally becomes more accurate the larger the longitudinal field.

Comparing the results from the two-kink approximation and ED, though they
agree on most of the important qualitative features, there are some
differences. This is due to the fact that we omit four-kink states, and
higher, in the two-kink approximation, while there is no such approximation in
the ED. For example, in $\chi_{xxxx}^{(3;1)}$, four-kink processes generate a
signal at $(\omega_{t},\omega_{\tau})\approx(E_{n},2E_{n})$, which is absent
in the two-kink approximation results but which we observe in the ED data.
Analogous to the PM phase, the four-kink states also introduce cancellation
terms in $\chi_{xxxx}^{(3;1,2)}$, regulating the apparent system size-
dependent behavior of R/NR and PP signals from ground state processes, leading
ultimately to $L$-independent signals.

There are also quantitative differences between the two-kink approximation and
ED results, especially in $\chi_{xxxx}^{(3;1,2)}$. The analytical predictions
tend to overestimate the intensities when compared to the ED data. This is
evident in, for example, the pump-probe signal in $\chi_{xxxx}^{(3;1)}$, and
the R and NR signals in $\chi_{xxxx}^{(3;1,2)}$. As in the PM phase discussed
previously, these differences stem from the proximity of poles of the
processes where $\ket{Q}$ is a four-kink state and $\ket{Q}$ is the ground
state. In ED, this proximity, coupled with a limited frequency resolution,
leads to an intensity reduction due to the cancellation of poles from the two
distinct processes. Indeed, time domain data (not shown) reveals an increasing
trend in the signal for $\chi_{xxxx}^{(3;1,2)}$ as time $t$ increases,
suggesting the presence of two peaks with opposite signs and slightly
different frequencies in these signals.

Finally, it’s worth noting that although the two-kink approximation and ED
exhibit broadly qualitative agreement for $\chi_{xxxx}^{(3;2)}$, the
individual $R_{a}$ that contribute to it, see Eqn. (2) (and Appendix A), are
actually quite different in the two methods. For example, the rephasing signal
in the two-kink approximation solely comes from $R_{2}$, while in the ED it
comes from a complex partial cancellation of $R_{1}$ and $R_{2}$. In some
cases, it may thus be useful to not only compare the non-linear
susceptibilites, but also the individual $R_{a}$, to better compare and
contrast different approaches. See Appendix F for more details on this point.

Figure 8: Finite temperature spectra of 1D-TFIM. A comparison between
analytical (a,b) and numerical (c) results of finite temperature with $L=500$
and $L=16$ is displayed. Parameters $J=0.7,h_{x}=0.3$ are used for
calculation. Figure 9: Temperature dependence of
$T\chi^{(3;2)}_{xxxx}(\omega_{t},\omega_{\tau})$ of the 1D-TFIM in a
longitudinal field from ED. The parameters $J=0.7,h_{x}=0.3$. (a-c) $h_{z}=0$,
i.e. exactly solvable case with inverse temperature $T=10.0,1.00,0.33$,
respectively. (d-f) $h_{z}=0.06$. (g-i) $h_{z}=0.4$.

##  V Finite temperature

Finally, we make use of one of the other advantages of ED and extend our study
to finite temperatures. As in the zero-temperature case, we first investigate
the exactly solvable case with $h_{z}=0$. Figure 8 shows $\chi^{(3;2)}_{xxxx}$
calculated for $L=500$ and $L=16$ at a temperature $T=1.43$. For this model,
one should note that this particular response function is rather special as
the only change with temperature is in the intensities of the signals. For
$L=500$, the overall peak intensity is suppressed as temperature increases,
but the rephasing R signal remains sharp in the two-dimensional frequency
space. For $L=16$, there are twice as many peaks in the signal, which is due
to the fact that at finite temperatures we also have a contribution from the
parity odd sector of the Jordan-Wigner transformed fermionic Hamiltonian of
the TFIM model (see Appendix C for details). Similar to the zero-temperature
case, we also observe a non-rephasing NR signal due to the finite system size,
and no other additional signal appears. Figure 8 (c) shows ED result for
$L=16$ at $T=1.43$. The quantitative agreement between the analytical and ED
results helps to validate our method for finite temperature.

Next, we investigate the finite-temperature 2DCS in the presence of the spinon
confining potential. The impact on $\chi^{(3;2)}$ is shown in Fig. 9 for the
zero, weakly, and strongly confined cases, with the data normalized by
$\beta=1/T$ so that they can be plotted with a common color scale. At a
sufficiently low temperature of $T=0.33$, which is below the specific heat
peaks for all three cases (see Fig. 11 in Appendix B), the intensities and
peak positions are close to the zero-temperature result. As temperature
increases, the signal intensities decrease for all values of $h_{z}$. For the
weakly confined case, $h_{z}=0.06$, the discreteness of the low-temperature
signal is smeared out as temperature increases, and the signal becomes
continuous. In particular, at $T=10.0$, the spectrum exhibits qualitative
similarities to the $h_{z}=0.00$ spectrum found in the thermodynamic limit,
suggesting the impact of confinement can only be resolved at low temperatures.
For $h_{z}=0.4$, the discreteness of the signal is preserved even at high
temperatures. We also observe additional peaks along the line
$\omega_{t}\approx 1$, which arise from the processes involving transitions
between two-kink states.

##  VI Discussion

The calculation of non-linear dynamical response functions for a quantum many-
body system is a formidable and challenging task. We explored here, on the one
hand, the utility of using ED to calculate non-linear susceptibilities in the
1D-TFIM, and, on the other hand, the impact of confinement on the unique
signatures of spinon fractionalization in the model. One of the key advantages
of ED, the ability to simulate to long times without any corresponding
increase in computational complexity, can be used to obtain high-resolution
2DCS spectra in frequency space. However, the key disadvantage, the relatively
small system sizes available with ED, generates additional issues with
interpreting the spectra due to potential finite size effects. We have
discussed in detail these issues in the context of the 1D-TFIM and its two
distinctive ground states. In both limits, ED can capture many of the
important qualitative features expected, though care must be taken in
distinguishing whether discrete excitations arise due to finite-size effects
or due to some underlying physical mechanism. Using this knowledge, we were
able to explore the impact of spinon confinement on the third-order
susceptibilities, and, in particular, on the characteristic spinon-echo
rephasing signal. Combining ED with an analytical two-kink approximation, we
were able to show how even a moderate longitudinal field can break up the
sharp anti-diagonal signal and induce a visible non-rephasing signal.

In terms of experimental relevance, there are a number of materials whose
dominant interactions can be written in the form of a 1D-TFIM. As an example,
in CoNb_2O_6, the FM chain of Co^2+ ions is considered to be a good
realization of the 1D-TFIM [38, 39, 27, 42]. In this material, the weak
interchain interaction effectively introduces a longitudinal field, and the
discrete spectrum of spinon bound states has been experimentally observed with
linear response probes [39, 38]. Therefore, with 2DCS, we would expect
qualitative features of the 2D spectrum computed here to also be observed in
the material. However, it would be necessary to investigate the effects of
additional terms, which are proposed to be relevant in the material, such as
XY and bond-dependent interactions, in order to gain a quantitative
understanding of the expected experimental spectrum.

Figure 10: Fourier transform of the non-linear magnetization $m_{\text{NL}}$
with increasing magnitude of incident THz pulses, $B_{0},B_{\tau}$, obtained
by ED. The parameters for the weakly confined case,
$J=0.7,h_{x}=0.3,h_{z}=0.06$ are used.

An important caveat, relevant for any experimental comparison, is that in
reality it is the non-linear magnetization $m^{\alpha}_{\text{NL}}$ that is
measured in experiment, as opposed to the individual susceptibilities. The
non-linear magnetization contains contributions from $\chi^{(2)}$,
$\chi^{(3;1)}$, and $\chi^{(3;2)}$, as well as higher order terms. The ratio
of the intensities of these terms is, roughly speaking, determined by the
strength of the THz field; $n$-th order terms are proportional to the $n$-th
power of the THz field. Therefore, if there is a finite $\chi^{(2)}$, it will
typically be the dominant contribution to the experimental signal. Figure 10
shows the Fourier transform of the non-linear magnetization
$m^{x}_{\text{NL}}$ obtained by ED with different strengths of the magnetic
field pulses. The parameters for the weakly confined case,
$J=0.7,h_{x}=0.3,h_{z}=0.06$ are used. Up to $B_{0}=B_{\tau}=0.01$, the
visible signal is dominated by the $\chi^{(2)}$ contribution. Only when
$B_{0}=B_{\tau}=0.1$ is the visibility of the third-order contribution
comparable to the second-order one, including the third-order pump-probe and
rephasing signals. Finally, at $B_{0}=B_{\tau}=1.0$, the field pulses can no
longer really be considered as perturbations, and the signal becomes dominated
by the pump-probe response.

Looking forward, there are a multitude of non-linear many-body phenomena that
can be further studied using ED and complementary analytical insights. For
example, in ordered phases, the 2DCS signatures of non-linear magnon decay and
lifetime effects, or the interplay between distinct excitation modes. On the
computational front, it would be useful to test the capabilities of ED in
studying 2D spin models. As an example, there are already known analytical
results for the non-linear response of the Kitaev honeycomb model [22, 26,
25], which can be used to compare and contrast with ED results and to diagnose
potential finite size effects. Taken together, it’s clear that there is a wide
world of exciting physics beyond the linear response regime still waiting to
be explored.

###### Acknowledgements.

We acknowledge partial funding from the DFG within Project-ID 277146847, SFB
1238 (projects C02, C03). This work has benefitted from multiple exchanges
during the 2023 KITP workshop “A New Spin on Quantum Magnets”, supported in
part by grant NSF PHY-1748958 to the Kavli Institute for Theoretical Physics
(KITP). The numerical simulations were performed on the JUWELS cluster at the
Forschungszentrum Juelich and the Noctua2 cluster at PC2 in Paderborn.

##  Appendix A Full expressions for the non-linear susceptibilities

We complement our technical discussion in the main manuscript here by
providing the full expressions for the non-linear susceptibilities,
$\chi^{(2)}$ and $\chi^{(3)}$. The second-order susceptibility can be
expressed as

$\chi^{(2)}_{\alpha\beta\gamma}(t_{2},t_{1})=-\frac{2}{N}\text{Re}[R_{1}-R_{2}],$
(14)

with

$\begin{split}R_{1}&=\left<M^{\alpha}(t_{2}+t_{1})M^{\beta}(t_{1})M^{\gamma}(0)\right>,\\\
&=\sum_{PQ}m^{\alpha}_{0Q}m^{\beta}_{QP}m^{\gamma}_{P0}e^{-\text{i}E_{P}t_{1}}e^{-\text{i}E_{Q}t_{2}},\\\
R_{2}&=\left<M^{\beta}(t_{1})M^{\alpha}(t_{2}+t_{1})M^{\gamma}(0)\right>,\\\
&=\sum_{PQ}m^{\beta}_{0Q}m^{\alpha}_{QP}m^{\gamma}_{P0}e^{-\text{i}E_{P}t_{1}}e^{-\text{i}(E_{P}-E_{Q})t_{2}}.\end{split}$
(15)

The third-order susceptibility $\chi^{(3)}$ can be expressed as

$\chi^{(3)}_{\alpha\beta\gamma\delta}(t_{3},t_{2},t_{1})=\frac{2}{N}\text{Im}[R_{1}+R_{2}+R_{3}+R_{4}],$
(16)

with

$\begin{split}R_{1}&=\left<M^{\gamma}(t_{1})M^{\beta}(t_{2}+t_{1})M^{\alpha}(t_{3}+t_{2}+t_{1})M^{\delta}(0)\right>\\\
&=\sum_{PQR}m^{\gamma}_{0R}m^{\beta}_{RQ}m^{\alpha}_{QP}m^{\delta}_{P0}\\\ &\
\ \ \ \times
e^{-\text{i}E_{P}t_{1}}e^{-\text{i}(E_{P}-E_{R})t_{2}}e^{-\text{i}(E_{P}-E_{Q})t_{3}},\\\
R_{2}&=\left<M^{\delta}(0)M^{\beta}(t_{2}+t_{1})M^{\alpha}(t_{3}+t_{2}+t_{1})M^{\gamma}(t_{1})\right>\\\
&=\sum_{PQR}m^{\delta}_{0R}m^{\beta}_{RQ}m^{\alpha}_{QP}m^{\gamma}_{P0}\\\ &\
\ \ \ \times
e^{\text{i}E_{R}t_{1}}e^{-\text{i}(E_{P}-E_{R})t_{2}}e^{-\text{i}(E_{P}-E_{Q})t_{3}}\\\
R_{3}&=\left<M^{\delta}(0)M^{\gamma}(t_{1})M^{\alpha}(t_{3}+t_{2}+t_{1})M^{\beta}(t_{2}+t_{1})\right>\\\
&=\sum_{PQR}m^{\delta}_{0R}m^{\gamma}_{RQ}m^{\alpha}_{QP}m^{\beta}_{P0}\\\ &\
\ \ \ \times
e^{\text{i}E_{R}t_{1}}e^{\text{i}E_{Q}t_{2}}e^{-\text{i}(E_{P}-E_{Q})t_{3}},\\\
R_{4}&=\left<M^{\alpha}(t_{3}+t_{2}+t_{1})M^{\beta}(t_{2}+t_{1})M^{\gamma}(t_{1})M^{\delta}(0)\right>\\\
&=\sum_{PQR}m^{\alpha}_{0R}m^{\beta}_{RQ}m^{\gamma}_{QP}m^{\delta}_{P0}\\\ &\
\ \ \ \times
e^{-\text{i}E_{P}t_{1}}e^{-\text{i}E_{Q}t_{2}}e^{-\text{i}E_{R}t_{3}}.\end{split}$
(17)

##  Appendix B Preparation of the thermal state

The thermal state $\left|\phi\right>$, which is used to calculate specific
heat and 2DCS, is obtained by the imaginary time evolution of the initial
state $\left|\phi_{0}\right>$. Instead of computing $\left|\phi\right>$
directly, we first compute $e^{\beta(l-H)/2}\ket{\phi_{0}}$, where $l$ is set
to larger than the largest eigenvalue of $H$. The expansion of the exponential
operator is given by:

$\begin{split}e^{\beta(l-H)/2}\ket{\phi_{0}}&=\sum_{k=0}^{\infty}\frac{(\beta/2)^{k}}{k!}(l-H)^{k}\ket{\phi_{0}}\\\
&=\sum_{k=0}^{\infty}\frac{(\beta/2)^{k}}{k!}\ket{k},\end{split}$ (18)

It is known that the relevant terms in the sum is localized in the range
$\left|k^{*}-k\right|=o(L)$, where the temperature of normalized
microcanonical thermal pure quantum state
$\ket{\phi_{k^{*}}}=\ket{k^{*}}/\|\ket{k^{*}}\|$ is close to $1/\beta$ [36,
37]. We choose the number of iterations as $k_{\text{max}}=2000$, which is
sufficiently large for the temperature range we are interested in in this
study. Then, multiplying scalar $e^{-\beta l/2}$ to this state gives
$\left|\phi\right>$. Similarly to the linear response [45], we calculate the
non-linear susceptibility from the time evolution of $\ket{\phi}$.

Figure. 11 shows the calculated specific heat $c_{V}$ with different values of
$h_{z}$. The temperatures used for the finite temperature 2DCS are chosen to
cover the different temperature regimes separated by the peaks of $c_{V}$.

Figure 11: Specific heat $c_{V}$ of 1D-TFIM in the longitudinal field from ED
with $L=16$. The parameters $J=0.7,h_{x}=0.3$, and several values of $h_{z}$
are considered.

##  Appendix C Zero/Finite temperature 2DCS in TFIM

Here, we provide a 2DCS in TFIM based on Ref. [19]. With the Jordan-Wigner
transformation,

$\begin{split}\sigma^{x}_{i}&=2c_{i}^{\dagger}c_{i}-1,\\\
\sigma^{z}_{i}&=(c_{i}^{\dagger}+c_{i})\exp\left(\text{i}\pi\sum_{j=1}^{i-1}c_{j}^{\dagger}c_{j}\right),\end{split}$
(19)

the Hamiltonian reads

$\begin{split}\mathcal{H}_{p}&=-J\sum_{i=1}^{L-1}(c_{i}^{\dagger}-c_{i})(c_{i+1}^{\dagger}+c_{i+1})\\\
+\
&(-1)^{p}J(c_{L}^{\dagger}-c_{L})(c_{1}^{\dagger}+c_{1})-h_{x}\sum_{i=1}^{L}(2n_{i}-1),\end{split}$
(20)

where $p$ is the parity of the number of fermions. Among $2^{L}$ Fock states
of each fermionic Hamiltonian, only $2^{L-1}$ states with the same parity as
$p$ are correct states of the spin Hamiltonian with the periodic boundary
condition. The ground state is known to lie in the parity even sector,
therefore the contribution from parity odd sector becomes finite only at the
finite temperature 2DCS. We introduce the momentum representation of the
fermionic operators as $c_{j}=(1/\sqrt{L})\sum_{k}e^{\text{i}kj}c_{k}$, where
$k=\pm(2n-1)\pi/L$ with $n=1,\cdots,L/2$ for parity even sector, and $k$
satisfies $k=2\pi n/L$ with $n=-L/2+1,\cdots,0,\cdots,L/2$ for parity odd
sector. Next, we introduce a pair Hamiltonian with momentums $k$ and $-k$ as

$\begin{split}\mathcal{H}_{k}&=(c_{k}^{\dagger},c_{-k})\left(\begin{array}[]{cc}\epsilon_{k}&\text{i}\Delta_{k}\\\
-\text{i}\Delta_{k}&-\epsilon_{k}\end{array}\right)\left(\begin{array}[]{c}c_{k}\\\
c_{-k}^{\dagger}\end{array}\right),\end{split}$ (21)

where $\epsilon_{k}=-2J\cos k-2h_{x}$ and $\Delta_{k}=-2J\sin k$, and a
Hamiltonian for $k=0$ and $k=\pi$ of the parity odd sector as

$\mathcal{H}_{k=0,\,\pi}=-2J(n_{0}-n_{\pi})-2h_{x}(n_{0}+n_{\pi}-1).$ (22)

Now, we can rewrite the Hamiltonian of the even sector as

$\mathcal{H}_{0}=\sum_{k>0}\mathcal{H}_{k},$ (23)

and the Hamiltonian of the odd sector as

$\mathcal{H}_{1}=\sum_{0<k<\pi}\mathcal{H}_{k}+\mathcal{H}_{k=0,\,\pi}.$ (24)

Performing the Bogoliubov transformation:

$\left(\begin{array}[]{c}c_{k}\\\
c_{-k}^{\dagger}\end{array}\right)=\left(\begin{array}[]{cc}\cos\frac{\theta_{k}}{2}&-\text{i}\sin\frac{\theta_{k}}{2}\\\
-\text{i}\sin\frac{\theta_{k}}{2}&\cos\frac{\theta_{k}}{2}\end{array}\right)\left(\begin{array}[]{c}\gamma_{k}\\\
\gamma_{-k}^{\dagger}\end{array}\right),$ (25)

where $\tan\theta_{k}=\Delta_{k}/\epsilon_{k}$, we obtain the diagonalized
Hamiltonian

$\mathcal{H}_{k}=\lambda_{k}(\gamma_{k}^{\dagger}\gamma_{k}-\gamma_{-k}\gamma_{-k}^{\dagger}),$
(26)

where $\lambda_{k}=\sqrt{\epsilon_{k}^{2}+\Delta_{k}^{2}}$.

Using the Anderson pseudo-spins:

$\begin{split}\tau_{k}^{x}&=\gamma_{-k}\gamma_{k}+\gamma_{k}^{\dagger}\gamma_{-k}^{\dagger},\\\
\tau_{k}^{y}&=\text{i}(\gamma_{-k}\gamma_{k}-\gamma_{k}^{\dagger}\gamma_{-k}^{\dagger}),\\\
\tau_{k}^{z}&=\gamma_{k}^{\dagger}\gamma_{k}-\gamma_{-k}\gamma_{-k}^{\dagger},\end{split}$
(27)

magnetization $M^{x}$ is given by

$\begin{split}M^{x}&=\frac{1}{2}\sum_{i=1}^{L}\sigma_{i}^{x}\\\
&=\sum_{0<k<\pi}\left(\cos\theta_{k}\tau_{k}^{z}+\sin\theta_{k}\tau_{k}^{y}\right),\\\
&\equiv\sum_{0<k<\pi}m_{k}^{x},\end{split}$ (28)

Note that we drop the contribution from $k=0$ and $k=\pi$ in the odd sector,
which is a constant of the motion and does not contribute to the 2DCS. In the
Heisenberg picture, $m_{k}^{x}$ becomes:

$m_{k}^{x}(t)=\cos\theta_{k}\tau_{k}^{z}+\sin\theta_{k}(\tau_{k}^{y}\cos(2\lambda_{k}t)+\tau_{k}^{x}\sin(2\lambda_{k}t)).$
(29)

Using the fact that each $m_{k}^{x}$ commutes with each other, the third-order
response is given by

$\begin{split}\chi^{(3)}_{xxxx}(t_{3},t_{2},t_{1})&=\frac{1}{L}\sum_{0<k<\pi}\left<[[[m_{k}^{x}(t_{3}+t_{2}+t_{1}),m_{k}^{x}(t_{2}+t_{1})],m_{k}^{x}(t_{1})],m_{k}^{x}(0)]\right>\\\
&=\frac{1}{L}\sum_{0<k<\pi}\left[A_{k}^{(1)}+A_{k}^{(2)}+A_{k}^{(3)}+A_{k}^{(4)}\right],\\\
\end{split}$ (30)

where

$\begin{split}A_{k}^{(1)}&=8\sin^{2}\theta_{k}\cos^{2}\theta_{k}\sin(2\lambda_{k}(t_{3}+t_{2}+t_{1}))\left<\tau_{k}^{z}\right>,\\\
A_{k}^{(2)}&=-8\sin^{2}\theta_{k}\cos^{2}\theta_{k}\sin(2\lambda_{k}(t_{2}+t_{1}))\left<\tau_{k}^{z}\right>,\\\
A_{k}^{(3)}&=4\sin^{4}\theta_{k}\sin(2\lambda_{k}(t_{3}+t_{1}))\left<\tau_{k}^{z}\right>,\\\
A_{k}^{(4)}&=4\sin^{4}\theta_{k}\sin(2\lambda_{k}(t_{3}-t_{1}))\left<\tau_{k}^{z}\right>,\end{split}$
(31)

Here we utilized the fact that
$\left<\tau_{k}^{x}\right>=\left<\tau_{k}^{y}\right>=0$. Zero temperature
limit Eq. (12) is obtained by $\left<\tau_{k}^{z}\right>=-1$, and consider the
even sector only, i.e., we take a sum for $k=(2n-1)\pi/L$ with
$n=1,\cdots,L/2$.

Finite temperature 2DCS is obtained by taking the thermal average of
$\left<\tau_{k}^{z}\right>$. The distribution function is given by

$Z=\text{Tr}\left(e^{-\beta\mathcal{H}_{\text{spin}}}\right)=\sum_{p}\text{Tr}\left(P_{p}e^{-\beta\mathcal{H}_{p}}\right)\equiv\sum_{p}Z_{p},$
(32)

where $P_{p}$ is the projection operator to the parity $p$ sector defined as

$P_{p}=\frac{1}{2}\left(1+(-1)^{p}e^{-\text{i}\pi N}\right).$ (33)

$Z_{0}$ can be computed as:

$Z_{0}=\frac{e^{-\beta
E_{0}}}{2}\left[\prod_{k}(1+e^{-\beta\lambda_{k}})+\prod_{k}(1-e^{-\beta\lambda_{k}})\right],$
(34)

where $E_{0}=-\sum_{k>0}\lambda_{k}$ is the ground state energy of the parity
even sector, and the product over $k$ runs the entire Brillouin zone.
Similarly, $Z_{1}$ can be computed as:

$Z_{1}=\frac{e^{-\beta
E_{1}}}{2}\left[\prod_{k}(1+e^{-\beta\lambda_{k}})-\prod_{k}(1-e^{-\beta\lambda_{k}})\right],$
(35)

Where $E_{1}=-\sum_{0<k<\pi}\lambda_{k}+2h_{x}$ is the energy of the state
obtained by annihilating the $k=0$ fermion from the odd-sector ground state.
For $k=0$ and $\pi$ we take $\lambda_{0}=-2J-2h_{x}$ and
$\lambda_{\pi}=2J-2h_{x}$, respectively. Finally, we obtain

$\begin{split}&\left<\tau_{k}^{z}\right>=-\frac{Z_{p;\bar{k}}\left(1-e^{-2\beta\lambda_{k}}\right)}{Z_{0}+Z_{1}}\end{split}$
(36)

where $p$ is the parity that is consistent with $k$, and $Z_{p;\bar{k}}$ is
obtained by omitting the contribution from $\pm k$ in the product over $k$.

##  Appendix D Phase untwisted 2DCS

In the main text, we discussed how the positive time constraint leads to a
distorted signal known as phase twisting, which can even eliminate the NR
signal at $\omega_{t}=\omega_{\tau}$. We demonstrate that this issue can be
resolved by incorporating negative $\tau$ values, which is feasible in
analytical calculations and ED simulations. Figure 12 presents the 2DCS
results obtained by including negative $\tau$. It is noteworthy that even for
a substantially large system size ($L=500$), the NR signal is observable.
Furthermore, the data for $L=24$ now appears as a discretized version of the
$L=500$ data. It is important to note that NR and R signals are not symmetric,
as predicted by Eq. 12. Therefore, phase untwisting cannot be simply achieved
by symmetrizing the experimental data around about $\omega_{\tau}$ [28].

Figure 12: Phase untwisted 2DCS obtained by including negative $\tau$ in the
calculation. The parameters $J=0.7,h_{x}=0.3,h_{z}=0.06$, and the system size
$L=500$ (a) and $L=24$ (b) are used.

##  Appendix E Two-kink model

Here, we provide the detailed calculation of 2DCS within the TK model
following the discussion in Refs. [44, 38, 29].
$\ket{j,l}\equiv\ket{\cdots\uparrow\uparrow\downarrow_{j}\cdots\downarrow_{(j+l-1)}\uparrow\uparrow\cdots}$
is the state with two kinks at $j$ and $j+l-1$. The projection to the subspace
spanned by $\ket{j,l}$ gives TK Hamiltonian which acts as:

$\begin{split}\mathcal{H}_{\text{TK}}\ket{j,l}&=(4J+2h_{z})\ket{j,l}\\\
&-h_{x}\left(\ket{j,l+1}+\ket{j,l-1}\right.\\\
&+\left.\ket{j+1,l-1}+\ket{j-1,l+1}\right).\end{split}$ (37)

Next, we introduce the momentum basis
$\ket{p,l}=(1/\sqrt{L})\sum_{j}e^{\text{i}pj}\ket{j,l}$. Here, normalization
factor $1/\sqrt{L}$ is introduced to make the qualitative comparison with ED
possible. As far as the zero-temperature 2DCS is concerned, we can restrict
ourselves to the zero-momentum sector. Using $\ket{l}\equiv\ket{p=0,l}$, we
can write the TK Hamiltonian as

$\begin{split}\mathcal{H}_{\text{TK}}\ket{l}&=(4J+2h_{z})\ket{l}\\\
&-2h_{x}(\ket{l+1}+\ket{l-1}).\end{split}$ (38)

Introducing the $n$-th eigen state
$\ket{\Phi_{n}}=\sum_{l}\psi_{n}(l)\ket{l}/\sqrt{\sum_{l=1}^{\infty}\left|\psi_{n}(l)\right|^{2}}$
and its eigen energy $E_{n}$ the Schrödinger equation becomes

$\mathcal{H}_{\text{TK}}\ket{\Psi_{n}}=E_{n}\ket{\Psi_{n}}.$ (39)

Using dimensionless paramters $\lambda_{n}=(E_{n}-4J)/(4h_{x})$ and
$\mu=h_{z}/(2h_{x})$, we obtain the following recursion relation:

$\begin{split}&(-\lambda_{n}+\mu l)\psi_{n}(l)\\\
&-(\psi_{n}(l+1)+\psi_{n}(l-1))/2=0,\end{split}$ (40)

with the boundary conditions $\psi_{n}(0)=\psi_{n}(+\infty)=0$. Eigen values
$\lambda_{n}$ can be obtained using the zeros $\nu_{n}$ of the Bessel function
$J_{\nu}(1/\mu)$ as a function of its order $\nu$ as

$\lambda_{n}=-\mu\nu_{n},$ (41)

Now, we consider processes where each intermediate state is a two-kink state.
For example, $\chi^{(3)}_{xxxx}$ is given by:

$\begin{split}&\chi^{(3)}_{xxxx}(t_{3},t_{2},t_{1})|_{\ket{Q}=\ket{\Phi_{n}}}\\\
&=\frac{2}{L}\sum_{nml}\sum_{p=1}^{4}\text{Im}\Big{[}\braket{0}{M^{x}}{\Phi_{n}}\braket{\Phi_{n}}{M^{x}}{\Phi_{m}}\\\
&\times\braket{\Phi_{m}}{M^{x}}{\Phi_{l}}\braket{\Phi_{l}}{M^{x}}{0}C_{p}(t_{3},t_{2},t_{1})\Big{]}\\\
&\equiv\sum_{nml}\sum_{p=1}^{4}\text{Im}\Big{[}A_{n,m,l}C_{p}(t_{3},t_{2},t_{1})\Big{]}\end{split}$
(42)

where $R_{p}$ defined in Eq. (17) is decomposed into the matrix elements
$A_{n,m,l}$ and the time-dependent part $C_{p}(t_{3},t_{2},t_{1})$. The factor
$(2/L)$ is incorporated into $A_{n,m,l}$. Here, the ferromagnetic ground state
$\ket{0}$ is defined as $\prod_{j=1}^{L}\ket{\uparrow}_{j}$, with an energy
$E_{0}=0$. Following the discussion outlined in Ref. [29], we obtain

$\braket{\Phi_{n}}{M^{x}}{0}=\frac{\sqrt{L}\psi_{n}(1)}{\sqrt{\sum_{l=1}^{\infty}\left|\psi_{n}(l)\right|^{2}}}$
(43)

and

$\begin{split}&\braket{\Phi_{n}}{M^{x}}{\Phi_{m}}\\\
&=-\frac{\psi_{n}(1)\psi_{m}(1)}{\sqrt{\sum_{l=1}^{\infty}\left|\psi_{n}(l)\right|^{2}}\sqrt{\sum_{l=1}^{\infty}\left|\psi_{m}(l)\right|^{2}}}\\\
&\times\left(\frac{1}{\lambda_{n}-\lambda_{m}-\mu}+\frac{1}{\lambda_{m}-\lambda_{n}-\mu}\right).\end{split}$
(44)

Therefore $A_{n,m,l}$ are obtained as

$\begin{split}&A_{n,m,l}=2\left(\frac{1}{\lambda_{n}-\lambda_{m}-\mu}+\frac{1}{\lambda_{m}-\lambda_{n}-\mu}\right)\\\
&\times\left(\frac{1}{\lambda_{m}-\lambda_{l}-\mu}+\frac{1}{\lambda_{l}-\lambda_{m}-\mu}\right)I_{n}I_{m}I_{l}\end{split}$
(45)

where the relative intensity of the $n$-th mode is calculated as:

$\begin{split}I_{n}&=\frac{\left|\psi_{n}(1)\right|^{2}}{\sum_{l=1}^{\infty}\left|\psi_{n}(l)\right|^{2}}\\\
&=2\mu\left.\left\\{\frac{\partial}{\partial\nu}\left[\frac{J_{\nu}(1/\mu)}{J_{\nu+1}(1/\mu)}\right]\right\\}^{-1}\right|_{\nu\rightarrow\nu_{n}},\end{split}$
(46)

which does not depend on the specific choice of $\psi_{n}(1)$. Note that
$\chi^{(3)}_{xxxx}|_{\ket{Q}=\ket{\Phi_{n}}}$ does not depend on $L$.

Figure 13: Separate calculation of the contributions from the processes
involving $\ket{Q}=\ket{\Phi_{n}}$, i.e., two-kink states, and
$\ket{Q}=\ket{0}$ to $\chi^{(3;1,2)}_{xxxx}$ within two-kink approximation.
The parameters $J=0.7,h_{x}=0.3,h_{z}=0.06$ (a) and $h_{z}=0.4$ (b) are used.

As can be seen in Fig. 13, the signal obtained by only considering two-kink
states as intermediate states does not reproduce the many of features obtained
by ED. In particular, the PP in $\chi^{(3;1)}_{xxxx}$ and R in
$\chi^{(3;2)}_{xxxx}$ are entirely missing. To obtain a non-rephasing signal,
we further consider the process that involves $\ket{Q}=\ket{0}$. Even though,
as discussed before, the contribution from such process is expected to be
mostly canceled out by the processes that involve four-kink states, we found
that the inclusion of this process gives a better agreement about the location
of the peaks with ED data.

For example, using Eq. (43), $\chi^{(3)}_{xxxx}(t,0,\tau)|_{\ket{Q}=\ket{0}}$
can be obtained as

$\begin{split}&\chi^{(3)}_{xxxx}(t,0,\tau)|_{\ket{Q}=\ket{0}}=2L\sum_{n,l}I_{n}I_{l}\\\
&\times\text{Im}\left[e^{-\text{i}E_{n}t}e^{-\text{i}E_{n}\tau}+2e^{\text{i}E_{l}\tau}e^{-\text{i}E_{n}t}+e^{-\text{i}E_{n}\tau}e^{-\text{i}E_{l}t}\right]\end{split}$
(47)

The second term of the left-hand side gives the signals at the fourth
quadrant, i.e. $(\omega_{t},\omega_{\tau})=(E_{n},-E_{l})$ with intensity
$2LI_{n}I_{l}$ (Fig. 13), proportional to $L$. This extensive terms should be
canceled out by the processes that involve four-kink states, making the
remaining signal $L$ independent.

##  Appendix F Separate calculation of each $R_{a}$

Figure 14: Separate calculation of each $R_{p}$ using ED (a, b) and two-kink
model (c). When $t_{1}=\tau,t_{2}=0,t_{3}=t$, the relation $R_{2}=R_{3}$
holds. Parameters $J=0.7,h_{x}=0.3$ are used.

In the main text, we numerically apply a magnetic field in a manner similar to
the MPS study in Ref. [30]. Alternatively, we can directly calculate the
imaginary part of $R_{a}$ using ED, which is more similar to the technique
used in Ref. [29]. Even though each $R_{a}$ is not accessible by experiment,
it can be useful to understand which processes contribute to each signal. Here
we provide an example in Fig. 14 where we calculate the imaginary part of each
$R_{a}$ using ED and the two-kink approximation.

Figure 14(a) shows the case of $h_{z}=0.0$, i.e. the exactly solvable case. In
addition to the peaks in $\chi^{(3;2)}_{xxxx}$, we observed additional signals
at $\omega_{t}=0$ and $\omega_{\tau}=0$. These contributions come from the
process involving the transition from $\ket{0}$ to $\ket{0}$ under the action
of $M_{x}$. $m^{x}_{00}$ is finite because of the finite value of $h_{x}$.
These contributions cancel each other out, giving zero intensity in
$\chi^{(3;2)}_{xxxx}$. Furthermore, in addition to the anti-diagonal rephasing
signal, we found off-diagonal rephasing signals of opposite sign in $R_{1}$
and $R_{2}$. This can be understood as follows: The rephasing signals in
$R_{2}$ come from the process where $\ket{Q}=\ket{0}$. Similar to what we saw
in the two-kink calculation, this results in a grid of peaks at
$(\omega_{t},\omega_{\tau})=(E_{k_{1}},-E_{k_{2}})$, where pairs of spinons
with momentum $(k_{1},-k_{1})$ and $(k_{2},-k_{2})$ are created in the
intermediate stage. On the other hand, the rephasing signals in $R_{1}$ come
from the process involving $\ket{Q}=\ket{k_{1},-k_{1},k_{2},-k_{2}}$, i.e.
four-kink states. Since two pairs of kinks are non-interacting, we can freely
create one on top of the other. An exception is when $k_{1}=k_{2}$, in which
case the creation of the second pair of kinks is forbidden due to the
fermionic nature of the spinons, which explains the absence of peaks at
$\omega_{t}=-\omega_{\tau}$ in $R_{1}$. After taking the sum of
$R_{1}+2R_{2}+R_{4}$, these off-diagonal parts of the rephasing signals are
exactly canceled out, and we obtain only a anti-diagonal rephasing signal,
where no cancellation occurs due to the absence of the rephasing signal in
$R_{1}$ at $\omega_{t}=-\omega_{\tau}$.

Figure 14(b) shows the case of $h_{z}=0.06$. Again we see a rephasing signal
in $R_{1}$ and $R_{2}$. The interacting nature of the confined spinons makes
the positions and intensities of the peaks different from each other, and the
resulting $\chi^{(3;2)}_{xxxx}$ signal has finite off-diagonal peaks. In the
two-kink approximation [Fig. 14(c)] we do not see any rephasing signals in
$R_{1}$, which is consistent with the fact that we have not included four-kink
states in the calculation. We also do not see a rather strong peak at
$\omega_{t}=0$ and $\omega_{\tau}=0$, again consistent with the fact that we
did not include the process involving the transition from $\ket{0}$ to
$\ket{0}$ under the action of $M_{x}$.

## References

  * Kampfrath  _et al._ [2011] T. Kampfrath, A. Sell, G. Klatt, A. Pashkin, S. Mährlein, T. Dekorsy, M. Wolf, M. Fiebig, A. Leitenstorfer, and R. Huber, Coherent terahertz control of antiferromagnetic spin waves, Nature Photonics 5, 31 (2011).
  * Morris  _et al._ [2014] C. M. Morris, R. Valdés Aguilar, A. Ghosh, S. M. Koohpayeh, J. Krizan, R. J. Cava, O. Tchernyshyov, T. M. McQueen, and N. P. Armitage, Hierarchy of Bound States in the One-Dimensional Ferromagnetic Ising Chain ${\mathrm{CoNb}}_{2}{\mathrm{O}}_{6}$ Investigated by High-Resolution Time-Domain Terahertz Spectroscopy, Phys. Rev. Lett. 112, 137403 (2014).
  * Zhang  _et al._ [2018] X. Zhang, F. Mahmood, M. Daum, Z. Dun, J. A. M. Paddison, N. J. Laurita, T. Hong, H. Zhou, N. P. Armitage, and M. Mourigal, Hierarchy of Exchange Interactions in the Triangular-Lattice Spin Liquid ${\mathrm{YbMgGaO}}_{4}$, Phys. Rev. X 8, 031001 (2018).
  * Bera  _et al._ [2021] A. Bera, S. Bera, S. Kalimuddin, S. Gayen, M. Kundu, B. Das, and M. Mondal, Review of recent progress on THz spectroscopy of quantum materials: superconductors, magnetic and topological materials, The European Physical Journal Special Topics 230, 4113 (2021).
  * Han  _et al._ [2012] T.-H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. Rodriguez-Rivera, C. Broholm, and Y. S. Lee, Fractionalized excitations in the spin-liquid state of a kagome-lattice antiferromagnet, Nature 492, 406 (2012).
  * Mourigal  _et al._ [2013] M. Mourigal, M. Enderle, A. Klöpperpieper, J.-S. Caux, A. Stunault, and H. M. Rønnow, Fractional spinon excitations in the quantum Heisenberg antiferromagnetic chain, Nature Physics 9, 435 (2013).
  * Wen  _et al._ [2019] J. Wen, S.-L. Yu, S. Li, W. Yu, and J.-X. Li, Experimental identification of quantum spin liquids, npj Quantum Materials 4, 12 (2019).
  * Banerjee  _et al._ [2016] A. Banerjee, C. A. Bridges, J.-Q. Yan, A. A. Aczel, L. Li, M. B. Stone, G. E. Granroth, M. D. Lumsden, Y. Yiu, J. Knolle, S. Bhattacharjee, D. L. Kovrizhin, R. Moessner, D. A. Tennant, D. G. Mandrus, and S. E. Nagler, Proximate Kitaev quantum spin liquid behaviour in a honeycomb magnet, Nature Materials 15, 733 (2016).
  * Zhu  _et al._ [2017] Z. Zhu, P. A. Maksimov, S. R. White, and A. L. Chernyshev, Disorder-Induced Mimicry of a Spin Liquid in ${\mathrm{YbMgGaO}}_{4}$, Phys. Rev. Lett. 119, 157201 (2017).
  * Shen  _et al._ [2016] Y. Shen, Y.-D. Li, H. Wo, Y. Li, S. Shen, B. Pan, Q. Wang, H. C. Walker, P. Steffens, M. Boehm, Y. Hao, D. L. Quintero-Castro, L. W. Harriger, M. D. Frontzek, L. Hao, S. Meng, Q. Zhang, G. Chen, and J. Zhao, Evidence for a spinon Fermi surface in a triangular-lattice quantum-spin-liquid candidate, Nature 540, 559 (2016).
  * Kimchi  _et al._ [2018] I. Kimchi, A. Nahum, and T. Senthil, Valence Bonds in Random Quantum Magnets: Theory and Application to ${\mathrm{YbMgGaO}}_{4}$, Phys. Rev. X 8, 031028 (2018).
  * Ma  _et al._ [2018] Z. Ma, J. Wang, Z.-Y. Dong, J. Zhang, S. Li, S.-H. Zheng, Y. Yu, W. Wang, L. Che, K. Ran, S. Bao, Z. Cai, P. Čermák, A. Schneidewind, S. Yano, J. S. Gardner, X. Lu, S.-L. Yu, J.-M. Liu, S. Li, J.-X. Li, and J. Wen, Spin-Glass Ground State in a Triangular-Lattice Compound ${\mathrm{YbZnGaO}}_{4}$, Phys. Rev. Lett. 120, 087201 (2018).
  * Lu  _et al._ [2017] J. Lu, X. Li, H. Y. Hwang, B. K. Ofori-Okai, T. Kurihara, T. Suemoto, and K. A. Nelson, Coherent Two-Dimensional Terahertz Magnetic Resonance Spectroscopy of Collective Spin Waves, Phys. Rev. Lett. 118, 207204 (2017).
  * Zhang  _et al._ [2023] Z. Zhang, F. Y. Gao, J. B. Curtis, Z.-J. Liu, Y.-C. Chien, A. von Hoegen, T. Kurihara, T. Suemoto, P. Narang, E. Baldini, and K. A. Nelson, Three-wave mixing of anharmonically coupled magnons (2023), arXiv:2301.12555 .
  * Zhang  _et al._ [2024] Z. Zhang, F. Y. Gao, Y.-C. Chien, Z.-J. Liu, J. B. Curtis, E. R. Sung, X. Ma, W. Ren, S. Cao, P. Narang, A. von Hoegen, E. Baldini, and K. A. Nelson, Terahertz-field-driven magnon upconversion in an antiferromagnet, Nature Physics (2024).
  * Khalil  _et al._ [2003] M. Khalil, N. Demirdöven, and A. Tokmakoff, Coherent 2D IR spectroscopy: Molecular structure and dynamics in solution, The Journal of Physical Chemistry A 107, 5258 (2003).
  * Cundiff  _et al._ [2012] S. T. Cundiff, A. D. Bristow, M. Siemens, H. Li, G. Moody, D. Karaiskaj, X. Dai, and T. Zhang, Optical 2-D Fourier Transform Spectroscopy of Excitons in Semiconductor Nanostructures, IEEE Journal of Selected Topics in Quantum Electronics 18, 318 (2012).
  * Fava  _et al._ [2023] M. Fava, S. Gopalakrishnan, R. Vasseur, F. H. L. Essler, and S. A. Parameswaran, Divergent nonlinear response from quasiparticle interactions (2023), arXiv:2208.09490 .
  * Wan and Armitage [2019] Y. Wan and N. P. Armitage, Resolving Continua of Fractional Excitations by Spinon Echo in THz 2D Coherent Spectroscopy, Phys. Rev. Lett. 122, 257401 (2019).
  * Li  _et al._ [2021] Z.-L. Li, M. Oshikawa, and Y. Wan, Photon Echo from Lensing of Fractional Excitations in Tomonaga-Luttinger Spin Liquid, Phys. Rev. X 11, 031035 (2021).
  * Potts  _et al._ [2023] M. Potts, R. Moessner, and O. Benton, Exploiting polarization dependence in two dimensional coherent spectroscopy: examples of Ce2Zr2O7 and Nd2Zr2O7 (2023), arXiv:2311.03200 .
  * Choi  _et al._ [2020] W. Choi, K. H. Lee, and Y. B. Kim, Theory of Two-Dimensional Nonlinear Spectroscopy for the Kitaev Spin Liquid, Phys. Rev. Lett. 124, 117205 (2020).
  * McGinley  _et al._ [2022] M. McGinley, M. Fava, and S. A. Parameswaran, Signatures of fractional statistics in nonlinear pump-probe spectroscopy (2022), arXiv:2210.16249 .
  * Krupnitska and Brenig [2023] O. Krupnitska and W. Brenig, Finite-temperature second harmonic generation in Kitaev magnets, Phys. Rev. B 108, 075120 (2023).
  * Negahdari and Langari [2023] M. K. Negahdari and A. Langari, Nonlinear response of the Kitaev honeycomb lattice model in a weak magnetic field, Phys. Rev. B 107, 134404 (2023).
  * Qiang  _et al._ [2023] Y. Qiang, V. L. Quito, T. V. Trevisan, and P. P. Orth, Probing Majorana wavefunctions in Kitaev honeycomb spin liquids with second-order two-dimensional spectroscopy (2023), arXiv:2301.11243 .
  * Fava  _et al._ [2020] M. Fava, R. Coldea, and S. A. Parameswaran, Glide symmetry breaking and Ising criticality in the quasi-1D magnet CoNb2O6, Proceedings of the National Academy of Sciences 117, 25219 (2020).
  * Hart and Nandkishore [2023] O. Hart and R. Nandkishore, Extracting spinon self-energies from two-dimensional coherent spectroscopy, Phys. Rev. B 107, 205143 (2023).
  * Sim  _et al._ [2023a] G. Sim, J. Knolle, and F. Pollmann, Nonlinear spectroscopy of bound states in perturbed Ising spin chains, Phys. Rev. B 107, L100404 (2023a).
  * Gao  _et al._ [2023] Q. Gao, Y. Liu, H. Liao, and Y. Wan, Two-dimensional coherent spectrum of interacting spinons from matrix product states, Phys. Rev. B 107, 165121 (2023).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011), january 2011 Special Issue.
  * Mukamel [1995] S. Mukamel,  _Principles of Nonlinear Optical Spectroscopy_ , Oxford series in optical and imaging sciences (Oxford University Press, 1995).
  * Nandkishore  _et al._ [2021] R. M. Nandkishore, W. Choi, and Y. B. Kim, Spectroscopic fingerprints of gapped quantum spin liquids, both conventional and fractonic, Phys. Rev. Res. 3, 013254 (2021).
  * Kuehn  _et al._ [2011] W. Kuehn, K. Reimann, M. Woerner, T. Elsaesser, and R. Hey, Two-Dimensional Terahertz Correlation Spectra of Electronic Excitations in Semiconductor Quantum Wells, The Journal of Physical Chemistry B 115, 5448 (2011).
  * Sidje [1998] R. B. Sidje, Expokit: a software package for computing matrix exponentials, ACM Trans. Math. Softw. 24, 130–156 (1998).
  * Sugiura and Shimizu [2012] S. Sugiura and A. Shimizu, Thermal Pure Quantum States at Finite Temperature, Phys. Rev. Lett. 108, 240401 (2012).
  * Sugiura and Shimizu [2013] S. Sugiura and A. Shimizu, Canonical Thermal Pure Quantum State, Phys. Rev. Lett. 111, 010401 (2013).
  * Rutkevich [2010] S. B. Rutkevich, On the weak confinement of kinks in the one-dimensional quantum ferromagnet CoNb2O6, Journal of Statistical Mechanics: Theory and Experiment 2010, P07015 (2010).
  * Coldea  _et al._ [2010] R. Coldea, D. A. Tennant, E. M. Wheeler, E. Wawrzynska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl, and K. Kiefer, Quantum Criticality in an Ising Chain: Experimental Evidence for Emergent E8 Symmetry, Science 327, 177 (2010).
  * Zhang and Song [2015] G. Zhang and Z. Song, Topological Characterization of Extended Quantum Ising Models, Phys. Rev. Lett. 115, 177204 (2015).
  * You  _et al._ [2016] W.-L. You, Y.-C. Qiu, and A. M. Oleś, Quantum phase transitions in a generalized compass chain with three-site interactions, Phys. Rev. B 93, 214417 (2016).
  * Morris  _et al._ [2021] C. M. Morris, N. Desai, J. Viirok, D. Hüvonen, U. Nagel, T. Rõõm, J. W. Krizan, R. J. Cava, T. M. McQueen, S. M. Koohpayeh, R. K. Kaul, and N. P. Armitage, Duality and domain wall dynamics in a twisted Kitaev chain, Nature Physics 17, 832 (2021).
  * Sim  _et al._ [2023b] G. Sim, F. Pollmann, and J. Knolle, Microscopic details of two-dimensional spectroscopy of one-dimensional quantum Ising magnets, Phys. Rev. B 108, 134423 (2023b).
  * Rutkevich [2008] S. Rutkevich, Energy spectrum of bound-spinons in the quantum Ising spin-chain ferromagnet, Journal of Statistical Physics 131, 917 (2008).
  * Iitaka and Ebisuzaki [2003] T. Iitaka and T. Ebisuzaki, Algorithm for Linear Response Functions at Finite Temperatures: Application to ESR Spectrum of $s=\frac{1}{2}$ Antiferromagnet Cu Benzoate, Phys. Rev. Lett. 90, 047203 (2003).