Quantum caustics and the hierarchy of light cones in quenched spin chains

We show that the light cone-like structures that form in spin chains after a quench are quantum caustics. Their natural description is in terms of catastrophe theory and this implies: 1) a hierarchy of light cone structures corresponding to the different catastrophes; 2) dressing by characteristic wave functions that obey scaling laws determined by the Arnol'd and Berry indices; 3) a network of vortex-antivortex pairs in space-time inside the cone. We illustrate the theory by giving explicit calculations for the transverse field Ising model and the XY model, finding fold catastrophes dressed by Airy functions and cusp catastrophes dressed by Pearcey functions; multisite correlation functions are described by higher catastrophes such as the hyperbolic umbilic. Furthermore, we find that the vortex pairs created inside the cone are sensitive to phase transitions in these spin models with their rate of production being determined by the dynamical critical exponent. More broadly, this work illustrates how catastrophe theory can be applied to singularities in quantum fields.


I. INTRODUCTION
According to Lieb and Robinson [1], there is a maximum speed v LR at which information can propagate in discrete quantum systems that obey the Schrödinger equation and have short range interactions. This is a powerful and generic statement because it implies that, despite the fact there is no intrinsic speed limit in the (non-relativisitic) Schrödinger equation, the response of these many-particle systems to a sudden quench should be in terms of a light cone-like time evolution of spatial correlations [2]. Physically, the "light cone" arises from the maximum group velocity of quasiparticles that are excited by the quench and that subsequently propagate through the sample [3]. Sophisticated methods of analysis have been applied to these quench problems including conformal field theory and tensor networks [3][4][5][6][7][8][9][10][11][12][13][14][15][16], and the theory has been tested in experiments on ultracold atoms [17][18][19] and ions [20,21] where quantum spin models [18,[22][23][24][25], the Bose-Hubbard (BH) model [26][27][28][29], 1D systems [30][31][32], and quantum walks on a lattice [33,34] can all be realized. The long coherence times of atomic systems make them particularly suited to studying such dynamics [35,36], and the ability to perform single-site manipulation and detection [37][38][39][40] has enabled unprecedented preparation and visualization of the relevant local observables.
In this paper we show that light cones in quenched spin chains are quantum caustics. These are quantum versions of wave focusing phenomena that occur widely in nature in the form of rainbows [41], ship wakes [42][43][44], tsunamis and tidal bores [45], and Cherenkov radiation [46] (including superfluid analogs [47][48][49]). In the geometric ray theory caustics occur where two or more rays coalesce, giving regions in space where the intensity diverges. By virtue of their singular nature, the natural mathematical description of caustics is via catastrophe theory which partitions them into a hierarchy of equivalence classes, each of which is structurally stable and has its own set of scaling relations [50][51][52]. To show specifically how this approach can be applied to spin chains we consider the exactly solvable 1D XY model [53,54], as well as the special case of the 1D transverse-field Ising model (TFIM) [55,56]. While both cases display light conelike behaviour, the more general XY model allows for an anisotropic coupling giving rise to a double cone [57,58]. Although we limit our calculations to these exactly solvable models, the structural stability of catastrophes (insensitivity to small perturbations) guarantees they must survive in the presence of weak non-integrability. This includes weak interactions between quasiparticles or disorder and therefore our results also apply to more general systems than just exactly solvable models.
Wave interference softens caustics and leads to structure on three scales [52]: at large scales we see divergent ray caustics, whereas at wavelength scales interference smoothes the divergences and dresses each caustic with a characteristic wave function which in the simplest case of two coalescing rays is the Airy function, and finally at subwavelength scales there are networks of vortexantivortex pairs. These robust features, including vortexantivortex networks, have been observed in optical fields [41], and more recently in electron microscopy [69]. They have also been discussed theoretically in the context of Bose-Einstein condensates [73,74] and various aspects seen experimentally in these systems [70][71][72]. Furthermore, the association between the Airy function (and its related kernels) and light cones has previously been noted by various authors [8,14,15,[59][60][61][62][63][64], and recent work has conjectured similar universal forms for wavefronts of outof-time-ordered correlators [65][66][67][68] by examining asymptotic limits of the Airy function. However, to the best of our knowledge the present paper is the first to study the hierarchy of universal wave functions that dress light cones, of which the Airy function is only the first, and also point out that light cones should generically contain networks of vortices which in the case of 1D chains appear as space-time vortices.
A fourth scale appears in quantum fields due to discretization of excitations leading to 'quantum catastro-
phes' [75][76][77][78][79][80] (rippling mirrors give analogous effects [81]). Going to the continuum (classical field) limit returns us to a wave catastrophe. As we shall show, light cones in spin chains have all the features of quantum catastrophes, including discretized versions of wave catastrophes and vortices which are regulated by the lattice constant. Although the cone itself is mildly affected by the presence of a quantum critical point (QCP) in the spin models we study, we find by contrast that the vortices are strongly affected and we use this feature to extract the dynamical critical scaling. The rest of this paper is organized as follows: In Sec. II we outline the relevant aspects of catastrophe theory, emphasizing the hierarchy of structures and their scaling properties. In Sec. III we show that light cones are in fact (quantum) caustics and hence their natural mathematical description is via catastrophe theory. In Sec. IV we introduce the XY and TFIM spin chains focusing on the quasiparticle dispersion relation which is the key ingredient we need to apply catastrophe theory. This program is implemented in Sec. V where we obtain the Airy and Pearcey functions for the wavefunctions dressing the fold and cusp catastrophes/cones in these models. In Sec. VI we verify the self-similar scaling properties of light cones that catastrophe theory predicts and in Sec. VII we describe how higher order catastrophes arise in the context of correlation functions. In Sec. VIII we identify and discuss the presence of vortex-antivortex pairs within light cones, while in Sec. IX we touch on the relevance of the theory to quench experiments, and in Sec. X we conclude with a discussion of the broader significance of the results. In order to make this paper self-contained we have included in appendices A-F the specifics of quantum spin chain diagonalization methods and various other details of our calculations.

II. GEOMETRIC AND WAVE CATASTROPHES
In what follows we will not need the full mathematical machinery behind catastrophe theory, but we will make use of a number of key results and for this reason we give a brief overview here. Our treatment is informal, but we emphasize that these results can be proved rigorously. The main idea can be stated simply: catastrophe theory classifies structurally stable singularities of functions and shows that such singularities can only take on certain characteristic shapes [50]. In up to four dimensions these are René Thom's seven elementary catastrophes which are listed in Table I.
Each catastrophe arises from two or more coalescing/bifurcating stationary points of its generating function Φ Q , the normal forms for which are given in the table. In the physical applications given in this paper Φ Q is the action functional and stationary points therefore correspond to classical paths or rays. From an optical/classical mechanics point of view a catastrophe is a caustic, i.e. the locus of points where the ray density diverges.
Thom's theorem states that the local behaviour of a function near coalescing stationary points can always be mapped by a smooth change of variables onto one of the catastrophes and in this sense catastrophes are universal. There is also a second sense in which catastrophes are universal: structural stability means stability against perturbations and thus catastrophes do not require special symmetry and hence occur generically in nature. Perturbations do not qualitatively change catastrophes and only quantitatively affect behaviour up to the strength of the perturbation.
The catastrophes in Table I are organized by the number n of state variables (their corank), and by the dimension Q of the control parameter space. Control space is the space where the function with its singularities actually lives. The control parameters C = {C 1 , C 2 , . . .} could be space and time coordinates as well as any other parameters. The state variables s = {s 1 , s 2 , . . . } characterize the rays. The simplest catastrophes (the cuspoids) have n = 1 and their generating functions are polynomials of the form with up to Q coalescing stationary points. The station-arity condition reads ∂Φ Q ∂s = 0 (2) and corresponds physically to Hamilton's principle of stationary action, while caustics arise from coalescing stationary points where the generating function is stationary to higher order [52] In the examples we provide in subsequent sections, we focus primarily on the fold and cusp catastrophes, as well as a discussion of the hyperbolic umbilic in the context of correlation functions. Folds and cusps are the only structurally stable singularities in the 2D (x, t) control plane where light cones in 1D chains live, while the higher catastrophes (although they may still exist in greater dimensions) can only be projected onto the plane by way of cusps and folds. This property is generic: catastrophes of higher order contain the lower ones [51]. The cusp is the meeting of two fold lines, the swallowtail contains two cusps, and so on. The wavefunctions, or wave catastrophes, associated with catastrophes can be obtained in a way analogous to Feynman path integrals by exponentiating the generating function and integrating over all paths, where λ plays the role of the wavenumber k or 1/ in quantum problems. In this form, the fact that the generating function plays the role of the physical action becomes clear. These functions are also known as diffraction integrals and many of their properties have been tabulated [44]. We emphasize that standard approximations such as the method of stationary phase where the integral over s is broken up into a sum of independent gaussian integrals around each of the stationary points are doomed to failure when the stationary points coalesce. One must instead keep the full form of Φ Q to get a result which is uniformly correct through the coalescence regions and this is precisely why diffraction integrals are crucial for treating bifurcation problems where solutions appear or disappear. The fold has a cubic action Φ 1 (s; C) = s 3 /3 + Cs, where in the case of a light cone in (1+1)-dimensions C = C(x, t). As the control parameter C is taken from positive values down through zero the cubic changes its form so as to describe two coalescing rays. The resulting wave catastrophe can be recognized as the integral form of the Airy function, In the absence of any special symmetry, two fold lines generically meet at cusps. In the region near the cusp 3) and which is shown as a black dashed line. The cusp is made of two fold lines which meet at the cusp tip at C1 = C2 = 0. There are three rays/waves inside the cusp and only one outside: two coalesce as we cross either of the fold lines, but all three coalesce at the cusp tip which is the most singular part of the classical caustic (a ray picture of the cusp can be seen in Fig. 2b in [80]). However, wave interference removes the classical singularities. The black dots show the locations of vortices: there is a line of vortices outside either edge of the cusp, and vortex-antivortex pairs inside. point the appropriate action is quartic and features two control parameters Φ 2 (s; C 1 , C 2 ) = s 4 /4 + C 2 s 2 /2 + C 1 s. This normal form, which formally resembles the Landau free energy for a continuous (2nd order) phase transition, describes the coalescence of up to three rays and results in a wave catastrophe known as the Pearcey function, which is a complex function of two variables. For our definitions/conventions for the Airy and Pearcey functions, see Eqns. (D5) and (C7), respectively. Plots of the absolute values |Ai(C)| and |Pe(C 1 , C 2 )| of the Airy and Pearcey functions are given in Fig. 1.
The fact that the Pearcey function is a two-dimensional complex function, with an amplitude and a phase at each point, allows for the possibility of vortices. This turns out to be the case: the black dots in Fig. 1(b) show the locations of vortices, or more precisely their cores. There is an ordered network of vortex-antivortex pairs inside the cusp and single rows of vortices lining the outer edges. These are subwavelength features that represent the finest layer of structure of a wave catastrophe. We find the vortices by densely covering the plane with loops around which we integrate the phase of the Pearcey function: loops that contain vortices give a ±2π phase change (the vortex cores also correspond to nodes of the Pearcey function, although in principle not all nodes need be vortices). An important feature of wave catastrophes is that they exhibit self-similar scaling. If the parameter λ is changed from λ to λ the wavefunctions will retain their functional forms but with rescaled coordinates, We can understand this scaling as follows: the overall amplitude scales as λ β Q , where β Q is known as the Arnol'd index. The distance between interference fringes is also rescaled, but generally the scale factor is different in each direction according to λ ςm , where ς m is the Berry index associated with coordinate C m . For the fold wave catastrophe β Ai = 1 6 and ς = 2 3 , and for the cusp wave catastrophe β Pe = 1 4 and ς = { 3 4 , 1 2 }. A complete list of Arnol'd and Berry indices for the seven elementary catastrophes is displayed in Table I.
The sets of Arnol'd and Berry indices accompanying the different catastrophes are reminiscent of the sets of critical exponents which define universality classes of equilibrium phase transitions. The underlying common cause of this similarity is the presence of singularities and singularities lead to universality. However, we emphasize that in the light cone case we study here, this universality occurs in non-equilibrium dynamics, and thus we have an example of universality in quantum dynamics.

III. LIGHT CONES AS QUANTUM CAUSTICS
Our approach to the light cone problem is based upon the idea that the build-up of correlations occurs through quasiparticle propagation [3]; this is known to be the case in a broad range of models including the BH, TFIM, and XY models. The Lieb-Robinson bound can then be expressed in terms of the maximal group velocity of quasiparticles [7,9] v LR = max k d k dk (8) where k is the dispersion relation for quasiparticles as a function of quasimomentum k. It can be seen immediately that this result is exactly equivalent to Eqns. (2) and (3) which give the conditions for a caustic (note that here we are implicitly considering real solutions to the caustic conditions; imaginary solutions correspond to phase velocity across the cone and are discussed in Appendix D. This aspect has also been discussed by Cevolani et al. in Ref. [16]). From this simple observation it follows that light cones are caustics and hence the results and insights of catastrophe theory can be applied to them. Let us focus on the case of a local quench where a single quasiparticle is created at position x = 0 in the middle of a spin chain (we briefly consider weakly nonlocal superpositions of multiple quasiparticles in Sec. IX, and also in Appendix E). Time evolving the state with the Hamiltonian H, the state vector at time t is where |0 b is the Bogoliubov quasiparticle groundstate and the operator b † x creates a quasiparticle at the site located at position x. For the remainder of the paper, we use the subscript 'b' to distinguish Fock states in the Bogoliubov basis from the Jordan-Wigner basis. Introducing the eigenstates |k of H we can write this as (see Appendix A for details) where N is the number of sites, and the phase θ(t) ≡ t/(2 ) k k is not observable but is included here for completeness. Projecting onto the position basis, the wavefunction Ψ(x n , t) ≡ x n |Ψ(t) on the n th lattice site is In these expressions n is an integer lying in the range {−(N − 1)/2, ..., (N − 1)/2}, and the separation between momenta in the sum is ∆k = 2π/(aN ).
In the continuum approximation (CA) the wavefunction corresponding to Eq.
where a = L/N is the lattice constant for a lattice of length L, and the quasimomentum k runs over the first Brillouin zone. A comparison of the exact (discrete) and CA wavefunctions is given in Fig. 8  and coalescence points (especially the latter). By Thom's theorem [50][51][52], we can therefore map Φ onto one of the normal forms Φ Q . However, although Thom's theorem guarantees that this can be done by smooth transformations, it does not tell us what these transformations actually are. Figuring out the mapping is part of the challenge in applying catastrophe theory to specific physical problems and it is to this task that we now turn.

IV. XY AND TFIM SPIN CHAINS
Let us consider a 1D XY model describing spins on a lattice interacting with a ferromagnetic coupling J, anisotropy parameter γ, and subject to an external field gJ. The Hamiltonian is where σ α i , α ∈ {x, y, z}, are Pauli operators. When γ = 1 this Hamiltonian reduces to that of the TFIM. The XY Hamiltonian can be diagonalized via the Jordan-Wigner transform followed by a Bogoliubov rotation, which maps spin operators to spinless fermions [82]. As shown in Appendix B, this leads to the free model is the annihilation (creation) operator for Bogoliubov modes with quasimomentum k and dispersion Thus, the phase/generating function in Eq. (12) takes the specific form (16) An exact numerical evaluation of the wavefunction given in Eq. (11) using the generating function Φ(k; x, t) for the XY model is plotted in Fig. 2. The fact that x n is discrete means that the light cone actually corresponds to a quantum catastrophe, for more discussion of quantum catastrophes in a spin context see Ref. [80]. However, in the semiclassical regime where N is large, the CA described by Eq. (13) works well. In this case Φ has the same functional form but with x and k taken as continuous variables, and the integral can be evaluated analytically in terms of Airy and Pearcey functions as will be explained in the next section.
Dividing Φ(k, x, t) as given in Eq. (16) by t we can identify three control parameters: (x/t, γ, g) [we reserve the energy scale J to play the role of k in Eq. (4)]. However, rays propagate in the 2D (x, t) plane rather than the full 3D control space and thus for generic values of the control parameters catastrophe theory predicts we should see folds and cusps. In fact, we find a double cone made of a cusp enclosed by two folds as shown in Fig. 2 (double cones occur both in spin systems and in coupled 1D gases [57,58]).
Mathematically speaking, the double cone arises because Eq. (16) has up to four stationary points within the first Brillouin zone, as shown by the green dots in the five overlays plotted in Fig. 2(b). Near the origin in Fig. 2 all four stationary points are present, but three are quasi-degenerate so Ψ is locally dominated by a Pearceylike function which gives the inner cone. As we cross the edges of the inner cone two stationary points annihilate (indicated by red stars in the overlays) leaving two rays which in turn annihilate at the edges of the outer cone so that locally it is dominated by an Airy function. Furthermore, the XY model has a QCP at g = 1 − γ 2 ; as the critical regime is approached the inner cone narrows and eventually collapses because the three inner stationary points in the generating function coalesce at this value of g. In the case of the TFIM (γ = 1) [4,83,84], Φ has only two stationary points and one finds a single cone with edges that are dressed by Airy functions. The insight from catastrophe theory is that the single cone is non-generic and only occurs due to the special symmetry of the Hamiltonian when γ=1.
Due to the presence of four stationary points, the careful reader might expect the XY model to show signatures of the swallowtail catastrophe. Indeed, this would generically be true, however it can be verified that the quadruple root coalescence do not occur for real k. It is the periodic dispersion relation of the model which keeps us from physically probing the highly-singular swallowtail point. The cusp and fold catastrophes that we observe here are however inherited from the part of the swallowtail which is physically permitted.

V. AIRY AND PEARCEY FUNCTIONS
Let us now demonstrate explicitly how the Airy and Pearcey catastrophe integrals emerge in the CA. Starting with the Pearcey integral, consider first the triple sta-tionary point coalescence responsible for the inner cone, which we have isolated in Fig. 2(c). One obvious difference between this wavefunction and the Pearcey function shown in Fig. 1 is that the cone boundary in the former is straight rather than the standard curved form of the cusp C 1 = ± 4C 3 2 /27. Physically, this is due to the free propagation of the fermionic quasiparticles. The required transformation to take us between physical coordinates and those of the standard curved cusp is similar to that used by Kaminski and Paris in Ref. [89]. In Appendix C we show that for our spin model it is and we have defined the Ising velocity, which is equal to v LR in the TFIM limit (in principle, v LR can be analytically solved for in closed form for general γ, however, the expression is complicated, and little physical insight is gained from writing it here).
To complete the diffraction integral we also need the integration variable s. This reads s = √ 2a (tΓ) 1 4 k and results in the Pearcey-like wavefunction Ψ Pe (C 1 , C 2 ; J) written out in Eq. (20). It rapidly tends to a true Pearcey function at longer times when S = √ 2π (tΓ) 1 4 1.
In order to display the close resemblance between Ψ Pe and the Pearcey function, we have plotted in Fig. 3 the wavefunction of the inner cone from Eq. (11) without expansions or approximations in terms of the transformed coordinates C 1 and C 2 . This can be compared with the actual Pearcey function plotted in Fig. 1. The only significant deviation is near C 2 = 0. Since the limit of integration S tends to 0 as t → 0, the cusp point itself becomes poorly defined, and we get a 'smearing' of the wavefunction as C 2 → 0. As a consequence, we cannot get a Pearcey function exactly at the origin, since the initial boundary condition requires the real-space wavefunction be entirely localized here. As we move away from the cusp point, however, the Pearcey function is indeed an excellent approximation to the true wavefunction.
As C 2 increases the Pearcey function can be approximated by two back-to-back Airy functions as the cusp evolves into two fold lines. Indeed, it is a general property of catastrophes that the higher ones evolve into the lower ones as we move away from the former's most singular points. This provides a rigorous explanation for why Airy functions, which are the simplest of the hierarchy of wave catastrophes, are commonly encountered in the asymptotics of light cones [8,14,15,[59][60][61][62][63][64].
To examine how the Airy function emerges in the CA we specialize to γ = 1 (TFIM Hamiltonian). We stress that the choice of γ does not affect the presence of the fold catastrophe (and thus Airy functions), only the simplicity of the subsequent calculations. To this end, note that for any g = 1 it can be readily checked that Φ(γ = 1) in Eq. (16) has only two stationary points as a function of k. We can therefore map onto the canonical fold generating function Φ 1 (s; C) by expanding Φ to third order in s. In the CA, and up to a global phase, we show in Appendix D that the correct control parameter in this case is The index j ∈ {1, 2} refers to cases g > 1, and g < 1, corresponding to above and below the QCP, respectively. The integration variable s j = (g 2−j t) When γ = 1 this process may be repeated around each fold catastrophe, including for any inner cones, and will result in the emergence of Airy functions with different definitions of the control parameter, C. For example, a particular limit of Eq. (21) has been conjectured to give a universal form for the wavefront of out-of-time-ordered correlators (OTOCs) [65][66][67][68]. According to catastrophe theory this is no surprise. Furthermore, closer to the 'brightest' parts of the OTOC the hierarchy of catastrophes allows for more elaborate structures beyond the Airy function.  23), since the initial condition precisely cancels the Arnol'd scaling to preserve particle number. Panel (d): The oscillation period, T , of Ψ(xn, t) for site x/a = 5 in the TFIM with g = 3; Eqns. (11) and (33) are plotted in blue circles and orange triangles, with blue-solid and orange-dashed trendlines, respectively (a range of 1 ≤ J/J ≤ 30 was chosen). Accounting for a geometric factor of sin[arctan (20)], we find the Berry index to be 0.654 ± 0.003 and 0.646 ± 0.009 for Ψ and ΨX, respectively.

VI. SCALING
The way the spin coupling strength J and the control parameters C appear in combination on the right hand sides of Eqns. (20) and (21) shows that light cones have nontrivial scaling properties: varying J is equivalent to rescaling the amplitude and coordinates. More specifically, increasing J causes the amplitude to increase at a rate determined by the Arnol'd index, and the interference patterns to oscillate more quickly in space and time at rates determined by the Berry index for each particular direction. The overall picture is that the fringes flow in towards the origin as J is increased and in the (singular) classical limit, which occurs when J → ∞, all wave structure is pulled into the origin. There are other choices we could have made for the scaling parameter since it need only fill the role of λ in Eq. (4): for the TFIM we could have alternatively chosen a or g, and in the case of the XY model we could also have chosen either of these or even γ. It is usually necessary to keep some physics constant during the scaling: we can keep the position of the classical ray caustics constant as J is varied by tuning a or g to keep v I unchanged.
Numerical verification of the catastrophe theory predictions for both the Arnol'd and Berry indices for the exact wavefunction Eq. (11) is presented in Fig. 4. Panels (a)-(c) show the scaling in the inner cone of the XY model: the fringe scaling is obtained by measuring the distance between peaks of the wavefunction along coordinates C 1 and C 2 as J is varied and match the Pearcey scaling given in Table I to within 1%. At first glance, it appears that panel (c) shows a contradiction between the expected amplitude scaling of the catastrophe integral and the wavefunction. However, a quick calculation involving the prefactor of the wavefunction which ensures that particle number is conserved shows that which exactly cancels the Arnol'd scaling. This is a peculiarity of our non-generic initial condition of starting with a completely localized initial state: when tracking a particular fringe, it will move towards the origin but this normalization factor means that its height does not scale with J. Panel (d) of Fig. 4 shows the predictions in the TFIM for the period T of oscillations near the caustic. Data is shown both for the exact wavefunction, given in Eq. (11), and also the 'spin-flip' state Ψ X , given in (33), which is easier to realize experimentally. Since the Berry index ς for the fold defines scaling perpendicular to the caustic, a geometric factor dependent on v LR must be applied. Numerical agreement to within 3% of Airy scaling given in Table I is found in both cases even for finite-sized systems at finite times.

VII. CORRELATION FUNCTIONS AND HIGHER ORDER CATASTROPHES
Rather than the probability distribution associated with the wavefunction itself, light cones are usually ob- served in correlation functions [17][18][19][20][21]. The equal time site-site correlation function is defined as Because Bogoliubov fermions are conserved, b † n (t) = b m (t) = 0, and the last term vanishes. The remaining piece is where we have used the state vector |Ψ(t) given in Eq. (10). Expressing all the operators in terms of quasimomentum (see Appendix A) we obtain In Fig. 5 we plot G(0, x m , t) on the upper half of the spin chain for two different values of g. It displays the same features as the wavefunction: a light cone, interference fringes, and vortices. In the CA, the equal time site-site correlation function becomes and expanding around the cone boundaries gives where C(x, t) is the same function of x and t as that given in Eq. (22).
Measurements and calculations (based on doublon and holon quasiparticles) on the BH model following a quench also find a product of two Airy functions for G(x n , x m , t) [8,17]. However, referring to Table I, generic dimension 3 singularities (i.e. two spatial coordinates x n and x m , as well as time t) of corank 2 (i.e. 2 integration variables, like in the two-site correlation function) are the elliptic umbilic, and hyperbolic umbilic catastrophes. The elliptic umbilic diffraction catastrophe has been studied by Berry, Nye and Wright [87] via the optics of a triangular water droplet lens, while the hyperbolic umbilic is a direct consequence of the primary coma aberration [86], and has been observed in matter waves using electron microscopy [69]. These catastrophes are generally more complicated than a squared Airy function, however, we note that in a certain plane the hyperbolic umbilic wave catastrophe does indeed reduce to the product of two Airy functions. More precisely, the hyperbolic umbilic wave catastrophe is given by [44] Ψ HU (x, y, z) = λ and when C 3 = 0 this reduces exactly to (30) Thus, both the XY model and the BH model give rise to a non-generic special case.
What physical quantity could the C 3 control parameter represent? Studying the form of Ψ HU given in Eq. (29) we note that C 3 controls the coupling between the s 1 and s 2 variables which in a spin chain correspond to the two quasimomenta k and k . For noninteracting quasiparticles, which is the case for the exactly solvable models considered in this paper, the two quasimomenta are uncoupled and thus C 3 is zero. Furthermore, the particular regime of the BH model where Refs. [8,17] obtained a product of Airy functions also corresponds to the free quasiparticle case. It is therefore clear that C 3 can be used to parameterize quasiparticle-quasiparticle scattering, and we predict that a model with interacting quasiparticles will give rise to light cones that sample hyperbolic umbilic wave catastrophes. This feature could be verified in an experiment where the strength of the coupling is varied for then the scaling along C 3 should go as ς 3 = 1/3. Other quantities, for example the spin-spin correlation function, Σ nm = σ x n σ x m − σ x n σ x m , may also be calculated exactly via the Jordan-Wigner and Bogoliubov transformations, and simplified using Wick's theorem. The functional forms of these quantities in the continuum approximation remain diffraction integrals, and thus will also display universal behaviour corresponding to catastrophes. . We define vortex density as being the total number of vortices that occur within a cone up to the time at which the light cone hits the edge of the system, taking care to normalize for different cone sizes at different values of g. Panel (b): Numerical determination of vortex pair creation times at a fixed point in space as g = gc = 1 is approached. In order to extrapolate to the critical point (inset), 30 data points (g,Jt/ ) are fitted to a quadratic and then differentiated. The resulting slopes are extrapolated to gc using a cubic and the intercept gives νz = 0.9999 ± 0.0004 (standard error on the fit). The range 0.02 ≤ |g − 1| ≤ 0.12 of g was chosen to optimize the proximity to the critical point along with data accuracy, since the wavefunction becomes highly oscillatory as g → 1. Numerical errors are smaller than the symbol sizes.

VIII. VORTICES AND CRITICALITY
As seen in Figs. 3 and 5, and also Fig. 8 in the appendices, we find that light cones contain lattices of vortexantivortex pairs. Vortices form the fine structure of wave catastrophes [86,[88][89][90], and in a continuum are zeros of Ψ where the phase χ ≡ ArgΨ is undefined (takes all values) and has the topological property where C is any closed path which contains a single vortex. On a discrete lattice we can still use such circuits to find vortices, but across lattice sites one must perform a sum instead of integrating, meaning that their spatial position is only known up to the lattice constant: in figures we place the vortices between lattice sites. Further-more, vortices on a lattice need not correspond to nodes or even phase singularities, but to points where the phase difference between adjacent sites is ±π (i.e. phase kinks or dark solitons). Thus, while phase interference regulates the amplitude divergence of ray caustics, the effect of a lattice is to regulate the phase singularities of wave theory. In recent work by some of the authors [80], the regularization of phase singularities by a lattice has been considered in Fock space. Whereas the classical light cone changes smoothly at the QCP [see, e.g., Eq. (19)], there is a sharp minimum in the vortex density, i.e. many vortex-antivortex pairs annihilate, see Fig. 6. In the CA, all vortices except those closest to the central axis annihilate at the QCP, while on a discrete lattice, more off-axis vortices survive but the same trend is observed. At a fixed point in space, the time at which a vortex is first detected increases as one approaches the critical point, becoming infinite in the CA. This diverging time scale τ is related to critical slowing and suggests a connection to the dynamical critical exponent, z. According to the scaling hypothesis of critical phenomena where ξ = |g − g c | −ν is the correlation length and ν is its equilibrium critical exponent. Fig. 5(b) plots τ as found from the wavefunction Eq. (11) as g is tuned to the QCP. By extrapolating the numerical data [ Fig. 5(b) inset] to the critical point we obtain νz = 1 and hence recover the known critical scaling for the 1D TFIM [91,92]. For purposes of clarity, we have only included the set of vortices which annihilate closest to the axis x = 0. Vortices which annihilate farther off-axis also display similar trends, which can be seen in Appendix F, along with further figures which help with visualization of this process. While a more complete understanding of the nature of the vortex-antivortex pairs within the light cone remains a subject of future work, we wish to highlight that their presence and scaling laws provide an interesting link between the predictions of catastrophe theory and universality (in and out of equilibrium). Due to the self-dual nature of the TFIM, qualitative behaviour for g > 1 is identical to that of the wavefunction below the transition with g → 1/g and t → gt.

IX. EXPERIMENTAL REALIZATION: SPIN FLIP STATE
The structural stability of catastrophes explains why they occur so frequently in nature. Apart from the examples given in the Introduction, they can also occur in disordered systems such as at the Anderson transition where an evanescent Airy function occurs [93], and it has also been shown that wave catastrophes have the property of self-healing after being disrupted [94]. There are, therefore, a broad range of initial conditions and spin models which will give rise to caustics in their dynamics.
So far we have used the initial condition of a localized single quasiparticle, as given in Eq. (9). This is a nongeneric initial condition and the reader may question how generic the resulting light cones really are. In fact, all our analysis is stable to perturbations around this initial condition. In particular, a state which is naturally generated in trapped ion experiments where individual ions can be addressed is a spin flip state which starts with all spins polarized in the x direction, except for the central spin, say, which is flipped [21], It is important to realize that physical spins are in general superpositions of multiple quasiparticles and vice versa. We elaborate upon the mathematical details of this point in Appendix E. What we find is that as long as the quench is not too close to the transition the number of quasiparticles created by a spin flip is close to one and hence we are perturbing around the single quasiparticle state given in Eq. (9). The evidence for this statement can be found in Figs. 4(d) and 6(a), which compare the results of using Ψ X with those of Ψ. We find that the scaling properties are essentially identical in the two cases whilst the behavior of the vortex density shows some finite differences but is qualitatively the same.

X. DISCUSSION AND CONCLUSIONS
Caustics are a natural phenomenon that can be seen by looking up in the sky on a rainy day. The primary bow of a rainbow is a fold caustic and careful observation reveals supernumerary arcs that are interference fringes described by the Airy function. This is the first in a hierarchy of caustics of increasing complexity whose underlying description is via catastrophe theory. This hierarchy has previously been explored in optics (particularly in the field of gravitational lensing [41]), thermodynamics [95,96], laser physics [97,98], hydrodynamics [43,45,99] and also cosmology [100,101]. By showing that light cones in many-body systems are also caustics, we are able to open the door to the application of a rigorous and unified mathematical framework for describing the dynamics of these systems following a quench.
The main conceptual result of this paper is that there is a hierarchy of light cone structures. They are stable against perturbations and dressed by characteristic wavefunctions that scale according to the sets of exponents given in Table I. The fold catastrophe and its attendant Airy function features in the TFIM, but breaking the symmetry of the TFIM leads us to the XY model and the second catastrophe, the cusp, which is dressed by the lesser-known Pearcey function. Choosing the spin coupling J as a tuning parameter, we show how the scaling exponents lead to non-trivial scaling of these wave catastrophes as J is varied.
The TFIM and XY models are exactly solvable and hence their quasiparticles are noninteracting. However, the defining feature of catastrophe theory is that it deals with structurally stable singularities and hence the light cone caustics we have described also occur in the presence of perturbations such as weak quasiparticle interactions. A related example of this is provided by the celebrated Kolmogorov-Arnold-Moser (KAM) theorem which shows that tori in the phase space of integrable systems are stable against nonintegrable perturbations. There is in fact a close connection between caustics and the quasiperiodic motion that arises in dynamical systems due to the existence of the tori [102].
Higher order catastrophes will become important in higher dimension spin lattices. Another way that higher order catastrophes become important is through n-body correlation functions. For the TFIM we find that the twosite equal time correlation function is described near the cone edge by the product of two Airy functions, which is, however, a special case of the hyperbolic umbilic catastrophe. We predict that adding quasiparticle interactions will lead to the full hyperbolic umbilic catastrophe.
On their finest scales, wave catastrophes contain vortex-antivortex pairs. We have seen that in the case of light cones in 1D spin chains these become vortexantivortex pairs in space-time. We note in passing that these are reminiscent of the Kosterlitz-Thouless transition that occurs in one space and one time dimension in the quenched 1D Bose-Hubbard model [103] and in quantum wires [104]. Being high energy features, we find that the vortices are strongly affected by critical slowing near a QCP, unlike the light cone itself which evolves smoothly. The vortices contain all the information about the QCP and can be used to extract the critical scaling behavior.
The fact that light cones are structurally stable and fall into distinct classes, each of which has its own set of scaling exponents, underlines that as a phenomenon they are an example of universality in out-of-equilibrium dynamics, somewhat akin to the universality classes of equilibrium phase transitions. The underlying reason for this universality in both cases is the presence of singularities, and the realization that light cones are caustics aids us in identifying and understanding their properties.

ACKNOWLEDGMENTS
We are grateful to Marc Cheneau for first pointing out to us the existence of Airy functions in light cones and to Laurent Sanchez-Palencia for discussions. We thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for funding.

Appendix A: Dynamics of a Bogoliubov Fermion
The spin models dealt with in this paper can be exactly diagonalized in terms of Bogoliubov fermions. Their Hamiltonians can therefore be written in the form where k is the dispersion relation and the operatorsb † k andb k create and annihilate, respectively, fermions with quasimomentum k. We shall denote the action of the creation operator on the Bogoliubov vacuum asb † k |0 b = |k b . These operators are related to their counterparts in position space via a discrete Fourier transform: where N is the number of sites/spins. Applying the time evolution operator to a single Bogoliubov fermion created at the center of the lattice we obtain the state vector: where θ(t) ≡ (t/2 ) k k . The corresponding spatial wavefunction is and inserting the standard result x|k = e ikx / √ N for the overlap gives If we allow ∆k = 2π/(aN ) to become very small (N >> 1) we can approximate the sum by the integral with generating function Φ = kx − t k . A comparison of the discrete and continuum cases for the TFIM is given in Fig. 8. In the semiclassical regime (N >> 1) both the sum and the integral are dominated by the points at which Φ is stationary. Along the caustic, however, a saddle-point approximation fails since we are at a degenerate stationary point.
The Hamiltonian for the XY model is, where σ α i , α ∈ {x, y, z} are the Pauli operators for the ith site. We will use the Jordan-Wigner (JW) transformation, followed by a Bogoliubov rotation, in order to diagonalize H. Following the conventions used by Dutta et al. in Ref. [92], the transformation to JW fermions is given by, We note that the JW fermions and Bogoliubov fermions have different vacuua; some more discussion of this point can be found in Appendix E. Next, we use a Fourier transform,c † k = j e ikxj c † j , and then rotate to Bogoliubov fermions viã along with the corresponding destruction operator and transformations for −k.
Here, x + J 2 y + 2h(J x + J y ) cos(ka) + 2J x J y cos(2ka) being a function of the parameters J x , J y , h, and a.
Next we introduce the anisotropy parameter γ so that we can write J x ≡ J(1 + γ)/2, J y ≡ J(1 − γ)/2 and let h ≡ gJ. We thereby arrive at the standard form of the Hamiltonian with dispersion k = 2J (cos(ka) + g) 2 + γ 2 sin 2 (ka). If we change our conventions in order to be consistent with Sachdev [82] we must rotate the Hamiltonian by taking σ x → σ x , σ y → σ y , and σ z → −σ z . Then we'll instead have Effectively this is like taking g → −g, allowing us to return to the standard form of the transverse-field Ising model in the γ → 1 limit, presented in the main text. The dispersion relation given in Eq. (B6) is plotted in Fig. 7. In this appendix we give more details of the calculations of the caustics and their wavefunctions that are presented in the main text. The XY model contains both fold and cusp catastrophes; we focus particularly on the cusp catastrophe and defer some of the treatment of the fold catastrophe to the next appendix (Appendix D) which is on the TFIM.
In the next section we describe how the triple coalescence of stationary points give rise to the Pearcey function which provides the inner cone in Fig. 2. The three stationary points which coalesce are those between the dashed lines in Fig. 7. For 0 < γ < 1 and 0 < g < 1, this coalescence occurs at k = 0, thus z = 1, and Eq. (C6) yields solutions g = 1 and g = 1 − γ 2 . The g = 1 solution is highly singular for nonzero anisotropy, while the solution g = 1 − γ 2 is the key for triple root coalescence.

Diffraction integral for the cusp wave catastrophe
Let us begin by defining the Pearcey function which is the canonical form of the wavefunction corresponding to the cusp catastrophe. The definition of the Pearcey function that we use is It features two parameters C 1 and C 2 and is generally a complex function. In fact, the common definition of the Pearcey function is the complex conjugate of (C7), however for our purposes the above definition is more convenient.
Since the coalescence of extrema in Φ occurs at k = 0, we expand to fourth-order, and factor out J/ which we will later use for scaling, where we have defined the following parameter, (C9) Note that the solution g = 1−γ 2 will kill off the quadratic piece.
We now rescale our integration variable, then our wavefunction locally takes the form, with definitions, and integration limit, Eq. (C11) shows that the wavefunction for the inner cone can locally be expressed as a diffraction integral which is generated by the cusp catastrophe Φ 2 = C 1 s + C 2 s 2 /2 + s 4 , and is thus directly related to the canonical Pearcey function when t is reasonably large and J/ = 1 (below we will see that we can choose any value of J/ and it will simply rescale the coordinates). Note, however, that the normalization of the wavefunction restricts the bounds of the integral as t → 0, and so no true cusp point can occur at the origin since Φ also vanishes there. Nevertheless, the region of integration is proportional to t 1/4 and so is larger than the separation between the stationary points as t → 0 since for any quartic equation of the form Φ 2 the position of the stationary points in the s coordinate is proportional to √ C 2 so that for any infinitesimal time dt the separation between them is proportional to only (dt) 1/2 . Thus it becomes imperative that we consider the effects of all three stationary points, giving rise to the Pearcey-like function described in Eq. (20).
Finally, in order to keep our expressions consistent for |g| < 1 and |g| > 1, we can instead factor out Jg/ overall. The above results are then identical up to a factor of 1/g, which can be absorbed into s and is irrelevant for the scaling. Thus, the expression v I given in Eq. (19) may be used in Eq. (C12) generally.
ds . (C17) These are the scaling factors for the cusp wave catastrophe as listed in Table I. As we tune J, it is convenient to keep the caustic in the same place. This is done by simultaneously tuning a such that the Ising velocity v I is constant.
Appendix D: Caustics in the Transverse-Field Ising Model As mentioned in the main text, the outer light cone in the XY model is dominated by its Airy-like behaviour because it arises from the coalescence of just two stationary points. Since this also occurs in the simpler TFIM (which is obtained by setting γ = 1) we focus on this case here. and which must be simultaneously fulfilled. Rearranging Eq. (D2), g(1 − cos 2 (ka)) = cos(ka)(g 2 − 2g cos(ka) + 1) (D3) leading to cos(ka) = g or cos(ka) = 1/g, as expected. Inputting this into Eq. (D1), along with sin(ka) = 1 − g 2 (or sin(ka) = 1 − 1/g 2 for g > 1), we can solve for the LR velocity, which is identical to the Ising velocity we defined in the previous section, Although the caustic lines are determined by the real solutions to Eq. (D3), there exist imaginary solutions for which the Lieb-Robinson velocity g designations are reversed. This seems to be responsible for lines of constant phase across the caustic (see Fig. 8). The presence of two separate speeds within the light cone is also demonstrated by Cevolani et al. in Ref. [16]. We term these imaginary solutions as 'imaginary caustics'.

Diffraction integral for the fold wave catastrophe
The canonical wave catastrophe corresponding to the fold catastrophe is the Airy function. The definition of the Airy function that we use is It features a single parameter C and is a real function if C is real. The stationary points of Φ coalesce when k = (1/a) arccos(g) for g < 1 and k = (1/a) arccos(1/g) for g > 1, respectively. Thus, for each of these cases, we will expand about these particular k values to third order and factoring out J/ overall, for g < 1 and, for g > 1. Of course, the expansion will only capture the behaviour of the wavefunction close to the light-cone, however this is our primary objective. Furthermore, we are guaranteed that (up to a smooth change of variables) this cubic form in particular is structurally stable and will capture the qualitative features of Φ. We now rescale our integration variables as Thus, for g < 1 and, Φ Ai (s; J) = J −2t g 2 − 1 + 2x v I arccos(1/g) for g > 1. Now we define the control variable as and limits, s Min = (gt) 1 3 (−π − arccos(g)) g < 1 t 1 3 (−π − arccos(1/g)) g > 1 (D16) s Max = (gt) 1 3 (π − arccos(g)) g < 1 t 1 3 (π − arccos(1/g)) g > 1 . (D17) Note that s Min < 0 and s Max > 0. Thus, if we assume long enough times, then it is reasonable to take these integration limits to plus and minus infinity. We now have a description of the wavefunction local to the light cone using a fold catastrophe integral, which in the limit of J/ → 1 will become the Airy integral. (E7) Next, we evolve in time using a more convenient representation of the time-evolution operator, (E8) Dropping the global phase factor, and projecting this state onto real space Ψ(x i , t) = x i |Ψ(t) using we arrive at, after a fair amount of algebra, FIG. 9. Graphic depicting vortex annihilation in the CA for the TFIM. Here, g = 1.75, and for clarity the lines x/a = 1.5 and x/a = 2.5 have been drawn. As g is tuned toward the transition, vortex-antivortex pairs (black dots) will approach one another and eventually annihilate at a particular point in space-time, denoted with an 'X'. It is the vortices which annihilate close to x/a = 1.5 that we refer to as 'primary' (red) and those which annihilate close to x/a = 2.5 we refer to as 'secondary' (blue, all annihilated in this image). In principle, there exist rows of vortices beyond these, but here we focus on those closer to the centre of the lattice and short times.

Appendix F: Vortex Scaling
Returning to our original initial condition of a single Bogoliubov fermion created at x = 0, we can identify space-time vortices in the time evolved system. The discrete (exact) and CA results are compared for the TFIM in Fig. 8 where the same general trend is observed in both cases: fewer vortices at the QCP at g = 1 than away from it at g = 0.5. The vortices that survive at the QCP are those near to the centre of the chain at x = 0, i.e. those closest to the position of the original excitation. In fact, in the CA only a single line of vortices on each side of the centre line survives.
The vortices are located by breaking the light cone up into small loops and integrating the phase of the wavefunction around each one. For a loop containing a single vortex, where the plus sign signifies a vortex and the minus sign an antivortex. For the discrete wavefunction the integral along the spatial part of the path C is replaced by a sum. If we track the positions of the vortices as g is varied we find that they flow in space-time in such a way that as the QCP is approached vortices and antivortices annihilate in pairs, each pair annihilating at a different point (x, t). This process is easier to follow in the CA than the discrete case because the discreteness in the lattice direction obscures the spatial location of vortices, so in this Appendix we specialize to the CA case (whereas the data presented in Fig. 6 in the main text are for the discrete The time at which vortex annihilation occurs along a particular set of vortices will diverge as we approach the QCP. Panel (b): Each consecutive vortex pair will annihilate at a point in space (x) which approaches the midpoint between two lattice sites. Thusx = 1.5a for the set of primary vortices, andx = 2.5a for the secondary vortices. case). In particular, Fig. 9 gives a pictorial representation of the annihilations occurring near the centre of the lattice for g = 1.75. We see that vortex-antivortex pairs converge on horizontal lines (i.e. spatial points) located at x/a = ±0.5, ±1.5, ±2.5, . . ..
The temporal behavior of the vortices can also be seen in Fig. 9. For values of g close the QCP the vortexantivortex pairs that occur at short times annihilate and so never occur, or, said another way, as g → 1 the creation time for vortex-antivortex pairs diverges, an example of critical slowing. Thus, there are two dimensions along which one can observe critical scaling: along t and along x, and the data for these two directions are shown in Fig. 10. It is clear from the way that the data falls onto straight lines on a log-log scale as g → 1 that the vortices display critical scaling. The figure shows two 'sets' of vortices, where each set annihilates within a small region of x/a at diverging time scales. The vortices we call primary vortices annihilate at positions approachingx = 1.5a, whilex = 2.5a for the secondary vortices. In the main text, we focus only on the primary vortices, since a greater number annihilate earlier in time and thus result in a less oscillatory integrand, allowing us to get closer to the transition while maintaining accuracy for a larger number of data points, but we see that the sec-ondary vortices obey the same scaling. The temporal scaling shown in Fig. 10(a) leads to a gradient of −1 and hence the relation νz = 1 as explained in the main text [see inset in Fig. 6(b)]. The spatial scaling is shown in Fig. 10(b) and leads to a gradient of 0.5.