Resolving continua of fractional excitations by spinon echo in THz 2D coherent spectroscopy

TL;DR

This paper demonstrates that terahertz 2D coherent spectroscopy can reveal sharp features and spinon echoes in fractionalized spin systems, providing new insights into their excitations and lifetimes beyond traditional methods.

Contribution

The study introduces the concept of spinon echo detection via 2D spectroscopy, enabling access to information on fractional excitations and their lifetimes in spin systems.

Findings

2D spectra show sharp features despite broad linear response continua.

Spinon echoes reveal the lifetime of two-spinon excited states.

Disorder and finite lifetime effects are distinguishable in 2D spectra.

Abstract

We show that the new technique of terahertz 2D coherent spectroscopy is capable of giving qualitatively new information about fractionalized spin systems. For the prototypical example of the transverse field Ising chain, we demonstrate theoretically that, despite the broad continuum of excitations in linear response, the 2D spectrum contains sharp features that are a coherent signature of a `spinon echo', which gives previously inaccessible information such as the lifetime of the two-spinon excited state. The effects of disorder and finite lifetime, which are practically indistinguishable in the linear optical or neutron response, manifest in dramatically different fashion in the 2D spectra. Our results may be directly applicable to model quasi-1D transverse field Ising chain systems such as CoNb$_2$O$_6$, but the concept can be applied to fractionalized spin systems in general.

Equations

$$H=-J(\sum _{n=1}^{L-1}{\sigma }_n^z{\sigma }_{n+1}^z+\eta {\sigma }_L^z{\sigma }_1^z)-h\sum _n{\sigma }_n^x,$$

$$M_{\mathrm{NL}}^{\alpha }(t,\tau )$$

$$+(A_0^{\beta }{)}^2{A}_{\tau }^{\beta }{\chi }_{\alpha \beta \beta \beta }^{(3)}(t,\tau +t,\tau +t)$$

$$+A_0^{\beta }(A_{\tau }^{\beta }{)}^2{\chi }_{\alpha \beta \beta \beta }^{(3)}(t,t,\tau +t)+O(A^4).$$

$${\chi }_{xx}^{(1)}(t)=\frac{2\theta (t)}{L}\sum _{k>0}{\sin }^2{\theta }_k{e}^{-t/T_2}\sin (2{\lambda }_{k}t),$$

$${\chi }_{xxx}^{(2)}(t,\tau +t)=\frac{4\theta (t)\theta (\tau )}{L}\sum _{k>0}{\sin }^2{\theta }_k\cos {\theta }_k$$

$$\times [e^{-t/T_1}\cos (2{\lambda }_k\tau )-e^{-(t+\tau )/T_2}\cos (2{\lambda }_k(\tau +t))].$$

$$A_k^{(1)}$$

$$\times e^{-(t_1+t_2+t_3)/T_2};$$

$$A_k^{(2)}$$

$$\times e^{-(t_1+t_2)/T_2}{e}^{-t_3/T_1};$$

$$A_k^{(3)}$$

$$\times e^{-(t_1+t_3)/T_2}{e}^{-t_2/T_1};$$

$$A_k^{(4)}$$

$$\times e^{-(t_1+t_3)/T_2}{e}^{-t_2/T_1}.$$

$$B^{\beta }(s)=A_0^{\beta }\delta (s)+A_{\tau }^{\beta }\delta (s-\tau ),$$

$$M_{\mathrm{AB}}^{\alpha }(t+\tau )$$

$$+{\int }_{-\infty }^{\infty }d{s}_1{\int }_{-\infty }^{\infty }d{s}_2{\int }_{-\infty }^{\infty }d{s}_3{\chi }_{\alpha \beta \beta \beta }^{(3)}(t+\tau -s_1,t+\tau -s_2,t+\tau -s_3)B^{\beta }(s_1)B^{\beta }(s_2)B^{\beta }(s_3)+\cdots .$$

$$M_{\mathrm{AB}}^{\alpha }(t+\tau )$$

$$+{\chi }_{\alpha \beta \beta }^{(2)}(t+\tau ,t+\tau )(A_0^{\beta }{)}^2+{\chi }_{\alpha \beta \beta }^{(2)}(t,t+\tau )A_{\tau }^{\beta }{A}_0^{\beta }+{\chi }_{\alpha \beta \beta }^{(2)}(t,t)(A_{\tau }^{\beta }{)}^2$$

$$+{\chi }_{\alpha \beta \beta \beta }^{(3)}(t+\tau ,t+\tau ,t+\tau )(A_0^{\beta }{)}^3+{\chi }_{\alpha \beta \beta \beta }^{(3)}(t,t+\tau ,t+\tau )A_{\tau }^{\beta }(A_0^{\beta }{)}^2$$

$$+{\chi }_{\alpha \beta \beta \beta }^{(3)}(t,t,t+\tau )(A_{\tau }^{\beta }{)}^2{A}_0^{\beta }+{\chi }_{\alpha \beta \beta }^{(3)}(t,t,t)(A_{\tau }^{\beta }{)}^3+\cdots .$$

$$M_{\mathrm{A}}^{\alpha }(t+\tau )$$

$$M_{\mathrm{B}}^{\alpha }(t+\tau )$$

$$M_{\mathrm{NL}}^{\alpha }(t+\tau )$$

$$={\chi }_{\alpha \beta \beta }^{(2)}(t,t+\tau )A_{\tau }^{\beta }{A}_0^{\beta }+{\chi }_{\alpha \beta \beta \beta }^{(3)}(t,t+\tau ,t+\tau )A_{\tau }^{\beta }(A_0^{\beta }{)}^2+{\chi }_{\alpha \beta \beta \beta }^{(3)}(t,t,t+\tau )(A_{\tau }^{\beta }{)}^2{A}_0^{\beta }+\cdots .$$

$$H=-J\sum _{n=1}^{L-1}{\sigma }_n^z{\sigma }_{n+1}^z-J{\sigma }_L^z{\sigma }_1^z-h\sum _n{\sigma }_n^x.$$

$${\sigma }_n^z=(c_n^{†}+c_n^{})\exp (i\pi \sum _{j=1}^{n-1}{c}_j^{†}{c}_j);{\sigma }_n^x=2c_n^{†}{c}_n^{}-1,$$

$$H=-J\sum _n(c_n^{†}-c_n^{})(c_{n+1}^{†}+c_{n+1}^{})-h\sum _n(c_n^{†}{c}_n^{}-c_n^{}{c}_n^{†}).$$

$$H=\sum _{k>0}(c_k^{†},c_{-k})\left(\begin{array}{cc}{ϵ}_k& i{\Delta }_k\\ -i{\Delta }_k& -{ϵ}_k\end{array}\right)\left(\begin{array}{c}{c}_k\\ c_{-k}^{†}\end{array}\right),$$

$$\left(\begin{array}{c}{c}_k\\ c_{-k}^{†}\end{array}\right)=\left(\begin{array}{cc}\cos \frac{{\theta }_k}{2}& -i\sin \frac{{\theta }_k}{2}\\ -i\sin \frac{{\theta }_k}{2}& \cos \frac{{\theta }_k}{2}\end{array}\right)\left(\begin{array}{c}{d}_k\\ d_{-k}^{†}\end{array}\right),$$

$$H=\sum _{k>0}{\lambda }_k(d_k^{†}{d}_k^{}-d_{-k}^{}{d}_{-k}^{†})=\sum _k{\lambda }_k(d_k^{†}{d}_k-\frac{1}{2}).$$

$$M^x=\frac{1}{2}\sum _n{\sigma }_n^x=\sum _{k>0}{m}_k^x,$$

$$m_k^x$$

$$=(d_k^{†},d_{-k})\left(\begin{array}{cc}\cos \frac{{\theta }_k}{2}& i\sin \frac{{\theta }_k}{2}\\ i\sin \frac{{\theta }_k}{2}& \cos \frac{{\theta }_k}{2}\end{array}\right)\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\left(\begin{array}{cc}\cos \frac{{\theta }_k}{2}& -i\sin \frac{{\theta }_k}{2}\\ -i\sin \frac{{\theta }_k}{2}& \cos \frac{{\theta }_k}{2}\end{array}\right)\left(\begin{array}{c}{d}_k\\ d_{-k}^{†}\end{array}\right)$$

$$=(d_k^{†},d_{-k})\left(\begin{array}{cc}\cos {\theta }_k& -i\sin {\theta }_k\\ i\sin {\theta }_k& -\cos {\theta }_k\end{array}\right)\left(\begin{array}{c}{d}_k\\ d_{-k}^{†}\end{array}\right).$$

$${\tau }_k^x\equiv d_{-k}{d}_k+d_k^{†}{d}_{-k}^{†},{\tau }_k^y\equiv i(d_{-k}{d}_k-d_k^{†}{d}_{-k}^{†}),{\tau }_k^z\equiv d_k^{†}{d}_k-d_{-k}{d}_{-k}^{†}.$$

$$H=\sum _{k>0}{\lambda }_k{\tau }_k^z,$$

$$m_k^x=\cos {\theta }_k{\tau }_k^z+\sin {\theta }_k{\tau }_k^y.$$

Full text

In the recent years wholly new classes of condensed matter systems have become of intense interest. Topological materials, quantum spin liquids, and strange metals are characterized by Berry phase effects, fractional excitations, and highly entangled ground states Balents (2010); Xiao et al. (2010); Savary and Balents (2016); Armitage et al. (2018). However, we can measure many of their correlations only imperfectly with existing tools. A promising direction is nonlinear response that has been used to characterize the symmetry of semiconductors Shen (1984) and magnets Fiebig et al. (2005), Berry phase in topological semimetals Sipe and Shkrebtii (2000); Morimoto and Nagaosa (2016); Moore and Orenstein (2010), and exotic ground states and excitations in correlated systems Zhao et al. (2016, 2015); Babujian et al. (2016).

For quantum spin liquids, one of their most remarkable properties is the emergence of fractional particles, known as spinons, that may be understood as carrying half a conventional spin degree of freedom. Spinons present a challenge for conventional spectroscopy as they must be excited in pairs. This typically leads to a broad continuum spectrum that represent a convolution of all possible ways that energy and momentum can be shared between two spinons. In conventional linear magnetic susceptibility $\chi^{(1)}(\omega)$ of a spin chain Morris et al. (2014) light excites a pair of spinons with opposite momenta (Fig. 1a). Each pair gives rise to a peak in the absorption spectrum $\mathrm{Im}\chi^{(1)}$ centered at the frequency $\omega=\lambda_{\mathbf{k}}+\lambda_{-\mathbf{k}}$, where $\lambda_{\mathbf{k}}$ is the dispersion relation of the spinon. As there are infinitely many such pairs, the absorption peaks congest the frequency axis, resulting in a broad continuum (Fig. 1b, top). While the broad continuum seen with terahertz (THz) optical spectroscopy and neutron scattering has reasonably been taken as evidence for spinons in spin chains Lake et al. (2005); Morris et al. (2014); Coldea et al. (2010), the situation is less straightforward in higher dimensional spin liquid candidates e.g. 2D Kitaev materials, herbertsmithite, and triangular lattices Han et al. (2012); Banerjee et al. (2017); Shen et al. (2016); Paddison et al. (2016); Zhang et al. (2018); Zhu et al. (2017). In such cases, the relative importance of finite lifetime and disorder and even fractionalization itself is unclear. In all cases, the intrinsic spectral properties of spinons such as the line width and shape are hidden in the continuum.

In this work, we show that the new technique of THz two-dimensional coherent spectroscopy (2DCS) Kuehn et al. (2011); Woerner et al. (2013) can provide qualitatively new information on the dynamical properties of spinons. We explore our ideas in the context of the simplest minimal model for fractionalization – the transverse field Ising chain (TFIC) – but the possibilities are more general. In the optical and radio frequency range Aue et al. (1976); Cundiff and Mukamel (2013); Hamm and Zanni (2011); Mukamel (1999) 2DCS is an established technique that probes nonlinear susceptibilities. Thanks to recent technical advances that enable table-top high intensity THz sources, it has been extended recently to the THz range to study graphene and quantum wells Kuehn et al. (2011); Woerner et al. (2013), molecular rotations Lu et al. (2016), and spin waves in the conventional magnet YFeO3 Lu et al. (2017). THz 2DCS uses two pulses in a collinear geometry to excite a system at one frequency and detect at another, thus producing a 2D spectrum. Applications of 2DCS include quantifying nonlinear couplings between excitations and – relevant to the present work – separating inhomogeneous and homogeneous broadening Cundiff and Mukamel (2013); Hamm and Zanni (2011); Mukamel (1999). A similar mechanism will allow for the resolution of the spinon continuum in the 2D spectrum (Fig. 1b, bottom), where congestion occurs along the 2D spectrum’s diagonal, but the intrinsic line width of each spinon pair is revealed by the spectral width along the anti-diagonal.

The TFIC Hamiltonian is Pfeuty (1970):

[TABLE]

Here $\sigma^{x,y,z}_{n}$ are Pauli matrices, $J>0$ is the ferromagnetic exchange, $h>0$ is the transverse field, and $L$ is the chain length. We shall use periodic ($\eta=1$) and open ($\eta=0$) boundary conditions as they suit our purposes. Macroscopic response functions are independent of such choices. This system admits a two-fold degenerate ferromagnetic (FM) ground state for $h<J$ and a single paramagnetic (PM) ground state for $h>J$. While strictly speaking the TFIC is not a spin liquid, the domain wall excitations of the FM phase are close analogues of spinons. Henceforth, we use “domain wall” and “spinon” interchangeably.

We consider a setup similar to that used in Ref. Lu et al. (2017). Two linearly polarized magnetic field pulses A and B arrive at the sample at time 0 and $\tau>0$ (Fig. 1a). The magnetization at a time $\tau+t$ along direction $\alpha$, $M^{\alpha}_{\mathrm{AB}}(\tau+t)$, is a convolution of applied field with the sample response Tim . The experiment is then repeated but with pulse A or B alone and the magnetization recorded as $M^{\alpha}_{\mathrm{A}}(\tau+t)$ and $M^{\alpha}_{\mathrm{B}}(\tau+t)$. The nonlinear signal is defined as $M^{\alpha}_{\mathrm{NL}}(t,\tau)=M^{\alpha}_{\mathrm{AB}}(\tau+t)-M^{\alpha}_{\mathrm{A}}(\tau+t)-M^{\alpha}_{\mathrm{B}}(\tau+t)$. The 2D spectrum is the Fourier transform (FT) of $M_{\mathrm{NL}}^{\alpha}(t,\tau)$ over the domain $t>0,\tau>0$.

The nonlinear magnetization $M_{\mathrm{NL}}^{\alpha}(t,\tau)$ is a direct measure of the second and/or third order magnetic susceptibilities. For simplicity, we model the magnetic field as two Dirac-$\delta$ pulses with the same polarization $\beta$, i.e. $B^{\beta}(s)=A^{\beta}_{0}\delta(s)+A^{\beta}_{\tau}\delta(s-\tau)$, where $s$ is time, and $A^{\beta}_{0,\tau}$ the pulse areas. In principle, the polarizations of pulse A and B can be different. The nonlinear signal (See Supplemental Material (SM) Sup ) is,

[TABLE]

Here, we have retained the dominant and sub-dominant contributions. The two $\chi^{(3)}$ terms encode different physical processes. In the first, pulse A couples to the sample at second order whereas pulse B couples at first order. In the second, the contributions of A and B are switched.

We are primarily interested in the spinons in the FM phase at zero temperature, and thus use the representative model parameters $h/(h+J)=0.3$ in the ensuing discussion. Since $\sigma^{x}_{n}$ excites spinon pairs, we focus on the polarization $\alpha=\beta=\hat{x}$. We calculate $\chi^{(2)}_{xxx}$ and $\chi^{(3)}_{xxxx}$ analytically through the following procedure (see SM Sup for details). We map Eq. (1) to free, fermionic spinons by using the Jordan-Wigner transformation Pfeuty (1970). Each pair of spinons with momenta $\pm k$ form a two-level system (TLS), whose ground (excited) state corresponds to the absence (presence) of the said pair. The energy level splitting is $2\lambda_{k}$, where $\lambda_{k}=2\sqrt{J^{2}+h^{2}+2Jh\cos k}$ is the spinon dispersion. As the TLSs formed by different spinon pairs are decoupled, the TFIC is equivalent to an ensemble of independent TLSs, thereby permitting a straightforward calculation of the nonlinear response Mukamel (1999).

In more realistic models, additions to Eq. (1) such as additional exchange interactions and spin-lattice couplings induce spinon interactions, which give effects such as spinon decay. By the TLS analogy, we incorporate these effects phenomenologically through a population time $T_{1}$ and decoherence time $T_{2}$ Shen (1984), which captures the essential physics while maintaining the analytic tractability Sup . We assume $T_{1,2}$ are $k$-independent for simplicity. The ideal TFIC then corresponds to $T_{1,2}=0$.

We begin with the linear susceptibility per site,

[TABLE]

where $\sin\theta_{k}=-2J\sin k/\lambda_{k}$ is the optical matrix element. The summation is over the positive half of the first Brillouin zone (1BZ). Using the above TLS picture, we interpret Eq. (3) as follows. The magnetic field pulse induces optical transitions in all TLS, producing a damped oscillatory signal with frequency $2\lambda_{k}$. The damping coefficient $1/T_{2}$ reflects the spinon decay. Since the frequency takes its value from a dense spectrum given by $k$, dephasing leads to an additional decay of $\chi_{xx}^{(1)}(t)$, which is difficult to distinguish from the intrinsic decay due to finite $T_{2}$. This difficulty of unravelling the dephasing and the intrinsic decay persists in the frequency domain. Here, each spinon pair contributes an absorption peak in $\mathrm{Im}\chi^{(1)}_{xx}$ centered at the energy $2\lambda_{k}$ with width $1/T_{2}$. As $k$ runs over the 1BZ, the peaks form a broad continuum, which disguises the intrinsic line width $1/T_{2}$. Comparing the continuum for $1/T_{2}=0$ (Fig. 2a) and $1/T_{2}=0.2(J+h)$ (Fig. 2f), the difference is merely quantitative.

We then turn to the lowest order nonlinear response:

[TABLE]

The first term on the right hand side of Eq. (4) is non-oscillatory in $t$. In the frequency domain, this gives rise to a peak centered at $\omega_{t}=0$, appearing as the streak along the $\omega_{\tau}$ axis shown in Fig. 2b. Increasing $1/T_{1}$ from 0 leads to broadening of the streak (Fig. 2g). Viewing $\omega_{\tau}$ as the pumping frequency and $\omega_{t}$ the detecting frequency, this streak is a THz rectification (TR) signal Lu et al. (2017). The second term of Eq. (4) is oscillatory in $t+\tau$. Yet, similar to Eq. (3), the dephasing leads to decay, which is further modulated by the intrinsic decay due to $T_{2}$. This results in a diffusive, barely discernible signal in the first frequency quadrant (Fig. 2b,g), which is similar to the non-rephasing (NR) signal usually found in $\chi^{(3)}$ Lu et al. (2017). See SM Sup for detailed discussion of these features.

Qualitatively different physics appears in $\chi^{(3)}_{xxxx}$. It is instructive to consider the more general form that corresponds to a three-pulse process (Fig. 3a): $\chi_{xxxx}^{(3)}(t_{3},t_{2}+t_{3},t_{1}+t_{2}+t_{3})=-(\theta(t_{1})\theta(t_{2})\theta(t_{3})/L)\sum_{k>0}A^{(1)}_{k}+A^{(2)}_{k}+A^{(3)}_{k}+A^{(4)}_{k}$, where

[TABLE]

$A^{(1\sim 4)}_{k}$ encode distinct evolution paths of the density matrix of the spinon pair with momenta $\pm k$ due to the THz pulses. While the forms of $A^{(1,2,3)}_{k}$ resemble that of $\chi^{(2)}_{xxx}$, $A^{(4)}_{k}$ is different in that $t_{1}$ and $t_{3}$ appear with opposite signs. Regardless the oscillation frequency $2\lambda_{k}$, the phase accumulated between the first and the second pulses ($t_{1}$) is cancelled after the third pulse at $t_{3}=t_{1}$. Said differently, the dephasing process during $t_{1}$ is countered by the rephasing process during $t_{3}$. This rephasing process is the incarnation of the photon echo in the context of spinon dynamics. Tracing $A^{(4)}_{k}$ back to its originating density matrix evolution sequence (Fig. 3a), we find the sequence is identical to the photon echo process from a TLS Kurnit et al. (1964); Mukamel (1999). Therefore, we term this process the “spinon echo”.

Photon echo and its analogues are a sensitive diagnostics of dissipation Kurnit et al. (1964); Mukamel (1999). Here, the rephasing signal from the spinon echo allows for a direct measurement of $T_{2}$. To see this, we return to the $\chi_{xxxx}^{(3)}$ measured in the two-pulse set up (Eq. (2)). $\chi_{xxxx}^{(3)}(t,t,\tau+t)$ corresponds to the limit $t_{1}\to\tau,t_{2}\to 0,t_{3}\to t$ (Fig. 3b). We may write $\sum_{k}A^{(4)}_{k}=f(t-\tau)\exp(-(t+\tau)/T_{2})$, where $f(t-\tau)$ comes from the sum of $\sin^{4}\theta_{k}\sin(2\lambda_{k}(t-\tau))$ and decreases as $|t-\tau|$ increases due to dephasing. Crucially, the arguments of $f$ and the $T_{2}$ term are orthogonal linear combinations of $t$ and $\tau$. The FT of $f$ is a broad continuum that depends on $\omega_{t}-\omega_{\tau}$, whereas the the FT of the $T_{2}$ term is a narrow Lorentzian function of $\omega_{t}+\omega_{\tau}$. The product of the two thus gives rise to a streak of rephasing signal in the imaginary part of the FT of $\chi_{xxxx}^{(3)}(t,t,\tau+t)$. The streak runs along the diagonal of the fourth quadrant, mirroring the energy range of spinon pairs. The width of the streak along the anti-diagonal is a direct measure of $1/T_{2}$: In the limit of $T_{2}\to 0$, the anti-diagonal width vanishes, reflecting the perfect phase cancellation in the spinon echo (Fig. 2d,e); With finite $T_{2}$, imperfect phase cancellation leads to a finite anti-diagonal width that scales with $1/T_{2}$ (Fig. 2f,g).

By contrast, $\chi_{xxxx}^{(3)}(t,\tau+t,\tau+t)$, corresponding to the limit $t_{1}\to 0,t_{2}\to\tau,t_{3}\to t$, does not contain a spinon echo (Fig. 3b). In the limit of $1/T_{1}\to 0$, $A^{(3,4)}_{k}$ are functions of $t_{3}=t$ only. In the frequency domain, this leads to a Dirac-$\delta$ peak on the $\omega_{\tau}=0$ line, which appears in the imaginary part as a streak along the $\omega_{t}$ axis (Fig. 2c). Taking $\omega_{\tau}$($\omega_{t}$) as the pumping (detecting) frequency, this may be interpreted as a pump-probe signal Lu et al. (2017). Increasing $1/T_{1}$ broadens the signal (Fig. 2h).

Both $\chi^{(3)}_{xxxx}$’s contain additional features that arise from $A^{(1,2,3)}_{k}$ terms in Eq. (5). Their FT contain a diffusive, weak NR signal in the first quadrant. They also contain a weak TR-like signal on the $\omega_{\tau}$ axis, which we discuss further in SM Sup .

To recap, the rephasing signal from the spinon echo process can directly reveal the $T_{2}$ time of spinon pairs. Crucially, in the absence of dissipation ($1/T_{1,2}=0$), the anti-diagonal width of the rephasing signal is zero. We now show that this feature is robust against quenched disorder. To this end, we set the transverse field $h_{n}$ and exchange constant $J_{n}$ to be site dependent, namely $h_{n}=ha_{n}$, $J_{n}=Jb_{n}$, where $a_{n},b_{n}$ are dimensionless random numbers drawn from a uniform distribution in the interval $[0.5,1.5]$. The linear response (Fig. 2k) shows only small changes comparing to the clean case. Since this model remains integrable, the spinons are still exact eigenstates, and therefore the anti-diagonal width of the rephasing signal remains resolution limited (Fig. 2n,o). Its strong sensitivity to dissipation, protected by the robustness against disorder, shows the utility of 2DCS.

In the FM phase, the $\sigma^{y}_{n}$ operators can also excite spinon pairs. We therefore expect the 2DCS spectrum with $\hat{y}$ polarization to be similar to $\hat{x}$. However, since the $\sigma^{y}_{n}$ is a non-local operator in the spinon basis, the analytic treatment made for $\hat{x}$ does not translate directly to $\hat{y}$. Nevertheless, as shown in the SM Sup , numerical calculation finds that the 2D spectra along $\hat{y}$ in the FM phase (Fig. S4) are qualitatively similar to that of $\hat{x}$. Note that in the PM phase the streak-like rephasing signal that is characteristic of fractional excitations is absent. $\chi^{(3)}_{yyyy}$ instead shows sharp isolated peaks Sup that are typical of nonlinear spin waves Lu et al. (2017); Babujian et al. (2016).

Using the TFIC as a prototypical example, we have demonstrated that THz 2DCS can resolve the spinon continuum and directly reveal the intrinsic line width of spinon pairs. We expect spinon echo to be a generic 2D spectral feature of models that host spinons. Provided that the spinons are coherent quasiparticles, the TLS picture naturally extends to higher dimensional spin liquids. In general, spinon echo will produce a rephasing streak qualitatively similar to that of the TFIC with finite $T_{1,2}$. In particular, the finite anti-diagonal width of the streak reflects the imperfect phase cancellation due to finite quasiparticle lifetime.

Our results may be applicable to CoNb2O6, which is the best known material example of a quasi-1D FM Ising chain Coldea et al. (2010); Morris et al. (2014); Viirok et al. (2018). CoNb2O6 orders at temperatures below $\sim 3$ K, but at slightly higher temperatures, the linear response is characterized by a broad lineshape characterized by superimposed 2- and 4-spinon continua, that hide information about spinon lineshapes. We expect that THz 2DCS can reveal the intrinsic spectral properties of spinons in this system. Experiments can be done in essentially the same fashion as previous THz 2DCS measurements on the conventional magnet YFeO3 Lu et al. (2017). Such experiments are underway. Analyzing theoretically the 2D spectra of more realistic material models will also prove fruitful. The spinon interactions present in these models will produce additional spectral features that are beyond the minimal model considered here.

With the information gained by establishing the technique on TFIC and its material realizations, we expect even richer information can be gained by applying THz 2DCS to higher dimensional materials that are suspected to harbor a spin liquid, but have only been characterized spectroscopically as having broad lineshapes such as 2D Kitaev magnets Banerjee et al. (2017), herbertsmithite Han et al. (2012), and triangular lattices Shen et al. (2016); Paddison et al. (2016); Zhang et al. (2018); Zhu et al. (2017). By direct analogy to the present results, we expect that one can measure the intrinsic lifetime of the multispinon excitations. Sharp anti-diagonal features may give direct evidence for fractionalized excitations and may be readily distinguished from highly damped conventional spin waves that could alternatively be present. Finally, we want to stress that our work is just an early step in understanding the utility of THz 2DCS for quantum materials. We believe important applications will be found in many systems including superconductors and topological materials.

Appendix A Expressing the 2D spectra in terms of nonlinear susceptibilities

In this section, we express the 2D spectra in terms of the nonlinear susceptibilities. In the setup considered in the main text, two linearly-polarized magnetic field pulses A and B arrive at the sample at time 0 and $\tau>0$. We assume the pulse profile takes the form of the Dirac-$\delta$ function for the sake of simplicity. The magnetic field pulse is given by:

[TABLE]

where $\beta$ labels the polarization, $s$ is time, and $A_{0,\tau}$ are the pulse areas.

The pulse field induces a response in the magnetization of the sample. The latter is related to the former through the linear and nonlinear magnetic susceptibilities:

[TABLE]

Here we have retained the leading and sub-leading nonlinear responses. $\alpha$ labels the spatial component of the magnetization. Subscript AB signifies the fact that the response mixes the effects due to both pulses. $t+\tau$ is the time at which the sample magnetization is measured/probed. $t>0$ is the delay between the second pulse (pulse B) and the time of measurement. Note that causality implies that $\chi^{(1)}_{\alpha\beta}(t+\tau-s_{1})=0$ if $s_{1}>t+\tau$. Similarly, $\chi^{(2)}_{\alpha\beta}(t+\tau-s_{1},t+\tau-s_{2})=0$ if $s_{1}>t+\tau$ or $s_{2}>s_{1}$; $\chi^{(3)}_{\alpha\beta}(t+\tau-s_{1},t+\tau-s_{2},t+\tau-s_{3})=0$ if $s_{1}>t+\tau$, or $s_{2}>s_{1}$, or $s_{3}>s_{2}$.

Substituting Eq. (6) in and using the causal properties of the susceptibilities, we find

[TABLE]

Setting $A^{\beta}_{\tau}\to 0$ ($A^{\beta}_{0}\to 0$) yields the response in the presence of the pulse A(B) alone:

[TABLE]

The nonlinear magnetization $M^{\alpha}_{NL}$ isolates the cross-correlation between the effects of the two pulses, which is defined as:

[TABLE]

This is the expression given in the main text.

Appendix B Calculating the 2D spectra along $\hat{x}$: without disorder

In this section, we provide the details of our calculation for the linear and nonlinear susceptibilities in the ferromagnetic phase of the quantum Ising chain without quenched disorder. In this specific case, we use periodic boundary condition for the spins. The Hamiltonian is

[TABLE]

B.1 Diagonalizing the Hamiltonian

The first step is to diagonalize the Hamiltonian Eq. (12) following the usual treatment for the TFIC Pfeuty (1970). To this end, we employ the Jordan-Wigner transformation:

[TABLE]

where $c^{\phantom{\dagger}}_{n}$ and $c^{\dagger}_{n}$ are the usual fermion annihilation/creation operators. The Hamiltonian reads:

[TABLE]

The $Z_{2}$ symmetry of the transverse field Ising model now becomes the fermion parity symmetry. It is known that the ground state lies in the parity even sector. This implies that the Jordan-Wigner fermions must obey anti-periodic boundary condition.

We switch to the momentum representation, $c_{n}=(1/\sqrt{L})\sum_{k}c_{k}\exp(ikn)$. The momentum $k=(2m+1)\pi/L$, $m\in Z$. Substituting in one gets,

[TABLE]

where $\epsilon_{k}=-2J\cos k-2h$, $\Delta_{k}=-2J\sin k$. The summation is restricted to the positive half of the first Brillouin zone. We perform the Bogoliubov transformation:

[TABLE]

where $\cos\theta_{k}=\epsilon_{k}/\lambda_{k}$, $\sin\theta_{k}=\Delta_{k}/\lambda_{k}$, $\lambda_{k}=\sqrt{\epsilon^{2}_{k}+\Delta^{2}_{k}}$. The Hamiltonian is now diagonalized in the spinon or $d_{k}$ basis:

[TABLE]

The ground state is annihilated by all $d_{k}$.

For later purposes, we may also express $M^{x}$, the total magnetization along $x$ direction, in the spinon basis:

[TABLE]

where

[TABLE]

It is instructive to use the Anderson pseudo-spins:

[TABLE]

The two-dimensional Hilbert space is spanned by the absence (pseudo-spin $\tau^{z}_{k}=-1$) and presence (pseudo-spin $\tau^{z}_{k}=1$) of the spinon pair with momenta $k,-k$. Using the Anderson pseudo-spins, the Hamiltonian reads:

[TABLE]

which describes an ensemble of independent two-level systems. The energy level splitting is $2\lambda_{k}$, which corresponds to the energy cost for exciting a pair of spinons out of vacuum. The $m^{x}_{k}$ operator reads:

[TABLE]

In particular, $m^{x}_{k}$ contains a term that switches the ground state and the excited state. In the Heisenberg picture:

[TABLE]

B.2 Computing the susceptibilities

We are now ready to calculate the various linear and nonlinear susceptibilities along $\hat{x}$. For now, we consider the standard transverse field Ising chain without any dissipation. Therefore, the spinons, being exact eigenstates of the Hamiltonian, possess infinite lifetime.

We first consider the linear susceptibility $\chi^{(1)}_{xx}(t)$. The starting point is the Kubo formula:

[TABLE]

Here, $\langle\cdots\rangle$ stands for the average in the ground state. The second equality follows from the fact that $m^{x}_{k}$ with different $k$ commute; the third equality follows from the Pauli algebra; the last equality follows from the property of the ground state that $\langle\tau^{x}_{k}\rangle=\langle\tau^{y}_{k}\rangle=0$, and $\langle\tau^{z}_{k}\rangle=-1$.

The second and third order nonlinear susceptibilities can be found in the same vein. The Kubo formula for the second order nonlinear susceptibility is given by,

[TABLE]

The nested commutators are computed step by step following the Pauli algebra:

[TABLE]

Since $\langle\tau^{x}_{k}\rangle=\langle\tau^{y}_{k}\rangle=0$, and $\langle\tau^{z}_{k}\rangle=-1$, we find the ground state expectation value of the nested commutators:

[TABLE]

Substituting the above back to the Kubo formula, we obtain:

[TABLE]

The Kubo formula for the third order nonlinear susceptibility is given by:

[TABLE]

Similar to the calculation for $\chi^{(2)}_{xxx}$, we compute the nested commutators step by step as:

[TABLE]

Here, we have omitted terms that are proportional to $\tau^{x}_{k}$ and $\tau^{y}_{k}$. These terms won’t contribute the ground state average, and, therefore, won’t appear in the final expression for the susceptibility. Using $\langle\tau^{x}_{k}\rangle=\langle\tau^{y}_{k}\rangle=0$, and $\langle\tau^{z}_{k}\rangle=-1$, we find:

[TABLE]

Substituting the above into the Kubo formula, we ultimately obtain:

[TABLE]

B.3 Incorporating dissipation effects

Here, we calculate the linear and nonlinear susceptibilities in the presence of spinon decay. Instead of starting from the first principle, we incorporate these effects phenomenologically, which is most conveniently achieved by using an approach akin to the optical Bloch equations Shen (1984). As we shall see, this phenomenological approach captures the essential features of dissipation and maintains the analytical tractability of the transverse field Ising chain.

To see how the optical Bloch equations naturally arise in this context, we note that each and every pair of spinons with momenta $\pm k$ effectively form a two-level system. Each two-level system is dynamically decoupled from the others. The state of the spinon pair $\pm k$ is specified by the $2\times 2$ density matrix $\rho_{k}$,

[TABLE]

where [math] and $1$ represent respectively the ground state and the excited state of the two-level system, or in terms of spinons, the absence and presence of the spinon pair. The time evolution of $\rho_{k}$ is canonically described by the optical Bloch equation (without the driving term). Its solution gives the explicit time evolution,

[TABLE]

Here, $2\lambda_{k}$ is the energy level splitting, or the energy cost for exciting the spinon pair $\pm k$ out of vacuum. $T_{1}$ is the time scale of spontaneous decay of the excited state, whereas $T_{2}$ is the decoherence time scale. The above may be rewritten as $\rho_{k}(t)=\mathcal{G}_{k}(t)\rho_{k}(0)$, where the super-operator $\mathcal{G}_{k}(t)$ propagates the density matrix. For the sake of simplicity, we assume $T_{1}$ and $T_{2}$ times are $k$-independent.

In what follows, we calculate the linear and nonlinear susceptibilities by taking an approach that is commonly used in the analysis of 2D coherent spectroscopy Mukamel (1999). Along the way, we will introduce the relevant terminology and tools.

We begin with the linear susceptibility $\chi^{(1)}_{xx}$. We expand the commutator that appear in the Kubo formula and obtain,

[TABLE]

Here, we take the following perspective: the correlation functions on the right hand side of the first equality represents a sequence of evolution for the density matrix $\rho_{k}$, known as the Liouville pathway in literature Mukamel (1999). $\mathcal{L}(m^{x}_{k})$ is the super-operator that means acting $m^{x}_{k}$ on $\rho_{k}$ from the left. Later, we shall use the super-operator $\mathcal{R}(m^{x}_{k})$, which means acting $m^{x}_{k}$ on $\rho_{k}$ from the right. $\mathcal{G}_{k}(t)$ is the aforementioned propagator of density matrix. $\rho_{k}(t)$ may be readily evaluated by using the definition of $m^{x}_{k}$ and $\mathcal{G}_{k}(t)$:

[TABLE]

Substituting it back to the expression for $\chi^{(1)}_{xx}$, we find,

[TABLE]

which is the result given in the main text.

It is convenient to use a bookkeeping device that graphically represents the various Liouville pathways that contribute to $\chi^{(1)}_{xx}$. Ladder diagrams, or so-called double-sided Feynman diagrams used commonly in 2D spectroscopy Mukamel (1999); Tomakoff (2011), fulfills this task. Here, we have slightly tailored them to suit our purposes. The left panel of Fig. 4 shows the ladder diagram that contributes $\chi^{(1)}_{xx}(t)$. We read the diagram as follows. Time flows upward. Each rung of the ladder (dashed line) represents a transition due to the action of $m_{k}$. The dot on the left corresponds to the action of $\mathcal{L}(m^{x}_{k})$, whereas the dot on the right corresponds to the action of $\mathcal{R}(m^{x}_{k})$. Note, however, that the last dot doesn’t correspond to any real transitions; it describes the final measurement process ($\mathrm{Tr}(m^{x}_{k}\rho_{k})$). Representing it as a “transition” is mere convention. The space between two rungs is the evolution of $\rho_{k}$ under $\mathcal{G}_{k}(t)$.

For the specific diagram for linear response given in the left panel of Fig. 4, the bottom of the ladder diagram is $|0\rangle\langle 0|$, corresponding to the pure ground state. Going upward, there is a transition to $|1\rangle\langle 0|$ due to the action of $m_{k}$. This picks up a matrix element $\langle 1|m^{x}_{k}|0\rangle$. This is followed by the evolution due to $\mathcal{G}_{k}(t)$, which picks up a time-evolution factor $\exp(-(1/T_{2}+2i\lambda_{k})t)$. Finally, the measurement yields another matrix element $\langle 0|m^{x}_{k}|1\rangle$. Collecting these, the total outcome is $|\langle 1|m^{x}_{k}|0\rangle|^{2}\exp(-(1/T_{2}+2i\lambda_{k})t)$, which is the term appearing in the summation of Eq. (37). It is important to bear in mind that other diagrams can in principle be drawn for $\chi^{(1)}_{xx}$. However, their contributions vanish due to cancellation.

We now consider the second-order nonlinear susceptibility $\chi^{(2)}_{xxx}$. A similar procedure allows us to relate the Kubo formula for $\chi^{(2)}_{xxx}(t,\tau+t)$ to the Liouville pathways. Out of a total of 9 possible pathways and their Hermitian conjugates, 5 diagrams have finite contribution to $\chi^{(2)}_{xxx}$, which are shown as ladder diagrams in the right panel of Fig. 4. Comparing to the diagram for $\chi^{(1)}_{xx}$, two new elements appear: First, the spontaneous decay of the excited state $|1\rangle\langle 1|$ to the ground state $|0\rangle\langle 0|$ is graphically represented as a red arrow. Second, the diagrams with odd number of dots on the right leg have a negative sign.

The diagrams in the right panel of Fig. 4 are organized in terms of their phase factors. The top row all contain a phase factor $\exp(2i\lambda_{k}\tau)$. The sum of these is:

[TABLE]

Note the additional minus sign from the pre-factor of the Kubo formula. This is the result given in the main text.

We finally turn to the third-order nonlinear susceptibility $\chi^{(3)}_{xxxx}(t_{3},t_{2}+t_{3},t_{1}+t_{2}+t_{3})$. As the complexity of the diagrams grow exponentially with the order of nonlinearity, we now have in total 42 Liouville pathways plus their Hermitian conjugates, out of which 16 independent pathways contribute. These are given in Fig. 5. The diagrams are organized in terms of their phase factors: from the top to the bottom, the phase factors of the diagrams in the first, second, third, and the last row are functions of $t_{1}+t_{2}$, $t_{1}+t_{2}+t_{3}$, $t_{1}+t_{3}$, and $t_{1}-t_{3}$ respectively. By now all rules of the ladder diagrammatic have been given. We now employ these rules to obtain the analytic results. The sum of the first row is given by:

[TABLE]

This is the result given in the main text.

We conclude this section by commenting that, when taking the limit of $T_{1,2}\to\infty$, the susceptibilities agree with the results given in Sec. B.2.

B.4 Performing the Fourier transform

With the analytic formula at hand, we are ready to compute the 2D coherent spectra. This would require performing a two-dimensional, one-sided Fourier transform with respect to the variables $t,\tau$ over the domain $t>0,\tau>0$. While it is possible to find the explicit expression directly, the resulted expressions are too complicated to be illuminating. We therefore choose a different route: we first calculate $\chi^{(2)}_{xxx}$ and $\chi^{(3)}_{xxxx}$ in the time domain using the aforementioned analytic expressions. We then perform a numerical fast Fourier transform (FFT) on the computed signal, thereby obtaining the results shown in the main text.

A couple of considerations come into choosing the time grid. Let $dt$ be the step width of the time grid. The maximal frequency that can be sampled by this grid is $\omega_{\mathrm{max}}=\pi/dt$. On the other hand, the maximal spinon pair energy $2\lambda_{k}$ is $4(J+h)$, and there is no signal above this frequency. In practice, we set $dt=\pi/(5(J+h))$, which corresponds to $\omega_{\mathrm{max}}=5(J+h)$ slightly above the maximal spinon pair energy.

The maximal number of steps of the grid, on the other hand, is limited by the system size $L$. The frequency spacing $\sim O(1/L)$, which sets a time scale $\sim O(L)$. Simulating over a time longer than this time scale will result in spurious peaks in the frequency domain due to the finite system size. In practice, we determine the maximal time steps by monitoring the real time signal. The finite size effect manifests itself as a “revival” of oscillations at late time. We limit the FFT to times before the onset of such a revival.

Fig. 6 shows the various nonlinear susceptibilities in real time. Here, $L=100$, $h/(J+h)=0.3$, and $1/T_{1}=1/T_{2}=0$. The diagonal feature in $\chi^{(3)}_{xxxx}(t,t,\tau+t)$ is the rephasing signal described in the main text. The “wavefronts” that propagate along the diagonal at the top right corner are due to the finite size effect. We therefore perform the 2D FFT over a smaller domain demarcated by the dashed box.

Appendix C Calculating the 2D spectra along $\hat{x}$: with quenched disorder

Here we provide the details of calculation for $\chi^{(1)}_{xx}$, $\chi^{(2)}_{xxx}$, and $\chi^{(3)}_{xxxx}$ in the presence of quenched disorder. As quenched disorder destroys the translational invariance, there is no strong reason for adhering to the periodic boundary condition. We thus choose open boundary condition. A major advantage of open boundary conditions is that it allows one to calculate the linear and nonlinear susceptibilities for other polarizations with only small modifications McCoy et al. (1971); Derzhko and Krokhmalskii (1997).

The Hamiltonian reads:

[TABLE]

Here, $J_{n}$ and $h_{n}$ are site dependent, random variables. We compute all physical observables for several realizations of disorder and then average over these realizations.

C.1 Diagonalizing the Hamiltonian

The diagonalization procedure is similar to the case without disorder. Here, it is more convenient to recast the Jordan-Wigner fermions $c_{n}$ and $c^{\dagger}_{n}$ in terms of even and odd Majorana fermion modes Kitaev (2001):

[TABLE]

In particular, $\{\alpha_{m},\alpha_{n}\}=2\delta_{mn}$, $\{\beta_{m},\beta_{n}\}=2\delta_{mn}$, and $\{\alpha_{m},\beta_{n}\}=0$. Furthermore, $\alpha_{n}$ is even under complex conjugation, whereas $\beta_{n}$ is odd. The Jordan-Wigner transformed Hamiltonian is then recasted as:

[TABLE]

Here, $\alpha,\beta$ are understood as the $L\times 1$ column vectors formed by $\alpha_{i}$ and $\beta_{i}$, respectively. $\mathrm{J}$ is a $L\times L$ bidiagonal matrix:

[TABLE]

We now seek the appropriate orthogonal transformation over $\alpha,\beta$ that brings the above Hamiltonian into diagonal form. To this end, we perform a singular value decomposition (SVD) on $\mathrm{J}$: $\mathrm{J}=\mathrm{U\Lambda V^{t}}$. Inserting the SVD into the above, we find,

[TABLE]

In the second line, we have defined new Majorana modes $\gamma\equiv\mathrm{V}^{t}\alpha$, and $\delta\equiv\mathrm{U}^{t}\beta$. They obey a similar set of anti-commutation relations: $\{\gamma_{m},\gamma_{n}\}=2\delta_{mn}$, $\{\delta_{m},\delta_{n}\}=2\delta_{mn}$, and $\{\gamma_{m},\delta_{n}\}=0$. Furthermore, $\gamma_{n}$ are even under complex conjugation and $\delta_{n}$ are odd. Now modes carrying different labels $n$ are independent. In the last equation, we have defined the complex fermions $d_{n}=(\gamma_{n}+i\delta_{n})/2$. Clearly, $d_{n}$ annihilates vacuum.

It is useful for latter purpose to define a $2L\times 2L$ correlation matrix $\mathrm{C}(t)$, whose matrix elements are defined as:

[TABLE]

To find its explicit value, we calculate a similar correlation matrix $\mathrm{K}(t)$ for $\gamma_{n}$ and $\delta_{n}$. The non-zero matrix elements of $\mathrm{K}(t)$ are:

[TABLE]

$\mathrm{C}$ is simply related to $\mathrm{K}$ through an orthogonal transformation:

[TABLE]

C.2 Computing the susceptibilities

We are now ready to calculate the susceptibilities for the polarization along $\hat{x}$ for the disordered case. We begin with the total magnetization:

[TABLE]

Note the leftmost $R$ sites and the rightmost $R$ sites of the chain are excluded from the summation. This is to avoid the effects from the zero energy edge modes. The choice of $R$ depends on the correlation length.

To compute the susceptibilities, we first relate them to multi-point spin-spin correlation functions:

[TABLE]

where $S^{x}_{jlmn}(s_{1},s_{2},s_{3},s_{4})=\langle\sigma^{x}_{j}(s_{1})\sigma^{x}_{l}(s_{2})\sigma^{x}_{m}(s_{3})\sigma^{x}_{n}(s_{4})\rangle$ is the four-point spin correlation function.

Thus the task of computing susceptibilities is reduced to computing spin correlation functions. To this end, we express $\sigma^{x}_{n}$ in terms of Majorana modes $\alpha_{n}$ and $\beta_{n}$, and then compute the average using the Wick theorem. The results are compactly expressed as matrix Pfaffians McCoy et al. (1971).

The two-point spin correlation function:

[TABLE]

In our calculation, we use a a chain of 40 sites, and set $R=5$ (Eq. 48). We first compute the susceptibilities in real time and then perform a numerical FFT. The choice of time grid is identical to Sec. B.4. We compute the matrix Pfaffian by using the numerical package developed by Wimmer Wimmer (2012).

Appendix D Calculating the 2D spectra along $\hat{y}$

D.1 Computing the susceptibilities

Here we provide the details for calculating the susceptibilities with polarization along $\hat{y}$. The procedure is largely in parallel with Section C: we first relate the susceptibilities to multi-point spin correlation functions, and then compute the spin correlation functions by going to the Jordan-Wigner fermion basis. Since this section concerns the susceptibilities along $\hat{y}$, we need to find the multi-point spin correlation functions $S^{y}_{jl}(s_{1},s_{2})$ and $S^{y}_{jlmn}(s_{1},s_{2},s_{3},s_{4})$. Note the second-order susceptibility $\chi^{(2)}_{yyy}$ vanishes by symmetry.

The spin operator $\sigma^{y}_{n}$ is non-local in the Jordan-Wigner fermion basis:

[TABLE]

Thus, the two-point spin correlation function becomes a highly non-local correlation function in the fermion basis. Nevertheless, it may be compactly expressed as matrix Pfaffian Derzhko and Krokhmalskii (1997):

[TABLE]

Apparently, computing $S^{y}_{jlmn}(s_{1},s_{2},s_{3},s_{4})$ is considerably more expensive comparing to $S^{x}_{jlmn}(s_{1},s_{2},s_{3},s_{4})$ due to the large size of the matrix. We benchmark this procedure on a short chain with $L=8$. We compute the susceptibilities by using two methods: (a) a direct numerical diagonalization of the spin Hamiltonian followed by the Kubo formula, and (b) the Pfaffian-based method. The absolute error is within $10^{-13}$ for all cases considered.

D.2 Results for ferromagnetic and paramagnetic phases

In this section, we present the 2D spectra in the ferromagnetic phase and the paramagnetic phase for the $y$ polarization. The calculation set up is identical to Sec. C. We shall focus on the case without disorder or spinon decay. We first consider the ferromagnetic phase. We set the model parameter to the representative value $h/(h+J)=0.3$. As the operator $\sigma^{y}_{n}$ flips the eigenvalues of $\sigma^{z}_{n}$ and thereby creating domain walls, we expect the spectra along $\hat{y}$ is qualitatively similar to $\hat{x}$. The calculation confirms this picture. The imaginary part of the linear susceptibility, $\mathrm{Im}\chi^{(1)}_{yy}$, shows a spinon continuum virtually identical to $\mathrm{Im}\chi^{(1)}_{xx}$ (Fig. 7(a)). The 2D spectra (Fig. 7(b)(c)) are also qualitatively similar to the spectra along $\hat{x}$. In particular, we observe a strong rephasing streak in the fourth quadrant of Fig. 7c. The streak runs along the diagonal direction, mirroring the energy range of the spinon pairs, whereas its width along the anti-diagonal direction is resolution limited, reflecting the infinite life time of the spinon pairs. The weak, box-shaped shadow is likely due to the fact that the $\sigma^{y}_{n}$ operator is non local in the spinon basis.

We next turn to the paramagnetic phase, where the spins are ordered in $\sigma^{z}_{n}$. We set the model parameters to the representative value $h/(J+h)=0.7$. For the pulse polarization along $\hat{y}$, the optically active excitation is the magnon with momentum $k=0$, and its energy $\omega_{k=0}=2(h-J)=0.8(J+h)$. We observe a sharp resonance at $\omega_{k=0}$ in $\mathrm{Im}\chi^{(1)}_{yy}$ consistent with the lack of fractionalized excitation in the paramagnetic regime. In the 2D spectra, the various peaks and their locations are expected for nonlinear spin waves Lu et al. (2017). In the Fourier transform of $\chi^{(3)}_{yyyy}(t,\tau+t,\tau+t)$, we observe a strong peak at the pumping frequency $\omega_{\tau}=0$ and the detection frequency $\omega_{t}=\omega_{k=0}$, known as the pump-probe signal, and a weaker peak at the pumping frequency $\omega_{\tau}=2\omega_{k}$ and the detection frequency $\omega_{t}=\omega_{k=0}$, known as the two-quanta signal Lu et al. (2017). On the other hand, the Fourier transform of $\chi^{(3)}_{yyyy}(t,t,\tau+t)$ contains a non-rephasing peak at $\omega_{\tau}=\omega_{k=0}$ and $\omega_{t}=\omega_{k=0}$, and a rephasing peak at $\omega_{\tau}=-\omega_{k=0}$ and $\omega_{t}=\omega_{k=0}$.

Appendix E Diagonal transitions and frequency vectors

Two-level systems are paradigmatic models for nonlinear optics Shen (1984). In a standard treatment such as in Ref. Shen (1984), it is assumed that the light can induce transitions between the ground state and the excited state. The “diagonal” transitions. i.e. from the ground state to the ground state, or from the excited state to the excited state, are usually not included in the model as they are often forbidden by symmetry. In this work however, the two-level systems that emerge from the transverse field Ising chain do not represent atomic orbitals but rather the spinon pairs. Diagonal transitions are allowed in this system. These diagonal transitions give rise to novel features in the 2D coherent spectra including a non-rephasing-like (NR-like) signal in $\chi^{(2)}$, and a terahertz-rectification-like (TR-like) signal in $\chi^{(3)}$, which should be absent in the 2D spectra of a canonical two-level system. These novel features may appear unusual at first glance. For instance, the TR is often thought of as a $\chi^{(2)}$ process. In this work, however, their signals seems to appear in the spectra associated with $\chi^{(3)}$ as well. As we shall see, the TR-like signals in $\chi^{(3)}$ spectra do not reflect true nonlinear rectification processes, but instead result from the aforementioned diagonal transitions. In other words, these TR-like signals are not true TR signals but merely resemble them in the frequency plane.

In this section, we elucidate the impact of the diagonal transitions on 2D coherent spectra by analyzing a toy model. We will also connect our results to the “frequency vector” scheme for understanding the 2D spectra [][.Notethepositionsofthefrequencyvectorsusedthereindifferfromoursduetodifferentparameterizationofthetimevariables.]kuehn2011nonlinear. The purpose of this section is pedagogical. It introduces the frequency vector scheme to readers that are unfamiliar with it, while, for readers with a background in nonlinear optics, it aims at showing how the aforementioned novel features in the 2D spectra can be understood by using a familiar language.

We first consider a standard two-level system. Written in spin language, the Hamiltonian is:

[TABLE]

Here, $\sigma^{z}=-1$ ($\sigma^{z}=1$) in the ground (excited) state. $\Omega>0$ is the energy level splitting.

Suppose we perform a THz 2D coherent spectroscopy experiment on this system. We use two linearly polarized, Dirac-$\delta$ pulses. The magnetic field is given by:

[TABLE]

Here, $s$ is the time, $A_{0,\tau}$ the pulse area, $\epsilon^{\alpha}$ the polarization vector, which we set to be:

[TABLE]

Here, $\theta$ is the angle between the polarization vector and the $\hat{z}$ axis. We measure the response parallel to $\epsilon^{\alpha}$. As our purpose is pedagogical, we will omit dissipation for simplicity.

Since the $\hat{z}$ component of the magnetic field pulse is non-vanishing when $\theta\neq\pi/2$, diagonal transitions are allowed, and consequently the 2D spectra should contain the novel features mentioned in the beginning of this section. Alternatively, we may think that mixed nonlinear response functions such as $\chi^{zyy}$ contribute to the 2D spectra. These two points of view are equally valid and complementary to each other. In what follows, we compute the 2D spectra of this system and relate these results to the frequency vector scheme.

E.1 Contributions to $\chi^{(2)}$

We first consider the $\chi^{(2)}$. To this end, we switch to the Heisenberg picture: $\sigma^{y}(t)=\cos(\Omega t)\sigma^{y}+\sin(\Omega t)\sigma^{x}$, and $\sigma^{z}(t)=\sigma^{z}$. The latter follows from the fact that $\sigma^{z}$ commutes with $H$. A straightforward calculation using the Kubo formula shows:

[TABLE]

We observe that $\chi^{(2)}=0$ when the diagonal transition is absent, i.e. $\theta=\pi/2$. $\chi^{(2)}_{zyy}$ is a function of $\tau$ only. In the frequency domain, this gives rise to a peak at $\omega_{\tau}=\Omega$, and $\omega_{t}=0$. This is the THz rectification (TR) signal. However, $\chi^{(2)}_{yzy}$ gives rise to a non-rephasing-like (NR-like) signal at $\omega_{\tau}=\omega_{t}=\Omega$. This is remarkable because usually NR signals occur in $\chi^{(3)}$.

We may rationalize these results by using a scheme for 2DCS known as frequency vectors. One may view it as a generalization of the usual frequency sum/difference picture discussed in nonlinear optics Shen (1984). First consider $\chi^{(2)}_{zyy}$ (Fig. 8, left). For this response function, both pulses are polarized along $\hat{y}$. Since $\sigma^{y}$ induces transitions between the ground state and the excited state, the first pulse will induce an oscillatory signal of the form $\exp[i\Omega(t+\tau)]$ and the second an oscillatory signal of the form $\exp(i\Omega t)$. These two oscillations are represented by two non-orthogonal vectors in the frequency plane spanned by $\omega_{\tau}$ and $\omega_{t}$ Kühn (2011): the first pulse (pulse A) gives a diagonal vector, whereas the second pulse (pulse B) gives a horizontal vector. The $\chi^{(2)}$ nonlinear coupling produce a difference vector (A$-$B) along the $\omega_{\tau}$ axis. This is precisely the TR signal we found by calculation. The sum vector (A$+$B) would produce a second-harmonic generation signal. However, it is absent in the two-level system as such a system can only emit radiation at fixed frequency $\Omega$. Note that we only consider the signals in the first and fourth quadrant of the frequency plane. Signals in the other two quadrants are related by reflection.

Somewhat different considerations apply to $\chi^{(2)}_{yzy}$ (Fig. 8). In this case, the first pulse (A) is polarized along $\hat{y}$, which gives a diagonal vector (A) in the frequency plane. Importantly, the second pulse (B) is polarized along $\hat{z}$, and such a pulse does not induce any oscillations. We may represent this as a null vector (B) in the frequency plane. The sum or difference of the two vectors (A$\pm$B) are thus the NR-like signal we found by calculation.

E.2 Contributions to $\chi^{(3)}$

Having understood $\chi^{(2)}$, we can now move to $\chi^{(3)}$. A straightforward calculation yields,

[TABLE]

The experimentally measured $\chi^{(3)}$’s are the two special limits of $\chi^{(3)}(t_{3},t_{2}+t_{3},t_{1}+t_{2}+t_{3})$. $\chi^{(3)}(t,t+\tau,t+\tau)$ corresponding to the limit $t_{1}\to 0,t_{2}\to\tau,t_{3}\to t$. The nonlinear response $\chi^{(3)}_{yyyy}(t,t+\tau,t+\tau)$ produces a pump-probe (PP) signal at $\omega_{\tau}=0$ and $\omega_{t}=\Omega$. This response function corresponds to a pulse sequence in which all three pulses are polarized along $\hat{y}$ (Fig. 9, left). In terms of the frequency vectors, this means the first two pulses gives a diagonal vector (A), and the third pulse gives a horizontal vector (B). The PP signal thus comes from the linear combination A$-$A$+$B. The other combinations, namely A$+$A$\pm$B, while in principle possible, are absent for the two-level system.

$\chi^{(3)}_{zyzy}(t,t+\tau,t+\tau)$ produces a TR-like signal at $\omega_{\tau}=\Omega$ and $\omega_{t}=0$. This response corresponds to a pulse sequence in which the first pulse is polarized along $\hat{y}$, the second along $\hat{z}$, and the third along $\hat{y}$ (Fig. 9, middle). In terms of frequency vectors, the first gives a diagonal vector (A1), the second a null vector (A2), and the third a horizontal vector (B). The TR-like signal comes from the combination A${}_{1}\pm$A${}_{2}-$B. The other combination A${}_{1}\pm$A${}_{2}+$B would give rise to a second-harmonic generation signal, which is absent for the two-level system.

$\chi^{(3)}_{yzzy}(t,t+\tau,t+\tau)$ produces a NR signal at $\omega_{\tau}=\omega_{t}=\Omega$. The pulse polarizations are successively along $\hat{y}$, $\hat{z}$, and $\hat{z}$ (Fig. 9, right). In the frequency plane, the first pulse gives a diagonal vector (A1), and the second and the last null vectors (A2 and B). Their combinations (A${}_{1}\pm$A${}_{2}\pm$B) has a unique outcome, namely the NR signal.

$\chi^{(3)}(t,t,t+\tau)$ correspond to the limit $t_{1}\to\tau,t_{2}\to 0,t_{3}\to t$. $\chi^{(3)}_{yyyy}(t,t,t+\tau)$ produces a NR signal and a rephasing signal; $\chi^{(3)}_{zyzy}(t,t,t+\tau)$, a TR-like signal; $\chi^{(3)}_{yzzy}(t,t,t+\tau)$, a NR signal. These signals can also be understood in terms of frequency vectors as shown in Fig. 10. For $\chi^{(3)}_{yyyy}(t,t,t+\tau)$, the first pulse (A) gives a diagonal vector, and the second and third a horizontal vector (B). The combination A$+$B$-$B yields the NR signal, whereas the other combination $-$A$+2$B gives the rephasing signal. For $\chi^{(3)}_{zyzy}(t,t,t+\tau)$, the first pulse (A) gives a diagonal vector, the second a null vector (B1), and the third a horizontal vector (B1). Their combinations A$\pm$B${}_{1}-$B2 yield the TR-like signal. Finally, for $\chi^{(3)}_{yzzy}(t,t,t+\tau)$, the first pulse gives a diagonal vector (A), and the second and third null vectors (B). Their combinations A$\pm$B${}_{1}\pm$B2 yield the NR signal.