Anomalous thermal relaxation and pump-probe spectroscopy of 2D topologically ordered systems

We study the behaviour of linear and nonlinear spectroscopic quantities in two-dimensional topologically ordered systems, which host anyonic excitations exhibiting fractional statistics. We highlight the role that braiding phases between anyons have on the dynamics of such quasiparticles, which as we show dictates the behaviour of both linear response coefficients at finite temperatures, as well as nonlinear pump-probe response coefficients. These quantities, which act as probes of temporal correlations in the system, are shown to obey distinctive universal forms at sufficiently long timescales. As well as providing an experimentally measurable fingerprint of anyonic statistics, the universal behaviour that we find also demonstrates anomalously fast thermal relaxation: correlation functions decay as a `squished exponential' $C(t) \sim \exp(-[t/\tau]^{3/2})$ at long times. We attribute this unusual asymptotic form to the nonlocal nature of interactions between anyons, which allows relaxation to occur much faster than in systems with quasiparticles interacting via local, non-statistical interactions. While our results apply to any Abelian or non-Abelian topological phase in two-dimensions, we discuss in particular the implications for candidate quantum spin liquid materials, wherein the relevant quantities can be measured using pre-existing time-resolved terahertz-domain spectroscopic techniques.


I. INTRODUCTION
Strongly correlated many-body systems in two spatial dimensions can host a remarkably rich variety of novel macroscopic quantum phenomena.Perhaps one of the most striking examples is the existence of emergent excitations that exhibit unconventional statistics-so-called 'anyons' [1,2].These quasiparticles are neither bosonic nor fermionic; rather, they possess nontrivial braiding statistics, meaning that the global wavefunction changes when one anyon moves along a path that encircles another.Remarkably, the wavefunction changes in the same way regardless of how far apart the anyons are throughout this process, which points to an effectively nonlocal interaction between excitations.This is only possible in systems whose ground states possess particular patterns of long-ranged entanglement; namely, in 2D topologically ordered phases [3,4].
Over the last several decades, a great deal of progress has been made in understanding the physics of anyons and the topological phases that host them.By now, there are a number of well-known phenomena that are established as being universal to 2D systems possessing excitations with fractional statistics: To name a few, ground state degeneracies appear on surfaces with nonzero genus [5]; quantum numbers can fractionalize [6,7]; and the entanglement entropy of large subregions contains a quantized topological contribution [8,9].These discoveries each provide important theoretical insight into the nature of topological order, and in some cases also serve as an experimental fingerprint of a given phase of matter.
In addition to the aforementioned properties, which pertain to equilibrium physics, one can also ask about the dynamics of systems with anyons.Besides transport measurements (which are challenging in systems with electrically neutral quasiparticles such as quantum spin liquids), the primary means of probing dynamics in solid-state systems is spectroscopy.Theoret-ical investigations into the behaviour of spectroscopic quantities in topologically ordered systems have begun comparatively recently, and for the most part the focus has been on linear spectroscopy, i.e. one analyses the signal using linear response theory.For instance, the spin structure factor in quantum spin liquids shows signatures of fractionalization [10][11][12][13][14][15][16][17], where excitations must be created in groups of at least two at a time.Similarly, it has been shown how fractional exclusion statistics (a consequence of anyonic statistics, generalizing Pauli's exclusion principle) can imprint themselves in absorption spectra [18].While these works provide useful insight into the nature of anyon creation and/or annihilation, there is only so much that can be learned about dynamics from linear response functions, which capture 'near-equilibrium' physics.
In this paper, we reveal universal dynamical phenomena associated with the braiding statistics of quasiparticles in topologically ordered systems, as witnessed by linear and nonlinear spectroscopic quantities.Our primary focus is on pumpprobe spectroscopy, where the system is perturbed by a series of two pulses, each of which excite quasiparticles.In particular, we highlight the significance of processes where anyons that were created at different times braid with one another-a possibility that does not arise in linear spectroscopy at zero temperature.As explained in a short paper that serves as a companion to this one [19], such processes dominate the latetime behaviour of the pump-probe response function, and the resulting signal takes a universal form [Eq. ( 4)], which constitutes an experimentally measurable signature of anyonic statistics.
One of our aims here is to present concrete calculations that support and generalize the results reported in Ref. [19], which were justified using more intuitive arguments, most of the time making reference to Z 2 quantum spin liquids.In brief, by considering the kinematics of those anyons generated by the sequence of pulses, we can compute the probability that their trajectories link in a way that leads to a nonzero braiding phase.Any such process gives a contribution to the pumpprobe response coefficient, and this is responsible for the universal form Eq. (4).Importantly, since anyons can braid without ever coming close to one another, the probability of braiding is always asymptotically higher than a scattering event due to short-ranged interactions between quasiparticles; therefore, as we shall argue, our result is robust against the inclusion of non-universal local interactions between excitations.
The pump-probe response coefficient is a particularly useful quantity in this context, since it allows one to isolate the effect that a single additional quasiparticles has on the motion of others.The insight we gain from studying pumpprobe spectroscopy is then applied to reveal salient features of linear response functions at finite temperature.Namely, we can consider the probability that anyons created by the timedependent perturbation braid with thermally activated quasiparticles.Again we find that these processes occur much more often than scattering does, which leads to an anomalously fast decay of the response function in the time domain: a 'squished exponential' form is seen C(t) ∼ exp(−[t/τ ] 3/2 ), for some temperature-dependent timescale τ [see Eq. (33)].This should be contrasted with the ordinary exponential decay that would be expected from local interactions.Since linear response functions serve as a quantifier of temporal correlations in the system, we conclude that topologically ordered systems exhibit much faster thermal relaxation that systems with quasiparticles having conventional statistics.
While we are not the first to study nonlinear spectroscopy in QSLs and other quantum magnets [20][21][22][23], previous works have focused on resolving the homogeneously broadened continuum of fractionalized excitations associated with the creation of multiple excitations, rather than detecting braiding statistics themselves.We also note that fractional exclusion statistics (a consequence of anyonic statistics, generalizing Pauli's exclusion principle) can imprint itself in absorption spectra even in the linear response regime, as discussed in Ref. [18].However, the the signal studied in this manuscript probes the braiding of excitations around one another, rather than the physics of their creation.Moreover, given that the pump-probe response coefficient involves the subtraction of two signals, one with a pump pulse and one without [see Eq. ( 3)], our approach has the advantage that the universal late-time behaviour can be disentangled from non-universal short-distance effects and background contributions, leading to a sharper signal.
Before embarking on any rigorous calculations, we begin our paper by specifying the systems and spectroscopic quantities that are to be studied in this work, and provide intuitive explanation of how the presence of anyonic excitations affects the signal measured in a pump-probe experiment.

A. Setup and Key Results
In this paper, we are concerned with gapped twodimensional systems, where quasiparticle excitations above the ground state can exhibit generalized statistics: indistinguishable particles can acquire exchange phases that interpolate between fermionic and bosonic, and mutual statistics can even be defined between distinguishable particles.We wish to study the dynamical response of these systems to external probes in regimes beyond linear response, and to understand how the mutual statistics of the lowest energy quasiparticles affects the relevant response coefficients.Primarily, we have in mind both mesoscopic systems in the quantum Hall regime and spin systems that are in (or proximate to) a spin liquid phase.
Our particular focus will be on the response of these systems to pulses of electromagnetic waves.In either of the aforementioned systems, the relevant energy scales correspond to a wavelength of light much greater than any realistic system size.Therefore, we restrict ourselves to external probes that are spatially homogeneous at zero wavevector k = 0 -that is, the operators to which the external electromagnetic fields couple are of the form where Â(⃗ r) is a Hermitian operator density.(On a lattice, the integral over space can be replaced by a sum over sites.)Later on, we comment on the possibility of accessing spatially resolved signatures either using inelastic neutron scattering rather than electron spin resonance, or moving to experimental platforms beyond solid state, e.g.ultracold atoms.
We mainly focus on a particular nonlinear response protocol known as pump-probe spectroscopy.Starting from the ground state (i.e. the quasiparticle vacuum) of the unperturbed Hamiltonian Ĥ0 , ρ 0 = |VAC⟩ ⟨VAC|, at time t = 0 the system is illuminated by a short, intense, 'pump' pulse of light which brings the state of the system out of equilibrium.Denoting the operator to which this pulse couples as Â0 , this results in an effectively instantaneous unitary rotation of ρ 0 for some constant κ controlling the strength of the pulse.After a time t 1 , a second 'probe' pulse is applied, whose purpose is to extract properties of the time-evolved nonequilibrium state.Deferring a proper treatment of the probe pulse and the relevant detection schemes to Section VI, for the time being we take it as given that the probe pulse allows one to extract the real part of the dynamical correlator ⟨ Â2 (t 1 + t 2 ) Â1 (t 1 )⟩ pert , where the expectation value ⟨ • ⟩ pert is taken with respect to the perturbed state in Eq. ( 2), and we work in the interaction picture with respect to Ĥ0 , i.e.Â1 (t 1 ) = e i Ĥ0t1 Â1 e −i Ĥ0t1 .The same experiment can be executed without the pump pulse and the results are subtracted to obtain a signal where the second expectation value is with respect to the original equilibrium state, which is independent of t 1 , and we have divided by the volume of the system L −2 such that χ PP is intensive.It is common practice in nonlinear spectroscopy to expand the signal in powers of κ; following standard nomenclature we write χ PP = ∞ n=1 κ n χ (n+1) PP . The coefficients χ (n+1) PP are nonlinear response functions of second order and higher.In particular, in the present setting the lowest order terms will turn out to be proportional to κ 2 , and therefore we write χ PP = κ 2 χ (3) PP that are valid at all times are prohibitively hard, and will depend on the details of the microscopic model in question.Nevertheless, here we argue that in systems where some excitations possess non-trivial braiding statistics, in the long-time limit t 1,2 → ∞ the response function χ (3) PP follows a universal behaviour.In particular, we will show that where χ (1) (t) = L −2 Tr Â2 (t) Â1 (0)ρ 0 denotes the linear response function, and c PP is a coefficient which depends on the details of the model and can generally be hard to explicitly compute.The relationship (4)-which is our main resultis a general feature of 2D systems whose excitations possess non-trivial braiding statistics, and therefore provides a powerful diagnostic tool to characterize fractional statistics using only pulses of light.The magnitude of the subleading term sets a timescale τ tr beyond which the transient effects represented by the o(t 2 ) term can be safely neglected and the ratio χ (3) PP /χ (1) takes its universal form = c PP t 3/2 2 ; this timescale will be characterised in later sections, see also Table I.
Most of the manuscript is dedicated to demonstrating the validity of Eq. ( 4), but first we find it instructive to review the following intuitive argument explaining this behaviour, which we reported in Ref. [19].In the following, and for most of our calculations, we will make explicit reference to systems where all anyons are Abelian, however the non-Abelian case can be treated in much the same way, as we show in Subsection IV E.
Due to their topological nature, quasiparticle excitations with non-trivial mutual statistics can only be created in multiplets of N > 1 particles by local operators.Let us focus on N = 2 for simplicity, and start by considering the behaviour of the unperturbed two-point function: the second term in Eq. (3).Since the expectation value is taken with respect to the quasiparticle vacuum, Â1 must create a quasiparticle pair at time t 1 and some position ⃗ r i , which will later be annihilated by Â2 at time t 1 + t 2 , position ⃗ r f [both ⃗ r i and ⃗ r f are to be integrated over according to Eq. ( 1)].Adopting a path integral formalism for this process, we must integrate over all possible trajectories of these particles ⃗ r 1 (t), ⃗ r 2 (t) for t ∈ [t 1 , t 1 + t 2 ], weighted by an appropriate action e iS[⃗ r1(t),⃗ r2(t)] ; these are drawn as blue lines in Fig. 1.
Supposing for now that the quasiparticles are free bosons S[⃗ r a (t)] = (m/2) dt (d⃗ r a /dt) 2 , then the amplitude can be evaluated exactly, and the result is proportional to e −2i∆t2 t −1 2 .The frequency of the oscillatory factor 2∆ is the energy required to excite two anyons relative to the quasiparticle vacuum, while the algebraic decay t −1 2 reflects the decreasing likelihood of finding two quasiparticles at the same point in space, which is necessary for them to be annihilated.
What changes when the pump pulse is applied beforehand?The post-pump state (2) contains additional quasiparticles, which we refer to as 'pump' quasiparticles, to distinguish them from the 'probe' excitations created by the probe pulse at time t 1 .In the absence of interactions (statistical or otherwise), the dynamics of the probe excitations are unchanged by the presence of these pump particles, and so the two terms in (3) exactly cancel.Now, suppose that the pump particles have non-trivial braiding statistics with respect to the probe particles.In this case, the action e iS[⃗ r1(t),⃗ r2(t)] must be multiplied by an extra statistical phase, equal to e 2πiα whenever a pump anyon passes through the spacetime loop formed by ⃗ r 1,2 (t) (see Fig. 1).Only trajectories that link in this way will contribute to χ PP , since the statistical phase prevents total cancellation of the two terms in (3); this is represented pictorially in the top right inset of Fig. 1.Therefore, to compute χ PP , we must integrate over ⃗ r 1,2 (t) as before, but now weighted by the probability that one of the excess pump anyons created by the pump pulse braids with the paths of the probe particles.
Working to leading order in κ, only a single pair of pump anyons will be created, with the quasiparticles being formed in wavepackets having opposite (crystal) momenta ⃗ k, − ⃗ k and being centred around some position ⃗ x i , which again is to be integrated over according to Eq. ( 1).These wavepackets propagate away from one another ballistically at their group velocities ±⃗ v = ± ⃗ ∇ k ϵ( ⃗ k), where ϵ( ⃗ k) is the single quasiparticle dispersion [24].The precise distribution of ⃗ k and the dispersion ϵ( ⃗ k) will depend on the microscopic model in question and details of Â0 , but this will not be relevant here; instead we can consider some fixed ⃗ v for now, and leave the averaging over ⃗ v at the end.Now we must integrate over ⃗ x i .Since the free action is independent of ⃗ x i , this gives a factor equal to the spatial area spanned by those initial positions for which the paths link (blue shaded region in Fig. 1).The component of ⃗ x i perpendicular to ⃗ v will be varied over a range of the order of the typical spatial separation of the two probe anyons ∼ |⃗ r 1 (t)−⃗ r 2 (t)|.By inspecting the free particle action, we see that for typical paths (those for which the phase does not oscillate too rapidly), this distance scales as ∼ t 2 /m in the long time limit.In the direction parallel to ⃗ v, a shift of ⃗ x i has the same effect as shifting the spacetime trajectory of the pump anyon upwards in the time direction (see Fig. 1).Therefore this component should be varied over a range ∼ |⃗ v|t 2 .Evidently, the space of initial positions ⃗ x i that yield linking trajectories has an area that asymptotically grows in time as It is this factor, coming from the integral over initial positions of the pump anyon, that leads to the universal form quoted in Eq. ( 4).Note that the average over ⃗ v does not have any 1. Schematic illustration of the processes contributing to the pump-probe response coefficient (3) in a (2 + 1)-dimensional spacetime, using a path integral picture.At time t = 0, the pump pulse generates a pair of pump anyons at position ⃗ xi, which in a semiclassical approximation propagate away from one another along trajectories with opposing velocities ±⃗ v (red lines).(We omit the backwards-time trajectory in this drawing, which brings these anyons back to their original position ⃗ xi; see Eq. ( 11).)A pair of probe anyons is created by the operator Â1 at time t1, position ⃗ ri, which are later annihilated by Â2 at time t1 + t2, position ⃗ r f .In a path integral formalism, the trajectories of the probe anyons are denoted ⃗ r1,2(t), and are drawn as blue lines.Statistical interactions between pump and probe anyons give rise to a phase e 2πiα whenever a pump anyon passes through the loop formed by the probe anyon trajectories.For a fixed ⃗ r1,2(t), we can integrate over all ⃗ xi such that the paths link.All other contributions cancel upon subtracting the terms in (3), as represented pictorially by the equation in the top right.For trajectories that contribute most to the path integral, the region of ⃗ xi satisfying this condition (light blue shaded region, dashed outline) has an area that scales as (see main text).This results in the asymptotic relation (4), valid in the limit of large t1,2.
bearing on the overall time-dependence; this simply controls the behaviour of the non-universal constant of proportionality c PP .
In the particular case we were considering, where N = 2 and there are no braiding statistics between the pairs of particles that are created at the same moment in time, we already saw that χ (1) up to an oscillatory phase factor, where the decay is due to the decreasing likelihood of anyon recombination.Hence, we have More generally, if anyons are created in multiplets of N > 2 particles, or if there are non-trivial statistics between particles in a given multiplet, then the recombination factor will be modified-see Sections II and IV D. Nevertheless, the t factor, which has a purely geometric origin, coming from the integral over ⃗ x i , remains the same.Thus, the relationship ( 4) is quite general.
While a number of assumptions have been made in this intuitive argument, these are not necessary for the relationship (4) to hold.Most notably, we have so far neglected nonstatistical interactions between quasiparticles, and assumed that the system is at exactly zero temperature.In Section IV, we will consider the effects of interactions and finite temperatures more quantitatively, but one can also understand the robustness of our result to such factors at the level of the above argument.Assuming that interactions are sufficiently shortranged (those decaying faster than ∼ |⃗ r 1 − ⃗ r 2 | −α at large separations, with α > 2 [25]), the presence of pump anyons can only appreciably affect the trajectories of the probe particles when the excitations are closer than some interaction radius r int .Using the same geometric approach as before, where one integrates over the initial coordinates of the pump particles keeping the probe anyons' trajectories fixed, the probability of these local scattering processes scales with the perimeter of the loop formed by ⃗ r 1,2 (t) [26].This gives a correction that is subleading compared to the long-ranged statistical interactions, where the relevant probability scales with the area (t At finite temperature, the presence of thermally excited quasiparticles (in addition to those created by the pump pulse) modifies the linear response coefficient χ (1) (t), since braiding between the trajectories of the probe anyons and the thermal excitations leads to an effective dephasing of the two-point correlator ⟨ Â2 (t) Â1 (0)⟩.However, the pump-probe response function will be modified in exactly the same way.While scattering between thermal and pump quasiparticles may alter the effective distribution of velocities, this only changes c PP , and so Eq.(4) continues to hold.This is shown explicitly later [Eqs.(33,34)].
This concludes our overview of the universal behaviour of the pump-probe response function.In summary, the-late time form of χ PP obeys a universal relationship Eq. ( 4), which can be understood as described above using a semiclassical picture.The structure of the remainder of our paper is as follows: To justify our intuitive arguments, in Section II we compute the main quantity of interest, namely the leading order contribution to χ PP (t 1 , t 2 ) [Eq. (3)], using an effective low-energy theory for a system with anyonic excitations.In Section III, we go beyond time-dependent perturbation theory to obtain the full response coefficient at all orders; doing so resolves an apparent paradox that the leading order contribution has an unphysical divergence in the long time limit.In Section IV, we discuss other effects that could not be included in our rigorous calculation, focusing on non-statistical interactions, finite temperatures, and non-Abelian statistics.To make connection between the low-energy theory used before and concrete microscopic models, in Section V we apply our results to the toric code model in a weak magnetic field, allowing us to connect phenomenological parameters with microscopic quantities.Finally, we discuss how the signal can be measured experimentally in Section VI, before concluding in Section VII.

II. CALCULATION OF NONLINEAR RESPONSE FUNCTION A. Effective low-energy theory
To begin a calculation of the pump-probe response coefficient, we will require a more detailed characterization of the operators Â0,1,2 appearing in Eqs.(2,3), which create and annihilate anyons, as well as a description of how anyons propagate once generated.For the systems we consider in this paper, the lowest-energy excitations are deconfined quasiparticles, which are separated from the ground state by a finite energy gap ∆ n > 0, where the label n is used to distinguish different quasiparticle species.Assuming translation invariance, we can specify a dispersion for each quasiparticle ϵ n (k).For the time being, we assume that the only interactions between anyons come through their braiding phases: the wavefunction acquires a phase of e 2πiα nn ′ when a particle of type n completes a loop that encircles a particle of type n ′ once in an anticlockwise direction.Later we will include the effect of additional short-range interactions, which do not modify the qualitative form of the response coefficients.
Our analysis applies to 2D topological phases in general, but it will often be helpful to make reference to a particular phase of matter as an example.For this purpose we consider the phase of matter in which the toric code lies [27,28].Systems in this universality class possess two types of excitations, known as electric and magnetic anyons (e and m respectively).While the electric-electric and magnetic-magnetic braiding phases are trivial α ee = α mm = 0, these particles are mutual semions with respect to one another α em = α me = 1/2.The toric code Hamiltonian is an exactly solvable model with these kind of excitations.At this fine-tuned point, anyons are motionless once created, meaning the dispersion is flat ϵ(k) = 0.However, perturbations that are weak compared to the excitation gap generically induce some dispersion, which endows these excitations with dynamics.In Section V, we will consider a specific perturbed toric code model, allowing us to relate our universal results to microscopic parameters.
In general, local operators can only excite quasiparticles in multiplets N := {n 1 , . . ., n N } that are statistically neutral with respect to all excitations when considered as a composite (i.e.N j=1 α nj n ′ ∈ Z for all n ′ ).For example, in the toric code the pairs {e, e} and {m, m} can be created locally, since braiding two electric anyons around a magnetic anyon gives a trivial phase of 2π.However, individual electric anyons {e} cannot be created locally, since they are not neutral with respect to the magnetic anyon.We can associate a threshold energy ∆ N = N j=1 ∆ nj to each valid N , which is the minimum energy required to create all the particles in the multiplet.For simplicity we will assume that different multiplets have threshold energies that are well-separated, although we expect that the existence of energetically degenerate multi-plets will not wash out the universal signal that we derive here.
The external probes we consider here will have frequencies that are close to these quasiparticle creation thresholds ∆ N .More formally, writing the microscopic light-matter coupling as a term in the Hamiltonian f (t) Âmicro , we take f (t) = ℜe −iω0t f 0 (t), where |ω 0 − ∆ N | ≪ ∆ N , and the function f 0 (t) varies on a timescale much longer than ω −1 0 .While the microscopic operator Âmicro could in principle connect the ground state to complicated states with a larger number of quasiparticles, these components oscillate quickly in the interaction picture, and hence can be ignored (provided one is interested in dynamics on timescales longer than ∆ −1 n ).After discarding these rapidly oscillating terms, the resulting Hamiltonian only contains operators Â0,1,2 that couple quasiparticle sectors differing by the creation/annihilation of the relevant multiplets.Furthermore, the discrepancy (ω 0 − ∆ N ) sets an amount of excess kinetic energy that the quasiparticles will have once created.We will assume that this energy is small enough such that the quasiparticle dispersions can be expanded to quadratic order about the band minimum (Anisotropy in the dispersion can also be accounted for in principle, however this will simply result in a rescaling of the pump-probe response function.)We make the above choices in order to progress with our analytical calculation, but we stress that the universal physics discussed here does not depend on the restrictions that we presently impose on the frequency profile of the pump pulse.Indeed, the response to a pulse with a broader range of frequencies will still include the signal we derive here, in addition to non-universal transient effects coming from other mechanisms.In this section, the only interactions between quasiparticles will be due to braiding phases only, and the effect of non-statistical interactions will be treated in Section IV A. Deferring a discussion of the effects of finite temperatures to Section IV B, we assume that the system is in its ground state ρ 0 = |VAC⟩ ⟨VAC| before any of the pulses have arrived, i.e. there are no quasiparticles present.Acting with one of the operators Â0,1,2 on the quasiparticle vacuum, we obtain a perturbed state |Ψ N ⟩ = Â0 |VAC⟩, where |Ψ N ⟩ is some translationally invariant wavefunction in the excitation sector with a single multiplet N .For the time being, we will add one additional restriction, namely that the particles within the multiplets N created by Â0,1,2 are statistically neutral with respect to one another.This does not preclude nontrivial braiding between excitations in different multiplets N , N ′ : For example, in the context of the toric code, we can consider N = {e, e} and N ′ = {m, m}.(It will be useful to keep this example in mind in the following.)The reason we make this assumption is that when particles within the set N possess mutual braiding phases, a short-distance regulator for the operator density A(⃗ r) appearing in Eq. ( 1) must be introduced, since such particles cannot be at the same point in space (otherwise the wavefunction would be ill-defined).There is some freedom in choosing this regulator, and non-universal features of the initial N -particle wavepacket may affect the subsequent dynamics.We address the more general case in Section IV D, but for now we can choose a simple form for |Ψ N ⟩, where the N pump particles begin in wavepackets localized at the same point in space, i.e.
In reality, anyons will not be perfectly pointlike but will have some characteristic size ξ that acts as an ultraviolet cutoff.We will eventually need to invoke this lengthscale to regularize divergent integrals in momentum space, but for now we can assume that anyons generated by each of the perturbing operators Â0,1,2 will be created and annihilated at the same point in space.
The post-pump state takes the form given in Eq. ( 7) for typical zero-momentum operators Â0 [Eq.(1)].However, one should bear in mind that in certain scenarios there may be selection rules that further constrain how the system is perturbed by the external pulses, besides those imposed by the fusion rules associated with the underlying topological order.For example, in a spin-half system with unbroken SU(2) spin-rotation invariance, the only translation invariant operators made up of single-site terms are the total magnetization operators M α := j Ŝα j , where Ŝα j is the spin operators for lattice site j along the quantization axis α; indeed, in electron spin resonance experiments this is the most natural operator to which light will couple.However, since M α generates the symmetry group, the ground state |VAC⟩ will be unperturbed by the pulse and no signal would be seen [29].In this specific case, one must either account for the small nonzero wavevector of the light pulse, or identify other microscopic operators to which light couples; for instance, the coupling operator describing Raman scattering at q = 0 is a spin bilinear, and hence not a symmetry generator [10,30].(For smaller symmetry groups, one can always find a polarization of light α such that excited states (7) are created by M α .)From hereon, we will assume that non-symmetry-generating coupling operators Â0,1,2 have been identified, for which the selection rules are not so stringent so as to prevent creation of anyons; thus Eq. ( 7) can be used.
Having specified the action of the operators Â0,1,2 within our low-energy effective description, we are now in a position to explicitly calculate response functions, starting with the simplest case, namely linear response.

B. Warm-up: Linear response
Before embarking on our calculation of the pump-probe response coefficient, it is useful to consider the behaviour of the second term in Eq. (3), i.e. the two-time correlation function in the absence of a pump pulse.This is effectively the linear response coefficient χ (1) (t) := ⟨ Â2 (t) Â1 (0)⟩.A common approach to calculating linear response quantities is to first calculate the Fourier transform of χ (1) (t) using a spectral representation.In the present setting, the excitations that can be created and annihilated by the operators Â1,2 [which have zero momentum; Eq. ( 1)] form a N -particle continuum spanned by plane wave states | ⃗ k 1 , . . ., ⃗ k N ⟩, subject to the condition

N j=1
⃗ k j = ⃗ 0 that is imposed due to conservation of momentum.The spectral density of these states exhibits nonanalytic behaviour at a frequency equal to the gap ∆ N .In the simplest case N = 2, a stepwise discontinuity appears, and this same kind of discontinuity will generically be present in the Fourier transformed linear response function.Transforming back to the time domain, this behaviour dictates that the late-time form of χ (1) (t) is proportional to e −i∆ N t t −1 , as quoted in Section I A.
Later, we will study the behaviour of the pump-probe response coefficient using a time-domain approach based on semiclassical trajectories.It is therefore worthwhile rederiving the above form using such a real-time picture.The effect of the operator Â1 (0) is to create a pair of quasiparticles in the state (7) at time t = 0.The wavefunction of the quasiparticles can be decomposed into wavepackets that have centre of mass position ⃗ x i and opposing momenta ⃗ k and − ⃗ k, where both ⃗ x i and ⃗ k are to be integrated over [31].In the semiclassical limit ℏ → 0, these quasiparticles propagate away from one another at their group velocity ⃗ v = ± ⃗ ∇ k ϵ( ⃗ k), and so their separation grows in time like 2|⃗ v( ⃗ k)|t.
If we modelled these wavepackets as perfectly pointike (as we would for classical particles), then we would not find any signal for large t, since the quasiparticles must be within some fixed distance of each other to be annihilated by the operator Â2 (0).However, quantum effects lead to a broadening of the profile of these wavepackets: they are not perfectly pointlike, but rather their width grows as ∼ ℏt/m (restoring ℏ for now).Consequently, at any given time t, quasiparticles with momenta that satisfy 2|⃗ v( ⃗ k)|t ≲ ℏt/m will have a finite amplitude of annihilation, and so contribute to χ (1) (t).Expanding ⃗ v( ⃗ k) ≈ ℏ ⃗ k/m for small ⃗ k, we see that the momenta giving a non-negligible amplitude have a magnitude ≲ ℏm/t, and such points occupy an area ∝ 1/t in 2D momentum space.If the quasiparticles within this multiplet are mutually bosonic, as in Eq. ( 7), then the integrand is approximately constant in this region, and we find χ (1) (t) ∼ 1/t as quoted before.
While we will assume trivial statistics within multiplets in the following pump-probe calculation, incidentally we can also use the above picture to understand the behaviour linear response coefficient when there are non-trivial exchange or braiding statistics between quasiparticles created at the same time.In this case we must be more careful in accounting for the matrix elements which controls the distribution of quasiparticle momenta created and annihilated by Â1,2 .If the particles are not mutual bosons, Pauli exclusion (or its generalization to anyons) prevents creation of two plane-wave states at the same momentum, and so the matrix element must vanish at ⃗ k = 0.For fermions, one readily finds ⟨ ⃗ k, − ⃗ k| Â1 |VAC⟩ ∼ | ⃗ k| at small | ⃗ k|, and a calculation analogous to that appearing in Ref. [18] generalizes this to | ⃗ k| α for anyons (subject to certain condi-tions on the structure of Â1 ; see Section IV D).This gives χ (1) (t) ∝ t −1−α , which is consistent with the results of Ref. [18].Additionally, if N > 2 mutually bosonic quasiparticles are created at the same time, then similar arguments can be used to show χ (1) (t) ∝ t −N +1 .The structure of matrix elements for N > 2 particles with non-trivial mutual statistics is more complicated; see Ref. [18].
Regardless of the intra-multiplet statistics, the key insight to take from the above is that to properly capture the latetime behaviour of the two-time correlator ⟨ Â2 (t) Â1 (0)⟩, we must account for quantum fluctuations about the semiclassical trajectories, i.e. the broadening of wavepackets as ∼ ℏt/m.

C. Pump-probe response
Now we turn to the full pump-probe response coefficient, Eq. ( 3), working perturbatively in the strength of the pump pulse κ.Here, we must distinguish the multiplet N created by the pump pulse via the operator Â0 from the multiplet N ′ created by the probe pulse operators Â1,2 .We assume that the pump and probe pulses have frequency profiles overlapping with the corresponding threshold energies ∆ N and ∆ N ′ , which may be different.Accordingly, we can again infer that each operator either creates or annihilates a multiplet, and so the leading order contributions come at second order in κ.Following standard naming conventions for nonlinear response coefficients [32], we define the perturbative response coefficient χ By Taylor expanding the exponentials in Eq. ( 2), we obtain These low-order contributions dominate the response in the limit of a weak pump pulse κ → 0, and so we will focus on them for now.However, it is important to bear in mind that the weak pulse limit does not commute with the long time limit t 2 → ∞, as will be clear once we derive the divergent growth χ [Eq. ( 5)].We will remedy this issue in Section III, where we calculate an expression for χ PP (t 1 , t 2 ) that includes contributions at all powers of κ, and thus remains valid as t 2 → ∞.
The quantity (9) describes a process where a multiplet N is created at time 0, followed by a multiplet N ′ at time t 1 , which is then annihilated at t 2 .This is precisely the process that was central to our intuitive argument in Section I A (see Fig. 1).Using the form of Â0 given above [Eq.( 7)], the response coefficient can be written as Tr[ζ Â2 (t 1 + t 2 ) Â1 (t 1 )], where we define ) Being unnormalized and not positive-definite, the operator ζ is not itself a valid density matrix; rather, it includes only the contributions to the pumped state (2) that are second order in κ.Nevertheless, it is helpful to think of the perturbative response coefficient as the expectation value of Â2 (t 1 + t 2 ) Â1 (t 1 ) with respect to a 'state' ζ, as one would if we were calculating the full response to all orders in κ.The contributions to this expectation value coming from each of the terms in (10) are represented pictorially in the top right inset of Fig. 1: in the first term, pump and probe anyons are both generated, whereas in the second term, the probe anyons are created on top of the vacuum, and the pump anyons only appear through the multiplicative factor ⟨Ψ N |Ψ N ⟩.
So far, we have not described how the statistical interactions between particles (specifically those between pump and probe excitations) can be included in our description.For this purpose, we find it useful to work in a path integral representation, which we now introduce.

D. Path integral representation of χ
(3) PP Using Eq. ( 7) and the local form of the operators Â1,2 [Eq.( 1)], we can express the response coefficient using a Feynman-Vernon functional integral for the dynamics [33], where both the forward and backward branches of the time evolution in (9) are expressed as a sum over paths.We consider all trajectories of the particles in N between times 0 and t 1 + t 2 , along with those of particles in N ′ between times t 1 and t 1 + t 2 .If the only interactions are statistical in nature, then the action can be written as a sum of the free particle actions S j [⃗ r j (t)] = (m j /2) dt ṙ2 plus a topological term Λ[{⃗ r(t)}] equal to the cumulative statistical phases associated with the braiding of probe anyons around pump anyons.An explicit formula for Λ[{⃗ r(t)}] will not be necessary, however we will later use the fact that Λ only depends on the relative coordinates between pump and probe anyons.Overall we have Here, ⃗ x + (t) and ⃗ x − (t) are the trajectories that describe the forward and backwards time evolution in (9), respectively.
(For clarity, we consistently use ⃗ x with appropriate subscripts to denote coordinates of pump anyons, and ⃗ r for probe anyons.)Note that no probe anyons are present on the backwards branch, and so the statistical phase Λ has no dependence on ⃗ x − j (t) and ⃗ r k (t).The above expression is an explicit representation of the processes illustrated in Fig. 1 (although the backwards trajectories are not drawn explicitly).The factor of (e iΛ −1) arises due to the subtraction of the two terms in ζ [Eq.(10)]; see the pictorial equation in the top right corner of Fig. 1.
Unfortunately, exact analytical calculations of the dynamics between times t 1 and t 1 + t 2 quickly become intractable as the number of particles increases.Even in the minimal case where |N | = |N ′ | = 2, the evaluation of the four-body path integral including the statistical interactions does not admit a closed-form solution.However, in the limit of long times t 1,2 , we can make two simplifications.Firstly, at sufficiently large t 1 we can consider just one of the pump anyons at a time.We make this approximation on the basis that in the long-time limit, the pump anyons will typically be separated by a large distance, and so the amplitude for the probe anyons braiding around more than one pump anyon is small.The result is that the statistical factor (e iΛ − 1) in ( 11) can be replaced by a sum where the new topological term Λj [{⃗ r k (t) − ⃗ x j (t)} k ] captures the statistical phase associated with the braiding of probe anyons around a single pump anyon j [34].
Our second simplification is to invoke a stationary phase approximation for the trajectories of the pump anyon j.To be specific, we decompose the path , plus a fluctuating part δ⃗ x + j (t), and the free part of the action then becomes m⃗ v 2 j (t As we argue in Appendix A, the dependence of the topological part of the action on δ⃗ x(t) can be neglected in the limit of large times, with relative corrections decaying at least as fast as O(t −1 2 ), i.e. we can take the trajectory of the pump anyon to be of constant velocity.The fluctuations δ⃗ x + j (t) can then be integrated over, along with the backwards trajectory ⃗ x − j (t) and its initial position ⃗ x − i , all of which can be expressed using the Feynman propagator.This leaves us with a manageable expression for the pump-probe response function Note that in the regime where the above applies, the pumpprobe coefficient has no t 1 -dependence.This is due to the constant-velocity nature of the pump anyon trajectories, meaning that any change of t 1 → t 1 + ∆t 1 can be thought of as equivalent to a rigid shift of ⃗ x cl,j (t) → ⃗ x cl,j (t) + ⃗ v∆t 1 .This is borne out in the above since the path integral over δ⃗ x(t) and ⃗ x − j (t) is proportional to (t 1 + t 2 ) −2 , which cancels with the factor of (t 1 + t 2 ) 2 that comes from the change of integration variables from ⃗ x f to ⃗ v. Additionally, the classical contribution to the action m⃗ v 2 (t 1 + t 2 )/2 cancels with the opposite phase coming from the backwards trajectory, which is why a factor of e imv 2 (t1+t2)/2 does not appear in (13).
Eq. ( 14) describes the propagator for N ′ probe particles moving from ⃗ 0 to ⃗ r f in the presence of a pump anyon whose trajectory is fixed, and given by ⃗ r cl (t) = ⃗ vt + ⃗ x i .Observe that we have made a semiclassical approximation for the path of the pump anyons and not the probe anyons.This is motivated by the insight gained from Section II B, where we saw that the behaviour of two-time correlation functions requires one to account for fluctuations of the relevant excitations; see also the discussion of Appendix A.
In Section II E, we will directly evaluate I(⃗ v, t 2 ), but for now it is helpful to briefly make connection with the arguments that we gave in Section I A to justify the scaling form (5). Evidently, the integral over ⃗ x i in ( 14) is precisely the integral that was responsible for the factor of t 3/2 2 in (5), and we can move it inside the path integral over ⃗ r k (t).Being a topological term, Λ only takes a finite number of distinct discrete values, and so we can split up the integral d 2 ⃗ x i (e i Λ −1) into patches where e i Λ takes different values, to get where c labels the distinct values Λc that the functional Λ can take, and A c is a functional of ⃗ r k (t) − ⃗ vt, equal to the (unsigned) area in the space of coordinates ⃗ x i that satisfy While we do not have a closed-form expression for A c , the intuitive arguments in Section I A indicate that this should scale as t 3/2 2 , and this will be backed up by our exact calculations.In fact, the scaling of χ PP (t 1 , t 2 ) can be seen fairly straightforwardly using the geometric interpretation offered by Eq. (15).First, note that since the only free parameters in this problem are ⃗ v, t 2 , and the quasiparticle masses m k , by dimension counting I(⃗ v, t) can only depend on velocity and time through the product |v| √ t 2 .Hence, the long-time limit is equivalent to the large-velocity limit.When we take |⃗ v| → ∞, the pump anyon will only ever be in the vicinity of the probe anyons for a short O(v −1 ) period of time.The winding number will then be entirely determined by the location of the probe anyons at this instant in time, which we call τ .In the case N ′ = 2, the trajectories contributing to the area A c in Eq. ( 15) are those where the ray traced by the fast pump anyon passes between two probe anyons at locations ⃗ r 1,2 (τ ).Thus, the component of ⃗ x i perpendicular to ⃗ v should be varied over a distance equal to |r 2,⊥ (τ Varying the component of ⃗ x i parallel to ⃗ v only changes the collision time τ , and so I(⃗ v, t) is given by the path integral of By evaluating the integral over trajectories ⃗ r k (τ ), one can show that this quantity is indeed proportional to t 1/2 2 , confirming Eq. ( 5).When N ′ > 2, an similar path integral describing the longtime limit of I(⃗ v, t 2 ) can be constructed, but the expression becomes more complicated.To determine I(⃗ v, t 2 ) in full generality, and to remove the need to rely on dimension-counting arguments, it is more convenient to return to the Schrödinger picture, wherein Eq. ( 14) can be computed exactly.
E. Evaluating Eq. (14) Our objective is now to evaluate the function I j (⃗ v, t) defined in (14).The action describes N ′ probe particles propagating in the presence of the pump anyon j, which moves along a fixed-velocity trajectory.Since the probe anyons are mutually non-interacting, we can consider the propagator for a single probe anyon where U k (t f ; t i ) is the unitary operator describing time evolution of particle k under the influence of the moving pump anyon from time t i to t f .In terms of this propagator, we have where is the propagator without the pump anyon.
Naturally, it is helpful to perform a Galilean boost to a frame moving with velocity ⃗ v relative to the laboratory frame.
We have where G is the propagator in the co-moving frame.In this frame, the pump anyon is static, and so we are free to place it at the origin.We will adopt polar coordinates (r, ϕ) with the x axis in the direction of ⃗ v.
A standard way to describe the effect of the pump anyon is to introduce an infinitesimally thin flux tube at the origin, whose strength is chosen such that an Aharonov-Bohm phase of 2πα jk is acquired every time particle k orbits around it.Any vector potential describing such a magnetic field will satisfy Γ d⃗ r • ⃗ A(⃗ r) = 2πα jk for any loop Γ circling the origin in an anticlockwise sense.It will be useful to start in the 'string gauge', where ⃗ A(⃗ r) is only on the negative yaxis, specifically ⃗ A(⃗ r) = Θ(−y)δ(x)x, where x is a unit vector in the direction along ⃗ v.We can then perform a gauge transformation ψ(r, ϕ) → e 2πiα jk ϕ ψ(r, ϕ), where we restrict ϕ ∈ (−π/2, 3π/2].This completely eliminates the vector potential at the expense of introducing twisted boundary conditions for all wavefunctions.In particular, wavefunctions can be assumed to be continuous functions of ϕ except at ϕ = 3π/2, where we have Since the statistical vector potential vanishes in the chosen gauge, the boosted Hamiltonian for particle k = 1, . . ., N (indexing the probe anyons) becomes To calculate the propagator for this Hamiltonian, we first have to construct all its eigenstates, subject to the boundary conditions imposed by anyonic statistics (18).This is most easily achieved by using a unitary transformation H ′ j = U † H boost U , where U = e imj xv shifts the momentum operator by m k v, which gives , where L = −i∂ ϕ is the angular momentum operator.The boundary condition (18) imposes that L must take values of ℓ − α k , where ℓ is an integer (we drop the label for the pump anyon j on all quantities for the time being).The radial part of the wavefunction must then satisfy Bessel's equation with constant (ℓ − α k ) 2 .The overall solution is with energy E q,ℓ = which, with the normalization given, form a complete set of states: The precise structure of these eigenstates stems from our assumption that the Hamiltonian in the boosted frame is rotationally invariant.This allows us to make analytical progress in the following, but we wish to highlight that the late-time form of the response coefficient will be qualitatively unchanged if rotational symmetry is broken, e.g.due to anisotropy in the dispersion ϵ n ( ⃗ k).We now have The integral over q can be evaluated using the standard integral The infinitesimal constant ensures that this expression remains valid in the limit t → 0, and we have used the identity I ν (ix + 0 + ) = e iνπ/2 J ν (+x) for real x > 0.
Our aim now is to evaluate this integral.As noted previously, on dimensional grounds I(v, t 2 ) can only depend on velocity and time through the combination v √ t 2 , and so the late-time limit can be understood by considering the behaviour as v → ∞.More quantitatively, the relevant dimensionless parameter in the problem is β := mv 2 t 2 , and so we expect the response to take its asymptotic form when β ≫ 1, i.e. t 2 ≫ τ tr , where τ tr is a timescale on the order of 1/mv 2 (see Table I).In this limit, the integrand becomes a rapidly oscillating function of ϕ i,f , which motivates a stationary phase approximation of these integrals.Points of stationary phase occur at ϕ i,f = 0, π, and of the four different combinations, the one that that gives a dominant contribution to I(v, t 2 ) is ϕ i = 0, ϕ f = π, i.e. the particles begin far along the positive x-axis and drift at a velocity v until they reach the negative x-axis.Details of the evaluation of this integral are given in Appendix B, the result of which gives where we have defined a topological quantity Performing the necessary integral over v (which should be cut off at large velocity ∼ 1/ξm to account for the finite spread of wavevectors created by the pulse), we obtain a pump-probe response coefficient χ PP (t 1 , t 2 ) that scales as t . A straightforward calculation gives the linear response coefficient χ (1) , and comparing the two we see agreement with the form originally stated in Eq. (4).

III. RESPONSE BEYOND PERTURBATION THEORY
As we showed in the previous sections, for the case N ′ = 2 the lowest order contributions to the pump-probe response coefficient χ (3) PP (t 1 , t 2 ) grow as t 1/2 2 in the limit t 2 → ∞.The fact that this quantity diverges at late times indicates that a perturbative expansion of the system's response to the external fields begins to fail.Specifically, if the late-time limit is taken while holding κ > 0 fixed, then higher order terms in Eq. ( 9) cannot be ignored, and the whole series must instead be resummed.In this section, we derive an expression for the full response of the system without relying on perturbation theory, using arguments that generalize those given above.The result, Eq. ( 30), remains valid in the long time limit for fixed κ.
When considering higher order contributions to χ PP (t 1 , t 2 ) [Eq. ( 3)], the main difference in our analysis is that we must consider the possibility that the pump pulse creates more than a single quasiparticle multiplet.Since the frequency of the pump pulse is tuned close to the threshold energy ∆ N (which we assume is not close to any other excitation threshold energy), the terms of order κ 2n will involve the creation of up to n copies of N .For the time being, we will continue to assume that excitations of the system interact with one another only through their statistical interactions, and that there are no nontrivial statistics among particles within either the pump or probe multiplet (generalizations of this scenario are addressed in Section IV).Therefore, if we work in the path integral formalism as in Section II D, for each trajectory of the probe anyons we can identify contributions where a particular number of pump anyons pass through the loop formed by ⃗ r 1,2 (t), and an appropriate statistical phase can be assigned to each contribution.Specifically, we can separate out processes where p pump anyons pass through the loop from one side, and p ′ from the opposite side, which yields a phase of e 2πiα(p−p ′ ) .Our task now is to determine, for each possible trajectory of the probe anyons ⃗ r 1,2 (t), the probability that the pump anyons follow paths such that this linking condition is satisfied.We denote this probability Q p,p ′ [⃗ r 1,2 (t)].The full non-perturbative response will then be given by the path integral over the probe anyon trajectories weighted by a factor of ⟨(e iΛ − 1) where again the subtraction of unity is due to the unperturbed correlator in (3).(The angled brackets ⟨ • ⟩ pr indicates that this averaging is being performed over the paths of the probe anyons.) Our previous perturbative calculation informs us that the probability to generate a single pump anyon that links with the loop in a particular sense is proportional to κ 2 A, where A is the area functional A[⃗ r(t) − ⃗ vt] integrated over velocities ⃗ v, which comes from integrating over the initial positions of the pump multiplet [see Eq. ( 15)].The pump pulse can produce many pump multiplets which are created and propagate approximately independently of one another (assuming their density is low enough), and so the probability that p particles link in a given sense will follow a Poisson distribution, with rate cκ 2 A for some constant c, i.e.
Applying the same logic to the paths that link in the opposite sense gives us an expression for Q p,p ′ [⃗ r 1,2 (t)].Thus, the trajectory of each probe anyon should be weighted by a factor Recalling that A is a functional of ⃗ r 1,2 (t), we must now perform the path integral over the trajectories of the probe anyons.Our previous arguments can be reapplied here, which tell us that for typical paths, A ∝ t 3/2 .The full response coefficient is now given by the same path integral expression as the linear response coefficient χ (1) (t 2 ), but with the additional weighting of (e iΛ − 1) pr , giving where the prefactor in the exponent is identified as the same constant c PP appearing in Eq. ( 4), to ensure agreement with our perturbative results upon expanding (30) to leading order in κ.Note that χ (1) (t 2 ) is bounded in the long-time limit, and so this nonperturbative expression for the pump-probe response coefficient does not diverge, in contrast to χ PP .Evidently, once short-time transient effects have decayed away, the ratio χ PP /χ (1) will depend on time only through a universal function of κ 2 t 3/2 2 , after choosing units where c PP = 1.The factor inside the square brackets in Eq. ( 30) is plotted in Fig. 2 for various values of κ.
The linear response coefficient itself is most easily evaluated in the case where N ′ = 2, and there are no non-trivial braiding phases between anyons created in the same multiplet (this was the case in the toric code example discussed in previous sections).There, one has χ (1) (t) ∝ t −1 , and hence the pump-probe response takes the form χ PP (t 1 , t 2 ) ∝ t −1 2 (e −cPPκ 2 t 3/2 2 − 1).This signal grows as √ t 2 for times much less than τ non−pert ∼ (c P P κ 2 ) −2/3 , after which nonperturbative effects become important.At late times, the pump anyons have such a strong effect that the phase coherence of the two-point function is completely lost, and hence the first term in (3) completely decays away.This leaves only the second term, which is the unperturbed correlation function, decaying as t −1 2 .Interestingly, even though the leading order perturbative response coefficient χ (3) (t 1 , t 2 ) does not diverge when N ′ > 2, our analysis shows that higher order terms, e.g.χ (5,7,...) , will always diverge for times beyond τ non−pert ; this can be understood as a consequence of the long-ranged nature of the interactions between anyons.
To summarise, the picture provided by these arguments is that the population of anyons produced by the pump pulse have the effect of dephasing the trajectories of the probe Late-time form of the ratio χPP(t1, t2)/χ (1) (t2), where χPP(t1, t2) is the full nonlinear response coefficient, including contributions at all orders in κ, Eq. ( 30).We use units where the nonuniversal constant cPP = 1, and vary κ 2 from 0.4 (blue) to 0.1 (orange) in steps of 0.1.Initially, the ratio of the response coefficients grow as t 3/2 2 (dashed line), in agreement with the perturbative expression derived in previous sections, see Eq. ( 4).After some timescale τnp ∝ κ −4/3 , nonperturbative effects become important, and we see a plateau of the ratio.
anyons through their mutual statistical interactions.This induces a relative suppression of the two-time correlator compared to its unperturbed value, which leads to a non-zero response coefficient (3).This interpretation will prove useful when we discuss the effects of thermally excited quasiparticles in Section IV B.

IV. ROBUSTNESS TO OTHER EFFECTS
In our calculation, we have made certain simplifications that allowed us to directly compute the pump-probe response coefficient.Here we consider processes and effects that were neglected above, and demonstrate that the qualitative form of the ratio χ PP /χ (1) remains universal in the long-time limit.Specifically, we will discuss the effects of short-ranged interactions (Sec.IV A), finite temperature (Sec.IV B), and nontrivial braiding statistics between within the multiplets that are created by each pulse (IV D).We also describe how our analysis can be generalised to systems with non-Abelian anyons in Section IV E. We will find that a number of timescales emerge from our analysis, which we summarise in Table I.

A. Short-ranged interactions
So far we have assumed that the only interactions between quasiparticles are through statistical braiding phases.However, if non-statistical interactions are present, as is the case generically, then the population of quasiparticles created by the pump pulse can influence the two-point correlator measured by the probe pulse through these interactions, and hence (v * σ) −1 e ∆/T Section IV B TABLE I. Summary of timescales that are relevant to pump-probe spectroscopy, when non-perturbative effects, short-ranged interactions, and finite temperatures are included.Here, T is the temperature, κ is the strength of the pump pulse [Eq.( 2)], v * is the maximum group velocity of quasiparticles, σ is the scattering cross section (having dimensions of length in 2D), and ∆ is the gap to excitations.The universal form of the finite-temperature linear response function, Eq. ( 33), can be seen when τ scat,th ≫ τ th ≫ τtr, which occurs at sufficiently low temperatures.Because thermal effects influence the pump-probe and linear response coefficients in the same way, their ratio remains unchanged; thus to see the universal form (30), we require only τscat,p ≫ τtr [see Eq. ( 31)], which occurs at sufficiently weak pump magnitude.
modify the response function (3).We argue that when interactions are sufficiently short-ranged, any such effect will be subleading compared to the contribution that we have identified above.
An intuitive way to see this is to employ the path integral perspective that we have used in the previous sections.The effects of short-ranged interactions are only felt by trajectories where a pump anyon comes within some characteristic radius r int of one of the probe anyons, and scatters off it.As before, we can integrate over the initial position of the pump anyons ⃗ x i , keeping all paths otherwise the same [this integral was responsible for the area functional A c in Eq. ( 15)].The range of ⃗ x i that result in paths where particles come within a distance r int of one another will scale with the perimeter of the probe anyon trajectories, rather than the area of the loop formed by them.The perimeter scales as t 2 (see Footnote [26]), which grows less quickly than the area ∼ t 3/2 2 ; hence interactions will only modify the subleading contributions to the response coefficient, represented by the term o(t 2 ) in Eq. ( 4).The above argument provides a relatively straightforward justification of why the late-time form of the perturbative response function χ (3) should not be altered by short-ranged interactions, but it is also useful to consider a more quantitative approach that does not rely on a perturbative expansion of χ PP .This is particularly important in light of the results of Section III, where we saw that non-perturbative effects can become important at late times.Looking at the ideal form Eq. ( 30), derived without non-statistical interactions, we see that the universal relationship will remain unchanged if the effects of local scattering between anyons occur on a timescale much longer than τ non−pert := (c PP κ 2 ) −2/3 .This scattering timescale is defined by the point at which the probability of a scattering event between a pump and probe anyon is order unity.This can be calculated in terms of a scattering crosssection σ, which in 2D is a length scale: Using standard scat-tering theory, we have τ scat,p = σv * λ pu , where λ pu is the density of anyons created by the pump pulse (which scales as κ −2 ), and v * is a typical velocity of the pump quasiparticles.Naturally, scattering between pump and probe anyons suppresses the two-time correlation function, and so we expect that the ratio χ PP /χ (1) will follow the form This modification to Eq. ( 30) makes no observable difference if τ scat,p ≫ τ non−pert , i.e. the statistical interactions alone fully compromise the phase coherence of the probe anyons before scattering processes have had time to take any effect.
In fact, as long as τ scat,p ≫ τ tr , then there will be an appropriate window of time in which the universal behaviour (30) can be seen: after transient effects have washed out, but before scattering effects have become appreciable.Note that this is always the case, independently of the system, if the pump pulse is weak enough, viz.κ is small enough.Alternatively, having weak interactions or small correlation lengths helps to satisfy this condition for larger values of κ.
To understand exactly what kinds of interactions count as sufficiently short-ranged, we can revisit the calculation that we described in Section II.Rotationally symmetric power-law interactions between pump and probe anyons can be included directly into the boosted Hamiltonian ( 19), and we suppose that at long distances these will decay as V (⃗ r j − ⃗ r k ) → V 0 |⃗ r j − ⃗ r k | −γ for some exponent γ.In this case, the angular part of the eigenstates (21) will remain unchanged, but the part of the Hamiltonian describing radial motion is now If γ > 2, then by applying dimensional analysis to the above differential operator, we can identify a crossover radius r int ∼ (V 0 /2m) 1/(γ−2) outside of which eigenstates are only weakly modified by the power-law interactions.(This length is not to be confused with the cross section σ, which would have to be computed via alternative means, e.g. through solving the appropriate Lippmann-Schwinger equation [36].)This radius is small for weak interactions, whereas the divergent contributions to the response coefficient are due to processes occurring at large distances r ≳ vt 2 .Hence, for γ > 2 these interactions will not qualitatively affect the late-time behaviour of the response function.
We do not directly address longer-ranged interactions γ ≤ 2 here, since in this case the assumption that quasiparticles separated by large distances propagate independently is not necessarily true.Indeed, there is no small length scale that can separate the regimes of small and large separation of quasiparticles, and so the key assumptions that went into our argument would be invalidated.It would be interesting to investigate such scenarios in future work, in particular in the context of the fractional quantum Hall effect, where anyons interact via long-ranged Coulomb forces γ = 1.

B. Finite temperature
Another assumption that has been made so far is that the system is in its ground state before the pump pulse arrives.In practice, with the system at finite temperature T , a population of thermally excited quasiparticles will be present, which themselves can affect the response of the system to external fields.In the regime T ≪ ∆, which we will focus on, the density of this population will be exponentially small ∼ e −∆/T , and so we can safely model the thermal excitations as a dilute gas of weakly interacting quasiparticles.
Firstly, let us neglect non-statistical interactions and, as a warm-up, consider the linear response coefficient, i.e. the twotime correlator χ lin (t) = ⟨ Â2 (t) Â1 (0)⟩.Focussing on the toric code for concreteness, as before we choose Â1,2 such that a pair of magnetic anyons are create at time 0 and annihilated at time t.The effect of thermal quasiparticles on χ lin (t) can then be understood using a picture analogous to that presented in Section III: For each trajectory of the magnetic anyons ⃗ r 1,2 (t), we can define a probability distribution for how many electric anyons pass through the loop (since e and m and mutual semions, we do not need to distinguish different linking orientations).The difference here is that the electric anyons are thermally activated, instead of being created out of the vacuum by the pump pulse, as before.
Thanks to the diluteness of the quasiparticle gas (the density λ th = d 2 k/(2π) 2 e −ϵ(k)/T scales as e −∆/T ≪ 1), the dynamics of the thermal electric anyons can be safely treated semiclassically [37].Accordingly, we describe the trajectories of the quasiparticles as straight lines with velocities independently distributed with probability density P (⃗ v), determined by the Boltzmann distribution.Since the electric anyons propagate independently, we can use the same logic as in Section III to argue that the probability of having p electric anyons linking with the loop formed by ⃗ r 1,2 (t) follows a Poisson distribution, and in this case the rate is given by s = λ th d 2 ⃗ v P (⃗ v)A[⃗ r 1,2 (t) − ⃗ vt], with A the area functional in (15), arising due to the integration over all initial positions of the electric anyons.Thus, each trajectory in the path integral over magnetic anyon trajectories ⃗ r 1,2 (t) should be weighted by a factor e −s p (−1) p s p /p! = e −2s , where s depends on ⃗ r 1,2 (t) through the area functional.
As usual, for typical trajectories the area functional scales as t 3/2 at late times, while the density follows an Arrhenius law λ th ∼ e −∆/T .Hence, comparing the finite-and zerotemperature response coefficients, we expect to find for some constant c.This allows us to define a new timescale τ th = (ce −∆/T ) −2/3 that describes how quickly the braiding phases between thermal and probe anyons degrades twopoint functions; see Table I.The prefactor may depend on how many anyons are created at a time by the probe pulse, among other factors, but will generally decay algebraically (as t −1 for the simple N ′ = 2 case considered in previous sections).
We see that at small finite temperatures, two-time correlation functions will decay via a characteristic 'squished exponen-tial' form e −(t/τ th ) 3/2 .Although this unusual form of broadening could in principle serve as a witness of nontrivial braiding even at linear response level, it is likely to be challenging to disentangle from other types of broadening, and as we will see there are constraints on the range of temperatures in which this decay mechanism will be the dominant one.This is why we propose measuring the pump-probe signal, where surplus anyons can be created in a controlled fashion using the pump pulse, and any background signals can be subtracted away according to Eq. ( 3).
With the above understood, we can determine the late-time behaviour of the pump-probe response coefficient χ PP (t 1 , t 2 ) at finite temperature by accounting for both thermal and pump-induced quasiparticles.The perturbed two-point function [the first term in Eq. ( 3)] is suppressed due to dephasing from both sources of quasiparticles, whereas the unperturbed correlator that is subtracted off has the same form as (33).The result is Comparing ( 33) and ( 34), we see that the universal form of the ratio (30) [which encompasses the perturbative result ( 4)] continues to hold at finite temperatures, since the linear and pump-probe response coefficients are modulated by the same decaying function.Of course, given the finite sensitivity of detectors in experiments, one wishes to work in a regime where τ th is large enough such that the individual signals χ PP and χ (1) do not become smaller than the experimental resolution before transient effects have worn off.Provided that temperatures can be lowered below ∆, this should be achievable thanks to the exponential dependence of τ th on 1/T .It is interesting to note parallels between these semiclassical arguments and an analogous derivation of the finitetemperature relaxational dynamics of the one-dimensional Ising chain in a transverse field, as studied in Ref. [37].In that context, quasiparticles are domain walls of separating domains of opposite magnetization, and so the two-time spin correlator C(t) = ⟨ Ẑj (t) Ẑj (0)⟩ ( Ẑj is a Pauli spin operator on some site j) acquires a phase of −1 each time a thermal excitation moves across site j.In the dilute-gas regime, when T is much less than the gap to excitations, C(t) is approximately equal its zero zero-temperature value multiplied by a decaying envelope ∼ e −t/τ that accounts for this dephasing due to thermal quasiparticles, which propagate with effectively random trajectories that are governed by the Boltzmann distribution.This multiplicative dephasing factor also arises in our results (33,34), with the difference that the mechanism of dephasing is non-local statistical interactions, rather than local scattering phases.This nonlocal mechanism gives rise to an envelope has with a different universal form: exp(−(t/τ th ) 3/2 ) instead of an ordinary exponential decay, for some timescale τ th ∝ e 2∆/3T .One additional effect that has not yet been accounted for is scattering between the probe anyons and the gas of thermal quasiparticles due to short-range non-statistical interactions.As we saw in the previous section, these scattering processes can lead to a further degradation of the phase coherence of the probe anyons, resulting in additional suppression of the twotime correlators.Assuming that the non-statistical interactions are short-ranged (decaying faster than an inverse square law, as in Section IV A), this will result in an ordinary exponential decay e −t2/τ scat,th , where in analogy to τ scat,th , the characteristic time is given by τ scat,th = (v * σλ th ) −1 , where again σ is the scattering cross-section, and λ th is the density of thermal quasiparticles.In the dilute gas regime, (low enough temperature and small enough κ), this envelope should affect the linear and pump-probe response coefficients equally, and hence the ratio χ PP /χ (1) should remain unchanged.Moreover, since the ratio τ scat,th /τ th grows as T is decreased, at sufficiently low temperatures we will have τ scat,th ≫ τ th , and hence the squished exponential form (33) will also be unaffected.
In addition, the combination of non-statistical interactions and finite temperatures provides a mechanism for the pump anyons to relax towards equilibrium, and this leads to a slow decay of the pump-probe signal with t 1 .The timescale for this to occur is again very slow due to the diluteness of the thermal excitations, on the order of τ scat,th , and hence it should be possible to find a suitable time delay t 1 that is large enough to see the asymptotic form of the response coefficient, but smaller than this thermalization timescale.
Finally, we remark on the possibility that the quasiparticles themselves may not be stable even at zero temperature, which occurs if the system in question is not actually in a topological phase, but only proximate to one, e.g. when anyons are weakly confined.In this case, the response coefficient will be altered nontrivially for times (t 1 + t 2 ) that exceed some cutoff, which is set by either the finite lifetime of quasiparticle excitations (which now remains finite even as T → 0), or the confinement lengthscale, whichever is reached first.(Note that this affects both the t 1 and t 2 dependence of χ PP , since the motion of pump anyons is also affected by such effects.)This cutoff diverges as one approaches the transition into the topological phase, and so if the system is proximate enough to a QSL, it will still be possible to observe the universal form described above.

C. Scattering from impurities
Realistic samples inevitably feature some amount of disorder.This can have two main effects for the dynamics of anyons: (a) impurities or defects can lead to elastic scattering of anyons, and, in certain cases, (b) disorder can generate and trap topological defects (see e.g.Ref. 38), which have nontrivial braiding properties with the dynamical anyons.In this subsection we discuss the consequences of these two impurity effects on the relaxation of linear and pump-probe response function.
a. Scattering effects.-Althoughimpurities in the sample are static, rather than mobile and dynamic, we can understand the effect of disorder at an approximate level in much the same way as scattering off thermally generated anyons: The impurities realise a short-ranged potential which is felt by the quasiparticles, and can scatter their momenta elastically.We can define an impurity scattering time τ imp = (v * σλ imp ) −1 , with λ imp the density of impurities and σ the impurity scattering cross-section.This gives us a typical time scale after which the momentum of a quasiparticle will be appreciably scattered.
Scattering of the probe anyons off impurities will degrade the amplitude for creation and re-annihilation, which will lead to a decay of the pump-probe signal.However, this effect is exactly reproduced in the linear response signal, and hence the ratio χ PP /χ (1) will remain unaffected.However, scattering of the pump anyons between times t 1 and t 1 + t 2 may modify the pump-probe signal in a way that is not counterbalanced by χ (1) .While a detailed calculation of the pumpprobe response coefficient in the presence of quenched disorder is beyond the scope of this work, we anticipate that these scattering events will make braiding between pump and probe anyons marginally less likely, since the straight-line trajectories shown in Fig. 1 will have to be modified.The universal signal we describe here will still be observable provided that the timescale τ imp is longer than the timescale for t 2 beyond which transient effects have subsided and the relation (4) becomes valid.Indeed, converting τ imp to a corresponding mean free path ℓ imp , we expect such a window of time to exist provided that disorder is not so strong such that ℓ imp ∼ a where a is the lattice spacing.This is certainly true in any 'weakdisorder' regime.
b. Braiding effects.-Theconsequences on χ (1) (t) of defects with nontrivial braiding can be understood along the lines of the argument provided in Subsec.IV B for thermal, i.e. dynamical, anyonic quasiparticles.However, since these topological defects are static, the average number of defects that braids with the anyon pair grows like ( t/m) 2 -to be compared with the vt × t/m when the thermal excitation have average velocity v -since √ t-spreading of the one-particle propagator is now the only contribution to braiding.Consequently, these effects produce a further exponential relaxation of χ (1) (t) scaling like exp(−t/τ an.imp. ) on top of the faster-than-exponential thermal suppression in Eq. (33).Therefore this extra contribution will be subleading for small concentration of impurities.
Instead, regarding the ratio χ PP /χ (1) , the braiding does not affect the pump anyons: the leading semiclassical contribution is obtained when their trajectories are the same in the forward and backward time evolution, so they cannot braid with the defects.Therefore, following the lines of the arguments in paragraph (a), we see that braiding effects do not impact the ratio χ PP /χ (1) .

D. Statistical interactions within multiplets
So far, we have considered response functions for perturbations that create multiplets of excitations within which all particles are mutually bosonic.An example that we regularly referred back to was the creation of a pair of electric anyons in the toric code, which have no non-trivial braiding or exchange statistics as a pair, despite being semionic with respect to magnetic excitations.Here we consider what happens if the multiplets created by the pump and/or probe pulses contain excitations that are not bosonic with respect to one another.One example of such a multiplet-again in the context of the toric code-is a pair of electric-magnetic (em) composite particles, which are fermionic with respect to one another.
As previously mentioned, an important consequence of non-bosonic statistics within a multiplet is that the constituent excitations cannot exist at the same point in space-a generalization of Pauli's exclusion principle.Thus, we cannot use wavefunctions of the form (7) as a sensible low-energy description of the state of the system immediately after the pulse.Since the wavefunction must vanish at points where particles coincide, one must invoke a regulator that specifies the limiting behaviour of |Ψ N ⟩ at small separations.
The effect of this generalized exclusion principle can already be seen in linear response functions, as was shown in Ref. [18].In brief, the authors of that work calculated the dynamical structure factor (the Fourier transform of a twotime correlator ⟨VAC| Â2 (t) Â2 (0)|VAC⟩) using a low-energy effective theory describing the dynamics of a pair of anyons between times 0 and t.Motivated by lattice models such as the toric code, the regularization of the post-pulse state Â2 |VAC⟩ that they chose was a rotationally symmetric wavefunction where the two anyons are separated by a finite exclusion radius a, i.e.
is the two-anyon state with centre of mass ⃗ R, and relative displacement (a, ϕ) in polar coordinates.Converting their frequency-space results into real time, the late-time behaviour of the correlator follows a power law t −1−α , where α is the statistical parameter as before.The same time-dependence can be shown to arise for any uniform state |Ψ N ⟩ where the initial distance between the two anyons does not exceed some fixed microscopic lengthscale a [39].With the exception of α = 0 (bosons), this clearly differs from the t −1 linear response behaviour that we argued for in Section I A, which is simply the amplitude for two free particles to recombine [the first factor in Eq. ( 5)].
While it is clear that individual response functions-linear or otherwise-will be modified by statistical interactions between multiplets, the central quantity in our work is the ratio of the pump-probe and linear response coefficients, which as we argue will continue to follow the universal form derived before Eqs.(4,30).Firstly, the effect of statistical interactions within the pump multiplet will only give rise to a quantitative modification of the distribution of quasiparticle velocities created by the pump pulse: once these quasiparticles are created, they will still propagate ballistically.This only leaves interactions within the probe multiplet.Even with these included, we can still use the path integral representation of the dynamics of probe anyons, described in Sections I A, II, which tells us that each trajectory of the probe anyons should be weighted by a factor of the area functional A c , equal to the size of the space of initial pump coordinates ⃗ x i that lead to non-trivial braiding [Eq.(15)].Crucially even when probe anyons are not mutually bosonic, as was the case considered before, for typical paths this area functional continues to follow the same latetime asymptotic form A c ∝ t 3/2 2 .Accordingly, we still ex-pect Eqs.(4,30) to hold, even though the individual response functions χ (1) , χ PP are modified.On the basis of the results of Ref. [18], in the case where two probe anyons are created at a time N ′ = 2, we expect to see the perturbative pump-probe response coefficient scaling as χ , where the braiding phase between probe anyons is given by 2πα pr .
The scaling of A c can be argued for solely using the dimension-counting arguments given at the end of Section II D, where the late-time limit is equated to the limit where the velocity of the pump anyon is taken to be large.At large velocities the area must scale linearly with v, and since the only velocity-independent length scale in the problem is t 2 /m, and the only dimensionless parameter is v √ mt 2 , this fixes We present a more concrete calculation that confirms this scaling of A c in Appendix C.
At the end of this section, we wish to highlight a difference between the results of Ref. [18], where the effects of particle statistics on linear response coefficients is studied, versus the effect we study in this paper, which shows up only beyond linear response.The former will be seen in systems that possesses fermionic excitations, which have nontrivial exchange statistics, but trivial braiding statistics.In contrast, the universal late time behaviour of the pump-probe response coefficient is a reflection of nontrivial braiding statistics: the phase e iΛ is determined by the linking of paths in spacetime, rather than an exchange of identical particles.Because of this, pump-probe spectroscopy serves as an identifier of topological excitations with braiding statistics, rather than just nontrivial exchange statistics, which arise in non-topological fermionic systems.

E. Non-Abelian statistics
Until now, we had only made explicit reference to systems with Abelian anyons, where the effect of braiding is to induce a complex phase in the wavefunction.However, our analysis also applies to topological phases whose excitations possess non-Abelian statistics.In such systems, excited states exhibit a topological degeneracy, meaning that an extra discrete quantum degree of freedom is required to fully specify the state of the system, in addition to the positions of the anyons [40].Braiding of excitations results in the application of a unitary rotation acting on this degenerate space.
These non-Abelian statistical interactions can be incorporated into a path integral language in a similar way to before.In place of the phase e iΛ in Eq. ( 11), we should instead substitute a matrix element of the unitary operator associated with the braid carried out by the trajectories ⃗ x + j (t), ⃗ r k (t).Specifically, we make the replacement (e iΛ − 1) where U is a functional of the trajectories, depending only on their braiding properties, and |χ i,f ⟩ are discrete wavefunctions in the discrete space, which are set by the specifics of the operators Â1,2 to which the probe pulse couples (see, e.g.Ref. [40], Sec.III C).
With the exception of this difference, all our arguments can be applied in exactly the same way as before.In particular, the decomposition of the path integral into topologically distinct contributions [Eq.(15)] still applies, just with the non-Abelian matrix element in place of the complex phase.The functionals A c depend only on the geometry of the trajectories, and the free part of the action is as before.Thus, the late-time form of the response coefficient should continue to obey the relationship (4).

V. APPLICATION TO PERTURBED TORIC CODE
In this section, we study a microscopic Hamiltonian that possesses anyonic excitations, which allows us to apply our general results to a more concrete setup.We are also able to relate the phenomenological parameters used in Section II (mass m, length scale ξ, etc.) to properties of the Hamiltonian.
The specific microscopic model that we consider is the toric code perturbed by a magnetic field.In the toric code, qubits are located at the edges j of a square lattice, which we describe using Pauli operators X j , Y j , Z j .The unperturbed Hamiltonian is a sum of four-body terms located at the vertices v and plaquettes p of the lattice [27,28] where the star operators Âv = j∈v Xe act on all edges around the vertex v, and the plaquette operators Bp = j∈p Ẑp act on all edges around a plaquette p.The ground state of Ĥ0 is the wavefunction stabilized by all star and plaquette operators, Starting from the ground state and acting with Ẑe on some edge creates a pair of excitations each of energy J Aone for each of the star operators Âv that act nontrivially on e and hence anticommute with Ẑj .Similarly, acting with Xj creates a pair of excitations on the two plaquettes shared by j, each with energy J B .These two types of excitation are referred to as electric (e) and magnetic (m) anyons respectively.The fact that they are semions with respect to one another can be seen by acting successively with operators Ẑj in a way that moves the electric particle around a path that encircles a magnetic particle (see Ref. [28] for details).Since these excited states are exact eigenstates of Ĥ0 , the anyons do not move once created in the absence of any external perturbation.The immobility of the excitations is reflected in the lack of dispersion in the spectrum of Ĥ0 : eigenstates come in highly degenerate multiplets with discrete energies n A J A + n B J B , where n A , n B are the number of electric and magnetic anyons, respectively.
To endow the anyonic excitations with a dispersion, we introduce a magnetic field, which for simplicity we place in the x-z plane.The full Hamiltonian that we consider in this section is We work in the limit J A,B ≫ h x,z .In this limit, we can neglect hybridization of eigenstates with different numbers of magnetic and electric anyons, and the main effect of the magnetic fields is to lift the degeneracy within each excitation number-sector.The x-magnetic field generates hopping of magnetic anyons in the dual lattice, and similarly the zmagnetic field allows electric anyons to hop in the original lattice.The dispersion of a single electric or magnetic anyon then becomes where (k x , k y ) is the quasimomentum, and a is the lattice spacing.We have written these dispersions relative to the band minima, which are at energies ∆ e = J A − 4h z and ∆ m = J B − 4h x (using the same notation for the threshold energies as in Section II A).
In a pump-probe experiment, the incoming pulses of light will naturally couple to the microscopic spins.We can choose the polarization of the incoming fields such that the pump pulse couples to the X-component of the spins, and the probe pulse couples to the Z-component.This way, assuming that the wavelength of the radiation is long compared to the sample size, the time-dependent fields experienced by the system are uniform in space We take the pump pulse to be a Gaussian wavepacket arriving at t = 0, centred around a frequency 2∆ e + δ 1 , where δ 1 is a detuning much smaller than J B , with a width of frequencies Due to its frequency profile, the pump pulse can only excite a pair of electric anyons, assuming J A and J B are separated by a gap larger than τ 1 .We can therefore write down the wavefunction of the system at times 0 where |j, j ′ ⟩ e,e is an excited state with electric anyons at lattice sites j, j ′ , and Ĥe,e is the Hamiltonian in the relevant excitation number-sector.Thanks to the lack of statistical interactions between these two particles, we can compute the time evolution by transforming to plane wave states | ⃗ k n ⟩ e = M −1 j e i ⃗ k•⃗ rj |j⟩ e and using the single-particle dispersion (37).Here, ⃗ r j is the real space coordinate for site j, M is the number of sites in the lattice, and the discrete set of wavevectors satisfying periodic boundary conditions are ⃗ k n = (2πn x /L, 2πn y /L), with n x , n y ∈ {−M/2 + 1, . . ., M/2}.
We then have where the two particle state | ⃗ k n , − ⃗ k n ⟩ e,e is the wavefunction of a pair of electric anyons in plane wave states with opposing quasimomenta ⃗ k n and − ⃗ k n , and we have defined f ( ⃗ k) := cos(k x a) + cos(k y a).(In performing the time evolution, we have neglected the effective hard-core constraint that two electric anyons cannot reside on the same vertex; however the effect of this is negligible in the regime of interest, as we will see.)The upper limit of the integral over t ′ can be extended to +∞ for times t ≫ τ 1 , which gives From this wavefunction we can read off the distribution of quasimomenta of the electric anyons created by the pump pulse.While various different hierarchies of energy scales can in principle be considered, for convenience we will work in a regime where δ 1 ≪ τ −1 1 ≪ h z , in which case this distribution is peaked near the bottom of the band, allowing us to expand (37) to quadratic order in k x,y .We can therefore consider quadratically dispersing electric anyons with isotropic mass The distribution of quasimomenta is then approximately proportional to e −ξ 4 e k 4 /2 , where the length scale is By transforming back to real space, we find that the wavefunction describes pairs of electric anyons in wavepackets of size ξ e centred around the same point.This provides a proper UV regularization of the wavefunction (7) that we employed previously.We observe that ξ e can be identified as the typical propagation length of the electric anyons over the time window τ 1 during which they are created.Note also that in the regime τ −1 1 ≪ h z , we know that ξ e is much greater than the lattice spacing, which allows us to approximate f ( ⃗ k) ≈ 2. This also justifies our choice to neglect the hard-core constraint on electric anyons in (41), since components of the wavefunction where two anyons are located at the same vertex are small.
The wavefunction (42) can be used in place of |Ψ N ⟩ in the operator ζ defined in Eq. (10).As in Section II, we will employ an approximation where we ignore the influence of the magnetic anyons generated by the probe pulse on the trajectories of the original electric anyons.Because of this, when the trace in ( 9) is taken, only contributions where the wavevector on the ket and bra parts of ζ coincide will survive.We then have + (terms annihilated by trace).
The probe pulse allows us to measure the two-time correlator appearing in Eq. (3) (see Section VI for details on how this is achieved).In our case, this pulse is polarized along the x axis, which means that the operators Â1,2 are simply j Xj .To isolate contributions coming from processes involving two magnetic anyons, the incoming waveform can be frequencymatched to the magnetic anyon pair threshold of 2∆ m , i.e. the pulse only contains frequency components near this energy.Because of this, we can again expand the magnetic anyon dispersion to quadratic order about the band minimum, and we identify the mass m m = (2h x a 2 ) −1 .While this assumption is useful for calculations, we expect to see the same qualitative results even if the range of frequencies is broader.
For each term in the sum in (45), we must compute the twotime correlator ⟨ Â2 (t 1 + t 2 ) Â1 (t 1 )⟩ of these magnetic anyons in the presence of electric anyons that propagate at the group velocity The frequency profile of the probe pulse ensures that Â1 excites a magnetic anyon pair which is de-excited by Â2 .The amplitude for this is precisely the propagator I(⃗ v n , t 2 ) that we computed in Section II E. Putting everything together, and using the normal- , the long-time limit of the perturbative pump-probe response coefficient becomes Using the expression (26), and restoring the original microscopic quantities using (43,44), we get where the factor of a −2 arises due to the normalization of χ by the volume L 2 , rather than the number of sites M .This calculation demonstrates how the universal t 1/2 2 divergence emerges starting within a specific microscopic model.
To derive this result, we have made certain assumptions about hierarchies of energy scales, namely that the fields h x,z should be weak enough such that hybridization between different anyon sectors is negligible, and that the pulse frequencies are close enough to threshold δ 1 ≪ h z .While deviations from these assumptions may affect the scaling of the prefactor, in general we expect the dependence on t 2 to be a universal feature of systems whose excitations possess non-trivial braiding statistics.

VI. EXPERIMENTAL CONSIDERATIONS
Having studied the behaviour of the pump-probe response coefficient in detail, we now provide a general discussion of the ingredients necessary to measure this quantity in experiment.For most of this section, our focus will be on putative solid-state realizations of quantum spin liquids, for which bulk probes are particularly useful.We comment on other settings later on, namely quantum Hall systems, ultracold atoms and Rydberg atom arrays.
The dynamics of spins in solid state systems typically occur on timescales of order ∼ 1 ps.As an example, in the candidate material α-RuCl 3 , for which there is evidence of a field-induced non-Abelian QSL phase [41][42][43][44], the magnetic couplings are estimated to be in the range 70-90 K [45], corresponding to a frequency of ∼ 1.5 THz.Recent technical advances have facilitated the generation of high-intensity THzdomain pulses with short time resolution [46,47], which have already been applied to study ultrafast magnetization dynamics in systems with spontaneous macroscopic spin ordering [48][49][50][51].Here, in analogy with standard pump-probe setups familiar from other kinds of nonlinear spectroscopy [32], we will describe a sequence of pulses which allows one to measure the particular response coefficient χ PP (t 1 , t 2 ) [Eq. ( 3)] in a candidate quantum spin liquid.In fact, this particular sequence has already been used in previous experiments, where the aim was to demonstrate coherent control of spin precessional motion [48].Thus, the effect we describe in this paper should be detectable using currently existing experimental techniques.
To be specific, we propose to first illuminate the sample with a short intense pump pulse whose frequency range overlaps with the creation threshold energy for a given quasiparticle multiplet (a pair of electric anyons, say).Since the wavelength of THz light is large, the incoming radiation couples directly to the total magnetization Mα , where the component α is set by the polarization of the magnetic field of the radiation (i.e.we neglect the momentum of the photons).After waiting for a time t 1 , a second weaker pulse is applied, which for now we model as infinitely short-lived, giving a magnetic field B pr (t) = B 0 δ(t − t 1 ) along a different direction β.This perturbation modifies the state of the electron spins at later times, and the resulting time-dependent magnetization Mγ (t) in turn leads to emission of radiation due to free induction decay (FID).The amplitude of the emitted FID radiation can be measured along a chosen polarization γ in a time-resolved fashion using e.g.electro-optic sampling [52], which allows one to infer the time-dependent magnetization ⟨ M (t 1 + t 2 )⟩.
We have already discussed the effect of the pump pulse in Sections I A and V: the state of the system immediately after the pulse can be described using the right hand side of (2), where Â0 includes components of the magnetization operator M that oscillate at frequencies within the frequency range of the pulse.As for the probe pulse, since this is weak and infinitesimally short-lived, we can expand to lowest order in B 0 .If ρ pert is the post-pump state, then immediately after the probe pulse the system is in the state ρpert −iB 0 [ M (t 1 ), ρpert ]+O(B 2 0 ) (we continue to work in the Heisenberg picture, where M = e i Ĥt M e −i Ĥt ).Then, the expectation value of the magnetization at time (t 1 + t 2 ) is given by pert (48) Therefore, by extracting the linear-in-B 0 part of the magnetization, and subtracting the same quantity without the probe pulse, we obtain the imaginary part of the desired response coefficient (3).
The above procedure is conceptually straightforward and achievable using currently available techniques.Nevertheless, it is also worth contemplating alternative setups that measure the response coefficient directly in the frequency domain, more akin to standard spectroscopic measurements.Rather than using electro-optic sampling to detect the emitted field, one can alternatively perform an absorption measurement, with the detector downstream of the probe pulse, such that the field being detected is a superposition of the probe pulse field and the FID signal E pr (t) + E FID (t).Using a spectrometer, the power spectrum I(ω) = |E pr (ω) + E FID (ω)| 2 can be obtained, and since the signal field is weak the signal will be found in the cross-term 2ℜ[E pr (ω) * E FID (ω)], since the quadratic term |E FID (ω)| 2 can be neglected.This measurement scheme constitutes an intrinsic heterodyne detection of the FID field, with the probe pulse serving as a local oscillator (see Ref. [32]).Due to the short probe pulse, E (ω) is approximately constant in ω, and so this gives us access to the one-sided Fourier transform of the imaginary part of the response coefficient χPP (t where t 1 is set by the time interval between the pump and probe pulses.Since the emitted FID field is π/2 out of phase with the magnetization [32], such an experiment would give us direct access to the imaginary part ℑ χPP (t 1 , ω), and the real part could be reconstructed using the Kramers-Kronig relations.
Given that there may be scenarios where the measured data is in the frequency domain, let us consider how the universal relationship between linear and pump-probe response coefficients manifests itself in Fourier space.Since Eq. ( 4) is valid in the limit of late times, we expect that the relationship will be most stark at frequencies that are close to the non-analytic points of χ(1) (ω).In particular, recall from Section II B that the imaginary part of χ(1) (ω)-which is proportional to the spectral function of the magnetization operator-exhibits non-analytic behaviour at the threshold frequency ∆ N ′ , the minimum energy required to create excitations above the quasiparticle vacuum.The nature of the edge singularity in χ(1) (ω) will determine the form of non-analytic behaviour seen in χ(3) PP (t 1 , ω) via the relationship (4).The simplest case, which applies to all the cases that we have studied in this work, is a power-law singularity, where the linear response coefficient in the time domain follows χ (1) (t) ∼ it −η e −i∆ N ′ t , TABLE II.Summary of the relationships between the linear and pump-probe response coefficients in the time-and frequencydomain; see Eqs. (4,51).Fourier transforms at frequency ω are taken with respect to the time t in linear response, and t2 in pump-probe response.We define δω := ω − ∆ N ′ as the frequency relative to the energy threshold for creation of excitations, and the function f (y) is the Fourier transform of exp(−|x| 3/2 ) with respect to x, and * denotes a convolution.These results are valid in the limit of long times t, t2, or sufficiently close to threshold, i.e. small δω, as appropriate.The exponent η depends on the particulars of how anyons are created and annihilated (see Section IV D), but the ratio of the linear and pump-probe response coefficients in both the frequency-and time-domain is universal.
where the exponent η depends on the number of anyons that can be created at a time by the probe pulse, and the statistical phases between them.In frequency space, this gives us where the above is expected to hold for |ω| sufficiently close to ∆ N ′ .Our time-domain results can be employed to determine the late-time form of the pump-probe response coefficient, which upon Fourier transforming gives We see a more drastic singularity in the pump-probe response coefficient by virtue of the fact that the ratio χ PP (t 1 , t 2 )/χ (1) (t 2 ) grows with t 2 .If χ (1) (ω) exhibits more complicated non-analytic behaviour (i.e.different from a power law), then one can instead use the convolution theorem to determine the corresponding form for χ(3) To properly capture the short-time behaviour of χ PP , before the universal signal (4) is dominant, the factor of |ω ′ | −5/2 should in principle be altered for values of ω ′ much larger than τ −1 tr .However this will not impact the qualitative form of χ(3) PP (t 1 , ω) near threshold.As discussed in Sections III and IV, at very long times the response functions may be modulated by a decaying envelope due to either non-perturbative effects or suppression due to scattering and/or finite temperatures.In the frequency domain, this results in a 'smoothing out' of any non-analytic behaviour over frequency scales on the order τ −1 , where τ is the appropriate timescale (see Table I).For example, at finite temperatures the linear and pump-probe response coefficients are modulated by a factor f (t/τ th ), where we define f (x) := exp(−|x| 3/2 ).Thus, χ(1) (ω) will be the convolution of the zero-temperature response function with f (ωτ th ), where f (y) is the Fourier transform of f (x), being a smooth function of y that peaks at y = 0 and has a width of order unity.In practice, given that such timescales are typically very large (at least for low temperatures and weak pump pulses), it is likely that the measurement apparatus will not be able to resolve these effects, and the formally divergent expressions given above can be used instead.These results are summarised in Table II.
We finally remark on some aspects of the generation of the pump anyons.Firstly, in pump-probe spectroscopy, the initial pulse is typically highly intense, with the aim to bring the system strongly out of equilibrium.In this case, assuming that the relevant matrix element for anyon generation [i.e. the coefficient of proportionality in Eq. ( 7)] is not small, then the density of pump anyons n pump will be fairly high.The pump anyon density can be related to the pulse strength factor κ discussed in Section III as n pump ∝ κ 2 , and hence the intensity of the pump pulse controls how long it takes for the nonperturbative regime to set in.We note that our analysis and prediction of the universal form (30) remains valid for finite pump densities, but begins to break down as n pump approaches 1/a 2 where a is the lattice spacing, i.e. one pump anyon per unit cell.This regime is unlikely to be reached in practice.
Another possibility is that the key physics described in this work might still be detectable even if the pump anyons are generated incoherently.Indeed, in Section IV B, we saw how thermally generated anyons can modify the linear response coefficient.In a sufficiently low temperature regime e −∆/T ≪ 1, such that the scattering time τ scat is much longer than the dephasing time τ th , we expect to see a characteristic linear response coefficient following Eq.(33).Thus, rather than using a pump pulse as a means of generating excess quasiparticles, an increase in temperature could be used.The temperature dependence of frequency-resolved THz absorption measurements at low temperatures could therefore also provide a signature of non-trivial braiding statistics.
We conclude this section by addressing other systems that may host topological phases with anyonic excitations.Twodimensional electron gases in the fractional quantum Hall regime can host Abelian and/or non-Abelian anyons [53]; however in practice performing spectroscopy on these systems may be challenging due to the presence of signals coming from other layers of the semiconductor heterostructure devices that are required to realise the electron gas.(In these systems, one can instead use novel device geometries to guide the motion of anyons through edge modes; this approach has recently been employed to detect braiding statistics [54,55].)Moreover, since anyons are charged and a magnetic field is present, our analysis would only remain valid up to a timescale set by the cyclotron frequency, see Footnote [24].Outside of the solid state, proposals have been put forward to realise topologically ordered phases in ultracold atomic gases [56][57][58][59][60], and more recently in arrays of Rydberg atoms in optical tweezers [61], which have since been implemented in Ref. [62].Light-based probes are natural in these settings, and thanks to the high levels of isolation from the environment and the lack of extraneous degrees of freedom, one expects to see clean spectroscopic signatures.Whether the signal we derive here can be seen in this context depends on the system sizes that can be reached, but with large enough samples, nonlinear spectroscopic probes could prove to be a useful probe of anyonic statistics, particularly in platforms where individual atoms cannot be addressed and measured in a spatially resolved way.

VII. CONCLUSION AND OUTLOOK
We have studied pump-probe spectroscopy of twodimensional systems that possess excitations with unconventional statistics.Our key result is a universal relationship dictating the late-time behaviour of the response coefficient.The origin of this behaviour can be intuitively understood using a path integral description for the dynamics of quasiparticles: The factor of t 3/2 2 in Eq. ( 4) arises when one calculates the probability that an anyon created by the pump pulse links with the trajectories of anyons created by the probe pulse, see Fig. 1.
After confirming this result through an explicit calculation of χ PP (t 1 , t 2 ), we considered the effects of non-statistical short-ranged interactions and finite temperatures, and argued that our result should remain valid even after these effects are included.Accordingly, the relationship between the linear and pump-probe response coefficients (4) serves as a robust fingerprint of anyonic statistics.While our rigorous calculations were performed using a low-energy effective theory for the dynamics of anyons, it is possible to make quantitative connections to specific microscopic models, as we demonstrated for the perturbed toric code.We finally discussed how the relevant signals can be measured using current THz-domain spectroscopic techniques.
Given that the experimental methods necessary to measure the relevant signal are already available, we anticipate that nonlinear spectroscopy could be used to obtain more information about the nature of magnetism in materials that are candidate quantum spin liquids.One of the most actively explored materials in this context is α-RuCl 3 , and neutron scattering and electron spin resonance data provide evidence that under certain applied magnetic fields this system is in or proximate to a QSL phase [44,[63][64][65][66].It would therefore be of great interest to investigate the behaviour of the pump-probe response coefficient in microscopic models that are thought to describe the spin dynamics in this material, as well as its close relatives [67].This would allow useful comparison with potential nonlinear spectroscopic experiments on this class of materials.Already, our results indicate that for such a proximate spin liquid, the pump-probe response coefficient should behave in the way discussed above, up to some characteristic timescale dictating the lifetime of quasiparticles, which should diverge close to the transition into a QSL.
In addition, our work suggests that universal relationships between linear and nonlinear response coefficients may arise in more general topologically ordered systems.For example, in three spatial dimensions excitations can be pointlike or looplike, and mutual statistics between particles and loops can be defined in analogy to the 2D case [68][69][70][71].Understanding how statistical phases between these excitations manifest themselves in nonlinear response will form an interesting direction for future work, which may prove to be useful in the search for topologically ordered materials in higher dimension e.g.Coulomb spin liquids [72].
⃗ x i,f , ⃗ r i,f in terms of the times t 1,2 themselves.We keep velocities rather than positions fixed because the upper limits of the velocity integrals will eventually be cut off by a nonuniversal UV scale v cutoff ∼ 1/ξm, where ξ is set by either the lattice constant or the size of the anyon wavepacket, and v cutoff should remain invariant under the scaling transformation.
Appendix B: Evaluation of Eq. (25) In this appendix we detail how the integral (25) is evaluated in the limit v → ∞ by means of a stationary phase approxi-mation.We start by transforming to dimensionless integration variables u i,f = M/2t 2 r i,f , and defining the dimensionless parameters As long as v ̸ = 0, the large-t limit of the above can be extracted by taking the limit β → ∞.If v = 0 then all t 2 -dependence drops out, and we obtain a contribution that is constant in t 2 .This contribution we ignore for now.From here on we assume v ̸ = 0, and take the long time limit via β → ∞; this is valid for t 2 ≫ 1/v 2 M .Additionally, since the statistical parameters α k are only defined modulo an integer, we can without loss of generality choose α k ∈ [0, 1).
Since β is large, the angular integrals can be evaluated using a stationary phase approximation, with stationary points at ϕ i,f = 0, π.This gives where The sums over σ i,f are for the different stationary points at ϕ i,f = 0, π, and we have split up the sums over ℓ into separate parts where ℓ k − α k is either positive or negative.Now we evaluate the sums using Ref. [73] Eq. 5. where µ can be any real value satisfying −1 < µ < 1 − α.
Now, we note that the integrand in (B2) is a fast-oscillating function of u f and u i , and hence will be dominated by contributions at large u f , u i ≳ β, where the Bessel functions oscillate equally quickly.Therefore, we can take the largez limit of (B3), which simplifies using the asymptotic form J µ (u) ≈ 2/πu cos(u − µπ/2 − π/4), valid for large real positive u: where terms that integrate quickly with u have been dropped.
We thus obtain (recalling that k α k is an integer and The above is now manifestly invariant under shifts of α k → α k + n, with n ∈ Z, as we would expect.Note that when becomes completely independent of α k , and so the difference of products in (B2) will vanish, leaving only the σ i = −σ f terms, whereupon we can set k A k cos(πα k ).We now have e iπσi/2 e −iσiβ(x f +xi) e −2iuiu f (B6) where the topological quantity Υ[{α k }] is given in Eq. ( 27) Now we make the transformation to variables X = u i + u f and x = u f − u i , giving Substituting the above into Eq.(B7), and restoring the original dimensionful quantities, we finally obtain Eq. (26).
along the direction perpendicular to ⃗ v.Here the angled brackets are a shorthand for the average over all paths ⃗ r j (t) that contribute to the two-time correlation function, i.e. ⟨C⟩ := ( Dr j (t)e iS C)/( Dr j (t)e iS ) for any functional C. (Note that this is not necessarily a real quantity, but we are interested in the typical magnitude of A c , and will therefore take an absolute value at the end.)The relation above was argued for in Subsec.II D only based on the ballistic trajectory of the pump anyon, and the argument carries over to the case where anyons within a multiple have non-trivial mutual statistics.
To compute the above, we will adopt the model of Ref. [18], where the local operator A(⃗ r) creates the two anyons at a microscopic distance a from one another.Specifically, |Ψ N ′ ⟩ = d 2 ⃗ R dϕ | ⃗ R, a, ϕ⟩ N ′ where | ⃗ R, a, ϕ⟩ N ′ is a two-particle state with centre-of-mass coordinate ⃗ R and relative separation (a, ϕ) in polar coordinates.Here a is a UV length scale, which is required to regularize the state of the quasiparticles created by the operators Â1,2 without violating the exclusion principle [previously given by Eq. ( 7)].By dimensional analysis, the final result will be proportional to t 3/2 2 multiplied by a function of the dimensionless ratio a m/t 2 .The behaviour of this function at small arguments will determine the latetime scaling behaviour of A c ; the following calculation will demonstrate that this function tends to a constant as a → 0, i.e. the lengthscale a falls out of the problem at late enough times.
By translation invariance, the centre-of-mass coordinate and the relative coordinate decouple, and the statistical phases depend only on the latter.In fact, the part of the Hamiltonian controlling the motion of the relative coordinate is precisely the same as the transformed Hamiltonian Ĥ′ k appearing in Section II E, which describes a single particle orbiting around a flux tube of strength 2πα at the origin.The eigenstates of this Hamiltonian are given in Eq. ( 21), and since the initial state |Ψ N ′ ⟩ is rotationally invariant we need only consider the zero angular momentum sector, ℓ = 0. We have To compute the necessary matrix elements, we require the expression for the eigenstates of Ĥ, Eq. ( 21), along with the standard integral given in Ref.Notice that the integral over τ is unchanged upon making the transformation τ → t 2 − τ ; we can thus change the upper limit to t 2 /2, and multiply the expression by 2. Defining the dimensionless parameter γ := ma 2 /t 2 , we now transform to new integration variables u = t 2 /τ , s = γ r/a, giving × e iγ −1 (γ 2 +s 2 )u 2 /(u−1) J α s u u − 1 J α (su) We are interested in the behaviour of this expression in the limit of small γ ≪ 1.In this limit, the integrand is a fastoscillating function of s, and so the integral will be dominated by contributions where s ≲ γ 1/2 .Since u lies in the interval [1,2], we can safely expand the second Bessel function for small arguments J α (su) ≈ (su/2) α /Γ(α + 1).(Note that the argument of the first Bessel function is large for u close to 1, and therefore should not be expanded.)The integral over s can be evaluated using another standard result [73]  ) is bounded as u → 1 for any a, the integral in the above converges to the constant π/16 as one takes the limit γ → 0. Thus, we can read the time dependence off as ⟨A c ⟩ χ (1) (t 2 ) ∝ t 1/2−α 2 . A simple calculation using Eq.(C4) gives the linear response coefficient as χ (1)