Divergent nonlinear response from quasiparticle interactions

We demonstrate that nonlinear response functions in many-body systems carry a sharp signature of interactions between gapped low-energy quasiparticles. Such interactions are challenging to deduce from linear response measurements. The signature takes the form of a divergent-in-time contribution to the response -- linear in time in the case when quasiparticles propagate ballistically -- that is absent for free bosonic excitations. We give an intuitive semiclassical picture of this singular behaviour, validated against exact results from a form-factor expansion of the Ising chain and tDMRG simulations in a non-integrable model -- the spin-1 AKLT chain. We comment on extensions of these results to more general settings, finite temperature, and higher dimensions.

The response of a quantum many-body system to external perturbations is central to experimentally extracting information about its properties.In typical settings, such as transport and scattering measurements, the system is weakly perturbed out of its equilibrium state, and the leading linear response [1,2] contribution can be related via the fluctuation-dissipation theorem to a twopoint equilibrium correlation function.Consequently, linear response functions are often relatively straightforward to interpret.For instance, at zero temperature (T = 0), as long as the external perturbation can create single quasiparticle (QP) excitations on top of the ground state, their dispersion can be read off directly from momentum-resolved linear-response data.However, the simplicity that lends power to linear response often limits its utility in more complex situations.For example, various distinct physical mechanisms can give rise to a broad frequency spectrum in linear response functions: e.g., selection rules requiring probe fields to excite multiple QPs, inhomogeneous broadening due to quenched disorder, and homogeneous broadening from QP decay.Differentiating among these using linear response data alone is a challenge.
Here, we offer such a qualitative perspective: a semiclassical theory of nonlinear response.We focus on the simplest non-trivial systems, whose low-energy spectrum consists of gapped QPs and where the perturbing field can excite single QPs.For ballistic QPs in one dimension (d = 1) and at T = 0 we find that the q = 0 third-order response functions diverge linearly in the time interval between distinct applications of the external field, with a strength set by the inter-QP scattering phase shifts, and a scaling function specific to the nonlinear protocol.This richness is to be contrasted with q = 0 linear response for such systems: a delta-function peak at the gap frequency, related to the stability of QPs, and is nonzero even for free bosonic QPs with no scattering.
Apart from enjoying the simplifying features of ballistic d = 1 QPs, the primary example we consider below -the transverse-field Ising chain -is also integrable.As is well known, it maps to a theory of free fermionic QPs, whose statistics enforce a scattering phase shift of −1 leading to singular nonlinear response.We benchmark semiclassical analysis for the Ising model against exact results using form-factor expansions, detailed in upcoming work [42].While the form-factor approach is applicable to a subset of integrable models, the semiclassical approach can be generalized more broadly.Within the same framework we can also treat non-integrable systems, as long as they feature stable, gapped QPs in some range of momentum.To confirm the validity of the semiclassical approach in this context we perform tDMRG simulations in the (non-integrable) AKLT spin-1 chain [43,44].Furthermore, with only slightly more work, the semiclassical approach can be generalized to treat [low] finite T .We also conjecture that many features persist in d > 1.
Relaxing the restriction to q = 0 to allow momentum resolution permits direct extraction of the QP scattering matrix by combining linear and nonlinear response data.Our work thus promises an intuitive way to compute and understand nonlinear responses in a variety of quantum many-body systems.
Setup.-As noted above, we initially focus on the transverse field Ising chain, with Hamiltonian H = −J In particular, we work in the paramagnetic phase (g > 1) and consider q = 0 perturbations coupling to the order parameter M = j σ z j (recall that only such 'Ising-odd' operators can excite single QPs above the ground state).Extensions to q = 0 are straightforward and will be reported in detail in [42].
The model is exactly solvable by means of a Jordan-Wigner transformation [45] that maps H to a quadratic fermion problem.This yields a QP dispersion relation (k) = 2J 1 + g 2 − 2g cos(k), with a gap ∆ = (0).However, since σ z maps to a non-local (string-like) operator, exact response functions involving M cannot be easily computed using Wick's theorem, and instead require a form-factor expansion using techniques developed in Refs [45][46][47][48][49][50][51][52][53].Similarly, any local spin operator that can create single fermionic QPs must be nonlocal in the fermion basis.Consequently, their nonlinear response is sharply distinct from that of fermion-local spin observables [4], that only excite even numbers of QPs.
We first study the pump-probe signal, which can be viewed as the difference in the linear response (measured between times t 1 , t 1 + t 2 ) computed in two states: the original QP vacuum |0 , and a 'pumped' state e −iµM (0) |0 obtained by perturbing the QP vacuum at t = 0 with a 'kick' of strength µ coupling to M [54].
In the response regime where the pump only weakly perturbs the system away from equilibrium, we can expand Ξ PP in µ and evaluate the resulting correlators in equilibrium: odd powers vanish by symmetry, so Ξ PP = µ 2 χ (3) PP + O(µ 4 ).The superscript denotes a thirdorder response, which we split into pieces divergent and convergent in time, χ is our focus.Only the connected correlator (denoted by the subscript C) contributes to the response, as required by causality.Both from the Heisenberg picture in (1) and the Schrödinger one (2), it is evident the correlators are not time-ordered.It is convenient to also distinguish the 'ket' and 'bra' sides of (2).Formally these correspond to forward and backward branches of the Keldysh time contour, which runs from t = 0 to t = t 1 +t 2 and back.At t = 0, the operator M acts on bra and ket sides, whereas at time t = t 1 it acts solely on the ket side.Both sides are then evolved up to time t = t 1 + t 2 , whereupon M is measured.This sequence is path-ordered on the Keldysh contour.We adopt the standard nomenclature where n th order response functions involve n external perturbations and n+1 operators; the extra operator corresponds to the measured observable.[Note that χ PP ∝ µ 2 , but appears at third order when expanding in terms of all external fields, as an extra perturbation is required to extract the linear response function in Eq. (1), cf.[54].] Long-time divergences and non-perturbative effects.-Ourmain result is that χ with a scaling function C PP whose numerical behaviour is shown in Fig. 2.This divergence can be given a simple semiclassical interpretation, involving ballistic propagation of quasiparticles and their scattering (Fig. 1).Before detailing the semiclassical analysis, we argue that we expect χ PP to diverge on general grounds.Even for arbitrarily small perturbation strength µ, Ξ PP will saturate to an O(1) value independent of µ at late times.This is most easily seen for sufficiently large t 1 , such that the system effectively thermalizes after the initial kick.The perturbed state e −iµM (0) |0 is then effectively at a finite (but small) T .[In the integrable case, it can be approximated by a generalized Gibbs ensemble (GGE), but this is not essential to the analysis.]The first line of (1) is then approximately a linear-response function in a finite-T state, which on general grounds is expected to decay exponentially with dephasing rate γ µ , e iµM (0) [M (t 1 + t 2 ), M (t 1 )]e −iµM (0) ∼ e −γµt2 , (4) intuitively due to stochastic scattering events with the QPs created by e −iµM (0) [55].Therefore, at long times, the effect of the perturbation become O(µ 0 ).A natural possibility is that γ µ ∝ µ 2 , suggesting that χ (3) PP diverges linearly in t 2 whenever we have a stable QP excitation at zero momentum.The full nonlinear response thus initially shows a linear divergence, probed in the response limit, which eventually saturates at late times t 1/µ 2 .
x t t Cartoon of processes contributing to the late-time divergence of the third-order response χ (3) Semiclassical picture.-Thescaling form (3) and the scaling function C PP can be quantitatively understood from a simple semiclassical picture inspired by the seminal approach of Refs.[55,56] (see also Refs.[57][58][59][60]).Our basic objects are wave-packet (WP) states |r, k of QPs, which we will think of as having approximately well-defined positions r and momenta k in the sense that the effects of the dispersion of the wave packets will be sub-leading.Multi-WP states |r; k are obtained as tensor products of single-WP states and by locality of the Hamiltonian (approximately) have simple dynamics as long as the WPs are spatially well separated.By construction n-WP states are in one-to-one correspondence with scattering states of n QPs.On the ket side, the action of the operator M (0) will after a sufficiently long time result in "intermediate" n-WP states (with n odd), where the individual WPs approximately originate from the same point x 0 (which is integrated over) and whose momenta sum to 0. Each WP (approximately) propagates ballistically with its group velocity v(k) = (k), i.e. e −iHt M (0) |0 is approximately a superposition of states of the form dx 0 |r; k with r specified by r j = x 0 +v(k j )t and k such that j k j ≈ 0, but otherwise arbitrary.We start by considering processes where all WPs proceed undisturbed to time t 1 + t 2 , whereupon they annihilate with the n WPs produced by M (0) on the bra side.Note that, in order to have a non-negligible overlap between the bra and the ket, the momenta k on the two sides must approximately coincide, as well as the position x 0 at which the WP shower is created (see Fig. 1).In the Ising model the amplitude associated with creating and annihilating a shower of n WPs with a given set of momenta is the so-called n-QP form factor on top of the ground state, whose precise value is unnecessary to proceed with the semiclassical calculation.[Note that in principle F n (k) and v(k) can be numerically computed for non-integrable systems using MPS [61].] In the processes of interest, M (t 1 ) creates a single q = 0 WP at x 1 , which is spatially well separated from the position at time t 1 of the n WPs produced by the action of M (0) on |0 .The WP created by M (t 1 ) is then annihilated by M (t 1 + t 2 ), giving rise to an amplitude |F 1 (0)| 2 e −i∆t2 .Naively, integrating over the arbitrary positions x 0 and x 1 generates a contribution to χ (3) PP that diverges proportionally to L in the thermodynamic limit.However, the contribution of processes in which the WP trajectories do not cross is simply equal to the product of two-point functions L −1 M (0)M (0) M (t 1 + t 2 )M (t 1 ) , and therefore the extensive spatial divergence cancels when we subtract the disconnected contributions.
In contrast, each time a pair of WP trajectories cross, the amplitude for the process picks up a factor of the scattering matrix S. In the Ising model, this simply encodes the fermionic statistics of QPs, i.e. S = −1.Therefore, a process involving an odd number of scattering events (cf.Fig. 1) does not cancel against disconnected contributions; we refer to processes which are connected only because of such events as scattering-connected.After subtracting the disconnected component we obtain a factor S − 1, evaluating to −2 in the Ising case.Finally, we must integrate over all possible x 1 that produce a scattering-connected process.[The integration over x 0 will produce a factor L, cancelling the L −1 in the definition of Ξ PP .]One can verify that, for a given set of momenta k of the WPs produced by the pump, the range of x 1 leading to scattering-connected processes has length v PP (t 1 + t 2 ), where v PP is a linear combination of the velocities of the various WPs [54].In this way we obtain a divergent contribution to χ Note that for g 1.1, C PP is dominated by C (3) PP , with higher values of n yielding numerically smaller corrections.This allows us to numerically estimate C PP , reported in Fig. 2 for representative values of g.
Other types of processes only give rise to subleading contributions at late times.Scattering-connected processes where M (t 1 ) creates more than one QP subsequently annihilated by M (t 1 +t 2 ) give contributions suppressed as t 2 → ∞.This follows since QPs created by M (t 1 ) spread ballistically from one another and therefore cannot be annihilated by a single local operator.For processes that are not simply scattering-connected, the position where all operators act is fixed by the ballistic propagation of the various QPs, therefore we expect their contributions to remain finite at late times.
Form-factor expansion.-Webolster the semiclassical result with an exact calculation of the late-time asymptotics of χ A key property of these matrix elements [62][63][64][65] is the presence of kinematic poles: as p i approaches k j the matrix element becomes singular, p|M |k ∼ 1 pi−kj .These contributions ultimately give rise to the late-time divergence in χ (3) PP;d [54].To benchmark the semiclassical picture, we consider its simplest non-zero contribution, from C (3) PP (sketched in Fig. 1).Counting the number of QPs before and after each operator M , we expect that C we see good agreement with the semiclassical expectation (see inset of Fig. 2).
We can also use the form-factor approach to compute the leading correction to (3), which scales as √ t 2 [66,67].
We conjecture that this can be interpreted as the effect of scattering-connected processes where M (0) and M (t 1 + t 2 ) create and annihilate a single QP, which can scatter with the QP being exchanged between A(t 1 ) and A(t 1 + t 2 ) when their trajectories are smeared by the broadening of QP wave-packets due to dispersion, omitted in the leading semiclassical analysis.
We have computed χ (3) PP;d (q 1 , q 2 ; t 1 , t 2 ) using tDMRG.The results for t 1 = 0, q 1 = π, and q 2 = 2π/3 are reported in Fig. 3.At sufficiently late times χ (3) PP;d is wellfitted by the functional form At 2 sin( (q 2 )t 2 −φ), which is consistent with our wave-packet analysis.Here A and φ are fitting parameters, while (q 2 ) is independently determined from the numerical computation of the two-point function.
Discussion.-We have shown that in the transversefield Ising model the nonlinear pump-probe response is characterised by a linear-in-time divergence.We will now argue that this is in fact a general feature of any interacting translationally-invariant many-body system that has a stable, gapped single-particle excitation.Such a behaviour is suggested by our analysis where we demonstrated that the late-time pump-probe signal Ξ PP can be re-expressed as the difference of two-point correlation functions at zero and finite temperatures and hence is O(1); the divergent response reconciles this result with the perturbative expansion in powers of the applied fields.We therefore expect this linear-in-time growth to be quite general and as we have argued, in many cases visible also on intermediate timescales at low finite temperature.We have developed a semiclassical picture of WP propagation and scattering that identifies the processes that give rise to the divergence.Our discussion has focused on the simple case of the Ising model; however, the form factor calculations in fact partially generalizes to other integrable theories as is shown in a forthcoming work [42].More importantly, the semiclassical arguments generalize to non-integrable models, and even to finite temperature [54].
An enticing possibility suggested by our work is the measurement of scattering matrices from third-order response functions.In χ (3) PP;d (q 1 , q 2 ; t 1 , t 2 ), if the n = 1 contribution is the dominant one -as can e.g.be achieved using a frequency-modulated pump with negligible amplitude to excite the system at energies greater than 2∆ -χ PP can be expresses in terms of the scattering matrix S(q 1 , q 2 ) and data that can be extracted from linear response, allowing us to read off S(|q|, 0) from the divergent piece of χ PP .[54] Finally, it would be interesting to understand if similar late-time divergences emerge in d > 1.The nonperturbative argument is evidently independent of dimension, suggesting that this is indeed the case.We leave a more detailed investigation of this intriguing possibility to future work.It would also be interesting to understand if these or similar mechanisms are responsible for the anomalously large nonlinear response recently observed in 2DCS experiments [12].

I. PUMP-PROBE SIGNAL AND ITS PERTURBATIVE EXPANSION
The pump-probe setup can be quite generally described by a time-dependent Hamiltonian For our purposes it is sufficient to consider where f (t) is a smooth function of time and we will take µ to be a very small parameter.Then we have At later times the quantum state of the system is then with T indicating that the exponential is time-ordered.Defining expectation values arXiv:2208.09490v2 [cond-mat.str-el]8 Feb 2023 where the response functions are given by Note that we have chosen to normalize every response function by 1/L, so that it remains finite in the L → ∞ limit.Finally, we imagine repeating the linear response measurement with µ = 0 and subtracting the two results to obtain This is the form of the pump-probe response given in eqn (1) of the main text.Choosing our initial state before the pump process to be the ground state |ψ(t = 0 − ) = |0 and then expanding Ξ PP (t, t ) to second order in µ we finally arrive at the definition of χ PP;d given in Eq. ( 2) Here {•, •} denotes an anticommutator and in the last line we have used that

II. LONG-TIME DIVERGENCES AND NON-PERTURBATIVE EFFECTS IN INTEGRABLE AND NON-INTEGRABLE SYSTEMS
In the main text, we argued that the late-time divergence of χ PP;d can be understood on general grounds, by noting that at late times the perturbation of strength µ has O(µ 0 ) effects.While in the main text we framed our arguments for non-integrable systems, in this section we further elaborate on the connection between long-time divergences and non-perturbative effects in integrable systems, where the argument can be presented in a more controlled fashion.
As in the main text, we work in the limit of large t 1 , when we can assume that the density matrix can be approximated by a generalized Gibbs ensemble (GGE), i.e. the density matrix takes the form ρ ∼ exp(− k λ k n k ), obtained by maximizing entropy subject to the conservation of the number of QPs n k in mode k.Here the λ k could in principle be fixed by imposing Tr(ρ n k ) = 0|e iµM n k e −iµM |0 .The first line of Eq. ( 1) is then a linear-response function on top of ρ.Given that ρ is a finite energy-density state, this is known to decay exponentially 1 , viz.
with γ µ determined by the set {λ k }.Thus Ξ PP becomes a difference of correlators at finite and zero energy density, and as the latter exhibits undamped oscillations in time we have Ξ PP ∼ O(µ 0 ).On the other hand, expanding Ξ PP in a power series in µ (at large, finite L) gives Ξ PP = n≥2 µ 2n g 2n (t 1 , t 2 ).This can only become O(µ 0 ) at late times t 1,2 if the limits of µ → 0 and t 1,2 → ∞ do not commute and the expansion is not uniformly convergent.A natural [Reproduced from Fig. 1 of the main text.]Cartoon of processes contributing to the late-time divergence of the third-order response χ P P ;d ∼ M (0)M (t1 + t2)M (t1)M (0) .A dashed line, corresponding to time t = t1 + t2 separates the bra (left) and ket (right) sides of the time evolution; time increases towards this line.Solid lines denote QP worldlines, and circles denote the action of the operator M .When two lines cross, the amplitude for the diagram is multiplied by S = −1.For a fixed x0, the red segments indicate the set of x1 giving rise to a scattering-connected contribution.The length of the red set is proportional to the overall timescale, leading to the linear divergence of χ possibility is to have an "infrared divergence" that needs to be summed to all orders in µ.In the Ising case we expect exponentiation of the dominant contribution in the late-time regime, so that These considerations suggest that χ PP diverges whenever we have a stable QP excitation at zero momentum, which gets damped at finite energy densities.Thus the nonlinear response initially shows a linear divergence, probed in the response limit, that eventually saturates at late times t 1/µ 2 .

III. EXPLICIT FORM OF vPP
Referring to the process depicted in Fig. 1 of the main text (also reproduce in Fig. S1 for the reader's convenience) the set of x 1 leading to scattering-connected processes has length v PP (t 1 + t 2 ), where v PP has the dimension of a velocity.We give here an explicit formula for v PP in terms of the velocities of the quasiparticles.Denoting with k 1 , . . ., k n the momenta of the WPs produced by the pump, v PP is obtained as follows: we define the set of positions

IV. EXTENSION TO TWO-DIMENSIONAL COHERENT SPECTROSCOPY (2DCS)
While the main text focuses on pump-probe spectroscopy, similar conclusions apply to 2DCS experiments, where the system is perturbed by a pair of field pulses at varying times.More explicitly, in Eq. ( 2) for 2DCS we have f (t) = δ(t − t 1 ) and µ is comparable to µ.By performing the experiment several times, one can measure We can then proceed to expand Ξ 2DCS perturbatively in µ and µ up to third order The starting point is the definition of single-QP wave-packet states in terms of asymptotic scattering states: Here is some arbitrary microscopic distance, which we assume to satisfy ξ.The evolution of a wave-packet state is simply given by up to corrections which can be neglected by choosing s.t.(k)(t 1 + t 2 ).Similarly, one can define wave-packet states in terms of asymptotic scattering states as Their time evolution will generally be complicated, however it simplifies in the asymptotic region, i.e. when |r i −r j | ∀ i, j.In this case, each r j coordinate will evolve ballistically with the group velocity v(k j ).Furthermore, whenever two wave-packets come in proximity of each other, they will scatter.In 1D and for systems with only one QP species, this will amount to the wavefunction acquiring an overall extra phase: the scattering matrix.

VII. SEMICLASSICAL PICTURE FOR NON-INTEGRABLE MODELS
We now detail how the wave-packet description above can be used to naturally extend the analysis in the main text to non-integrable models.
The starting point is characterizing the state M |0 .In non-integrable systems this state will generally be complicated since M can create multiple excitations in the neighborhood of the same point.After some time t, however, we expect asymptotic QP states to emerge.In this case we can write with r specified by r j = x 0 + v(k j )t.Importantly, for a fixed k, the integral over y is negligible, unless j k j −1 .
Here F denotes the probability amplitude to produce a given wave-packet state, and is generally challenging to compute from first principles.Note that in this subsection, the " " sign between two states should not be intended in the sense that the difference between the two has a small norm, but rather that we are isolating a component of the state that will ultimately give the leading contribution to χ PP;d .
The validity of the equation above can be argued as follows.In fact, (23) is not the most general form for a WP state, which would instead be for some function F n and with r j = x j +v(k j )t.[In describing the time dependence we have assumed that no scattering process is taking place.This is justified at sufficiently long times as the n WPs are moving away from each other.]The relation above is to be understood in the sense that both sides coincide when projected onto scattering states of quasi-particles, and is an approximation in non-integrable models as it neglects components of the state not captured by the QP description.Nonetheless, given that the LHS is a low-energy state, we expect this approximation to yield accurate results for most models.Furthermore, the initial positions x j must necessarily be close: i.e. |x i − x j | ξ A for some microscopic length ξ A .However, since we are interested in divergences originating at large lengthscales and timescales, we can choose ξ A , so that we can approximate x j x 0 for all j and F n (k, x j ) F n (k, x 0 ).Finally, the fact that the left-hand side is invariant under translation by any length a requires the x 0 −dependence of F n (k, x 0 + a)e ia j kj = F n (k, x 0 ), implying that F n (k, x 0 ) is of the form e −ix0 j kj F n (k).
Finally, for the semiclassical analysis of the divergence we also need to use that Here the first term amounts to adding a zero-momentum WP at position r , which is possible if |r − r j | for all j.The second line accounts for processes where M annihilates one of the previously-present WPs, instead.Finally, the • • • denote more complicated outcomes, e.g., where more then one WP is created or annihilated.
We can now discuss in detail the bra and ket side of Eq. ( 2).The bra side can be simply expanded as in ( 23) with r 4 = x 0 + v(k j )(t 1 + t 2 ).The ket side can be built progressively.|ψ 1 = U (t 1 )M |0 is given by (23).We are then interested in the scenario where M produces an extra WP at q = 0 on top of |ψ 1 : The states can then be evolved semiclassically up to time t 2 + t 1 .Along this evolution, any of the WPs produced at t = 0 crosses the trajectory of the q = 0 WP produced at t = t 1 , the overall wavefunction acquires a phase given by the product of the scattering matrices for each crossing.We denote this phase as S(x 0 , x 1 ; k).Finally the ket side is given by |ψ 3 = M U (t 2 ) |ψ 2 .Here, as described in the main text, we consider the case where M annihilates the q = 0 WP, thus producing √ n! dx 1 dx 0 d n k (2π) n−1 e −ix0 j kj F n (k)e −i j (kj )(t1+t2) |r; k S(x 0 , x 1 ; k) + . . ., (28) with r j = x 0 + v(k j )(t 1 + t 2 ).
Finally, we can take the overlap between the bra and ket ψ 4 |ψ 3 .Comparing the expressions for the two we see that only the terms with k − k −1 and |x 0 − x 0 | give a significant overlap.This gives an overall view of how to extend the semiclassical picture to the non-integrable case.The crucial point, however, remains that for any values of the other integration variables s.t.j k j ≈ 0 we have dx 0 (S(x 0 , x 1 ; k) − 1) ∝ (t 1 + t 2 ), (29) meaning that upon rescaling both times by some factor λ, the integral is also rescaled by λ.These considerations therefore lead us to the scaling form (3) also in non-integrable theories.A more detailed analysis will be presented in Ref. 8.

VIII. DETAILS OF THE TDMRG SIMULATIONS
To compute χ PP;d numerically we proceed as follows.We consider a finite system of size L and compute its ground state |0 through DMRG.[Alternatively, for the AKLT chain the Matrix-Product-State (MPS) form of the ground state is explicitly known.]We then compute the two states |ψ 2 = e iq1L/2 U (t 2 )M (q 2 )U (t 1 )S z L/2 |0 (31) by using tDMRG for the time evolution and by encoding M (q) as a Matrix Product Operator (MPO) and using standard algorithms to compress MPO-MPS products.Finally, we compute χ [ ψ 1 |M (−q 2 )|ψ 2 ] − C2 (q 1 , 0) C 2 (−q 2 , q 2 , t 2 ) + − C * 2 (q 2 , q 1 , t 1 + t 2 ) C2 (−q 2 , t 1 ) + C * 2 (−q 2 , q 1 , t 1 ) C2 (q 2 , t 1 + t 2 ) ( where C2 (k, t) = e iq1L/2 0|M (−k, t)S z L/2 (0)|0 (33) Here E 0 denotes the ground-state energy.We note that while we have used translational invariance to obtain the above form, the numerical computations is carried out for an open chain.This is justified for sufficiently large system sizes.The asymmetry between the C2 and C 2 is to guarantee that the disconnected component of ( 32) is properly subtracted.[The terms in the second line are 0 in a system with periodic boundary conditions if q 1 = q 2 , but not necessarily in a systems with open boundary conditions.Nonetheless, they are O(L 0 t 0 1,2 ), so they do not play an important role in our analysis.] In our numerical calculation (Fig. 3) the system size is L = 200, the maximum bond-dimension is D = 400, and the Trotter step δt = 0.01, and we checked that our results are converged in bond dimension and Trotter step.Simulations have been performed using the ITensor library 9,10 .
Finally, we note that other techniques to compute nonlinear response functions using MPS techniques have been recently developed 11,12 .
diverges at late times, line, corresponding to time t = t1 +t2 separates the bra (left) and ket (right) sides of the time evolution; time increases towards this line.Solid lines denote QP worldlines, and circles denote the action of the operator M .When two lines cross, the amplitude for the diagram is multiplied by S = −1.For a fixed x0, the red segments indicate the set of x1 giving rise to a scattering-connected contribution.The length of the red set is proportional to the overall timescale, leading to the linear divergence in (3).
PP,d (t 1 + t 2 , t 1 ) + χ Cartoon of low temperature semiclassical processes.The thermal state is approximated as an ensemble of WP states representing thermal QPs, which follow the same trajectories on bra and ket sides (dashed red lines).As they propagate, the thermal QPs can scatter against those produces by M (0) (orange dots) and by M (t1) (red dot).Each orange dot on the bra side is canceled by its partner on the ket side, and hence they do not affect the amplitude.This reduces the problem in essence to that solved in Ref. 2.