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

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, twodimensional 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 * yuan.wan@iphy.ac.cn response χ (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 χ (3) response ( Fig. 1(a)). Let τ be the time delay between the second and first pulses, t 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 ≈ τ 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 τ [1,2]. Had the dynamics been unitary, the quantum mechanical phase accumulated in τ would be completely removed when t = τ , resulting in a perfect rephasing. This perfect rephasing process would produce a photon echo signal at t = τ regardless the value of τ . 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 τ is a manifestation of the decoherence time (T 2 time), whereas the decay of the signal as a function of t 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 ⊥ 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, χ +−−+ and its complex conjugate χ (3) −++− show photon echo, which appears as a peak on the t axis at t ≈ τ . 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 − τ rather than t and τ both, resembling the perfectly rephasing photon echo in a few-body system. This implies that the signal measured at a given value of t − τ is independent of τ . Moreover, the echo signal depends weakly on the waiting time t w , and saturates when t w → ∞.
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 leftmoving 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 = τ , 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 τ or t 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 RGirrelevant 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 τ or t w exceeds the lifetime of the wave packets. This is manifest as the decay of the echo signal as a function of τ or t 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 τ 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 = τ , decays as a stretched exponential ∼ exp(−C(τ /τ disp ) 1/2 ), where C is a numerical constant. Meanwhile, the echo signal shows weak dependence on the waiting time t w and saturates when t w → ∞. We attribute the lack of t w dependence to the absence of thermalization in the toy model -The decay of the photon echo as a function of t 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 χ (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 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 τ , the second and the third by tw, 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 ≈ τ .
in Sec. V. In Sec. VI, we point out a few interesting open problems.

II. SETUP
We consider the S = 1/2 XXZ spin chain: j labels the lattice sites. S ± j , S z are the S = 1/2 spin operators. J ⊥ is the exchange constant in the spin xy plane. We shall consider both "ferromagnetic" (J ⊥ < 0) and "antiferromagnetic" (J ⊥ > 0) chains. J z is the exchange constant in the spin z axis. We include the Zeeman term due to an external field B z (Fig. 1). Throughout this work, we use the natural units with = k B = µ = 1 where µ 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: χ 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, τ , τ + t w , where τ, t 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: Here, χ  +−−+ 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 χ +−−+ (t, t w , τ ) by holding t w constant and scanning τ and t. We obtain the two-dimensional spectra χ (3) (ω t , t w , ω τ ) by performing a two-dimensional one-sided Fourier transform of χ (3) (t, t w , τ ) over the do-main t > 0 and τ > 0. Note alternative protocols for visualizing the nonlinear response exist [4,6]. Ours is closely related to that of Refs. 5 and 7.

III. BOSONIZATION
In this section, we compute the nonlinear response of Eq. (1) by using bosonization. We show that the nonlinear susceptibility χ (3) +−−+ exhibits photon echo and characterize its features.

A. 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: Here, u is the speed of sound. K is the Luttinger parameter. They may be computed from the microscopic model parameters J ⊥ , J z , B by using the Bethe ansätz [20,21]. φ and θ are boson fields with the compactification conditions: They obey the non-local commutation relation: where Θ(·) is the Heaviside step function. Note there is freedom in choosing the commutation relation between φ and θ. 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: 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. (6) the spatially staggered (∝ (−) 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 ⊥ < 0) and the antiferromagnetic chain (J ⊥ > 0) are exactly mapped to each other by a π-rotation of spins about the z-axis at every other site, i.e. S ± j → (−1) j S ± 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. (6).

B. Four-point response function
The experimentally measured signal is related to the spacetime-dependent nonlinear magnetic susceptibility χ (2)). Its Kubo formula [22] reads: 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. χ +−−+ measures the system's response at the spacetime origin due to successive perturbations at From Eq. (6), 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: where 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 We also impose the charge neutrality condition j m j = j n j = 0 to ensure G R does not vanish in the thermodynamic limit.
Eq. (8) 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 ] = −2sinh [A,B] 2 e i(A+B) , where A, B are arbitrary linear combinations of φ and θ. This permits a straightforward evaluation of the nested commutators in Eq. (8). We then compute the thermal average by using exp(iA) = exp(− A 2 /2), where A is an arbitrary linear combination of θ and φ. We obtain: The above is the main result of this subsection. Here, α ij and C ij are defined for the vertex operator V i and V j . α ij comes from the commutator of vertex operators: We have used light cone coordinates is the sign function. l j and r j are real parameters related to m j , n j through: Meanwhile, where 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. (11) is causal and under what conditions. In Appendix B, we show that Eq. (11) is causal provided that the vertex operators are local, i.e. Eq. (10) holds for all V j . It is straightforward to check that the vertex operators that appear in the expression of S ± (x) (Eq. (6)) indeed fulfill this condition. Thus, χ +−−+ is causal as expected.
The next step is to find χ 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 χ +−−+ valid when t, t w , τ are large. To this end, we use the following approximation for the function C ij (Eq. (12c)): We have omitted a cutoff dependent prefactor. Eq. (14) 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 = ±ut ij . The latter is short distance physics and should not affect the asymptotic behavior. The validity of Eq. 14 will be further verified a posteriori by numerical integration.
With the approximation Eq. (14), 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 τ increases, which we discard as they are uninteresting for our purpose. The second group of terms show signature of photon echoas we fix τ and scan t, they exhibit a maximum near t ≈ τ . Retaining these terms, we find: The above is the key result of this subsection. Eq. (15) suggests that the photon echo is perfect, i.e. it is a function of t − τ . The time scale of the photon echo signal is set by 2K/(πT ). Moreover, it is independent of the waiting time t 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. (15) by performing the integration Eq. (2) numerically. We use a short distance cutoff π T /u = 1, which smoothens the algebraic singularities and discontinuities that would otherwise appear in the integrand in the limit → 0 + . With finite , 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 πT R/u = 2 with the relative error O(10 −3 ).
+−−+ for representative Luttinger parameters K = 0.7, K = 1, and K = 1.5. πT t w = 1 and πT τ = 20, 40. For all cases, the nonlinear response shows a surge near t ≈ τ , exhibiting the clear signature of photon echo. Note the maximum is not exactly located at t = τ but fairly close to it. The data for different val-ues of τ overlay within numerical error when plotted as a function of t − τ . This demonstrates the photon echo is independent of τ for large τ .
Numerical integration also indicates that the value of χ (3) +−−+ measured at t = τ does not depend on t w within numerical error in the asymptotic regime πT τ 1, which is in agreement with Eq. (15) We further examine the asymptotic behavior of χ (15) suggests the product would be ∝ τ − t when t < τ , and ∝ −(t − τ ) 3 when t > τ . The insets of Fig. 2 show that its behavior is in excellent agreement with Eq. (15).
Having analyzed the photon echo signal for fixed value of τ , we now scan τ and present the nonlinear response as a function of both τ and t. Fig. 3(a) shows χ (3) +−−+ for K = 1 and πT t w = 1. The photon echo appears as a bright feature at the diagonal of the (t, τ ) 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 limitedthe photon echo signal in the time domain ( Fig. 3(a)) is independent of t + τ at late time, and hence its Fourier transform ∼ δ(ω t + ω τ ).
To recapitulate, the χ +−−+ of the ferromagnetic chain is real and independent of the magnetization density. It exhibits perfect photon echo, which depends on t − τ rather than t and τ both.

D. AFM chain
We turn to the antiferromagnetic (J ⊥ > 0) case in this subsection. We write the spin operator as: where we have defined vertex operators a and b for later convenience. Inserting the above into Eq. (7) yields: We have defined a set of response functions with the form of Eq. (8). For instance, G R a † aaa † is defined as: The other response functions are defined in the same vein. The expression for these response function can be read off from Eq. (11) by plugging in appropriate values of l 1,2,3 and r 1,2,3 : l = √ K +1/(2 √ K) and r = − √ K +1/(2 √ 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 † aab † ) as they vanish in the thermodynamic limit. The calculation of χ +−−+ parallels the ferromagnetic case (Sec. III C). Note, however, χ +−−+ is now complex and depends on the magnetic field through the magnetization density m. We find only G a † aaa † and G b † bbb † contribute to photon echo. Using the approximation (Eq. (14)), we obtain after lengthy calculation: Here, we have defined a parameter ∆ = 2K + 1/(2K). The above is the key result of this subsection.
Eq. (19) shows that the χ +−−+ 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/((∆ − 2)πT ). In particular, in the Heisenberg limit where K → 1/2 (∆ → 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. (19) by numerical integration. Fig. 4 show χ To test the asymptotic behavior of χ +−−+ , we multiply its complex modulus with exp(π(∆ − 2)|t − τ |). Eq. (19) suggests the product shows linear behavior for t > τ and approaches a constant for t < τ . The insets of Fig. 4 show good agreement with Eq. (19). For K = 0.7 and πT τ = 50, the data deviate slightly from the constant behavior; we think this occurs because τ is not sufficiently large to suppress the non-asymptotic contributions.
We then present the dependence of the photon echo signal on both τ and t. Fig. 5(a)&(b) show the real and imaginary parts of χ The photon echo appears as the bright feature persists along the diagonal direction (t = τ ). Performing Fourier transform over t and τ 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 ω t = 2πmu and ω τ = −2πmu as suggested by Eq. (19). 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 χ

IV. SPINON LENSING PICTURE
In the previous section, we have shown that the nonlinear magnetic susceptibility χ +−−+ 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θ). Let us consider the time evolution after the operator exp(−iθ(0, 0)) at t = x = 0 on the initial state |Ψ(0) at t = 0 which is assumed to be an energy eigenstate. At time t, the state is given by where we used the assumption that Ψ(0) is an energy eigenstate and dropped the overall phase factor. We can decompose the field θ to chiral components as where The equation of motion implies that each chiral components depends only on the corresponding light cone coordinate, as in Eq. (21). Using this, we can translate the time dependence into the position dependence, so that where θ L,R (∓ut) is now understood as a static operator at locations ∓ut, up to the overall phase factor which includes the one coming from the commutation relation between θ R and θ L . Now, recall the equal-time commutation relation: Namely, exp(−iθ L,R (y)) creates a kink of step π/2 in φfield at the location y. Since this step corresponds to the magnetization density (1/2) δ(x − y), this kink can be interpreted as a spinon [17]. Eq. (23) then implies that the application of the operator exp(−iθ) 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θ L,R ) can be viewed as spinon creation operators of the corresponding chirality. This justifies the following simple visual picture: Applying exp(−iθ(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 ±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) 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 antispinons) 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 ± 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θ(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. (25) maximizes the first factor, which is formally divergent in the scaling limit (the cutoff → 0 + ) but bounded to be a finite constant due to the finite in a realistic system. However, x + 1 = 0 implies x − 1 = 2ut 1 , leading to the exponential suppression of the other factor in Eq. (25). 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. (13) as a sum over various four-point correlation functions by expanding the nested commutators: where R i is related to R i by reversing the order of operators in the product. Following Ref. 1 and 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. 25), here, we can catch and annihilate all the spinons created earlier by choosing (Fig. 6c): 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θ)! 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. (27). 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. (27) 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 ± ut 3 and two at 0, and antispinon creation operators at locations −x 2 ± ut 2 and −x 1 ± ut 1 (Fig. 6b). Under the condition Eq. (27), 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 χ (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. (27). Fig. 6d shows the behavior of χ (3) +−−+ as a function of the detection position in a representative lensing configuration. Here we choose πT t = πT τ = 10 and πT t w = 10. It shows that the response is maximum at the "focal point" (0, 0). As we have used a finite short distance cutoff 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 to suppress the effect of the regularizer. Now, the optical nonlinear response χ  17)). Among these, G R a † aaa † and G R b † bbb † 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 χ 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 temper-atures. 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φ) [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 τ or t w exceeds the life time of these excitations. This, in turn, will be manifest as the decay of the photon echo signal as τ or t w increases, in close analogy with the dissipationinduced photon echo decay in few-body systems [1,2]. Straightforward dimensional analysis suggests the decay rate vanishes as T 1+η as T → 0 where η is the dimension of the umklapp term.
On the other hand, there also exist a different family of RG-irrelevant perturbations, such as (∇ 2 φ) 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 τ 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
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: Here, a is the lattice constant. θ n resides on the lattice site labeled by n, whereas φ n+1/2 resides on the midpoint of the lattice link connecting the site n and n + 1. Their commutation relation is given by: [φ n+1/2 , θ n ] = −iπ Θ(n − n + 1/2). It is sufficient for our purpose to consider the ferromagnetic chain where S ± n = exp(∓iθ n ). We set the Luttinger parameter K = 1.
The phonon modes are exact eigenstates of Eq. (28) owing to its quadratic form. The phonon dispersion relation is now nonlinear: ω q = 2u/a | sin(qa/2)|. Therefore, Eq. (28) 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 χ +−−+ as a function of t and τ at temperature T = 1.2u/a. The waiting time t 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, τ ) 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 disper-sion near q = 0, we obtain ω q = u(|q| − a 2 |q| 3 /24 + · · · ). 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 τ disp : Here, ξ 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 τ disp . Our data support this hypothesis. Fig. 7(c) shows the semi-log plot of χ (3) +−−+ at t = τ as a function of (τ /τ 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(τ /τ disp ) ν ), where the exponent ν = 1/2 and C is a numeric constant. The decay rate 1/τ disp ∼ T 3 , consistent with the dimension η = 2 of the higher derivative term (∇ 2 φ) 2 . The origin of the stretching exponent ν = 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 τ , we turn to its t w dependence. Fig. 7(e) shows χ (3) +−−+ at t = τ as a function of t w for two representative values of τ . The echo signal initially decrease as t 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 χ +−−+ saturates to finite value as t w → ∞ at the lensing configuration Eq. (27). In a few-body system, the decay of photon echo as a function of t w reflects the thermal relaxation of the optical population. Here, as the phonon modes are exact eigenstates of the Hamiltonian Eq. (28), the phonon population cannot relax. We therefore heuristically attribute the saturation to the absence of thermal relaxation in Eq. (28).
To recapitulate, by a numerical experiment, we have shown that the photon echo signal decays as a function of τ in the presence of wave packet dispersion. The decay is controlled by a dispersion time scale τ disp . Meanwhile, the echo signal saturates as t w → ∞, 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 χ −++− 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, Cs 2 CoCl 4 is thought to be a realization of the easy-plane antiferromagnetic (J ⊥ > 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 χ (3) response with linear polarization (e.g. χ (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 τ increases and eventually disappears when τ far exceeds the dispersion time scale. However, the echo signal does not disappear when the waiting time t w → ∞. This points toward different uses of the pulse delays τ and t 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 treat-ment [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 = 0 due to spinon interactions [37,38]. In light of these results, we expect that the q = 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 = 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.
In this appendix, we discuss the different choices for the commutation relation between φ and θ 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, φ and ∇θ form a pair of canonical conjugate variables: The above does not uniquely fix the commutator between φ and θ. Several popular choices exist in literature. In this work, we use: We term this choice the "Heaviside prescription" for latter convenience. A closely related variant is [φ(x), θ(y)] = iπΘ(y − x) [15]. Since the two variants are rather similar, we shall focus the prescription Eq. (A2). Another popular choice is [16]: 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 noninteracting fermions (Luttinger parameter K = 1), the left and right chiral boson fields are dynamically coupled: The commutation relation of φ L , and likewise φ R , are identical for both prescriptions. However, the commutator between φ L and φ R depends on the prescription: 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: anti-commute thanks to the commutator between φ L and φ R . By contrast, in the signum prescription, Klein factors are needed to ensure the anti-commutation relation: where η L,R are the Klein factors obeying the Clifford algebra: η 2 L = η 2 R = 1, and η L η R = −η R η 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 = 1) does not spoil the commutator between the left and right chiral boson fields. In this case, the two dynamically decoupled fields read: and we find Eq. (A5) 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: where c n is the fermion annihilation operator on lattice site n. In particular, We now bosonize the string operator following the standard procedure [15,16]: In the Heaviside prescription, the boundary term φ(−∞) commutes with θ(x). We thus may regard it as a number and omit it: By contrast, in the signum prescription, we must keep the boundary term as [φ(−∞), θ(x)] = 0: In the signum prescription, the boundary term φ(−∞) is needed to reproduce the commutation relation between the fermion field operator and the string operator (Eq. A9). Had we dropped the boundary term, we would have found: which disagrees with Eq. A9. 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 ⊥ < 0) without loss of generality. Following Refs. 15 and 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: .
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: using the signum prescription. This is incorrect because these two operators must commute. For the uniform components, we have: which is independent of the prescription, and M − j ∼ e iθ cos(2φ) (H.) η L e i(θ+2φ−φ(−∞)) + η R e i(θ−2φ+φ(−∞)) (S.) .

(A17)
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. (A13), (A14), and (A17) 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 fourpoint response function: 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. (A15)). As a result, G R = 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 θ and φ 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.
Combining the above results, we conclude that the response function G R (Eq. (8)) is causal. We note that one could also verify the causality of G R by using its explicit expression Eq. (11).

Appendix C: Antiferromagnetic Heisenberg chain
In this appendix, we compute χ The second equality follows from the fact that the Heisenberg interaction, being SU (2) symmetric, commutes with M ± . Solving the above, we find: Substituting the above into the Kubo formula for χ +−−+ , and using the spin commutation relation, we obtain: 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 = 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 ω = 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.