# Photon echo from lensing of fractional excitations in Tomonaga-Luttinger
spin liquid

Zi-Long Li  Institute of Physics, Chinese Academy of Sciences, Beijing 100190,
China  University of Chinese Academy of Sciences, Beijing 100049, China
Masaki Oshikawa  Institute for Solid State Physics, University of Tokyo,
Kashiwa, Chiba 277-8581, Japan  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba
277-8583, Japan  Yuan Wan<EMAIL_ADDRESS>Institute of Physics, Chinese
Academy of Sciences, Beijing 100190, China  University of Chinese Academy of
Sciences, Beijing 100049, China  Songshan Lake Materials Laboratory, Dongguan,
Guangdong 523808, China

###### Abstract

We study theoretically the nonlinear optical response of Tomonaga-Luttinger
spin liquid in the context of terahertz (THz) two-dimensional coherent
spectroscopy (2DCS). Using the gapless phase of the XXZ-type spin chain as an
example, we show that its third-order nonlinear magnetic susceptibilities
$\chi^{(3)}_{+--+}$ and $\chi^{(3)}_{-++-}$ exhibit photon echo, where $\pm$
refers to the left/right-hand circular polarization with respect to the
$S^{z}$ axis. The photon echo arises from a “lensing” phenomenon in which the
wave packets of fractional excitations move apart and then come back toward
each other, amounting to a refocusing of the excitations’ world lines.
Renormalization group irrelevant corrections to the fixed point Hamiltonian
result in dispersion and/or damping of the wave packets, which can be
sensitively detected by lensing and consequently the photon echo. Our results
thus unveil the strength of THz-2DCS in probing the dynamical properties of
the collective excitations in a prototypical gapless many-body system.

##  I Introduction

Progress in condensed matter physics is intimately connected to the
development of new spectroscopic techniques. Among the many emerging
spectroscopies, _two-dimensional coherent spectroscopy_ (2DCS) stands out as a
promising tool for investigating strongly correlated systems. The 2DCS use
multiple coherent electromagnetic waves to probe the nonlinear optical
properties of a sample, thereby producing a two-dimensional spectrum that
visualizes the sample’s nonlinear response as a function of the frequencies of
probing electromagnetic waves [1, 2].

Comparing to the more familiar one-dimensional spectroscopy that probes linear
optical properties, the 2DCS reveals not only the optical excitations in a
sample but also their relationship. In the infrared and higher frequency
range, its ability to diagnose the interplay between optical excitations has
been widely leveraged by chemists to unravel the structure of complex
molecules and map out the kinetic pathways of chemical reactions [1, 2, 3].
The advent of terahertz (THz) 2DCS now puts this technique in the right energy
window to study many-body phenomena. The THz 2DCS has offered new experimental
insights into quantum wells [4], antiferromagnets [5], and electronic glasses
[6]. On the theory front, it has recently been suggested that the THz 2DCS can
resolve the spectral continua formed by optical excitations in several clean
and disordered many-body system and characterize their dynamical properties,
which would be challenging to accomplish with linear spectroscopy, if at all
[7, 8, 9, 10, 11].

A main strength of the 2DCS lies in the _photon echo_ [12] signal from the
third-order nonlinear optical response $\chi^{(3)}$. The photon echo is an
optical analogue of the spin echo in nuclear magnetic resonance (NMR) [13].
Schematically, one may observe the photon echo by exciting the system with
three successive short optical pulses and detecting the resulted $\chi^{(3)}$
response (Fig. 1(a)). Let $\tau$ be the time delay between the second and
first pulses, $t_{\mathrm{w}}$ (the waiting time) be the delay between the
third and second, and $t$ be the delay between the time of detection and the
last pulse. The photon echo signal appears as a surge of the nonlinear
response at $t\approx\tau$ in close analogy with the nuclear magnetic
resonance (NMR) spin echo (Fig. 1(b)).

Similar to its NMR cousin, the photon echo signal can diagnose dissipation in
a few-body system by effecting a time-reversal operation — the system’s
evolution during the time delay $t$ reverses the evolution occurred during the
time delay $\tau$ [1, 2]. Had the dynamics been unitary, the quantum
mechanical phase accumulated in $\tau$ would be completely removed when
$t=\tau$, resulting in a perfect rephasing. This perfect rephasing process
would produce a photon echo signal at $t=\tau$ regardless the value of $\tau$.
Deviation from the perfect rephasing is thus a direct measure of dissipation
in the few-body system. Specifically, the decay of echo signal with increasing
$\tau$ is a manifestation of the decoherence time ($T_{2}$ time), whereas the
decay of the signal as a function of $t_{\mathrm{w}}$ probes the population
time ($T_{1}$ time).

This unique ability of diagnosing dissipation makes one naturally wonder if
the photon echo could also find its success in strongly correlated many-body
systems. Although the latest theoretical inquiries have suggested interesting
applications of photon echo in gapped systems [7, 8, 10] and disordered
systems [9], less attention is paid to clean, gapless many-body systems.
Adapting this technique to a gapless many-body setting poses new theoretical
challenges. The standard framework for analyzing the photon echo uses the
language of energy levels [1, 2], which is best suited for few-body systems
with discrete energy spectra. While it is possible to analyze the zero
temperature optical response of a gapped system by truncating the Fock space
to a subspace containing finite number of excitations [14] and thereby making
use of the established framework, such truncation is not permitted in general
for gapless systems. Important questions such as the existence of photon echo
in gapless strongly-correlated systems, its underlying mechanism, and its
features require theoretical investigation.

In this work, we address these questions by studying the nonlinear optical
response of a prototypical gapless strongly-correlated system, namely the
Tomonaga-Luttinger spin liquid (henceforth “Luttinger spin liquid” for short)
[15, 16]. For concreteness, we consider the Luttinger spin liquid hosted by
the XXZ spin chain, which possesses a global $U(1)$ spin rotational symmetry
with respect to the spin $z$ axis. We consider exclusively the nonlinear
magnetic response $\perp z$ and decompose the electromagnetic wave
polarization in the left-handed ($+$) and right-handed ($-$) basis.

Using the bosonization Hamiltonian at the renormalization group (RG) fixed
point, we find that, among the six symmetry-allowed third-order magnetic
susceptibilities, $\chi^{(3)}_{+--+}$ and its complex conjugate
$\chi^{(3)}_{-++-}$ show photon echo, which appears as a peak on the $t$ axis
at $t\approx\tau$. The echo signal possesses a universal asymptotic form,
which we obtain analytically and verify numerically. Crucially, the photon
echo is “perfect” in the sense that the signal is a function of $t-\tau$
rather than $t$ and $\tau$ both, resembling the perfectly rephasing photon
echo in a few-body system. This implies that the signal measured at a given
value of $t-\tau$ is independent of $\tau$. Moreover, the echo signal depends
weakly on the waiting time $t_{\mathrm{w}}$, and saturates when
$t_{\mathrm{w}}\to\infty$.

This perfect photon echo, although resembles the one due to the prefect
rephasing process in few-body systems, comes as a pleasant surprise — The
rephasing process is understood as a result of the optical transitions between
discrete energy levels. Here, the rephasing picture does not directly apply as
the present system has a continuous energy spectrum. Instead, we trace its
origin back to a unique “lensing” phenomenon of fractional excitations in
Luttinger spin liquids (Fig. 6c): The first THz pulse creates two wave packets
of fractional excitations with opposite chirality [17]. The second pulse
converts the left-moving wave packet into a right-moving one, whereas the
third pulse converts the right-moving wave packet into a left-moving one.
These two wave packets then meet each other later at $t=\tau$, thereby
producing an echo. For the fixed-point Hamiltonian, the phonon modes are exact
eigenstates of the Hamiltonian, and their dispersion relation is linear.
Consequently, the wave packets of fractional excitations can propagate through
the system indefinitely without decay or dispersion. This naturally explains
the perfect photon echo, namely the echo signal does not decay as either
$\tau$ or $t_{\mathrm{w}}$ increases.

Our lensing picture immediately suggests that the photon echo is a sensitive
diagnostic to the RG-irrelevant perturbations to the fixed point Hamiltonian.
These RG-irrelevant corrections give only minor corrections to the most
physical quantities at low temperature, and consequently they are often
elusive to experimental probes. In the 2DCS, these corrections manifest
themselves in the decay of the photon echo thanks to lensing.

In the XXZ spin chain, the RG-irrelevant corrections include the umklapp terms
and higher order gradient terms [15, 16]. The umklapp terms give rise to the
dissipation of phonon modes and therefore the decay of the wave packets at
finite temperature. As a result, the lensing becomes unattainable when the
pulse delay $\tau$ or $t_{\mathrm{w}}$ exceeds the lifetime of the wave
packets. This is manifest as the decay of the echo signal as a function of
$\tau$ or $t_{\mathrm{w}}$, which is analogous to the dissipation-induced
decay of the phonon echo in few-body systems mentioned above.

The higher order gradient terms, on the other hand, may result in the decay of
the photon echo through a different mechanism. By adding these terms, one may
introduce a small curvature to the dispersion relation of the phonon modes
while keeping them as the exact eigenstates. The curvature results in the
dispersion of the wave packets. We expect the lensing to be ineffective beyond
a time scale $\tau_{\mathrm{disp}}$, at which point the width of the wave
packet is comparable with the correlation length.

This dispersion-induced decay of the photon echo is distinct from the
dissipation-induced decay and finds no immediate analogue in few-body systems.
We study this decay mechanism on a toy model, namely the harmonic chain, which
is a lattice discretization of the fixed point Hamiltonian. The photon echo,
when measured at $t=\tau$, decays as a stretched exponential
$\sim\exp(-C(\tau/\tau_{\mathrm{disp}})^{1/2})$, where $C$ is a numerical
constant. Meanwhile, the echo signal shows weak dependence on the waiting time
$t_{\mathrm{w}}$ and saturates when $t_{\mathrm{w}}\to\infty$. We attribute
the lack of $t_{\mathrm{w}}$ dependence to the absence of thermalization in
the toy model — The decay of the photon echo as a function of $t_{\mathrm{w}}$
reflects the population time. Since the population of the phonon modes cannot
relax, its population time is effectively infinity.

To summarize, our analysis shows that the $\chi^{(3)}$ photon echo from the
Luttinger spin liquid is a sensitive diagnostic of the RG-irrelevant
perturbations to the fixed point Hamiltonian, which are difficult to detect
with linear optical spectroscopy. It also uncovers a dispersion-induced photon
echo decay mechanism unique to many-body systems. Conceptually, the lensing of
fractional excitations is a convenient picture for understanding the photon
echo in the Luttinger spin liquid. The lensing picture extends the phase
interference picture, commonly invoked for the photon echo in few-body systems
[1, 2], from the time domain to the spacetime domain.

The rest of this work is organized as follows: In Sec. II, we describe the
problem setup. We present the bosonization analysis in Sec.. III and the
lensing picture in Sec. IV. We investigate the dispersion-induced photon echo
decay in Sec. V. In Sec. VI, we point out a few interesting open problems.

##  II Setup

Figure 1: (a) The Faraday configuration considered in this work. A magnetic
field $B$ is applied in the $z$ axis. Three short electromagnetic pulses with
circular polarizations pass through the $S=1/2$ spin chain. The propagation
direction is parallel with the spin $z$ axis. The first pulse is right-handed,
whereas the second and the third are left-handed. We denote the time delay
between the first and the second pulse by $\tau$, the second and the third by
$t_{\mathrm{w}}$, and the third pulse and the time of detection by $t$. (b)
The photon echo appears as the surge of nonlinear optical response at the time
$t\approx\tau$.

We consider the $S=1/2$ XXZ spin chain:

$\displaystyle
H=\sum_{j}\frac{J_{\perp}}{2}(S^{+}_{j}S^{-}_{j+1}+h.c.)+J_{z}S^{z}_{j}S^{z}_{j+1}-BS^{z}_{j}.$
(1)

$j$ labels the lattice sites. $S^{\pm}_{j},S^{z}$ are the $S=1/2$ spin
operators. $J_{\perp}$ is the exchange constant in the spin $xy$ plane. We
shall consider both “ferromagnetic” ($J_{\perp}<0$) and “antiferromagnetic”
($J_{\perp}>0$) chains. $J_{z}$ is the exchange constant in the spin $z$ axis.
We include the Zeeman term due to an external field $B\parallel z$ (Fig. 1).
Throughout this work, we use the natural units with $\hbar=k_{B}=\mu=1$ where
$\mu$ is the magnetic moment carried by the spin. We may extend Eq. (1) by
including additional terms so long as they preserve the symmetries. Our
analysis is applicable to the Luttinger spin liquid phase of Eq. (1) and its
extensions.

The 2DCS measures a sample’s nonlinear optical response [1, 2]. The
electromagnetic wave interacts with an insulating spin system such as Eq. (1)
primarily by the Zeeman coupling. Therefore, the nonlinear optical response of
Eq. (1) is chiefly due to its leading-order nonlinear magnetic susceptibility.
Note the photon echo that arises from this nonlinear magnetic response closely
resembles the NMR spin echo in that both are induced by a sequence of magnetic
field pulses, and the former may be justifiably called spin echo as well [5].
However, different from the canonical NMR spin echo set up [13], 2DCS does not
require precise control of field pulse area. Here, we adhere to the term
“photon echo” to emphasize this difference from the NMR spin echo.

The Hamiltonian Eq. (1) possesses a $U(1)$ spin rotational symmetry with
respect to the $z$ axis. Since the total magnetization in $z$ commutes with
$H$ and does not evolve in the Heisenberg picture, it is natural to consider
the magnetic response in the $xy$ plane. Experimentally, this corresponds to
the _Faraday geometry_ where the propagation direction of the electromagnetic
wave is parallel/anti-parallel to the external field (Fig. 1a) [18].

The $U(1)$ symmetry of the Hamiltonian Eq. (1) forbids second-order in-plane
nonlinear magnetic susceptibilities. The same symmetry allows six third-order
in-plane nonlinear magnetic susceptibilities, out of which three are
independent: $\chi^{(3)}_{+-+-}$, $\chi^{(3)}_{+--+}$, and
$\chi^{(3)}_{++--}$, where $+(-)$ refers to the left (right)-handed circular
polarization of the electromagnetic wave. The other three susceptibilities,
namely $\chi^{(3)}_{-+-+}$, $\chi^{(3)}_{-++-}$, and $\chi^{(3)}_{--++}$, are
related to the former three by complex conjugation and thus offer no new
information. Throughout this work, we focus on $\chi^{(3)}_{+--+}$, which
exhibits the photon echo, and defer the discussion on the other two
susceptibilities to Sec. IV.

The THz 2DCS typically measures the nonlinear response in the time domain [4].
Following the discussion in Sec. I, we consider a three-pulse set up (Fig. 1):
three circularly polarized electromagnetic pulses arrive at the sample
successively at time $0$, $\tau$, $\tau+t_{\mathrm{w}}$, where
$\tau,t_{\mathrm{w}}>0$ are the delay between successive pulses. The signal
from the sample detected at a later time $t>0$ after the last pulse contains
contributions from both linear and nonlinear responses. Rerunning the
experiment with individual pulses, one may subtract off the linear response
and thereby isolate the nonlinear response.

The nonlinear signal is a convolution of the pulse profile and the third-order
magnetic susceptibility:

$\chi^{(3)}_{+--+}(t,t_{\mathrm{w}},\tau)=\iiint^{\infty}_{-\infty}dx_{1}dx_{2}dx_{3}\\\
\widetilde{\chi}^{(3)}_{+--+}(t,x_{1};\,t+t_{\mathrm{w}},x_{2};\,t+t_{\mathrm{w}}+\tau,x_{3}).$
(2)

Here, $\chi^{(3)}_{+--+}$ is the optical susceptibility that depends only on
time, while $\widetilde{\chi}^{(3)}_{+--+}$ is the _spacetime-dependent_
susceptibility. Note the time parametrization of the latter quantity follows
the standard convention for nonlinear susceptibility [19], whereas the former
does not. We use the symbol with and without the tilde to emphasize these
differences.

We visualize $\chi^{(3)}_{+--+}(t,t_{\mathrm{w}},\tau)$ by holding
$t_{\mathrm{w}}$ constant and scanning $\tau$ and $t$. We obtain the two-
dimensional spectra $\chi^{(3)}(\omega_{t},t_{\mathrm{w}},\omega_{\tau})$ by
performing a two-dimensional one-sided Fourier transform of
$\chi^{(3)}(t,t_{\mathrm{w}},\tau)$ over the domain $t>0$ and $\tau>0$. Note
alternative protocols for visualizing the nonlinear response exist [4, 6].
Ours is closely related to that of Refs. 5, 7.

##  III Bosonization

In this section, we compute the nonlinear response of Eq. (1) by using
bosonization. We show that the nonlinear susceptibility $\chi^{(3)}_{+--+}$
exhibits photon echo and characterize its features.

###  III.1 Bosonization essentials

The bosonization of $S=1/2$ XXZ chain is standard [15, 16]. We briefly review
the results to establish notations. Upon bosonization, the Hamiltonian at the
RG fixed point reads:

$\displaystyle
H=\frac{u}{2\pi}\int(\frac{1}{K}(\nabla\phi)^{2}+K(\nabla\theta)^{2})dx.$ (3)

Here, $u$ is the speed of sound. $K$ is the Luttinger parameter. They may be
computed from the microscopic model parameters $J_{\perp},J_{z},B$ by using
the Bethe ansätz [20, 21]. $\phi$ and $\theta$ are boson fields with the
compactification conditions:

$\displaystyle\phi\sim\phi+\pi,\quad\theta\sim\theta+2\pi.$ (4)

They obey the non-local commutation relation:

$\displaystyle[\phi(x),\theta(y)]=-i\pi\Theta(x-y),$ (5)

where $\Theta(\cdot)$ is the Heaviside step function. Note there is freedom in
choosing the commutation relation between $\phi$ and $\theta$. We discuss the
difference between the different commutation relation prescriptions and the
associated subtleties in Appendix A.

Up to a cutoff dependent prefactor, the spin operators assume the following
form:

$\displaystyle
S^{-}_{j}\approx\left\\{\begin{array}[]{cc}\exp(i\theta(x))&(J_{\perp}<0)\\\
\exp(i\theta(x))\cos(2\phi(x)-2\pi mx)&(J_{\perp}>0)\\\ \end{array}\right.,$
(8)

where $m$ is the magnetization density, and $x$ is the spatial coordinate of
site $j$. The external field $B$ is subsumed into the expression of spin
operator through $m$.

Note we have omitted in Eq. (8) the spatially staggered ($\propto(-)^{j}$)
component, which does not contribute to the optical response as the wavelength
of the THz probe is typically much larger than the lattice spacing. Although
the ferromagnetic chain ($J_{\perp}<0$) and the antiferromagnetic chain
($J_{\perp}>0$) are exactly mapped to each other by a $\pi$-rotation of spins
about the $z$-axis at every other site, i.e.
$S^{\pm}_{j}\to(-1)^{j}S^{\pm}_{j}$, such mapping also exchanges the uniform
and staggered components of the spin operator. Since we only consider the
uniform component, it is necessary to distinguish the two cases for our
purpose as shown in Eq. (8).

###  III.2 Four-point response function

The experimentally measured signal is related to the spacetime-dependent
nonlinear magnetic susceptibility $\widetilde{\chi}^{(3)}_{+--+}$ (Eq. (2)).
Its Kubo formula [22] reads:

$\widetilde{\chi}^{(3)}_{+--+}(1,2,3)=i^{3}\Theta(t_{1})\Theta(t_{2}-t_{1})\Theta(t_{3}-t_{2})\\\
\times\langle[[[S^{+}(0),S^{-}(-1)],S^{-}(-2)],S^{+}(-3)]\rangle.$ (9)

$0,1,2,3$ are shorthand notations for spacetime coordinates $(0,0)$,
$(t_{1},x_{1})$, $(t_{2},x_{2})$, and $(t_{3},x_{3})$, respectively.
$\widetilde{\chi}^{(3)}_{+--+}$ measures the system’s response at the
spacetime origin due to successive perturbations at $(-t_{3},-x_{3})$,
$(-t_{2},-x_{2})$, and $(-t_{1},-x_{1})$.

From Eq. (8), we see that the spin operators are linear combinations of vertex
operators. It will be convenient to seek a general expression for the four-
point response function with the following form:

$\displaystyle
G^{R}(1,2,3)=i^{3}\Theta(t_{1})\Theta(t_{2}-t_{1})\Theta(t_{3}-t_{2})$
$\displaystyle\times\langle[[[V_{0}(0),V_{1}(-1)],V_{2}(-2)],V_{3}(-3)]\rangle,$
(10)

where

$\displaystyle V_{j}=e^{i(m_{j}\theta+2n_{j}\phi)},$ (11)

are the vertex operators. We focus on the _local_ vertex operators, i.e. those
preserve the boson compactification conditions (4). This imposes the condition
on the coefficients

$\displaystyle m_{j},n_{j}\in\mathbb{Z}.$ (12)

We also impose the charge neutrality condition $\sum_{j}m_{j}=\sum_{j}n_{j}=0$
to ensure $G^{R}$ does not vanish in the thermodynamic limit.

Eq. (10) can be calculated by using the established technique [16]. Here we
only sketch the key steps. Using the Baker-Campbell-Hausdorff formula, we find
the commutator $[e^{iA},e^{iB}]=-2\mathrm{sinh}\frac{[A,B]}{2}\,e^{i(A+B)}$,
where $A,B$ are arbitrary linear combinations of $\phi$ and $\theta$. This
permits a straightforward evaluation of the nested commutators in Eq. (10). We
then compute the thermal average by using
$\langle\exp(iA)\rangle=\exp(-\langle A^{2}\rangle/2)$, where $A$ is an
arbitrary linear combination of $\theta$ and $\phi$. We obtain:

$G^{R}(1,2,3)=-8\Theta(t_{1})\Theta(t_{2}-t_{1})\Theta(t_{3}-t_{2})\\\
\times\sin(\alpha_{10})\sin(\alpha_{20}+\alpha_{21})\sin(\alpha_{30}+\alpha_{31}+\alpha_{32})\\\
\times C_{10}C_{20}C_{21}C_{30}C_{31}C_{32}.$ (13)

The above is the main result of this subsection. Here, $\alpha_{ij}$ and
$C_{ij}$ are defined for the vertex operator $V_{i}$ and $V_{j}$.
$\alpha_{ij}$ comes from the commutator of vertex operators:

$\alpha_{ij}=\frac{\pi}{2}(l_{i}l_{j}\mathrm{Sgn}(x^{+}_{ij})+r_{i}r_{j}\mathrm{Sgn}(x^{-}_{ij})\\\
+l_{i}r_{j}-r_{i}l_{j}).$ (14a) We have used light cone coordinates
$x^{\pm}\equiv ut\pm x$. $x^{\pm}_{ij}\equiv x^{\pm}_{i}-x^{\pm}_{j}$.
$\mathrm{Sgn}(\cdot)$ is the sign function. $l_{j}$ and $r_{j}$ are real
parameters related to $m_{j},n_{j}$ through: $\displaystyle
l_{j}=\frac{m_{j}}{2\sqrt{K}}+\sqrt{K}n_{j};\quad
r_{j}=\frac{m_{j}}{2\sqrt{K}}-\sqrt{K}n_{j}.$ (14b) Meanwhile, $\displaystyle
C_{ij}=\left|\frac{\sinh(\pi Tx^{+}_{ij}/u)}{\pi
T\epsilon/u}\right|^{l_{i}l_{j}}\left|\frac{\sinh(\pi Tx^{-}_{ij}/u)}{\pi
T\epsilon/u}\right|^{r_{i}r_{j}},$ (14c)

where $\epsilon$ is the short-distance cutoff. $T$ is the temperature.

Causality is an important property shared by all experimentally accessible
response functions. For a relativistic system described by the fixed point
Hamiltonian Eq. (3), the response vanishes if the perturbations are outside
the past light cone of the detection event. It is then natural to ask if Eq.
(13) is causal and under what conditions. In Appendix B, we show that Eq. (13)
is causal provided that the vertex operators are _local_ , i.e. Eq. (12) holds
for all $V_{j}$. It is straightforward to check that the vertex operators that
appear in the expression of $S^{\pm}(x)$ (Eq. (8)) indeed fulfill this
condition. Thus, $\widetilde{\chi}^{(3)}_{+--+}$ is causal as expected.

###  III.3 FM chain

Figure 2: (a) Nonlinear magnetic susceptibility $\chi^{(3)}_{+--+}$ of a
ferromagnetic chain as a function of $\pi T(t-\tau)$ at $\pi T\tau=20$ (blue)
and $\pi T\tau=40$ (cyan). The waiting time $\pi Tt_{\mathrm{w}}=1$. Luttinger
parameter $K=0.7$. Data for different $\tau$ are scaled by the _same_ factor
such that the maximum of $\pi T\tau=20$ data is 1. Upper right inset: the
$t>\tau$ part of the data, multiplied by $\exp(\gamma|t-\tau|)$ with
$\gamma=\pi T/(2K)$, plotted against $(\pi T(t-\tau))^{3}$. Lower left inset:
the $t<\tau$ part of the data, multiplied by $\exp(\gamma|t-\tau|)$, plotted
against $\pi T(t-\tau)$. (b) The same as (a) but for $K=1$. (c) The same as
(a) but for $K=1.5$.  Figure 3: (a) Nonlinear magnetic susceptibility
$\chi^{(3)}_{+--+}$ of a ferromagnetic chain as function of $\pi Tt$ and $\pi
T\tau$. The waiting time $\pi Tt_{\mathrm{w}}=1$. The Luttinger parameter
$K=1$. The data are real numbers and have been rescaled such that the maximum
is 1. (b)(c) The real and imaginary parts of the two-dimensional spectrum,
obtained by Fourier transforming the data of (a). Only the first and fourth
quadrants are shown. The data in the other two quadrants are obtained by
complex conjugation, i.e. the real part (b) is invariant under the
transformation $\omega_{\tau}\to-\omega_{\tau},\omega_{t}\to-\omega_{t}$,
whereas the imaginary part (c) changes sign.

In this subsection, we consider the ferromagnetic ($J_{\perp}<0$) chain.
Plugging Eq. (8) into the Kubo formula (Eq. (9)) yields:

$\widetilde{\chi}^{(3)}_{+--+}(1,2,3)\overset{(J_{\perp}<0)}{=}i^{3}\Theta(t_{1})\Theta(t_{2}-t_{1})\Theta(t_{3}-t_{2})\\\
\times\langle[[[e^{-i\theta(0)},e^{i\theta(-1)}],e^{i\theta(-2)}],e^{-i\theta(-3)}]\rangle.$
(15)

The above has the form of Eq. (10). We may read off its explicit expression
from Eq. (13) by setting $l=r=1/\sqrt{4K}$ for $\exp(i\theta)$ and
$l=r=-1/\sqrt{4K}$ for $\exp(-i\theta)$. An immediate consequence of Eq. (13)
is that $\chi^{(3)}_{+--+}$ is strictly real and independent of the
magnetization density $m$.

The next step is to find $\chi^{(3)}_{+--+}$ by integrating
$\widetilde{\chi}^{(3)}_{+--+}$ over spatial coordinates (Eq. (2)). Given the
complex structure of the integrand, the integral is unlikely to admit a
simple, closed form. We instead seek the asymptotic form of
$\chi^{(3)}_{+--+}$ valid when $t,t_{w},\tau$ are large.

To this end, we use the following approximation for the function $C_{ij}$ (Eq.
(14c)):

$\displaystyle C_{ij}(x_{ij},t_{ij})$ $\displaystyle\approx\exp\left(\frac{\pi
T}{u}(l_{i}l_{j}|x^{+}_{ij}|+r_{i}r_{j}|x^{-}_{ij}|)\right).$ (16)

We have omitted a cutoff dependent prefactor. Eq. (16) captures the
exponential falling-off/growth of $C_{ij}$ away from the light cone but
neglects the algebraic singularity near the light cone $x_{ij}=\pm ut_{ij}$.
The latter is short distance physics and should not affect the asymptotic
behavior. The validity of Eq. 16 will be further verified _a posteriori_ by
numerical integration.

With the approximation Eq. (16), the integral Eq. (2) is now elementary. After
lengthy calculations, we find the integration produces two groups of terms:
The first group of terms simply decrease as $t$ or $\tau$ increases, which we
discard as they are uninteresting for our purpose. The second group of terms
show signature of photon echo — as we fix $\tau$ and scan $t$, they exhibit a
maximum near $t\approx\tau$. Retaining these terms, we find:

$\displaystyle\chi^{(3)}_{+--+}\overset{(J_{\perp}<0)}{\sim}\left\\{\begin{array}[]{cc}(\tau-t)\,e^{-\frac{\pi
T(\tau-t)}{2K}}&(t\ll\tau)\\\ -(t-\tau)^{3}e^{-\frac{\pi
T(t-\tau)}{2K}}&(t\gg\tau)\end{array}\right..$ (19)

The above is the key result of this subsection.

Eq. (19) suggests that the photon echo is perfect, i.e. it is a function of
$t-\tau$. The time scale of the photon echo signal is set by $2K/(\pi T)$.
Moreover, it is independent of the waiting time $t_{\mathrm{w}}$ in the
asymptotic regime. As we have discussed in Sec. I, the photon echo in a few-
body system is understood as the result of the rephasing process, which in
turn builds on transitions between discrete energy levels. Since the energy
spectrum of the Luttinger spin liquid is continuous, the rephasing picture
does not apply. The physics behind the photon echo in the present system will
be discussed in detail in Sec. IV.

We test the validity of Eq. (19) by performing the integration Eq. (2)
numerically. We use a short distance cutoff $\pi\epsilon T/u=1$, which
smoothens the algebraic singularities and discontinuities that would otherwise
appear in the integrand in the limit $\epsilon\to 0^{+}$. With finite
$\epsilon$, the integrand decreases rapidly outside the light cone of $(0,0)$.
We thus limit the domain of integration to the light cone plus a small
interval of size $R$ beyond the light cone. We set $\pi TR/u=2$ with the
relative error $O(10^{-3})$.

Fig. 2 shows $\chi^{(3)}_{+--+}$ for representative Luttinger parameters
$K=0.7$, $K=1$, and $K=1.5$. $\pi Tt_{\mathrm{w}}=1$ and $\pi T\tau=20,40$.
For all cases, the nonlinear response shows a surge near $t\approx\tau$,
exhibiting the clear signature of photon echo. Note the maximum is _not_
exactly located at $t=\tau$ but fairly close to it. The data for different
values of $\tau$ overlay within numerical error when plotted as a function of
$t-\tau$. This demonstrates the photon echo is independent of $\tau$ for large
$\tau$.

Numerical integration also indicates that the value of $\chi^{(3)}_{+--+}$
measured at $t=\tau$ does not depend on $t_{\mathrm{w}}$ within numerical
error in the asymptotic regime $\pi T\tau\gg 1$, which is in agreement with
Eq. (19)

We further examine the asymptotic behavior of $\chi^{(3)}_{+--+}$ by
multiplying it with $\exp(\pi T|t-\tau|/(2K))$. Eq. (19) suggests the product
would be $\propto\tau-t$ when $t<\tau$, and $\propto-(t-\tau)^{3}$ when
$t>\tau$. The insets of Fig. 2 show that its behavior is in excellent
agreement with Eq. (19).

Having analyzed the photon echo signal for fixed value of $\tau$, we now scan
$\tau$ and present the nonlinear response as a function of both $\tau$ and
$t$. Fig. 3(a) shows $\chi^{(3)}_{+--+}$ for $K=1$ and $\pi
Tt_{\mathrm{w}}=1$. The photon echo appears as a bright feature at the
diagonal of the $(t,\tau)$ plane. This feature persists along the diagonal
direction, highlighting the fact that the photon echo is perfect.

Performing the FFT of the time domain data yields the two-dimensional spectrum
(Fig. 3(b)&(c)). In the frequency domain, the photon echo manifest itself as a
pair of highly anisotropic peaks in the second/fourth quadrants, symmetrically
distributed with respect to the origin. The peak width along the diagonal
direction of the second/fourth quadrant scales with $T$. The peak width along
the anti-diagonal direction is resolution limited — the photon echo signal in
the time domain (Fig. 3(a)) is independent of $t+\tau$ at late time, and hence
its Fourier transform $\sim\delta(\omega_{t}+\omega_{\tau})$.

To recapitulate, the $\chi^{(3)}_{+--+}$ of the ferromagnetic chain is real
and independent of the magnetization density. It exhibits perfect photon echo,
which depends on $t-\tau$ rather than $t$ and $\tau$ both.

###  III.4 AFM chain

We turn to the antiferromagnetic ($J_{\perp}>0$) case in this subsection. We
write the spin operator as:

$\displaystyle S^{-}(x)$ $\displaystyle=e^{-2\pi
imx}e^{i(\theta(x)+2\phi(x))}+e^{2\pi imx}e^{i(\theta(x)-2\phi(x))}$
$\displaystyle=e^{-2\pi imx}a(x)+e^{2\pi imx}b(x),\quad(J_{\perp}>0)$ (20)

where we have defined vertex operators $a$ and $b$ for later convenience.
Inserting the above into Eq. (9) yields:

$\displaystyle\widetilde{\chi}^{(3)}_{+--+}(1,2,3)\overset{(J_{\perp}>0)}{=}e^{2\pi
im(x_{1}+x_{2}-x_{3})}G^{R}_{a^{\dagger}aaa^{\dagger}}(1,2,3)$
$\displaystyle+e^{-2\pi
im(x_{1}+x_{2}-x_{3})}G^{R}_{b^{\dagger}bbb^{\dagger}}(1,2,3)$
$\displaystyle+e^{2\pi
im(x_{1}-x_{2}+x_{3})}G^{R}_{a^{\dagger}abb^{\dagger}}(1,2,3)$
$\displaystyle+e^{-2\pi
im(x_{1}-x_{2}+x_{3})}G^{R}_{b^{\dagger}baa^{\dagger}}(1,2,3)$
$\displaystyle+e^{2\pi
im(x_{1}-x_{2}-x_{3})}G^{R}_{b^{\dagger}aba^{\dagger}}(1,2,3)$
$\displaystyle+e^{-2\pi
im(x_{1}-x_{2}-x_{3})}G^{R}_{a^{\dagger}bab^{\dagger}}(1,2,3).$ (21)

We have defined a set of response functions with the form of Eq. (10). For
instance, $G^{R}_{a^{\dagger}aaa^{\dagger}}$ is defined as:

$\displaystyle
G^{R}_{a^{\dagger}aaa^{\dagger}}(1,2,3)=i^{3}\theta(t_{1})\theta(t_{2}-t_{1})\theta(t_{3}-t_{2})\langle[[[a^{\dagger}(0),a(-1)],a(-2)],a^{\dagger}(-3)]\rangle.$
(22)

The other response functions are defined in the same vein. The expression for
these response function can be read off from Eq. (13) by plugging in
appropriate values of $l_{1,2,3}$ and $r_{1,2,3}$: $l=\sqrt{K}+1/(2\sqrt{K})$
and $r=-\sqrt{K}+1/(2\sqrt{K})$ for vertex operator $a$; for vertex operator
$b$, the value of $l$ and $r$ are switched. We have dropped the response
functions that violate the charge neutrality condition (e.g.
$G^{R}_{a^{\dagger}aab^{\dagger}}$) as they vanish in the thermodynamic limit.

The calculation of $\chi^{(3)}_{+--+}$ parallels the ferromagnetic case (Sec.
III.3). Note, however, $\chi^{(3)}_{+--+}$ is now complex and depends on the
magnetic field through the magnetization density $m$. We find only
$G_{a^{\dagger}aaa^{\dagger}}$ and $G_{b^{\dagger}bbb^{\dagger}}$ contribute
to photon echo. Using the approximation (Eq. (16)), we obtain after lengthy
calculation:

$\displaystyle\chi^{(3)}_{+--+}$
$\displaystyle\overset{(J_{\perp}>0)}{\sim}e^{-2\pi imu(t-\tau)}$
$\displaystyle\times\left\\{\begin{array}[]{cc}e^{-\pi(\Delta-2)T(\tau-t)}&(t\ll\tau)\\\
(t-\tau)e^{-\pi(\Delta-2)T(t-\tau)}&(t\gg\tau)\end{array}\right..$ (25)

Here, we have defined a parameter $\Delta=2K+1/(2K)$. The above is the key
result of this subsection.

Eq. (25) shows that the $\chi^{(3)}_{+--+}$ of antiferromagnetic chain
exhibits perfect photon echo similar to the ferromagnetic chain. Different
from the ferromagnetic chain, the signal now shows oscillations with the
frequency set by the magnetization density $m$, which, in turn, depends on the
magnetic field $B$. The time scale of the signal is $1/((\Delta-2)\pi T)$. In
particular, in the Heisenberg limit where $K\to 1/2$ ($\Delta\to 2$), the time
scale diverges — this reflects the Larmor precession of the total
magnetization, which we elaborate on in Appendix C.

We assess the validity of Eq. (25) by numerical integration. Fig. 4 show
$\chi^{(3)}_{+--+}$ for representative Luttinger parameters $K=0.7$, $K=1$,
and $K=1.5$. $\pi Tt_{\mathrm{w}}=1$. The echo is clearly visible from the
data. The data for different values of $\tau$ overlay when plotted as a
function of $t-\tau$, demonstrating that the echo is perfect. Numerical
integration suggests $\chi^{(3)}_{+--+}$ is independent of $t_{\mathrm{w}}$
within numerical error when $t$ and $\tau$ are large.

To test the asymptotic behavior of $\chi^{(3)}_{+--+}$, we multiply its
complex modulus with $\exp(\pi(\Delta-2)|t-\tau|)$. Eq. (25) suggests the
product shows linear behavior for $t>\tau$ and approaches a constant for
$t<\tau$. The insets of Fig. 4 show good agreement with Eq. (25). For $K=0.7$
and $\pi T\tau=50$, the data deviate slightly from the constant behavior; we
think this occurs because $\tau$ is not sufficiently large to suppress the
non-asymptotic contributions.

We then present the dependence of the photon echo signal on both $\tau$ and
$t$. Fig. 5(a)&(b) show the real and imaginary parts of $\chi^{(3)}_{+--+}$ as
a function of $t$ and $\tau$ for the Luttinger parameter $K=1$. The
magnetization density $2mu/T=1.15$. The waiting time $\pi Tt_{\mathrm{w}}=1$.
The photon echo appears as the bright feature persists along the diagonal
direction ($t=\tau$). Performing Fourier transform over $t$ and $\tau$
produces the two-dimensional spectrum (Fig. 5(c)&(d)). The photon echo appears
in the frequency domain as a highly anisotropic peak in the fourth quadrant.
The peak is approximately located at $\omega_{t}=2\pi mu$ and
$\omega_{\tau}=-2\pi mu$ as suggested by Eq. (25). Its width along the
diagonal of the fourth quadrant scales with $T$, whereas its width along the
anti-diagonal direction is resolution limited.

To recapitulate, the $\chi^{(3)}_{+--+}$ of the antiferromagnetic chain shows
clear signature of perfect photon echo. The echo signal is oscillatory with
the frequency set by magnetization density.

Figure 4: (a) Nonlinear magnetic susceptibility $\chi^{(3)}_{+--+}$ of an
antiferromagnetic chain as a function of $\pi T(t-\tau)$ for $\pi T\tau=50$
(blue) and $\pi T\tau=80$ (cyan). The waiting time $\pi Tt_{\mathrm{w}}=1$.
Luttinger parameter $K=0.7$. The data for different $\tau$ are rescaled by the
_same_ factor such that the maximum of the complex modulus of the $\pi
T\tau=50$ data is 1. Upper right inset shows the complex modulus of the
$t>\tau$ part of the data, multiplied by $\exp(\gamma|t-\tau|)$ with
$\gamma=\pi(\Delta-2)T$. Lower left inset the complex modulus of the $t<\tau$
part of the data, multiplied by $\exp(\gamma|t-\tau|)$ (b) The same as (a) but
for $K=1$. (c) The same as (a) but for $K=1.5$ Figure 5: (a)(b) The real and
imaginary parts of the nonlinear magnetic susceptibility $\chi^{(3)}_{+--+}$
of an antiferromagnetic chain as function of $\pi Tt$ and $\pi T\tau$. The
Luttinger parameter $K=1$. Magnetization density $2mu/T=1.15$. The waiting
time $\pi Tt_{\mathrm{w}}=1$. The data are scaled such that the maximum of the
absolute value of the data is 1. (c)(d) The real and imaginary parts of the
two-dimensional spectrum.

##  IV Spinon Lensing Picture

Figure 6: (a) The two-point correlation function (Eq. (31)) can be understood
in terms of a spinon creation/annihilation process. The spin raising operator
at $(-t_{1},-x_{1})$ creates a pair of spinons with left and right chirality,
which then propagate with velocities $\pm u$ (blue solid lines). The spin
lowering operator tries to annihilate the spinon pair at $(0,0)$. (b) The
four-point correlation function Eq. (32b) corresponds to the spinon creation
at $(-t_{3},-x_{3})$ and $(0,0)$ (blue solid lines), and antispinon creation
at $(-t_{2},-x_{2})$ and $(-t_{3},-x_{3})$ (blue dashed lines). (c) In the
lensing configuration (Eq. (33)), the left-moving spinon created at
$(-t_{3},-x_{3})$ is captured by the spin lowering operator at
$(-t_{2},-x_{2})$ and converted to a right-moving antispinon. Likewise, the
right-moving spinon is captured by the spin lowering operator at
$(-t_{1},-x_{1})$ and converted to a left-moving antispinon. The antispinons
then meet at $(0,0)$ and are annihilated by the spin raising operator. (d)
Behavior of $\widetilde{\chi}^{(3)}_{+--+}$ in the shaded region in (c), which
shows its magnitude reaches the maximum near the “focal point” $(0,0)$. Here,
$\pi Tt=\pi T\tau=10$, and $\pi Tt_{\mathrm{w}}=10$.

In the previous section, we have shown that the nonlinear magnetic
susceptibility $\chi^{(3)}_{+--+}$ of the Luttinger spin liquid exhibits
photon echo that resembles the perfectly rephasing echo in a few-body system
such as a single spin. In this section, we provide an intuitive picture that
clarifies its origin in the many-body system under consideration, namely
Luttinger spin liquids. We illustrate our picture on the ferromagnetic spin
chain and then generalize it to the antiferromagnetic chain, and discuss its
various features and implications.

To set the stage, we review the effect of the bosonized spin raising operator
$\exp(-i\theta)$. Let us consider the time evolution after the operator
$\exp(-i\theta(0,0))$ at $t=x=0$ on the initial state $|\Psi(0)\rangle$ at
$t=0$ which is assumed to be an energy eigenstate. At time $t$, the state is
given by

$\displaystyle|\Psi(t)\rangle$
$\displaystyle=e^{-i\mathcal{H}t}e^{-i\theta(0,0)}|\Psi(0)\rangle$
$\displaystyle=e^{-i\mathcal{H}t}e^{-i\theta(-t,0)}e^{+i\mathcal{H}t}e^{-i\mathcal{H}t}|\Psi(0)\rangle$
$\displaystyle\sim e^{-i\theta(-t,0)}|\Psi(0)\rangle,$ (26)

where we used the assumption that $\Psi(0)\rangle$ is an energy eigenstate and
dropped the overall phase factor. We can decompose the field $\theta$ to
chiral components as

$\displaystyle\theta(t,x)=$
$\displaystyle\theta_{L}(x^{+})+\theta_{R}(x^{-}),$ (27)

where

$\displaystyle\theta_{L}=\frac{1}{2}\left(\theta+\frac{\phi}{K}\right);\quad\theta_{R}=\frac{1}{2}\left(\theta-\frac{\phi}{K}\right).$
(28)

The equation of motion implies that each chiral components depends only on the
corresponding light cone coordinate, as in Eq. (27). Using this, we can
translate the time dependence into the position dependence, so that

$\displaystyle|\Psi(t)\rangle\sim$ $\displaystyle
e^{-i\theta_{L}(-ut)}e^{-i\theta_{R}(ut)}|\Psi(0)\rangle,$ (29)

where $\theta_{L,R}(\mp ut)$ is now understood as a static operator at
locations $\mp ut$, up to the overall phase factor which includes the one
coming from the commutation relation between $\theta_{R}$ and $\theta_{L}$.

Now, recall the equal-time commutation relation:

$\displaystyle\phi(x)e^{-i\theta_{L,R}(y)}=e^{-i\theta_{L,R}(y)}\left(\phi(x)-\frac{\pi}{2}\Theta{(x-y)}\right).$
(30)

Namely, $\exp(-i\theta_{L,R}(y))$ creates a kink of step $\pi/2$ in
$\phi$-field at the location $y$. Since this step corresponds to the
magnetization density $(1/2)\,\delta(x-y)$, this kink can be interpreted as a
spinon [17]. Eq. (29) then implies that the application of the operator
$\exp(-i\theta)$ at $t=0$ is equivalent to creation of a right-moving spinon
at $x=+ut$ and a left-moving spinon at $x=-ut$ at time $t$, when the other
operations are applied after the time $t$. The chiral vertex operators
$\exp(-i\theta_{L,R})$ can be viewed as spinon creation operators of the
corresponding chirality. This justifies the following simple visual picture:
Applying $\exp(-i\theta(0,0))$ creates a pair of the right-moving and left-
moving spinons at $t=x=0$. Then these spinons propagate with the constant
velocity $\pm u$.

Real-time correlation functions of the vertex operators can be qualitatively
understood in terms of this visual picture. As the simplest example, consider
the two-point correlation function (Fig. 6a)

$\displaystyle\langle S^{-}(0)S^{+}(-1)\rangle\sim\langle
e^{i\theta(0,0)}e^{-i\theta(-t_{1},-x_{1})}\rangle$
$\displaystyle\quad\sim\langle
e^{i\theta_{L}(0)}e^{i\theta_{R}(0)}e^{-i\theta_{L}(-x_{1}-ut_{1})}e^{-i\theta_{R}(-x_{1}+ut_{1})}\rangle$
$\displaystyle\quad\sim\left|\frac{\sinh(\pi Tx^{+}_{1})}{\pi
T\epsilon}\right|^{-\frac{1}{4K}}\left|\frac{\sinh(\pi Tx^{-}_{1})}{\pi
T\epsilon}\right|^{-\frac{1}{4K}}.$ (31)

Here, we set $t_{1}>0$ and omit an overall phase factor. This correlation
function corresponds to the creation of a pair of spinons at
$(-t_{1},-x_{1})$, and the annihilation of a pair of spinons (or equivalently
the creation of a pair of anti-spinons) at the origin. From the equation of
motion (the second line), this correlation function can be interpreted as an
expectation value of spinon creation operators at the split locations
$-x_{1}\pm ut_{1}$ and _two_ spinon annihilation operators at $0$. For a
generic choice of $t_{1},x_{1}$, these three locations are different. As a
result, the expectation value, although does not vanish completely, is
exponentially suppressed with respect to the correlation length $u/T$ at
finite temperature $T$.

If we wish to maximize the correlation function, we should try to “catch” the
spinons created at $(-t_{1},-x_{1})$ at the later time and annihilate them. If
we can manage to annihilate all the spinons created earlier, the correlation
function would be reduced to the equal-time correlation function of the
identity operator, which does not decay as a function of time $t$ at all!
Unfortunately, it is impossible to annihilate both spinons created at
$(-t_{1},-x_{1})$ by the single operator $\exp(i\theta(0,0))$, since spinons
split and are located on different positions at later times.

We note it is possible to annihilate one of them, for example the left-moving
spinon, by choosing the location of the second operator on the left side of
the light cone ($x_{1}=-ut_{1}$, or $x^{+}_{1}=0$). Setting $x^{+}_{1}=0$ in
Eq. (31) maximizes the first factor, which is formally divergent in the
scaling limit (the cutoff $\epsilon\to 0^{+}$) but bounded to be a finite
constant due to the finite $\epsilon$ in a realistic system. However,
$x^{+}_{1}=0$ implies $x^{-}_{1}=2ut_{1}$, leading to the exponential
suppression of the other factor in Eq. (31). Thus the two-point correlation
function is exponentially suppressed for large $t_{1}$, for any choice of the
relative spatial location $x_{1}$.

The situation is quite different for four-point correlation functions and
corresponding response functions. It is convenient to rewrite the response
function Eq. (15) as a sum over various four-point correlation functions by
expanding the nested commutators:

$\displaystyle\widetilde{\chi}^{(3)}_{+--+}(1,2,3)$ $\displaystyle\sim
R_{1}+R_{2}+R_{3}+R_{4}$
$\displaystyle-R^{\prime}_{1}-R^{\prime}_{2}-R^{\prime}_{3}-R^{\prime}_{4},$
(32a) where $\displaystyle R_{1}$ $\displaystyle=\langle
e^{-i\theta(0)}e^{i\theta(-1)}e^{i\theta(-2)}e^{-i\theta(-3)}\rangle;$ (32b)
$\displaystyle R_{2}$ $\displaystyle=\langle
e^{-i\theta(-3)}e^{i\theta(-1)}e^{-i\theta(0)}e^{i\theta(-2)}\rangle;$ (32c)
$\displaystyle R_{3}$ $\displaystyle=\langle
e^{-i\theta(-3)}e^{i\theta(-2)}e^{-i\theta(0)}e^{i\theta(-1)}\rangle;$ (32d)
$\displaystyle R_{4}$ $\displaystyle=\langle
e^{i\theta(-2)}e^{i\theta(-1)}e^{-i\theta(0)}e^{-i\theta(-3)}\rangle.$ (32e)

$R^{\prime}_{i}$ is related to $R_{i}$ by reversing the order of operators in
the product. Following Ref. 1, 2, each contribution can be identified as a
Liouville pathway.

We find that each four-point correlation function can be made free from
exponential suppression for an appropriate choice of locations of the
operators at given set of times. As exchanging the order of operators in
$R_{i}$ only results in a phase factor, it is sufficient for our purpose to
consider the correlation function $R_{1}$ (Fig. 6b). The spinon pairs are
created at $(-t_{3},-x_{3})$ and $(0,0)$, and annihilated at $(-t_{2},-x_{2})$
and $(-t_{1},-x_{1})$. Unlike the case with the two-point correlation function
(Eq. 31), here, we can catch and annihilate all the spinons created earlier by
choosing (Fig. 6c):

$\displaystyle x_{1}=u\tau;\quad x_{2}=-u(\tau+t_{\mathrm{w}});\quad
x_{3}=-ut_{\mathrm{w}}.$ (33)

For this configuration, the operator at $(-t_{2},-x_{2})$ annihilates the
left-moving spinon created at $(-t_{3},-x_{3})$ and creates a right-moving
antispinon. The operator at $(-t_{1},-x_{1})$ annihilates the right-moving
spinon created at $(-t_{3},-x_{3})$ and creates a left-moving antispinon. The
two antispinons moving to the opposite directions finally meet at the origin
at time $0$, and are annihilated by the operator $\exp(i\theta)$! As a
consequence, the magnitude of the four-point correlation function reaches its
maximum and does not decay even in the long-time limit as long as the spatial
coordinates are chosen according to Eq. (33). We call this phenomenon, which
was absent in the two-point correlation function, the _spinon lensing_ , that
is, it is possible to focus the two (anti)spinons to the same point at time
$0$, by placing the operators judiciously at earlier times. We also note that
the spatial mirror reflection of Eq. (33) is equally valid for lensing.

The above reasoning can be made precise by applying the equation motion to
reduce the $4$-point correlation function to an equal-time correlation
function at time $0$ of multiple operators — spinon creation operators at
locations $-x_{3}\pm ut_{3}$ and two at $0$, and antispinon creation operators
at locations $-x_{2}\pm ut_{2}$ and $-x_{1}\pm ut_{1}$ (Fig. 6b). Under the
condition Eq. (33), the locations of the spinon creation operators match with
those of the spinon annihilation operators, and consequently the correlation
function $R_{1}$ reaches the maximum.

The same phenomena occur in the other four-point correlation functions.
Summing them up, the nonlinear response $\widetilde{\chi}^{(3)}_{+--+}$ as a
function of spacetime coordinates $(t_{1},x_{1})$, $(t_{2},x_{2})$, and
$(t_{3},x_{3})$, reaches a maximum at the lensing configuration Eq. (33). Fig.
6d shows the behavior of $\chi^{(3)}_{+--+}$ as a function of the detection
position in a representative lensing configuration. Here we choose $\pi Tt=\pi
T\tau=10$ and $\pi Tt_{\mathrm{w}}=10$. It shows that the response is maximum
at the “focal point” $(0,0)$. As we have used a finite short distance cutoff
$\epsilon$ to regularize the algebraic singularity at the light cone, we have
slightly shifted $x_{1}$ and $x_{2}$ away from the light cone by a small
distance $\epsilon$ to suppress the effect of the regularizer.

Now, the optical nonlinear response $\chi^{(3)}_{+--+}$ is related to
$\widetilde{\chi}^{(3)}_{+--+}$ by spatial integration. When the time delays
$t=\tau$, the spatial integration is dominated by the neighborhood of the
lensing configuration Eq. (33) and its mirror reflection, and, as a result,
remains non-vanishing even in the long-time limit
$t=\tau\to\infty,t_{\mathrm{w}}\to\infty$. This explains the origin of the
prefect photon echo observed in the $\chi^{(3)}_{+--+}$ of the ferromagnetic
chain.

Turning to the antiferromagnetic chain, $\widetilde{\chi}^{(3)}_{+--+}$ is now
a linear combination of various four-point response functions (Eq. (21)).
Among these, $G^{R}_{a^{\dagger}aaa^{\dagger}}$ and
$G^{R}_{b^{\dagger}bbb^{\dagger}}$ support the lensing phenonenon similar to
Fig. 6 with the only difference being the type of fractional excitations
created/annihilated by the spin operators. In the antiferromagnetic chain,
$S^{+}$ creates a pair of spinons plus a pair of Laughlin quasiparticle and
quasihole [17]. Note the created spinon and the Laughlin quasiparticle (and
similarly the created antispinon and the Laughlin quasihole) are superimposed
on each other, and therefore the world lines shown in Fig. 6 should be
interpreted as the world line of the composite object. Our numerical
calculations show that, indeed, only these two response functions produce
photon echo whereas the others response functions do not.

In Sec. II, we have pointed out that the nonlinear susceptibilities
$\chi^{(3)}_{+-+-}$ and $\chi^{(3)}_{++--}$ do not exhibit photon echo. This
can now be easily understood in terms of lensing. For these two
susceptibilities, it is impossible to arrange the spin operators in such a way
that all the created fractional excitations get annihilated at a later time.
Our calculations indeed confirm this.

Our discussion so far is based on the ideal Luttinger spin liquid at the RG
fixed point. In general, there are RG-irrelevant perturbations to the fixed
point Hamiltonian Eq. (3), but they only give sub-leading corrections to most
of the physical quantities at low temperatures. This also makes these RG-
irrelevant corrections difficult to detect in experiments. Nevertheless, in
the 2DCS in Luttinger spin liquid, the RG-irrelevant perturbations can produce
pronounced effects because the lensing of fractional excitations relies on two
features of the fixed point Hamiltonian: First, the phonon modes are the exact
eigenstates of the Hamiltonian. Second, the phonon dispersion relation is
exactly linear. These features ensure that a wave packet of fractional
excitation can propagate indefinitely through the system without dissipation
or dispersion. RG-irrelevant corrections to the fixed point Hamiltonian would
prevent the indefinite propagation of the wave packet, and, therefore, could
be sensitively detected by the suppression of lensing.

In the XXZ spin chain, the RG-irrelevant perturbations include the umklapp
terms such as $\cos(4\phi)$ [15, 16], which result in the damping of phonon
modes at finite temperature. Consequently, the wave packets of fractional
excitations acquire finite life time. The lensing phenomenon shown in Fig. 6
is suppressed when $\tau$ or $t_{\mathrm{w}}$ exceeds the life time of these
excitations. This, in turn, will be manifest as the decay of the photon echo
signal as $\tau$ or $t_{\mathrm{w}}$ increases, in close analogy with the
dissipation-induced photon echo decay in few-body systems [1, 2].
Straightforward dimensional analysis suggests the decay rate vanishes as
$T^{1+\eta}$ as $T\to 0$ where $\eta$ is the dimension of the umklapp term.

On the other hand, there also exist a different family of RG-irrelevant
perturbations, such as $(\nabla^{2}\phi)^{2}$, which produce a small curvature
in the phonon dispersion relation while keeping them as the exact eigenstates
of the Hamiltonian. As a result, the wave packet disperses as it propagates,
but there is no true dissipation. The lensing picture suggests that the
dispersion of wave packet would also result in the decay of the photon echo.
For instance, the wave packet of the left-moving spinon created by $S^{+}(-3)$
would be quite dispersed when $\tau$ is large. It could not completely
annihilate with the left-moving antispinon created by $S^{-}(-2)$ as the
latter’s wave packet is still sharp. This would spoil lensing. Thus, the
lensing picture suggests a dispersion-induced decay mechanism for the photon
echo in the Luttinger spin liquid. In the next section, we examine this
mechanism in more detail.

##  V Dispersion-induced decay

Figure 7: (a) The nonlinear magnetic susceptibility $\chi^{(3)}_{+--+}$ as a
function of $t$ and $\tau$, computed from the harmonic chain model Eq. (34).
The waiting time $t_{\mathrm{w}}=4a/u$. The temperature $T=1.2u/a$. The data
are scaled such that the maximum is 1. (b)(c) The real and imaginary parts of
the two-dimensional spectrum obtained by Fourier-transforming the data in
panel (a). Only the first and the fourth quadrants are shown. The other two
quadrants are obtained by complex conjugation. (d) $\chi^{(3)}_{+--+}$,
measured at $t=\tau$, as a function of $(a^{2}Tt/u^{2})^{1/2}$ for various
temperatures. $t_{\mathrm{w}}=4a/u$. (e) $\chi^{(3)}_{+--+}$, measured at
$t=\tau$, as a function of waiting time $t_{\mathrm{w}}$ for two
representative values of $\tau$. The temperature $T=0.8u/a$.

The spinon lensing picture suggests that the dispersion of the wave packets
could also lead to the decay of photon echo. In this section, we put this idea
to test by a numerical experiment. We consider the harmonic chain, which is a
discretization of the fixed point Hamiltonian:

$\displaystyle H=\frac{u}{2\pi
a}\sum_{n}\frac{(\phi_{n+\frac{1}{2}}-\phi_{n-\frac{1}{2}})^{2}}{K}+K(\theta_{n+1}-\theta_{n})^{2}.$
(34)

Here, $a$ is the lattice constant. $\theta_{n}$ resides on the lattice site
labeled by $n$, whereas $\phi_{n+1/2}$ resides on the midpoint of the lattice
link connecting the site $n$ and $n+1$. Their commutation relation is given
by: $[\phi_{n+1/2},\theta_{n^{\prime}}]=-i\pi\,\Theta(n-n^{\prime}+1/2)$. It
is sufficient for our purpose to consider the ferromagnetic chain where
$S^{\pm}_{n}=\exp(\mp i\theta_{n})$. We set the Luttinger parameter $K=1$.

The phonon modes are exact eigenstates of Eq. (34) owing to its quadratic
form. The phonon dispersion relation is now nonlinear:
$\omega_{q}=2u/a\,|\sin(qa/2)|$. Therefore, Eq. (34) represents an idealized
model system to study the dispersion-induced decay without the complications
of dissipation effects. By contrast, in a microscopic spin lattice model, both
effects are present and difficult to disentangle.

Fig. 7(a) shows $\chi^{(3)}_{+--+}$ as a function of $t$ and $\tau$ at
temperature $T=1.2u/a$. The waiting time $t_{\mathrm{w}}=4a/u$. Note the data
are strictly real as the case in continuum (Fig. 3). We find a clear signature
of photon echo running along the diagonal of the $(t,\tau)$ plane. However,
the echo signal decays slowly along the diagonal direction owing to the
dispersion effect. In the frequency domain, this decay manifests itself as a
slight broadening of the photon echo peaks along the anti-diagonal direction
of the frequency plane (Fig. 7(b)&(c)).

We now investigate the temperature dependence of the dispersion-induced decay.
Expanding the phonon dispersion near $q=0$, we obtain
$\omega_{q}=u(|q|-a^{2}|q|^{3}/24+\cdots)$. Therefore, the width of the wave
packet grows as $(a^{2}ut)^{1/3}$ as it propagates. At temperature $T$, this
defines a dispersion time scale $\tau_{\mathrm{disp}}$:

$\displaystyle(a^{2}u\tau_{\mathrm{disp}})^{1/3}\sim\xi\sim\frac{u}{T};\Rightarrow\tau_{\mathrm{disp}}\sim\frac{u^{2}}{a^{2}T^{3}}.$
(35)

Here, $\xi$ is the spin correlation length. Beyond this time scale, the wave
packet is essentially indistinguishable from thermal fluctuations. Thus, we
hypothesize that the decay of the echo signal is controlled by
$\tau_{\mathrm{disp}}$.

Our data support this hypothesis. Fig. 7(c) shows the semi-log plot of
$\chi^{(3)}_{+--+}$ at $t=\tau$ as a function of
$(\tau/\tau_{\mathrm{disp}})^{1/2}$. The data for different $T$ approximately
collapse to a straight line. This suggests the echo signal decays as a
stretched exponential, $\exp(-C(\tau/\tau_{\mathrm{disp}})^{\nu})$, where the
exponent $\nu=1/2$ and $C$ is a numeric constant. The decay rate
$1/\tau_{\mathrm{disp}}\sim T^{3}$, consistent with the dimension $\eta=2$ of
the higher derivative term $(\nabla^{2}\phi)^{2}$. The origin of the
stretching exponent $\nu=1/2$ is unclear at the moment but likely associated
with the asymptotic behavior of Airy functions.

Having investigated the decay of the echo as a function of $\tau$, we turn to
its $t_{\mathrm{w}}$ dependence. Fig. 7(e) shows $\chi^{(3)}_{+--+}$ at
$t=\tau$ as a function of $t_{\mathrm{w}}$ for two representative values of
$\tau$. The echo signal initially decrease as $t_{\mathrm{w}}$ increases,
reflecting the effect of dispersion, but eventually saturate to a finite value
This is also consistent with the lensing picture. A straightforward
calculation shows that the spacetime-dependent susceptibility
$\widetilde{\chi}_{+--+}$ saturates to finite value as
$t_{\mathrm{w}}\to\infty$ at the lensing configuration Eq. (33). In a few-body
system, the decay of photon echo as a function of $t_{\mathrm{w}}$ reflects
the thermal relaxation of the optical population. Here, as the phonon modes
are exact eigenstates of the Hamiltonian Eq. (34), the phonon population
cannot relax. We therefore heuristically attribute the saturation to the
absence of thermal relaxation in Eq. (34).

To recapitulate, by a numerical experiment, we have shown that the photon echo
signal decays as a function of $\tau$ in the presence of wave packet
dispersion. The decay is controlled by a dispersion time scale
$\tau_{\mathrm{disp}}$. Meanwhile, the echo signal saturates as
$t_{\mathrm{w}}\to\infty$, which we attribute to the lack of thermalization in
the model system.

##  VI Discussion

In this work, we have shown that the nonlinear magnetic susceptibility
$\chi^{(3)}_{+--+}$ and its complex conjugate $\chi^{(3)}_{-++-}$ of the
Luttinger spin liquid exhibit photon echo that resembles the perfectly
rephasing photon echo in a few-body problem. However, the rephasing picture
does not directly apply to the present system in that its energy spectrum is
continuous. Instead, the echo signal arises from the lensing of the fractional
excitations, and its decay as a function of the pulse delays is a sensitive
diagnostic for the dissipation and/or dispersion.

The photon echo can be directly measured by THz 2DCS on spin chain materials
that host Luttinger spin liquids. For example, Cs2CoCl4 is thought to be a
realization of the easy-plane antiferromagnetic ($J_{\perp}>J_{z}>0$) $S=1/2$
XXZ chain [23]. Although our analysis is carried out in the circular
polarization basis, one could simplify the experimental set up by measuring
the $\chi^{(3)}$ response with linear polarization (e.g. $\chi^{(3)}_{xxxx}$)
since the linear polarization basis and the circular polarization basis are
related by a linear transformation.

We find the lensing of fractional excitations to be a convenient picture for
understanding the photon echo in the Luttinger spin liquid. A crucial feature
of the lensing is the refocusing of the wave packets world lines, reminiscent
of the refocusing of quantum phase accumulation in the NMR spin echo or photon
echo in few-body systems. However, the lensing is unique to many-body system
in that it entails the propagation of wave packets. It could be viewed as a
conceptual extension of the more familiar interference picture [1, 2] commonly
used in the study of photon echo in few-body systems from the time domain to
the spacetime domain.

An earlier work on the photon echo of quantum Ising chain uses the time-domain
interference picture [7], which is made possible by a mathematical mapping
that relates the nonlinear response of a many-body system to that of an
ensemble of independent few-body systems. Although it might be possible to
adapt this methodology for the Luttinger spin liquid, we think it will be much
less illuminating given the complex structure of the optical matrix elements.

An interesting consequence of the lensing picture is that the dispersion of
the wave packet alone could lead to the decay of photon echo signal. This
decay mechanism finds no immediate analogue in few-body systems. In the
presence of dispersion but no dissipation, the echo signal decays as $\tau$
increases and eventually disappears when $\tau$ far exceeds the dispersion
time scale. However, the echo signal does not disappear when the waiting time
$t_{\mathrm{w}}\to\infty$. This points toward different uses of the pulse
delays $\tau$ and $t_{\mathrm{w}}$ — the former could be used as a dial to
monitor both dispersion and dissipation effects, whereas the latter is
sensitive to dissipation alone.

Our work opens a few directions worth further investigation. First of all, our
discussion on the umklapp terms, or more broadly dissipation effects, has been
qualitative. It will be useful to examine their impact on the photon echo
quantitatively by employing the quantum kinetic theory [24, 25]. Secondly, it
would be interesting to compare our predictions with lattice model
calculations. Our preliminary results on the XY spin chain [26, 27, 28, 29,
30] show good agreement with the bosonization predictions and will be
published elsewhere. Thirdly, our analysis focuses on the gapless phase of the
XXZ type spin chain. It would be interesting to explore the nonlinear photon
echo in the gapped phase by using a semiclassical treatment [31, 32, 33], or
in the vicinity of the Heisenberg point where the proximate $SU(2)$ symmetry
might gives rise to new universal features [34, 35]. We also note an
interesting recent work where the generalized hydrodynamics is employed to
compute some nonlinear responses of integrable systems including the XXZ spin
chain [36]. Finally, even though we exclusively consider the spatially uniform
(momentum $q=0$) magnetic response in this work, our calculations can be
generalized to finite $q$. It has been shown that, in a magnetized
antiferromagnetic Heisenberg chain, the magnetic resonance mode acquires a
nonlinear dispersion at $q\neq 0$ due to spinon interactions [37, 38]. In
light of these results, we expect that the $q\neq 0$ nonlinear magnetic
resonance may exhibit rich physics. Although optical measurements usually
probes the $q=0$ response, the presence of Dzyaloshinskii-Moriya (DM)
interaction in a spin chain may allow for optical measurement of the $q\neq 0$
responses. The DM interaction may be eliminated by a gauge transformation at
the expense of a momentum boost [39]. Hence, the $q=0$ response of the spin
chain with a uniform DM interaction is equivalent to the finite $q$ response
of the corresponding spin chain without the DM interaction [40, 41].

Looking beyond the Luttinger spin liquids, we think our method and the
physical picture could be applicable to other one-dimensional critical systems
and potentially higher dimensional systems. Perhaps the most experimentally
relevant are gapless systems with charged excitations, where the charge
degrees of freedom could directly couple to the electric component of the THz
pulse and therefore produce stronger nonlinear response. In short, we believe
future investigation on the nonlinear response of many-body systems will
uncover far richer dynamical phenomena and offer deeper insight into these
systems.

###### Acknowledgements.

YW and MO thank the hospitality of the Kavli Institute for Theoretical
Physics, University of California, Santa Barbara, where a part of this work
was performed. MO also thanks Shunsuke Furukawa, Shunsuke C. Furuya, and
Masahiro Sato for useful comments. This work was supported in part by the
National Natural Science Foundation of China (Grant No. 11974396), the
Strategic Priority Research Program of the Chinese Academy of Sciences (Grant
No. XDB33020300), Japan Society for the Promotion of Science (Grant Nos.
JP18H03686 and JP19H01808), Japan Science and Technology Agency CREST (Grant
No. JPMJCR19T2), and the National Science Foundation (Grant No. NSF
PHY-1748958). A part of the numerical calculations were carried out on
Tianhe-1A at the National Supercomputing Center in Tianjin, China.

##  Appendix A Commutation relation prescriptions

In this appendix, we discuss the different choices for the commutation
relation between $\phi$ and $\theta$ fields, and the subtle issues [42]
associated with these choices. Although elementary, these issues are pertinent
to the calculation of higher order response functions.

In bosonization, $\phi$ and $\nabla\theta$ form a pair of canonical conjugate
variables:

$\displaystyle[\phi(x),\nabla\theta(y)]=i\pi\delta(x-y).$ (36)

The above does not uniquely fix the commutator between $\phi$ and $\theta$.
Several popular choices exist in literature. In this work, we use:

$\displaystyle[\phi(x),\theta(y)]=-i\pi\Theta(x-y).$ (37)

We term this choice the “Heaviside prescription” for latter convenience. A
closely related variant is $[\phi(x),\theta(y)]=i\pi\Theta(y-x)$ [15]. Since
the two variants are rather similar, we shall focus the prescription Eq. (37).
Another popular choice is [16]:

$\displaystyle[\phi(x),\theta(y)]=-i\frac{\pi}{2}\mathrm{Sgn}(x-y).$ (38)

We term this choice the “signum prescription”.

In what follows, we show that the different prescriptions result in different
bosonization dictionaries. We first consider the fermion field operator. For
non-interacting fermions (Luttinger parameter $K=1$), the left and right
chiral boson fields are dynamically coupled:

$\displaystyle\phi_{L}=\theta+\phi,\quad\phi_{R}=\theta-\phi.$ (39)

The commutation relation of $\phi_{L}$, and likewise $\phi_{R}$, are identical
for both prescriptions. However, the commutator between $\phi_{L}$ and
$\phi_{R}$ depends on the prescription:

$\displaystyle[\phi_{L}(x),\phi_{R}(y)]=\left\\{\begin{array}[]{cc}-i\pi&(\mathrm{H.})\\\
0&(\mathrm{S.})\end{array}\right..$ (42)

In other words, boson fields with different chiralities do not commute
(commute) in the Heaviside (signum) prescription. As a result, in the
Heaviside prescription, the left and right chiral fermions:

$\displaystyle\psi_{L}\sim e^{i\phi_{L}};\quad\psi_{R}\sim
e^{i\phi_{R}},\quad(\mathrm{H.})$ (43a) anti-commute thanks to the commutator
between $\phi_{L}$ and $\phi_{R}$. By contrast, in the signum prescription,
Klein factors are needed to ensure the anti-commutation relation:
$\displaystyle\psi_{L}\sim\eta_{L}e^{i\phi_{L}};\quad\psi_{R}\sim\eta_{R}e^{i\phi_{R}},\quad(\mathrm{S.}),$
(43b)

where $\eta_{L,R}$ are the Klein factors obeying the Clifford algebra:
$\eta^{2}_{L}=\eta^{2}_{R}=1$, and $\eta_{L}\eta_{R}=-\eta_{R}\eta_{L}$.

It is important to bear in mind that choosing different prescriptions do not
change the dynamics of the boson fields. Furthermore, including fermion
interactions ($K\neq 1$) does not spoil the commutator between the left and
right chiral boson fields. In this case, the two dynamically decoupled fields
read:

$\displaystyle\phi_{L}=\sqrt{K}\theta+\frac{\phi}{\sqrt{K}},\quad\phi_{R}=\sqrt{K}\theta-\frac{\phi}{\sqrt{K}},$
(44)

and we find Eq. (42) holds regardless the value of $K$.

In the next, we consider the Jordan-Wigner string operator, which plays a
crucial role in the bosonization of the XXZ spin chain:

$\displaystyle\mathcal{S}_{j}=\cos(\pi\sum_{n<j}c^{\dagger}_{n}c^{\phantom{\dagger}}_{n}),$
(45)

where $c_{n}$ is the fermion annihilation operator on lattice site $n$. In
particular,

$\displaystyle
c_{i}\mathcal{S}_{j}=(-)^{\Theta(j-i-1/2)}\mathcal{S}_{j}c_{i}.$ (46)

We now bosonize the string operator following the standard procedure [15, 16]:

$\displaystyle\mathcal{S}_{j}$
$\displaystyle\approx\cos(k_{F}x-\int^{x}_{-\infty}\nabla\phi(y)dy)$
$\displaystyle=\cos(k_{F}x-\phi(x)+\phi(-\infty)).$ (47)

In the Heaviside prescription, the boundary term $\phi(-\infty)$ commutes with
$\theta(x)$. We thus may regard it as a number and omit it:

$\displaystyle\mathcal{S}_{j}=\cos(k_{F}x-\phi(x)).\quad(\mathrm{H.}).$ (48a)
By contrast, in the signum prescription, we must keep the boundary term as
$[\phi(-\infty),\theta(x)]\neq 0$:
$\displaystyle\mathcal{S}_{j}=\cos(k_{F}x-\phi(x)+\phi(-\infty)).\quad(\mathrm{S.}).$
(48b)

In the signum prescription, the boundary term $\phi(-\infty)$ is needed to
reproduce the commutation relation between the fermion field operator and the
string operator (Eq. 46). Had we dropped the boundary term, we would have
found:

$\displaystyle
c_{i}\mathcal{S}_{j}\overset{!}{=}e^{-i\pi/2\mathrm{Sgn}(x-y)}\mathcal{S}_{j}c_{i},$
(49)

which disagrees with Eq. 46.

Finally, we are ready to bosonize the spin operators of the XXZ chain using
the bosonization expressions of the fermion operator and the string operator.
We assume the chain is antiferromagnetic ($J_{\perp}<0$) without loss of
generality. Following Refs. 15, 16 but keeping a close eye on the Klein
factors and the boundary term, we find the staggered components of the spin
operators are given by:

$\displaystyle
N^{z}_{j}\sim\left\\{\begin{array}[]{cc}\cos(2\phi)&(\mathrm{H.})\\\
\eta_{L}\eta_{R}e^{-2i\phi}+\mathrm{H.C.}&(\mathrm{S.})\end{array}\right.,$
(52) $\displaystyle
N^{-}_{j}\sim\left\\{\begin{array}[]{cc}e^{i\theta}&(\mathrm{H.})\\\
\eta_{L}e^{i(\theta+\phi(-\infty))}+\eta_{R}e^{i(\theta-\phi(-\infty))}&(\mathrm{S.})\end{array}\right..$
(55)

We stress that the Klein factors and the boundary term are essential in
reproducing the correct commutation relations. For instance, had we dropped
the Klein factors and the boundary term, we would have obtained:

$\displaystyle N^{z}(x)N^{-}(y)\overset{!}{=}-N^{-}(y)N^{z}(x),$ (56)

using the signum prescription. This is incorrect because these two operators
must commute.

For the uniform components, we have:

$\displaystyle M^{z}_{j}\sim-\frac{1}{\pi}\nabla\phi,$ (57)

which is independent of the prescription, and

$\displaystyle
M^{-}_{j}\sim\left\\{\begin{array}[]{cc}e^{i\theta}\cos(2\phi)&(\mathrm{H.})\\\
\eta_{L}e^{i(\theta+2\phi-\phi(-\infty))}+\eta_{R}e^{i(\theta-2\phi+\phi(-\infty))}&(\mathrm{S.})\end{array}\right..$
(60)

We see that the bosonization dictionary takes a much simpler form in the
Heaviside prescription. In fact, the Heaviside prescription has been used in
the literature [43, 44] when the precise bosonization formulae are needed. On
the other hand, the full bosonization formulae of the spin operators in the
$S=1/2$ chain in the signum prescription given in Eqs. (52), (55), and (60)
are new to the best of our knowledge. For the signum prescription, one might
hope that we could omit the Klein factors and the boundary term in
calculations and still get correct results. This is indeed the case when
calculating two-point functions, which is a fortunate coincidence. In fact,
dropping the Klein factors and the boundary term from the bosonization
dictionary can lead to incorrect results in the signum prescription when
calculating higher order response functions.

To illustrate this point, consider the following four-point response function:

$G^{R}(1,2,3)=i^{3}\Theta(t_{1})\Theta(t_{2}-t_{1})\Theta(t_{3}-t_{2})\\\
\times\langle[[[N^{x}(0),N^{z}(-1)],N^{x}(-2)],N^{z}(-3)]\rangle.$ (61)

Had we dropped the Klein factors and the boundary term in the Signum
prescription, we would have found $N^{x}$ and $N^{z}$ anti-commute when they
are spacelike separated (Eq. (56)). As a result, $G^{R}\neq 0$ even when the
point of detection is outside the light cone of the perturbation, thereby
violating the causality principle.

To recapitulate, different choices for the commutation relation of $\theta$
and $\phi$ lead to different bosonization dictionaries. The bosonization
dictionary in the signum prescription includes Klein factors and the boundary
term, which cannot be omitted in calculating the high-order response
functions.

##  Appendix B Causality of the response function

Figure 8: (a)(b)(c) illustrates three possible cases in which the spacetime
points $(-t_{1},-x_{1})$ (red dot), $(-t_{2},-x_{2})$ (green dot),
$(-t_{3},-x_{3})$ (blue dot) lie outside of the past light cone (dashed line)
of the origin (black dot). In any of the three cases, the response function
$G^{R}$ (Eq. (10)) vanishes. (d) Causality requires that all three points to
lie inside the past light cone of the origin to induce a non-zero response.

In this appendix, we prove the causality of the response function $G^{R}$ (Eq.
(10)).

To set the stage, we recall that causality in relativistic quantum field
theory requires that two observables commute if their separation is spacelike:

$\displaystyle[O_{1}(t_{1},x_{1}),O_{2}(t_{2},x_{2})]=0,\quad(u^{2}t_{12}^{2}-x^{2}_{12}<0)$
(62)

where $O_{1,2}$ are arbitrary observables.

We first prove the following general statement: $G^{R}$ is causal provided
that all the operators that appear in Eq. (10) fulfil the condition Eq. (62).
Specifically, we need to show $G^{R}=0$ if any of the three spacetime points
$(-t_{1},-x_{1})$, $(-t_{2},-x_{2})$, $(-t_{3},-x_{3})$ are outside the light
cone of $(0,0)$. We prove this statement by exhaustion:

  1. (i)

If $(-t_{1},-x_{1})$ is outside the light cone of $(0,0)$ (Fig. 8(a)), then
Eq. (62) implies $[V_{0},V_{1}]=0$, and consequently $G^{R}=0$.

  2. (ii)

Now suppose $(-t_{1},-x_{1})$ is inside the light cone of $(0,0)$. If
$(-t_{2},-x_{2})$ is outside the light cone of $(0,0)$, it must also be
outside the light cone of $(-t_{1},-x_{1})$ (Fig. 8(b)). Thus,
$[V_{0},V_{2}]=[V_{1},V_{2}]=0$ by virtue of Eq. (62). It then follows that
$[[V_{0},V_{1}],V_{2}]=0$, and consequently $G^{R}=0$.

  3. (iii)

Now suppose both $(-t_{1},-x_{1})$ and $(-t_{2},-x_{2})$ are inside the light
cone of $(0,0)$. If $(-t_{3},-x_{3})$ is outside the light cone of $(0,0)$, it
must also be outside the light cone of $(-t_{1},-x_{1})$ and $(-t_{2},-x_{2})$
(Fig. 8(c)). Thus, $[V_{0},V_{3}]=[V_{1},V_{3}]=[V_{2},V_{3}]=0$ thanks to Eq.
(62). It follows that $[[[V_{0},V_{1}],V_{2}],V_{3}]=0$, and thus $G^{R}=0$.

This completes our proof.

In the next, we turn to the Tomonaga-Luttinger liquid theory and show that the
local vertex operators:

$\displaystyle V=e^{im\theta+2in\phi},\quad(m,n\in\mathbb{Z})$ (63)

indeed fulfill the condition Eq. (62). To this end, we consider two vertex
operators $V_{1}$ and $V_{2}$. Suppose they are spacelike separated. We choose
a reference frame in which they are synchronous and compute their equal time
commutator:

$\displaystyle V_{1}V_{2}$
$\displaystyle=e^{i(m_{1}\theta(x_{1})+2n_{1}\phi_{(}x_{1}))}e^{i(m_{2}\theta(x_{2})+2n_{2}\phi_{(}x_{2}))}$
$\displaystyle=e^{2\pi
i(n_{1}m_{2}\Theta(x_{1}-x_{2})-m_{1}n_{2}\Theta(x_{2}-x_{1}))}V_{2}V_{1}$
$\displaystyle=V_{2}V_{1}.$ (64)

The second equality follows from the Baker-Campbell-Hausdorff formula; the
third equality follows from the fact $m_{1,2},n_{1,2}\in\mathbb{Z}$. The above
immediately implies:

$\displaystyle[V_{1}(x_{1}),V_{2}(x_{2})]=0.$ (65)

Thus, we have verified that the local vertex operators fulfill the causality
condition Eq. (62).

Combining the above results, we conclude that the response function $G^{R}$
(Eq. (10)) is causal. We note that one could also verify the causality of
$G^{R}$ by using its explicit expression Eq. (13).

##  Appendix C Antiferromagnetic Heisenberg chain

In this appendix, we compute $\chi^{(3)}_{+--+}$ of the antiferromagnetic
Heisenberg chain ($J_{z}=J_{\perp}>0$ in Eq. (1)) by using the equation of
motion of the lattice model [34, 35].

The equation of motion for the total magnetization $M^{\pm}=M^{x}\pm iM^{y}$
reads:

$\displaystyle i\frac{\partial M^{\pm}}{\partial
t}=[M^{\pm},H]=[M^{\pm},-BM^{z}]=\pm BM^{\pm}.$ (66)

The second equality follows from the fact that the Heisenberg interaction,
being $SU(2)$ symmetric, commutes with $M^{\pm}$. Solving the above, we find:

$\displaystyle M^{\pm}(t)=e^{\mp iBt}M^{\pm}.$ (67)

Substituting the above into the Kubo formula for $\chi^{(3)}_{+--+}$, and
using the spin commutation relation, we obtain:

$\displaystyle\chi^{(3)}_{+--+}=-4im\Theta(t)\Theta(t_{\mathrm{w}})\Theta(\tau)e^{-iB(t-\tau)},$
(68)

where $m$ is the magnetization density. Therefore, for the Heisenberg chain,
the nonlinear magnetic susceptibility shows oscillation with constant
magnitude, which reflects the Larmor precession of the total magnetization
vector. We stress that this simplicity is unique to the spatially uniform
(momentum $q=0$) response function. By contrast, the $q\neq 0$ magnetic
response functions may exhibit rich physics that is beyond the Larmor
precession. For instance, the dispersion relation of the magnetic resonance
mode, emanating from the Larmor frequency $\omega=B$ at $q=0$, acquires a
curvature thanks to the spinon interactions [37, 38]. These effects are not
visible in the $q=0$ response function.

## References

  * Mukamel [1995] S. Mukamel,  _Principles of Nonlinear Optical Spectroscopy_ (Oxford University Press, Oxford, 1995).
  * Hamm and Zanni [2011] P. Hamm and M. Zanni,  _Concepts and Methods of 2D Infrared Spectroscopy_ (Cambridge University Press, Cambridge, 2011).
  * Cundiff and Mukamel [2013] S. Cundiff and S. Mukamel, Optical multidimensional coherent spectroscopy, Physics Today 66, 44 (2013).
  * Woerner  _et al._ [2013] M. Woerner, W. Kuehn, P. Bowlan, K. Reimann, and T. Elsaesser, Ultrafast two-dimensional terahertz spectroscopy of elementary excitations in solids, New Journal of Physics 15, 025039 (2013).
  * 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).
  * Mahmood  _et al._ [2020] F. Mahmood, D. Chaudhuri, S. Gopalakrishnan, R. Nandkishore, and N. P. Armitage, Observation of a marginal Fermi glass using THz 2D coherent spectroscopy (2020), arXiv:2005.10822 [cond-mat.str-el] .
  * 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).
  * 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).
  * Parameswaran and Gopalakrishnan [2020] S. A. Parameswaran and S. Gopalakrishnan, Asymptotically exact theory for nonlinear spectroscopy of random quantum magnets, Phys. Rev. Lett. 125, 237601 (2020).
  * Nandkishore  _et al._ [2020] R. M. Nandkishore, W. Choi, and Y. B. Kim, Spectroscopic fingerprints of gapped quantum spin liquids, both conventional and fractonic (2020), arXiv:2010.07947 [cond-mat.str-el] .
  * Phuc and Trung [2020] N. T. Phuc and P. Q. Trung, Direct and ultrafast probing of quantum many-body interaction and mott-insulator transition through coherent two-dimensional spectroscopy (2020), arXiv:2009.08598 [cond-mat.str-el] .
  * Kurnit  _et al._ [1964] N. A. Kurnit, I. D. Abella, and S. R. Hartmann, Observation of a Photon Echo, Phys. Rev. Lett. 13, 567 (1964).
  * Hahn [1950] E. L. Hahn, Spin Echoes, Phys. Rev. 80, 580 (1950).
  * Babujian  _et al._ [2016] H. M. Babujian, M. Karowski, and A. M. Tsvelik, Probing strong correlations with light scattering: Example of the quantum Ising model, Phys. Rev. B 94, 155156 (2016).
  * Affleck [1988] I. Affleck, Field theory methods and quantum critcal phenomena, in  _Les Houches Summer School in Theoretical Physics: Fields, Strings, Critical Phenomena_ (1988) pp. 563–640.
  * Giamarchi [2003] T. Giamarchi,  _Quantum Physics in One Dimension_ (Oxford University Press, Oxford, 2003).
  * Pham  _et al._ [2000] K.-V. Pham, M. Gabay, and P. Lederer, Fractional excitations in the luttinger liquid, Phys. Rev. B 61, 16397 (2000).
  * Zvezdin and Kotov [1997] A. Zvezdin and V. Kotov,  _Modern Magnetooptics and Magnetooptical Materials_ (Taylor & Francis, New York, 1997).
  * Boyd [2008] R. Boyd,  _Nonlinear Optics_ , 3rd ed. (Academic Press, Cambridge, Massachusetts, 2008).
  * Luther and Peschel [1975] A. Luther and I. Peschel, Calculation of critical exponents in two dimensions from quantum field theory in one dimension, Phys. Rev. B 12, 3908 (1975).
  * Haldane [1980] F. D. M. Haldane, General Relation of Correlation Exponents and Spectral Properties of One-Dimensional Fermi Systems: Application to the Anisotropic $S=\frac{1}{2}$ Heisenberg Chain, Phys. Rev. Lett. 45, 1358 (1980).
  * Kubo [1957] R. Kubo, Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems, Journal of the Physical Society of Japan 12, 570 (1957).
  * Breunig  _et al._ [2013] O. Breunig, M. Garst, E. Sela, B. Buldmann, P. Becker, L. Bohatý, R. Müller, and T. Lorenz, Spin-$\frac{1}{2}$ $XXZ$ Chain System ${\mathrm{Cs}}_{2}{\mathrm{CoCl}}_{4}$ in a Transverse Magnetic Field, Phys. Rev. Lett. 111, 187202 (2013).
  * Tavora and Mitra [2013] M. Tavora and A. Mitra, Quench dynamics of one-dimensional bosons in a commensurate periodic potential: A quantum kinetic equation approach, Phys. Rev. B 88, 115144 (2013).
  * Tavora  _et al._ [2014] M. Tavora, A. Rosch, and A. Mitra, Quench dynamics of one-dimensional interacting bosons in a disordered potential: Elastic dephasing and critical speeding-up of thermalization, Phys. Rev. Lett. 113, 010601 (2014).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Two soluble models of an antiferromagnetic chain, Annals of Physics 16, 407 (1961).
  * McCoy  _et al._ [1971] B. M. McCoy, E. Barouch, and D. B. Abraham, Statistical mechanics of the $\mathrm{XY}$ model. iv. time-dependent spin-correlation functions, Phys. Rev. A 4, 2331 (1971).
  * Derzhko and Krokhmalskii [1998] O. Derzhko and T. Krokhmalskii, Numerical Approach for the Study of the Spin-$1/2$ XY Chains Dynamic Properties, Physica Status Solidi (B) 208, 221 (1998).
  * Derzhko  _et al._ [2000] O. Derzhko, T. Krokhmalskii, and J. Stolze, Dynamics of the spin $1/2$ isotropic XY chain in a transverse field, Journal of Physics A: Mathematical and General 33, 3063 (2000).
  * Maeda and Oshikawa [2003] Y. Maeda and M. Oshikawa, Numerical analysis of electron-spin resonance in the spin-$\frac{1}{2}$ $\mathrm{XY}$ model, Phys. Rev. B 67, 224424 (2003).
  * Sachdev and Young [1997] S. Sachdev and A. P. Young, Low temperature relaxational dynamics of the ising chain in a transverse field, Phys. Rev. Lett. 78, 2220 (1997).
  * Damle and Sachdev [1998] K. Damle and S. Sachdev, Spin dynamics and transport in gapped one-dimensional heisenberg antiferromagnets at nonzero temperatures, Phys. Rev. B 57, 8307 (1998).
  * Damle and Sachdev [2005] K. Damle and S. Sachdev, Universal relaxational dynamics of gapped one-dimensional models in the quantum sine-gordon universality class, Phys. Rev. Lett. 95, 187201 (2005).
  * Oshikawa and Affleck [1999] M. Oshikawa and I. Affleck, Low-Temperature Electron Spin Resonance Theory for Half-Integer Spin Antiferromagnetic Chains, Phys. Rev. Lett. 82, 5136 (1999).
  * Oshikawa and Affleck [2002] M. Oshikawa and I. Affleck, Electron spin resonance in $S=\frac{1}{2}$ antiferromagnetic chains, Phys. Rev. B 65, 134410 (2002).
  * Fava  _et al._ [2021] M. Fava, S. Biswas, S. Gopalakrishnan, R. Vasseur, and S. A. Parameswaran, Hydrodynamic non-linear response of interacting integrable systems (2021), arXiv:2103.06899 [cond-mat.str-el] .
  * Keselman  _et al._ [2020] A. Keselman, L. Balents, and O. A. Starykh, Dynamical Signatures of Quasiparticle Interactions in Quantum Spin Chains, Phys. Rev. Lett. 125, 187201 (2020).
  * Kohno [2009] M. Kohno, Dynamically Dominant Excitations of String Solutions in the Spin-$1/2$ Antiferromagnetic Heisenberg Chain in a Magnetic Field, Phys. Rev. Lett. 102, 037203 (2009).
  * Oshikawa and Affleck [1997] M. Oshikawa and I. Affleck, Field-Induced Gap in $\mathit{S}=1/2$ Antiferromagnetic Chains, Phys. Rev. Lett. 79, 2883 (1997).
  * Bocquet  _et al._ [2001] M. Bocquet, F. H. L. Essler, A. M. Tsvelik, and A. O. Gogolin, Finite-temperature dynamical magnetic susceptibility of quasi-one-dimensional frustrated spin-$\frac{1}{2}$ Heisenberg antiferromagnets, Phys. Rev. B 64, 094425 (2001).
  * Gangadharaiah  _et al._ [2008] S. Gangadharaiah, J. Sun, and O. A. Starykh, Spin-orbital effects in magnetized quantum wires and spin chains, Phys. Rev. B 78, 054436 (2008).
  * von Delft and Schoeller [1998] J. von Delft and H. Schoeller, Bosonization for beginners — refermionization for experts, Annalen der Physik 7, 225 (1998).
  * Hikihara and Furusaki [2004] T. Hikihara and A. Furusaki, Correlation amplitudes for the spin-$\frac{1}{2}$ $\mathrm{XXZ}$ chain in a magnetic field, Phys. Rev. B 69, 064427 (2004).
  * Teo and Kane [2014] J. C. Y. Teo and C. L. Kane, From Luttinger liquid to non-Abelian quantum Hall states, Phys. Rev. B 89, 085101 (2014).