Spectroscopic fingerprints of gapped quantum spin liquids, both conventional and fractonic

We explain how gapped quantum spin liquids, both conventional and 'fractonic', may be unambiguously diagnosed experimentally using the technique of multidimensional coherent spectroscopy. 'Conventional' gapped quantum spin liquids (e.g. $Z_2$ spin liquid) do not have clear signatures in linear response, but do have clear fingerprints in non-linear response, accessible through the already existing experimental technique of two dimensional coherent spectroscopy. Type I fracton phases (e.g. X-cube) are (surprisingly) even easier to distinguish, with strongly suggestive features even in linear response, and unambiguous signatures in non-linear response. Type II fracton systems, like Haah's code, are most subtle, and may require consideration of high order non-linear response for unambiguous diagnosis.

Quantum spin liquids are phases of quantum matter which feature fractionalized excitations and emergent deconfined gauge fields (for reviews, see [1,2]).Recently an exotic 'fractonic' variant of quantum spin liquids has been proposed, which additionally host emergent excitations with restricted mobility (for reviews, see [3,4]).The most stable and theoretically best understood examples of spin liquids, both conventional and fractonic, have the ground state manifold separated from the rest of the spectrum by an energy gap.The Z 2 spin liquid provides a paradigmatic example of a conventional gapped quantum spin liquid, while the X-cube model and Haah code provide paradigmatic examples of gapped fractonic spin liquids, of type I and type II respectively.While great strides have been made in the theoretical understanding of gapped quantum spin liquids, these phases have not yet been unambiguously observed in any experimental system.Part of the challenge here is that it is hard to find unambiguous experimental diagnostics for these exotic phases, which are accessible using currently available experimental techniques.This is in contrast to gapless spin liquids (both conventional and fractonic), which do have clean signatures e.g. in the form of 'pinch points' in the dynamical structure factor, which may be probed via neutron scattering [5,6].
In this paper we explain how gapped quantum spin liquids, both conventional and fractonic, may be unambiguously diagnosed using already existing spectroscopic tools.Gapped conventional spin liquids do not have clear diagnostics in linear response, but do have unambiguous fingerprints in non-linear response.These may be identified through the technique of two dimensional coherent spectrsocopy (2DCS), originally pioneered in the context of nuclear magnetic resonance and physical chemistry [7][8][9][10], and recently applied also to solid state systems [11][12][13][14].Gapped fracton phases of type-I are actually easier to identify, with strongly suggestive features even in linear response, and unambiguous fingerprints in non-linear response.Fracton phases of type II are most subtle, and may require consideration of high order non-linear response for unambiguous diagnosis.
The rest of this paper is structured as follows.In Section I we introduce the paradigmatic models that we use as examples of gapped spin liquid, type I fracton and type II fracton phases respectively.In Section II we explain the properties of the respective phases in linear response, and point out how type I fracton phases (and only type I fracton phases) have clear signatures therein.In Sec.III we explain how to calculate the non-linear susceptibility measured in a 2DCS experiment.In Sec IV we explain the key signatures of each of the phases considered within 2DCS.We conclude in section V with a discussion of outlook and implications.Technical details are relegated to the appendices.

I. PHASES OF INTEREST
We are interested in three types of phases: conventional gapped spin liquids, type I fractons, and type II fractons.We consider a paradigmatic example of each, as well as a 'control' example which is not a spin liquid.While the model Hamiltonians we write down are for the most part fine tuned to an exactly solvable point, our interest will be specifically in those features of the phase that are robust to perturbations, and which thus do not rely on being at the fine tuned point.The solvable Hamiltonians that we present therefore simply represent 'fixed point Hamiltonians' that characterize the respective phases.
As a non-spin liquid control example, we consider the quantum Ising model in two or three spatial dimensions.This is defined on a two (three) dimensional square (cubic) lattice with spin-1/2 variables living on the vertices.The Hamiltonian is where •, • denotes nearest neighbors, Z and X are Pauli opreators, and λ is a parameter controlling strength of perturbations.Here (...) denotes arbitrary local perturbations that fail to commute with the ZZ term but which arXiv:2010.07947v1 [cond-mat.str-el]15 Oct 2020 respect Ising symmetry (the simplest such perturbation being a transverse field, k X k ).Note that in two (three) dimensions, this is an interacting model in its own right.We work with periodic boundary conditions.We work mostly in the ferromagnetic phase, λ 1, when the ground state spontaneously breaks Z → −Z symmetry, and the elementary excitations are spin flips.This is not a spin liquid phase and serves as our control example.
As a paradigmatic example of a gapped conventional spin liquid, we consider the perturbed two-dimensional toric code [15].This is defined on the two-dimensional square lattice with spin-1/2 variables living on the links.The Hamiltonian is where v indicates a sum over vertices, A v is a product of four X operators on the four links connected to a vertex v, p indicates a sum over square plaquettes, B p is a product of four Z operators around a plaquette, (...) denotes arbitrary local perturbations, and λ 1.The exactly solvable point is λ = 0. We work with periodic boundary conditions and the ground state which have zero flux penetrating the great circles of the torus.At the exactly solvable point, it has eigenvalue +1 under every A and B operators.
As a paradigmatic example of a type I fracton phase, we consider the (perturbed) X-cube model [16,17].This is defined on a three-dimensional cubic lattice with spin-1/2 variables living on the links.The Hamiltonian is where v denotes sum over vertices, A ab v is a product over four Z operators on links connected to vertex v lying in the ab plane, c denotes sum over cubes, B c is a product over twelve X operators framing an elementary cube, λ 1, and (...) denotes arbitrary local perturbations.The exactly solvable point is λ = 0. We work with periodic boundary conditions and the ground state which have zero flux penetrating the great circles of the torus.At the exactly solvable point it has eigenvalue +1 under every A, B. The elementary excitations consist of lineons, planeons, and fractons.Lineons are created at the end of a string of X operators and are one dimensional particles (i.e.able to move in only one direction, since the string can change length but cannot change direction without creating additional lineons); planeons are created at the end of a string of Z operators and are two dimensional particles (i.e.able to move in a two dimensional plane).Fractons are created at the corners of a rectangular membrane of Z operators and are totally immobile under local perturbations.All these properties are explained at length in [3,17].
As a paradigmatic example of a type II fracton phase, we consider the (perturbed) Haah code [18].This is de-fined on a three dimensional cubic lattice with two spin-1/2 variables on every vertex.The Hamiltonian is where A is a particular product of X type Pauli's and identities around a cube, and B is likewise with X → Z, Λ ≈ 1, λ 1, and (...) denotes arbitrary local perturbtions.The exactly solvable point is λ = 0.The ground state is unique with open boundary conditions (eigenvalue +1 under every A and B at the solvable point).The elementary excitations are fractons (totally immobile excitations).There are no subdimensional particles, but a composite of four fractons in a tetrahedral arrangement is locally creatable (and therefore mobile).
It should be emphasized that while these model Hamiltonians appear exceedingly baroque, nonetheless the phases are robust.The value of these model Hamiltonians is that they allow for exact calculations.While we will present some exact results obtained at the solvable points (λ = 0), our emphasis in what follows will be on universal features that should be present throughout the phase, and not just at the solvable points.Also, while we have chosen periodic boundary conditions for analytic convenience, experiments of course will have open boundary conditions.However, since the non-linear response we are discussing does not connect distinct topological sectors, the choice of boundary conditions should not make a difference.

II. LINEAR RESPONSE
The quantities of interest in linear response spectroscopy are correlators M z (t)M z (0) and M x (t)M x (0) and M y (t)M y (0) , where M z = Z measures the magnetization in the z direction and M x and M y are defined analogously.All of these are measurable with the suitable polarization of the incident light.The expectation values are taken on the ground state.Inserting a resolution of the identity in terms of eigenstates of the Hamiltonian, these correlators take the form Now lets evaluate these for each of our models of interest, working close to the solvable points.

A. Ising model
At the exactly solvable point, M z (t)M z (0) will only have a Fourier component at zero frequency, because Z operators will not produce any transitions between states.However X and Y operators will produce spin flips.|α will thus be a state with a single flipped spin.The energy cost at the exactly solvable point will be z, where z is the coordination number of the lattice, producing a sharp δ function peak at frequency ω = z.Away from the exactly solvable point the spins will acquire dispersion E(k).We assume that because of the large (effectively infinite) speed of light, while the light can supply energy it cannot supply momentum, so the transitions involved can only be to states with zero total momentum.For the Ising model |α has to be the k = 0 spin flip state, in which case we will observe a sharp feature in linear response at the corresponding frequency.This sharp feature is indicative of the existence of a sharp excitation at zero momentum.This assumes that the energy of a single k = 0 spin flip does not overlap the multi-excitation continuum, so there is no decay channel available for the single spin flip state.This assumption should be safe as long as we are away from the critical point.

B. Two dimensional toric code
Lets start by examining the M z correlator.A Z operator will anticommute with two A operators.Thus, |α will be a state with two 'spinons.'At the exactly solvable point it will have energy 2 and a sharp signal in the spectrum at frequency 2. If we perturb about the exactly solvable point the spinons will acquire dispersion E(k).Let us assume as before that the light does not supply appreciable momentum, so |α has to be a state with zero total momentum.Nonetheless the two spinons can have momentum ±k and so the state can have energy 2E(k) [19].So away from the exactly solvable point here will be no sharp signature in linear response, but only a broad continuum.The behavior of the M x correlator is analogous with the only difference an altered dispersion relation.Finally, a Y operator will create two A spinons and two B spinons.At the solvable point there will be a sharp signal at frequency 2 + 2J.Away from the solvable point it will blur into a continuum.
The observation of a broad continuum in linear response (with any polarization) indicates that Z polarized light does not create a single sharp excitation at zero momentum.However, it cannot distinguish between the case where it creates a pair of sharp quasiparticles (as is the case here, at least in the frequency range where the two spinon state does not overlap the multi-spinon continuum), and the case where there is no sharp excitation because e.g.there are no good quasiparticles.We will see in Section IV how 2DCS can disambiguate between these possibilities.

C. X-cube
A Z operator will anticommute with the four B operators containing that link.Let us suppose it is an x directed link for specificity.It will thus create a multiplet of four cube excitations (fractons).At the exactly solvable point this will produce a sharp signal at frequency 4K.Now lets consider the generic case, where we perturb about the solvable point.As before we restrict to states |α with total momentum zero.Now pairs of two fractons are planeons (mobile excitations restricted to move in either the xy plane or the xz plane, depending on how we group the fractons).These will acquire dispersion upon turning on perturbations so there will be a range of possible energies for |α .As before there will be nothing sharp in linear response.Now lets consider the M x correlator.An X type operator will anticommute with four A type operators (two for each vertex sharing that link).The excitation at each vertex is a lineon -i.e. when we turn on perturbations each lineon can move freely, but only in the direction of the link.Nevertheless, the lineons will acquire dispersion.Naively we might think that this means there will not be any sharp signal in linear response, but this is too fast -actually there is a sharp signal here.Namely, the lineons, being one dimensional particles, will have a one dimensional Van Hove singularity in their density of states, which will diverge near the gap edge as 1/ √ E. If we assume that the energy of the lineons is minimized at k = 0 (which is a highly plausible assumption) then this should be readily observable in linear response spectroscopy, and should provide a clean signature of the existence of one dimensional excitations in this (otherwise three dimensional) system.
Finally, a Y operator will anticommute with four B operators and four A operators, and will thus excite four 'fractons' (which could be grouped into two planeons), as well as four vertex operators (grouped into two lineons).Again the lineons will have a crisp signature through their one dimensional density of states, and so the M y correlator will also have clear signatures of the existence of one dimensional particles.
In all the above, we have implicitly worked in the infinitesimal neighborhood of the exactly solvable point where e.g. an X operator only creates a single pair of lineons.Away from the solvable point there will also be some admixture of four lineon creation and higher multiparticle processes.(One way to see this is to perform a Schrieffer Wolff transformation [20] to remove the perturbation, and to note that the new 'A' operators have some admixture of highly multiparticle operators in the original basis).These 'multiparticle absorbtion' channels should be added on to the absorbtion spectrum discussed above, and will give rise to additional multiparticle continua.However, the existence of multiparticle continua will not alter the fact that there is a one dimensional Van Hove singularity, which should be clearly detectable in absorbtion experiments.(In addition, the threshold frequency for the multiparticle continua will generically be higher than the threshold frequency for the two lineon continuum, so the above analysis will be strictly accurate in the frequency range where the light supplies enough energy to create two lineons, but not to create more than two excitations).
Thus, type I gapped fracton phases which support one dimensional (i.e.lineon) excitations will have a clear signature in linear response, in the form of a one dimensional Van Hove singularity in the absorption spectrum of an otherwise three dimensional material.The only challenge (for crisply identifying this signature as coming from lineons) is to show that the signature does not come because e.g. the material has the structure of a system of weakly coupled one dimensional wires.This should be straightforward to do in combination with some alternative experiment e.g.transport.In transport experiments (performed at non-zero temperature, so that there is transport) a system of coupled wires would have a clear preferred axis, whereas a fracton system with lineons would not, since there are lineons that can move along any lattice axis.An alternative (purely spectroscopic) way to tell is that there will not be any signature of 'one dimensional' physics in the M z correlator.

D. Haah code
Now X and Z correlators behave the same way.Acting with either X or Z creates four excitations in a 'tetrahedral' arrangement (four fractons), with energy cost 4 (or 4Λ) at the exactly solvable point.These cannot be grouped into mobile excitations (away from the solvable point), so we continue to have a sharp 'delta function' signature at frequency 4 even away from the exactly solvable point.This remains the case even if we e.g.turn on a magnetic field and sweep the angle.
The existence of a sharp robust signature in linear response is tempting to identify as diagnostic of a type II fracton phase.However we should bear in mind that something similar could show up if we had e.g.well isolated two level systems in the problem, or indeed if local fields created sharp individual excitations (as in the Ising model) which are then constrained to have zero momentum.As such, the observation of a sharp robust signal in linear response alone is not sufficient to conclude that one is dealing with a type II fracton phase.

E. Summary of linear response
To conclude: type I fracton phases (X-cube) do have a crisp diagnostic in linear response, in the form of a one dimensional Van Hove singularity in the absorbtion spectrum which is present only for certain polarizations of the incident light.This is a signature that the phase contains fractionalized 'lineon' excitations i.e. excitations that can only move in one dimension.In contrast, neither conventional gapped spin liquids (toric code) nor type II fracton phases (Haah code) have crisp signatures in linear response.To diagnose these phases we will need to turn to 2DCS and non-linear response.

III. 2DCS AND NONLINEAR RESPONSE
We consider a 2DCS experiment where two δ function pulses of magnetic fields are applied at time 0 and τ .Two magnetic fields are linearly polarized along α and β direction, respectively: The pulse induced magnetization M γ (t + τ ) is measured at later time t + τ : where N spin is the total number of spins on a lattice.There is no sum over repeated indices.The canonical 2DCS experiment extracts the nonlinear response by subtracting off the signal observed in the presence of either pulse alone.This leaves us with Now let us consider two possible experiments.Experiment I has (α, β, γ) = (x, x, x).Experiment II has (α, β, γ) = (z, z, z).Other combinations of polarizations are of course possible, and may be interesting to consider, but these two are sufficient to provide an unambiguous diagnostic of spin liquids, both conventional and frac-tonic.If we consider the models on a square and cubic lattice, the second order susceptibility χ γαβ (t, t + τ ) = 0 for those two experiments because we need even number of X or Z operators to pair up creation and annihilation of excitations.Then, the leading contributions to the nonlinear signals are the third order susceptibilities, χ γααβ (t, t + τ, t + τ ) and χ γαββ (t, t, t + τ ).They are related to four-point correlation functions via the generalized Kubo formula [13]: where Θ is Heaviside step function (equal to one when the argument is larger than zero), and the operators have been time evolved (Heisenberg picture) with respect to the system Hamiltonian.Note that M β is always applied at time τ , M α at time 0, and M γ at time t + τ .Let us define where in the last line we have reverted to the Schrodinger picture, and E µ,ν,λ are the energy difference between the excited state and the ground state.Now we can write We always have t γ = t + τ , t β = τ and t α = 0. Let S abcd be just the matrix element part of R abcd (i.e.without the phases in Eq. ( 11)).Then we can write S γβαα e −i(Eµt+Eν τ ) − S βγαα e i((Eµ−Eν )t−Eν τ ) − 2S αγβα e i((Eµ−Eν )t+(Eµ−E λ )τ ) +2S αβγα e i((Eν −E λ )t+(Eµ−E λ )τ ) + S ααγβ e i((Eν −E λ )t+Eν τ ) − S ααβγ e i(E λ t+Eν τ ) , where we have set = 1, and the expressions are now in the Schrodinger picture.Here S abcd is the implicit function of µ, ν, and λ.We note that the experimentally observed signal will be χ γααβ (t, t+τ, t+τ )B 2 α B β +χ γαββ (t, t, t+τ )B α B 2 β .Thus both third-order susceptibilities could, in principle, be extracted by individually varying the strength of the two applied pulses.We will thus generally present the two third-order susceptibilities separately, even though what will be observed will of course be the sum of the two.
There is one subtlety to note here.In what follows we will generally be looking at the response in frequency space, and the Fourier transform of the terms in square brackets is the quantity of interest.However, because the susceptibility is multiplied by step functions, the observable signal will consist of the 'interesting' signal (the terms in square brackets) convolved with (δ(ω t ) + Pi/ω t )(δ(ω τ ) + Pi/ω τ ) = δ(ω t )δ(ω τ ) − P 1 ωtωτ + i(δ(ω t )P 1 ωτ + δ(ω τ )P 1 ωt ), where P denotes principal value.How this effects the observable signal will be noted below, where appropriate.
We start by working in two spatial dimensions.In this case, at leading order close to the solvable point, the experiment II with zzz polarization does not produce any transitions or lead to any nonlinear signal.So we can focus on the experiment I with xxx polarization.There is only one matrix element to be evaluated.The states |µ and |λ both contain a single flipped spin with zero momentum, with energy E µ = E λ = 4 + δ 0 , where δ 0 denotes corrections to the energy coming from perturbations [21].Meanwhile |ν contains either zero or two flipped spins which may be adjacent but do not need to be.Thus we have E ν = 0, E ν = 8 + δ k or E ν = 6 + δ , where the correction δ k is the kinetic energy contributions coming from configurations where the two flipped spins have momentum ±k, and δ −2δ 0 is the 'binding energy' for two adjacent spin flips (with zero momentum).In the former case we need to sum over k.Note that because E ν is dispersive, the signal will be broadened in the frequency directions corresponding to E ν .However it will be sharp in all other directions, modulo the broadening from convolution discussed above.We further note that the pathway with E ν = 0 just corresponds to linear response done twice, and cannot contribute to the nonlinear response (See Appendix for detailed justification, but one can also simply note that the contribution coming from such a pathway would scale as the system size, whereas the nonlinear susceptibility must of course be independent of the system size.).We can therefore focus on the terms with E ν = 8 + δ k or E ν = 6 + δ .
Let (ω t , ω τ ) be the frequencies conjugate to t and τ , respectively.Then the channel with E ν = 6 + δ will give rise to sharp (delta function) signals at frequencies ω t,τ equal to 2+δ −δ 0 , 4+δ 0 , or 6+δ +δ 0 .Thus, signals will appear with a spacing in frequency space equal to roughly half the linear response gap.It is important to bear in mind that this occurs in a clearly non-spin-liquid system, so the appearance of signals at a fraction of the linear response gap is not in itself a signal of fractionalization.
There will also be contributions from intermediate states with E ν = 8 + δ k , These will be broadened along one direction (that corresponding to E ν but will be sharp in the other direction.This will give rise to streaks in the ω t direction, with and without offset 4 in the ω τ direction, to a streak in the ω τ direction with offset 4 in the ω t direction, and also to a diagonal streak with offset 4 in the ω t direction -see Fig. 1 for an illustration.This is our 'control' example of a non-spin liquid.Features that appear herein cannot be viewed as spin liquid signatures.
Meanwhile in three dimensions, the possible energies are shifted to E µ = 6 = E λ and E ν = {0, 10 + δ , 12 + δ k } and the signatures are analogous to the two dimensional case.
The actual results in Fig. 1 are produced via a direct calculation from Eq.14, 15 in two dimensions, including matrix elements.Subfigures (a, b) corresponds to the exactly solvable point λ = 0.In subfigures (c, d) we have assumed that the excited states are plane waves, and have calculated matrix elements accordingly.Excited states containing two flipped spins are a symmetrized product of two plane waves, as appropriate for bosonic excitations.For a detailed discussion of the calculation, see the Appendices.
We now have to discuss an annoying subtlety associated with the fact that the experimentally observed signal gets convolved with δ(ω t )δ(ω τ ) − P 1 ωtωτ + i(δ(ω t )P 1 ωτ + δ(ω τ )P 1 ωt ), because the Fourier transforms only run over the positive time axes.This produces a weak 1/ω broadening in both directions (from convolution with the second term), and a somewhat stronger 1/ω broadening along both axes (from convolution with the third and fourth terms).This 'broadening from convolution' can make it hard to see the more physical broadening from dispersion.Within our present approximations, the physical part of the spectrum is pure imaginary in Fourier space (see appendices), and given that the strongest part of the broadening from convolution comes with an extra factor of i, it may simply be removed by taking the imaginary part of the non-linear susceptibility -this is done in Fig. 1, and will be done throughout for all the models we study.However, once decay of the excitations is re-introduced, whether through coupling to extraneous degrees of freedom such as phonons, or through decay into multi-excitation sectors (at frequencies overlapping the multi-magnon continuum), then the 'physical' part of the signal will in general become complex, and the 'broadening from convolution' problem will become unavoidable.However, the 1/ω broadening from convolution still leaves the intensity of the signal sharply peaked where it would have been, so if one simply applies a high-pass filter on intensity then this may suffice to deal with the problem.

B. Toric code
For the toric code (our paradigmatic example of a conventional gapped spin liquid), both polarization configurations yield interesting results.A Z operator applied to any link creates a pair of vertex excitations residing on the two vertices adjacent to that link.(More generally, the end points of strings of Z operators produce vertex excitations).Meanwhile, X operators applied to a link create a pair of plaquette excitations on the two plaquettes containing that link.More generally, strings of X operators create plaquette excitations are the ends.Both vertex and plaquette excitations are mobile but are not locally creatable (the locally creatable things are pairs of vertex or plaquette excitations).We will restrict the analysis to the sector where each Z or X operator locally creates only two plaquette or vertex excitations and will ignore mixing with the 'many excitation' sector.'Multiexcitation' channels will make additional contributions to the signal, but as we will see unambiguous diagnostics appear already at leading order (and will be there regard-

(a)
< l a t e x i t s h a 1 _ b a s e 6 4 = " m r U less of what additional signatures appear from multiexcitation pathways).As usual, neglect of multi-excitation pathways should also be safe at frequencies below the multi-excitation gap.We will allow all excitations to have dispersion (i.e.we will implicitly perturb about the solvable point).
Here |λ and |µ contain a pair of plaquette excitations with net momentum zero.We have E µ = 2J + 2δ k and E λ = 2J + 2δ k , where the plaquette excitations have momenta ±k and ±k respectively, and a corresponding kinetic energy δ k (δ k ).Meanwhile, |ν could contain any of: zero plaquette excitations (E ν = 0), two plaquette excitations (E ν = E µ = E λ ), or four plaquette excitations (E ν = E µ + E λ ).All momenta have to be summed over.As before, the channel with E ν = 0 corresponds to linear response done twice and cannot make a (properly extensive) contribution to the non-linear response, and will therefore be ignored.(In the appendix we explicitly show how this channel cancels to give zero contribution to the non-linear response).We will therefore focus on the channels with H e n Y 9 5 6 4 q T z x z B H z i f P 1 y M j U U = < / l a t e x i t > 4w < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 q o v 7 / t C z v 7 one direction, after allowing dispersion.
We begin by considering the sequence with E ν = E µ = E λ .This channel corresponds to creation of a string with visons (plaquette excitations) at the two ends, followed by elongation of the string, followed by annihilation of the string, with the fact that all states have the same energy being a consequence of the string having zero line tension.Since zero line tension for strings is equivalent to deconfinement, we can reasonably expect signals from this channel to contain signatures of deconfinement.Now the first and last terms in Eq. (14,15) give rise to diagonal stripes (outside the two vison gap), whereas the second and fifth terms in Eq.14), together with the second, third, fourth and fifth terms in Eq.15 collectively give rise to a sharp signal along the ω τ axis in the two dimensional Fourier transform (again, above the two vison gap).Both these features (the sharp streak along the ω τ axis, and the streak along the diagonals with no offset are signatures of deconfinement, at least in conjunction with the lack of a sharp signal in linear response.The crucial aspect here is that the streaks have no offset -while the signal is only present outside the two vison gap, the extrapolation of the streak goes through the origin.This is a consequence of the fact that, as mentioned, a 'string' hosting visons at the two ends can change length with no change in energy -a signature of deconfinement.Another way to view it is that each X operator creates a pair of visons, but the visons are their own antiparticles so the application of an X operator to a two vison state can leave us in a two vison state.That is, the first pulse creates a pair of visons with momentum ±k, the second creates another pair of visons with momentum ±k (one of which annihilates one of the visons in the original pair), and the third pulse then returns us to the ground state.Again, the fact that a local operator creates a pair of sharply defined Z 2 charged excitations is a signature of fractionalization.In contrast, in the two dimensional Ising model the equivalent streaks had an offset, indicating that the second X operator necessarily changed the energy of the state.In the one dimensional Ising model, equivalent signals do arise [12] -but then the one dimensional Ising model also has deconfined fractionalized excitations (the domain walls).
In Fig. 2, we present computations of the signal including appropriate matrix elements in 14, 15, both at and away from the solvable point.The wavefunctions are taken in real space at the exactly solvable point, and are assumed to be symmetrized products of plane waves away from the solvable point.Full details of the calculation are presented in the appendices.As we see, the features present in the explicit calculation are precisely those expected by inspection of 14, 15, and knowledge of the toric code.Now there is the channel with E ν = E µ + E λ , where E µ and E λ may be individually varied.Naively, this produces a contribution that scales like the square of the volume, so to leading order the contributions from this channel must cancel.However, this may leave a nontrivial extensive piece behind -after all this channel is not just linear response done twice.Explicit calculation of this channel is impractically tedious away from the exactly solvable point -the states |µ and |λ are a symmetrized product of two plane waves, whereas the state |ν is a symmetrized product of four plane waves (4! terms in the sum), and altogether there are 2 * 4! * 2 = 96 possible terms to evaluate and sum.However, there is also no good reason to believe this channel should have clear diagnostics of deconfinement, and we have already found clear signatures of deconfinement from the channel with E µ = E ν = E λ .The only concern might be that the channel with E ν = E ν + E λ might somehow cancel the clear signatures we found coming from the channel with E µ = E ν = E λ .We have verified that this does not happen at the exactly solvable point by explicit calculation (see Appendix, and also Fig. 2(a)).Accordingly, we believe it is safe to ignore the contributions from this channel -there may be additional contributions to the non-linear response from it, but they do not cancel the contributions from the E µ = E ν = E λ channel, and those already contain clear signatures of deconfinement.
Thus, the combination of linear response and 2DCS data can provide clean fingerprints of deconfined fractionalized excitations.If linear response does not show any sharp features this could be because local fields create multiplets of dispersive fractionalized excitations (as in the toric code), or it could just mean that they create excitations that are not good quasiparticles.However, if 2DCS also shows a sharp stripe along the ω τ axis and along the diagonal, then this indicates that local fields in fact create a pair of quasiparticles which are their own antiparticles, such that a second pulse can leave us in the 'two quasiparticle' sector.This is the case for deconfined fractionalized excitations that live at the ends of strings with zero line tension -the second pulse just elongates the string -but it is not the case if there are non-fractionalized excitations which are just not good quasiparticles.
The analysis of the experiment for fields in the z direction proceeds analogously -fields create vertex excitations instead of plaquette excitations, but we obtain identical results.
Thus, from linear response, we can learn whether local fields can create sharply defined individual quasiparticles.If the answer is no, then using 2DCS we can determine if local fields actually create pairs of deconfined fractionalized quasiparticles.

C. X-cube
We now consider 2DCS on the X-cube phase.Now an X operator applied to a link creates a pair of vertex excitations (lineons) on the two vertices adjacent to that link.These lineons can move freely, but only in one dimension (the dimension of the link).
A Z operator applied to a link creates a quartet of cube excitations (fractons) on the four cubes sharing that link.These fractons are individually immobile, but can move in pairs as planeons (i.e. can move freely in the plane orthogonal to the long axis of the planeon).Now |µ and |λ both contain a single pair of lineons with zero total momentum (and with associated Van Hove singularity in the DOS).Meanwhile, |ν can contain either zero lineons, a single pair of lineons (in which case E µ = E ν = E λ , or two pairs of lineons.The lineons are deconfined fractionalized excitations which are their own antiparticles, much like the visons in the toric code.The analysis parallels the analysis of the toric code, with the one exception that the DOS is that of a one dimensional system.In particular, we will have the same diagnostics of deconfinement as in the toric code.These include streaks along the ω τ axis and along the diagonal, with no offset.This is illustrated in Fig. 3. Full calculations are in the appendices -we have performed exact calculations at the solvable point, and away from the solvable point we have calculated the signal from the channel with E µ = E ν = E λ , which is the channel expected to give rise to signatures of deconfinement.We have also verified that the contribution from the intermediate channel with E ν = E µ + E λ does not cancel the signal of interest at the solvable point, where this extra (more complicated) channel can be analytically treated.Thus, linear response on X-cube can tell that a local field creates one dimensional excitations (lineons), through DOS.2DCS can then verify that the lineons are deconfined.Now µ and λ both contain a pair of planeons with zero net momentum.Since the planeons are mobile there is an associated continuum of energies.Meanwhile options for ν include zero planeons (E ν = 0), one pair of planeons Thus far the analysis parallels the toric code case, with planeons taking the place of visons.In particular, there will be a diagonal streak in the two dimensional Fourier transform, with no offset, as well as a streak along the ω τ axis with no offset, as a signature of deconfined fractionalized excitations which are their own antiparticles.However, there is one additional intermediate option, whereby ν contains one pair of planeons and also two fractons (which each have zero momentum), such that E µ = E λ but E ν = E µ + 2K.Lets explore the consequences of this additional channel, bearing in mind that E µ can take a continuum of values (outside the two planeon gap), and that we are only interested in features that are sharp in at least one direction.
The first and sixth terms in Eq.14 will give rise to diagonal streaks with offset 2K in the ω t direction, outside the two planeon gap, while the second and fifth terms will give rise to streaks along ω τ with offset 2K in the ω t direction outside the two planeon gap.Meanwhile the first and sixth terms in Eq.15 will again produce diagonal streaks with no offset, while the middle four terms in Eq.15 generate streaks in the ω τ direction with offset   2K in the ω t direction outside the planeon gap.Thus, the main 'new' consequence is offset stripes.Now, offset stripes also appeared in the TFIM, but there the offset was equal to the linear response gap.Here the offset is equal to half the linear response gap.Altogether this is quite informative.The lack of any sharp features in linear response (besides the gap) tells us that local fields do not create individual local sharp quasiparticles.The appearance of vertical streaks along axis in the two dimensional Fourier transform tell us that local fields do create deconfined pairs of quasiparticles (planeons).The additional appearance of sharp streaks at an offset a fraction of the linear response gap tells us that the planeons can further fractionalize, but also that the objects generated by planeon fractionalization are non-dispersive (fractons).Altogether 2DCS provides a crisp diagnostic for planeons, fractons, and lineons in the X-cube phase.
In the Appendix, we have provided explicit calculations for the X-cube model in the z polarization.Calculations are exact at the solvable point, whereas away from the solvable point they include the channels with E ν = E µ = E λ and E ν = E µ + 2K, which produce the key signals of interest.Figure 4 plots the results.

D. Haah code
We will consider here only the xxx experimental geometry -the zzz geometry has analogous behavior.Now both |µ and |λ are produced by acting on the ground state with a singe X operator, and hence contain a 'tetrahedron' of four fractons, with energy E µ = E ν = 4.This cannot be subdivided into two (or more) mobile excitations, so given the constraint that |µ and |λ have zero momentum, it follows that they also have a sharply defined energy, rather than a continuum.What about |ν ?This can contain either: zero fractons (E ν = 0)), six fractons (E ν = 6 + δ ), or eight fractons (E ν = 8 + 2E k ).
In the last we have taken account of the fact that composites of four fractons are locally creatable and hence mobile objects, which can carry momentum.
Let us compare to the Ising model in three dimensions (i.e. the same spatial dimensionality as the Haah code).In the Ising model, we have E µ = E ν = 6 (equal to the linear response gap) and E ν = 0 or E ν = 12+2E k (2 spin flips with equal and opposite momentum), or E ν = 10+δ (bound state of two adjacent spin flips).While there is a quantitative difference to the Haah code, there is not a qualitative difference -each time we act with an X operator we create a single mobile excitation (a single spin flip in the Ising model, a tetrad of excited cubes in the Haah code), and two adjacent mobile excitations can form a bound state.While the quantitative differences are suggestive that what is going on in the Haah code may be different, there is not therefore an unambiguous diagnostic in 2DCS.
It should however be possible to diagnose the Haah code in higher order non-linear response.A key property of the Haah code is that acting with an appropriate fractal membrane of X operators only produces four cube excitations at the 'corners' [17].Thus, while acting with a single X operator produces a state with energy equal to the linear response gap, acting with four X operators in a tetrahedral arrangement also produces a state with the same energy gap.In contrast, in the Ising model there is no way to produce a state with energy equal to the linear response gap by acting with four operators.If we go to sufficiently high order to access this intermediate state,      then we should be able to qualitatively distinguish the Haah code from the trivial case of the Ising model.However, this will require consideration of at least the seventh order non-linear response, which is beyond the scope of the present paper.Experimentally also, measurement of a seventh order non-linear response may be challengingalthough it could certainly be accomplished in principle, e.g. with a seven pulse experiment.Unambiguous diag-nosis of the Haah code is therefore a more challenging task than identification of its conventional spin liquid or type-I fracton counterparts.
V. CONCLUSIONS Thus, we have identified purely spectroscopic fingerprints of both conventional gapped spin liquids and gapped fractonic phases using the technique of 2DCS spectroscopy.The easiest to diagnose are type I fracton phases with lineon (one dimensional) excitations -these can be diagnosed even in linear response through the existence of one dimensional Van Hove singularities in the absorbtion spectrum.If we turn to non-linear response and the 2DCS spectrum, then type I fracton phases have further fingerprints of deconfinement, of the existence of lineon excitations, and also of the existence of totally immobile fracton excitations.Conventional gapped spin liquids are intermediate in subtlety to detect.In linear response there are no sharp features besides the bulk gap.However in non-linear response and the 2DCS spectrum, there appear sharp signatures of the existence of deconfined quasiparticles.Finally, type II fracton phases (such as the Haah code) are most subtle to diagnose.Linear response has no clear signatures, and even 2DCS gives only quantitative (but not qualitative) distinctions from trivial possibilities.An unambiguous diagnosis of the Haah code likely requires a consideration of high order non-linear response (seventh order susceptibility should suffice), which may be challenging to access experimentally.
A number of future directions present themselves for consideration.For instance, thus far we have considered only the simplest version of the experiment, where each of the pulses (and the observed signal) has the same polarization.Crossed polarizations would be interesting to consider in future work, and could yield additional information beyond the 'single polarization' experiments discussed herein.Additionally, the calculations presented herein have been for a system prepared in the ground state.Extension to low temperature Gibbs states would also be a natural problem for future work.Perhaps most significantly, we have thus far ignored dissipative line broadening.However, dissipation will inevitably be present, whether extrinsic (due to coupling to e.g.phonons), or intrinsic (due to decay to the multiexcitation continuum).This will produce line broadening, and a quantitative understanding thereof could be a fruitful endeavour.Indeed, it seems plausible [22] that a careful analysis of the temperature dependence of line broadening might itself provide a clean diagnostic for otherwise difficult to detect phases like the Haah code.The ability of 2DCS to distinguish between energy relaxation (T 1 time) and dephasing (T 2 time) may also be useful in this regard.
We have thus far also assumed that we are dealing with clean systems, whereas realistic experimental systems will inevitably be disordered.Disorder is not expected to be important for the signals discussed herein, but would certainly be important for any analysis of dissipative lineshapes, where the ability of 2DCS to distinguish 'intrinsic' line broadening from inhomogenous broadening could be particularly useful.Incorporating disorder into the analysis would thus also be a fruitful project for the future.
Furthermore, thus far we have identified sharp signatures of deconfinement, lineons, fractons etc, but while these serve as crisp diagnostics for a fractionalized phase (or a fracton phase), they do not establish which fractionalized (or fractonic) phase we are dealing with.For example, the diagnostics we identified would not be able to distinguish between the Z 2 spin liquid represented by the toric code and that represented by the doubled semion model [23], nor between the X-cube phase and the semionic X-cube phase [24].Identifying diagnostics able to unambiguously identify which phase within each class we were dealing with would also be a fruitful topic for future work.Similarly, diagnostics capable of identifying symmetry enriched phases, akin to [25] would also be worth identifying.
Finally, thus far we have considered only the third order non-linear response, probed in a two-pulse experiment with 2DCS.There are other possibilities.For instance, one could consider the third order response probed through a three pulse experiment, and analyzed via the three dimensional Fourier transform, or we could consider fifth (or higher) order nonlinear susceptibilities.Indeed we have argued that the seventh order non-linear response should be particularly interesting, as it likely offers a clean diagnostic for the Haah code.All these would no doubt yield additional information, at the cost of complicating the necessary experiment.Regardless, it appears clear that the multidimensional spectroscopy technique places a powerful new tool at our disposal, which can be used to identify exotic phases that would be difficult to diagnose via conventional techniques.
By plugging in the energy cost into Eqs.(S3) and (S4), we can calculate the third order susceptibilities χ (a/b) x (t, τ ) for Ising model on a square lattice: The terms proportional to L 2 are precisely cancelled, which is consistent with the nonextensive definition of the nonlinear susceptibilities.Then Under generic perturbations, a single spin-flip excitation (magnon) becomes dispersive.Since the contribution from perturbative process involving |ν = |0 is eventually cancelled, we only need to consider nonlinear dynamics of two neighboring spin flip composite (ε ν = 6 + δ) and dispersive motion of two separated magnons (ε ν = 8 + ε q1 + ε q2 ), where δ and ε q1,2 are perturbative corrections to the static energy cost.

B. 2D Toric code
Magnetic fields along x and z-axis excite e (vertex excitation, i.e., bosonic charge) and m (plaquette excitation, i.e., Z 2 flux) particles of the two-dimensional toric code, respectively.y-polarized fields excite both e and m particles.On a square lattice, dynamics of the fluxes (χ xxxx ) and the bosonic charges (χ zzzz ) are the same.So let's focus on the nonlinear response of the bosonic charges under two z-polarized incident pulses: Im R (a/b) (µ, ν, λ; t, τ ) 0|Z l |µ µ|Z l |ν ν|Z l |λ λ|Z l |0 (S47) Im R (a/b) (µ, ν, λ; t, τ ) 0|Z l |µ µ|Z l |ν ν|Z l |λ λ|Z l |0 (S48) where denotes the four edges of a plaquette, and 2 and 4 denote the sets of all two-particle states and four-particle states, respectively.The last two sums, Eqs.(S49) and (S50), are the consequence of the ground stating being a symmetric superposition of all possible closed loop.

Nonlinear dynamics of lineons
The third-order susceptibilities χ xxxx (t, t + τ, t + τ ) ≡ χ x (t, τ ) are quantifying nonlinear response of the lineons.Without any perturbations, the susceptibilities can be exactly calculated with the resolution of identity in the real space basis: Like 2D toric code, the terms proportional to the system size L 3 are vanishing [Eq.(S93) = 0].With a single lineon excitation energy ε l = 2, the third-order susceptibilities are In the presence of weak perturbations, the nonlinear susceptibilities can be derived by following the similar calculations done for the toric code.Again we focus on the case where |ν ∈ 2: R (a/b) (p, p, r; t, τ ) cos 2 p(cos p + cos r) cos r + R (a/b) (p, r, r; t, τ ) cos p(cos p + cos r) cos 2 r (S98) p,q,r R (a/b) (p, q, r; t, τ ) cos p(cos p + cos q)(cos q + cos r) cos r, where R (a/b) (p, q, r; t, τ ) are calculated from Eqs. (S82) and (S83) with the one-dimensional dispersion for lineons, ε l (p) = 2 + w cos p.This reproduces ν ∈ 2 contributions of Eq.(S95) and (S96) when w = 0.
t e x i t s h a 1 _ b a s e 6 4 = " m r U x V X w U b J U 2 4 f 2 n + e f C e / U j B c A = " > A A A B 8 3 i c b V D L S g M x F L 3 x W e u r 6 t J N s A h 1 U 2 a k o O 6 K b l x W s A / o D C W T Z t r Q T G Z I M k I Z + h t u X C j i 1 p 9 x 5 9 + Y a W e h r Q c C h 3 P u 5 Z 6 c I B F c G 8 f 5 R m v r G 5 t b 2 6 W d 8 u 7 e / s F h 5 e i 4 o + N U U d a m s Y h V L y C a C S 5 Z 2 3 A j W C 9 R j E S B Y N 1 g c p f 7 3 S e m N I / l o 5 k m z I / I S P K 1 p 9 x 5 9 + Y a W e h r Q c C h 3 P u 5 Z 6 c I B F c G 8 f 5 R m v r G 5 t b 2 6 W d 8 u 7 e / s F h 5 e i 4 o + N U U d a m s Y h V L y C a C S 5 Z 2 3 A j W C 9 R j E S B Y N 1 g c p f 7 3 S e m N I / l o 5 k m z I / I S P K Q U 2 K s 5 H k R M e M g z G r B x W x Q q T p 1 Z w 6 8 S t y C V K F A a 1 D 5 8 o Y x T S M m D R V E 6 7 7 r J M b P i D K c C j Y r e 6 l m C a E T M m J 9 S y W J m P a z e e Y Z P r f K E I e x s k 8 a P F d / b 2 s n 8 A f o 8 w e s A Z F 3 < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 Y Z Z t j H A J b I C 3 D k F i T 6 3 6 8 d U A B 0 = " > A A A B 8 3 i c b V D L S g M x F L 3 x W e u r 6 t J N s A h 1 U 2 a k o O 6 K b l x W s A / o D C W T Z t r Q T G Z I M k I Z + h t u X C j i 1 p 9 x 5 9 + Y a W e h r Q c C h 3 P u 5 Z 6 c I B F c G 8 f 5 R m v r G 5 t b 2 6 W d 8 u 7 e / s F h 5 e i 4 o + N U U d a m s Y h V L y C a C S 5 Z 2 3 A j W C 9 R j E S B Y N 1 g c p f 7 3 S e m N I / l o 5 k m z I / I S P K FIG. 2. Imaginary part of 2D Fourier transform Ft,τ of the third-order susceptibilities χxxxx for 2D toric code (a, b) without perturbations and (c) with generic perturbations.(a) Im Ft,τ [χxxxx(t, t + τ, t + τ )] and (b) Im Ft,τ [χxxxx(t, t, t + τ )] of the toric code show sharp pointlike signals due to plain motion of a pair of two plaquettes (marked with green circles, Eµ = Eν = E λ process in Sec.IV B) and quantum processes involving two pairs of plaquettes (marked with red circles, Eν = Eµ + E λ process in Sec.IV B).The signals at (ωτ , ωt) = (±2, 0) are the distinct signatures of the model having deconfined excitations because the signals suggest vanishing string tension between two plaquettes of the toric code.(c) Im Ft,τ [χxxxx(t, t + τ, t + τ )] of the perturbed toric code.Im Ft,τ [χxxxx(t, t, t + τ )] shows indistinguishably equal nonlinear spectrum.The nonlinear signals due to the process involving only two plaquettes (Eµ = Eν = E λ ) are shown.Unlike Ising model, presence of the dispersive deconfined excitations yield clear spread of the signals by the full bandwidth of the single plaquette dispersion 4w = 0.8.Green stars mark locations of the signals without the perturbations.

2 .
(α, β, γ) = (z, z, z) t e x i t s h a 1 _ b a s e 6 4 = " m r U x V X w U b J U 2 4 f 2 n + e f C e / U j B c A = " > A A A B 8 3 i c b V D L S g M x F L 3 x W e u r 6 t J N s A h 1 U 2 a k o O 6 K b l x W s A / o D C W T Z t r Q T G Z I M k I Z + h t u X C j i 1 p 9 x 5 9 + Y a W e h r Q c C h 3 P u 5 Z 6 c I B F c G 8 f 5 R m v r G 5 t b 2 6 W d 8 u 7 e / s F h 5 e i 4 o + N U U d a m s Y h V L y C a C S 5 Z 2 3 A j W C 9 R j E S B Y N 1 g c p f 7 3 S e m N I / l o 5 k m z I / I S P K Q U 2 K s 5 H k R M e M g z G r k Y j a o V J 2 6 M w d e J W 5 B q l C g N a h 8 e c O Y p h G T h g q < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 Y Z Z t j H A J b I C 3 D k F i T 6 3 6 8 d U A B 0 = " > A A A B 8 3 i c b V D L S g M x F L 3 x W e u r 6 t J N s A h 1 U 2 a k o O 6 K b l x W s A / o D C W T Z t r Q T G Z I M k I Z + h t u X C j i 1 p 9 x 5 9 + Y a W e h r Q c C h 3 P u 5 Z 6 c I B F c G 8 f 5 R m v r G 5 t b 2 6 W d 8 u 7 e / s F h 5 e i 4 o + N U U d a m s Y h V L y C a C S 5 Z 2 3 A j W C 9 R j E S B Y N 1 g c p f 7 3 S e m N I / l o 5 k m z I / I S P K Q U 2 K s 5 H k R M e M g z G r 0 Y j a o V J 2 6 M w d e J W 5 B q l C g N a h 8 e c O Y p h G T h g q FIG. 3. Imaginary part of 2D Fourier transform Ft,τ of the third-order susceptibilities χxxxx for X-cube model (a, b) without perturbations and (c) with generic perturbations.Qualitatively, the Fourier-transformed χxxxx of X-cube model is similar to those of the toric code.(a) Im Ft,τ [χxxxx(t, t + τ, t + τ )] and (b) Im Ft,τ [χxxxx(t, t, t + τ )] of X-cube model show sharp pointlike signals suggesting the deconfinement, i.e., a pair of two lineons move with zero string tension (marked with green circles, Eµ = Eν = E λ process in Sec.IV C 1). Red circled signals are due to quantum processes involving four lineons (Eν = Eµ + E λ process in Sec.IV C 1). (c) Im Ft,τ [χxxxx(t, t + τ, t + τ )] of the perturbed X-cube model.Im Ft,τ [χxxxx(t, t, t + τ )] is almost the same.Only the signals relevant to Eµ = Eν = E λ process are shown.Similar to the toric code, dispersive lineons exhibit clear spread of the signals by the full bandwidth of a lineon dispersion 2w = 1.Green stars mark locations of the peaks without the perturbations.
t e x i t s h a 1 _ b a s e 6 4 = " V Q I 0 k M k a u 0 6 Y d X 7 x 1 r a W i B Q U N U I = " > A A A B 8 3 i c b V D L S g M x F L 2 p r 1 p fV Z d u g k W o m z I j B X V X c O O y g n 1 A Z y i Z N N O G Z j J D k h H K 0 N 9 w 4 0 I R t / 6 M O / / G T D s L b T 0 Q O J x z L / f k B I n g 2 j j O N y p t b G 5 t 7 5 R 3 K 3 v 7 B 4 d H 1 e O T r o 5 T R V m H x i J W / Y B o J r h k H c O N Y P 1 E M R I F g v W C 6 V 3 u 9 5 6 Y 0 j y W j 2 a W M D 8 i Y 8 l D T o m x k u d F x E y C M K u T y / m w W n M a z g J 4 n b g F q U G B 9 r D 6 5 Y 1 i m k Z M G i q I 1 g P X S Y y f E W U 4 F W x e 8 V L N E k K n Z M w G l k o S M e 1 n i 8 x z f G G V E Q 5 j Z Z 8 0 e K H + 3 s h I p P U s C u x k n l G v e r n 4 n z d I T X j j Z 1 w m q W G S L g + F q c A m x n k B e M Q V o 0 b M L C F U c Z s V 0 w l R h B p b U 8 W W 4 K 5 + e Z 1 0 r x p u s 3 H 7 0 K y 1 W k U d Z T i D c 6 i D C 9 f Q g n t o Q w c o J P A M r / C G U v S C 3 t H H cr S E i p 1 T + A P 0 + Q O p 4 Z F 0 < / l a t e x i t > 8w < l a t e x i t s h a 1 _ b a s e 6 4 = " e E C l v M M + L F w g C w 3 y S o 0 o n O g s U z Y = " > A A A B 6 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e x K w H g L e v E Y x T w g W c L s p D c Z M j u 7 z M w q I e Q P v H h Q x K t / 5 M 2 / c Z L s Q R M L G o q q b r q 7 g k R w b V z 3 2 8 m t r W 9 s b u W 3 C z u 7 e / s H x c O j p o 5 T x b D B Y h G r d k A 1 C i 6 x Y b g R 2 E 4 U 0 i g Q 2 A p G N z O / 9 Y h K 8 1 g + m H G C f k Q H k o e c U W O l + + p T r 1 h y y + 4 c Z J V 4 G S l B h n q v + N X t x y y N U B o m q N Y d z 0 2 M P 6 H K c C Z w W u i m G h P K R n S A H U s l j V D 7 k / m l U 3 J m l T 4 J Y 2 V L G j J X f 0 9 M a K T 1 O A p s Z 0 T N U C 9 7 M / E / r 5 O a s O p P u E x S g 5 I t F o W p I C Y m s 7 d J n y t k R o w t o U x x e y t h Q 6 o o M z a c g g 3 B W 3 5 5 l T Q v y l 6 l f H V X K d W u s z j y c A K n c A 4 e X E I N b q E O D W A Q w j O 8 w p s z c l 6 c d + d j 0 Z p z s p l j + A P n 8 w d i o I 1 J < / l a t e x i t > 8w < l a t e x i t s h a 1 _ b a s e 6 4 = " e E C l v M M + L F w g C w 3 y S o 0 o n O g s r 5 O a s O p P u E x S g 5 I t F o W p I C Y m s 7 d J n y t k R o w t o U x x e y t h Q 6 o o M z a c g g 3 B W 3 5 5 l T Q v y l 6 l f H V X K d W u s z j y c A K n c A 4 e X E I N b q E O D W A Q w j O 8 w p s z c l 6 c d + d j 0 Z p z s p l j + A P n 8 w d i o I 1 J < / l a t e x i t > 8w < l a t e x i t s h a 1 _ b a s e 6 4 = " e E C l v M M + L F w g C w 3 y S o 0 o n O g s U z Y = " > A A A B 6 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e x K w H g L e v E Y x T w g W c L s p D c Z M j u 7 z M w q I e Q P v H h Q x K t / 5 M 2 / c Z L s Q R M L G o q q b r q 7 g k R w b V z 3 2 8 m t r W 9 s b u W 3 C z u 7 e / s H x c O j p o 5 T x b D B Y h G r d k A 1 C i 6 x Y b g R 2 E 4 U 0 i g Q 2 A p G N z O / 9 Y h K 8 1 g + m H G C f k Q H k o e c U W O l + + p T r 1 h y y + 4 c Z J V 4 G S l B h n q v + N X t x y y N U B o m q N Y d z 0 2 M P 6 H K c C Z w W u i m G h P K R n S A H U s l j V D 7 k / m l U 3 J m l T 4 J Y 2 V L G j J X f 0 9 M a K T 1 O A p s Z 0 T N U C 9 7 M / E / r 5 O a s O p P u E x S g 5 I t F o W p I C Y m s 7 d J n y t k R o w t o U x x e y t h Q 6 o o M z a c g g 3 B W 3 5 5 l T Q v y l 6 l f H V X K d W u s z j y c A K n c A 4 e X E I N b q E O D W A Q w j O 8 w p s z c l 6 c d + d j 0 Z p z s p l j + A P n 8 w d i o I 1 J < / l a t e x i t > 8w < l a t e x i t s h a 1 _ b a s e 6 4 = " e E C l v M M + L F w g C w 3 y S o 0 o n O g s r 5 O a s O p P u E x S g 5 I t F o W p I C Y m s 7 d J n y t k R o w t o U x x e y t h Q 6 o o M z a c g g 3 B W 3 5 5 l T Q v y l 6 l f H V X K d W u s z j y c A K n c A 4 e X E I N b q E O D W A Q w j O 8 w p s z c l 6 c d + d j 0 Z p z s p l j + A P n 8 w d i o I 1 J < / l a t e x i t > (b)< l a t e x i t s h a 1 _ b a s e 6 4 = " S 6 z v 1 j t K s g 4 z J J D Z N 0 e S h M B T I x y g t A I 6 4 Y N W J m C a a K 2 6 y I T r D C 1 N i a K r Y E b / X L 6 6 R 7 1 f C a j d u H Z q 3 V K u o o w x m c Q x 0 8 u I Y W 3 E M b O k A h g W d 4 h T c n d V 6 c d + d j O V p y i p 1 T + A P n 8 w e r Z 5 F 1 < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 Y Z Z t j H A J b I C 3 D k F i T 6 3 6 8 d U A B 0 = " > A A A B 8 3 i c b V D L S g M x F L 3 x W e u r 6 t J N s A h 1 U 2 a k o O 6 K b l x W s A / o D C W T Z t r Q T G Z I M k I Z + h t u X C j i 1 p 9 x 5 9 + Y a W e h r Q c C h 3 P u 5 Z 6 c I B F c G 8 f 5 R m v r G 5 t b 2 6 W d 8 u 7 e / s F h 5 e i 4 o + N U U d a m s Y h V L y C a C S 5 Z 2 3 A j W C 9 R j E S B Y N 1 g c p f 7 3 S e m N I / l o 5 k m z I / I S P K L b t b 6 u 0 s b m 1 v V P e r e z t H x w e V Y 9 P e i p O J a F d E v N Y D n y s K G e C d j X T n A 4 S S X H k c 9 r 3 p 3 e 5 3 3 + i U r F Y P O p Z Q r 0 I j w U L G c H a S K 4 b Y T 3 x w 6 w e X M 5 H 1 Z r d s B d A 6 8 Q p FIG. 4. Imaginary part of 2D Fourier transform Ft,τ of the third-order susceptibilities χzzzz for X-cube model (K = 1) (a, b) without perturbations and (c, d) with generic perturbations.(a) Im Ft,τ [χzzzz(t, t + τ, t + τ )] and (b) Im Ft,τ [χzzzz(t, t, t + τ )] of X-cube model are qualitatively similar to those of Ising model except the signals marked with green circles.The green circled signals are coming from free motion of a pair of planons having restricted mobility within two-dimensional planes (Eµ = Eν = E λ process in Sec.IV C 2).Again, the signals at (ωτ , ωt) = (±4, 0) suggest deconfined excitations of X-cube model.Red circled signals are due to quantum processes involving four planons (Eν = Eµ + E λ process in Sec.IV C 2), and blue circled signals are originating from the process with two planons and two immobile fractons (Eν = Eµ + 2K = E λ + 2K process in Sec.IV C 2).(c) Im Ft,τ [χzzzz(t, t + τ, t + τ )] and (d) Im Ft,τ [χzzzz(t, t, t + τ )] of the perturbed X-cube model.The nonlinear signals from the dispersive two-planon (Eµ = Eν = E λ ) and two-planon-two-fracton (Eν = Eµ + 2K = E λ + 2K) processes are shown.Recall that the blue circled signals of Ising model (FIG. 1) are coming from a single zero momentum excitation of the two-magnon bound state.Hence, they do not show significant spreading in the presence of perturbations.However, the blue circled signals of X-cube model are originating from nonlinear dynamics of two finite momentum carrying planons scattering with two immobile fractons.Thus, we can see clear spread of the signals by the full bandwidth of the single planon dispersion 4w = 0.8.Blue and green stars mark locations of the signals without the perturbations.
t e x i t s h a 1 _ b a s e 6 4 = " m r U x V X w U b J U 2 4 f 2 n + e f C e / U j B c A = " > A A A B 8 3 i c b V D L S g M x F L 3 x W e u r 6 t J N s A h 1 U 2 a k o O 6 K b l x W s A / o D C W T Z t r Q T G Z I M k I Z + h t u X C j i 1 p 9 x 5 9 + Y a W e h r Q c C h 3 P u 5 Z 6 c I B F c G 8 f 5 R m v r G 5 t b 2 6 W d 8 u 7 e / s F h 5 e i 4 o + N U U d a m s Y h V L y C a C S 5 Z 2 3 A j W C 9 R j E S B Y N 1 g c p f 7 3 S e m N I / l o 5 k m z I / I S P K Q U 2 K s 5 H k R M e M g z G r k Y j a o V J 2 6 M w d e J W 5 B q l C g N a h 8 e c O Y p h G T h g q < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " g D J T Q U Q x D 9 + n 5 6 J d 2 S o x P m q O j B 0 = " > A A A B 8 3 i c b V D L S g M x F L 3 x W e u r 6 t J N s A h 1 U 2 a k o O 6 K b l x W s A / o D C W T Z t r Q T G Z I M k I Z + h t u X C j i 1 p 9 x 5 9 + Y a W e h r Q c C h 3 P u 5 Z 6 c I B F c G 8 f 5 R m v r G 5 t b 2 6 W d 8 u 7 e / s F h 5 e i 4 o + N U U d a m s Y h V L y C a C S 5 Z 2 3 A j W C 9 R j E S B Y N 1 g c p f 7 3 S e m N I / l o 5 k m z I / I S P K Q U 2 K s 5 H k R M e M g z G r B x W x Q q T p 1 Z w 6 8 S t y C V K F A a 1 D 5 8 o Y x T S M m D R V E 6 7 7 r J M b P i D K c C j Y r e 6 l m C a E T M m J 9 S y W J m P a z e e Y Z P r f K E I e x s k 8 a P F d / b 2 s n 8 A f o 8 w e s A Z F 3 < / l a t e x i t > (c) < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 Y Z Z t j H A J b I C 3 D k F i T 6 3 6 8 d U A B 0 = " > A A A B 8 3 i c b V D L S g M x F L 3 x W e u r 6 t J N s A h 1 U 2 a k o O 6 K b l x W s A / o D C W T Z t r Q T G Z I M k I Z + h t u X C j i 1 p 9 x 5 9 + Y a W e h r Q c C h 3 P u 5 Z 6 c I B F c G 8 f 5 R m v r G 5 t b 2 6 W d 8 u 7 e / s F h 5 e i 4 o + N U U d a m s Y h V L y C a C S 5 Z 2 3 A j W C 9 R j E S B Y N 1 g c p f 7 3 S e m N I / l o 5 k m z I / I S P K Q U 2 K s 5 H k R M e M g z G r 0 Y j a o V J 2 6 M w d e J W 5 B q l C g N a h 8 e c O Y p h G T h g q

5 F 4 <
FIG. S1.(a) After Z acts on the red link, different fracton configurations are available depending on which link is flipped subsequently.(b) If a blue link l ∈ B is flipped, a planon (face-sharing two-cube composite) is displaced by a unit distance.When (c) an orange link l ∈ O or (f) a green link l ∈ G is flipped, six fractons are excited.(d, e, g, h) These fractons can be grouped into two mobile planons (blue face-sharing cubes) and two immobile fractons (green cubes).While the two planons in (g) have restricted mobility within two parallel planes, the planons in (d, e, h) are constrained to two planes perpendicular to each other.If the secondly flipped link does not belong to the set of nearby links B ∪ O ∪ G, then total eight fractons are excited by the two spin flips.