Dynamical quantum phase transitions from random matrix theory

We uncover a novel dynamical quantum phase transition, using random matrix theory and its associated notion of planar limit. We study it for the isotropic XY Heisenberg spin chain. For this, we probe its real-time dynamics through the Loschmidt echo. This leads to the study of a random matrix ensemble with a complex weight, whose analysis requires novel technical considerations, that we develop. We obtain three main results: 1) There is a third order phase transition at a rescaled critical time, that we determine. 2) The third order phase transition persists away from the thermodynamic limit. 3) For times below the critical value, the difference between the thermodynamic limit and a finite chain decreases exponentially with the system size. All these results depend in a rich manner on the parity of the number of flipped spins of the quantum state conforming the fidelity.

New phases have concomitant new types of phase transitions.Among these, dynamical quantum phase transitions (DQPTs)  (for reviews [48][49][50][51]) are still far from being well understood, partly because numerical simulations of time evolution in interacting quantum systems are extremely challenging, and partly because there are very few examples for which one can rigorously prove the existence of a DQPT.The first problem starts to be resolved thanks to the outstanding advance of quantum simulators, either digital or analogue, which can accurately simulate the long time evolution of some target Hamiltonian.Indeed, quantum simulations of DQPT have been already reported for instance in [52][53][54][55][56][57][58][59].
In this work, we contribute to the second problem, by finding and describing a new type of DQPT which differs significantly from all previous examples and, at the same time, can be analyzed rigorously using Random Matrix Theory (RMT).As it is standard in the context of DQPTs, we will consider the Loschmidt echo and the associated dynamical free energy [60][61][62].We will show that, for the isotropic XY model, the dynamical free energy undergoes a third order phase transition in the spirit of, but different than, the celebrated Gross-Witten-Wadia (GWW) phase transition in lattice gauge theory [63][64][65].We will also show that the transition can be already observed in rather small system sizes.

SETUP AND MAIN RESULTS
A probe of dynamical phase transitions is the Loschmidt echo, or fidelity [62].A system of L qubits is prepared in the initial state |ψ 0 and evolved forward in time through a Hamiltonian H.One is usually interested in the Loschmidt amplitude and the corresponding probability, the Loschmidt echo: A DQPT is signalled by the non-analytic behaviour of the dynamical free energy f (t) ∝ − ln L(t).
In this paper we will focus on the Loschmidt echo defined by the isotropic XY Hamiltonian 1 with periodic boundary conditions (PBC) and initial state 1 Also referred to as XX model and well-known to be equivalent to free fermions [66].Therefore, L denotes the size of the chain and N is the number of flipped spins in the state.The corresponding Loschmidt amplitude is the echo is L N (t) = |G N (t)| 2 , and we define the dynamical free energy as where t and τ are related through scaling: Notice that f N (τ ) is normalized by N 2 , not by the total number L of qubits, hence it significantly differs from previous approaches [15,67].The Loschmidt return rate − ln L N (t) is not an extensive quantity here and remains finite in the thermodynamic limit (see Methods).We will define f (τ ) as the result of taking first the thermodynamic limit L → ∞ and then the limit N → ∞, when we consider only even values of N to compute the limit.
Our main results are: Main result 1.There is a third order dynamical phase transition in f (τ ) occurring at τ cr ≈ 0.33137.
Main result 2. The third order phase transition persists in the thermodynamic limit L, N → ∞ with t = N τ and constant ratio L/N > 1.2.
Main result 3.For τ < τ cr the difference in f N (τ ) between considering the thermodynamic limit and a chain of length L decreases exponentially with L−N if L > 2N .
Moreover, we prove that the phase transition is robust against the presence of impurities in the preparation of the initial state (4).
We stress that, whilst main results 1 & 2 require N → ∞ taking only even values, main result 3 also holds for the odd case.This will be shown in detail below.
All our main results are mathematically rigorous, and are proven using RMT.Dealing with real-time evolution requires the development of new techniques, which constitutes the main technical contribution of this paper.
The connection between the XY model and RMT is well understood, but so far it mostly dealt with static quantities, like thermal states [68][69][70][71][72][73].The Loschmidt echo (and generalizations thereof) has of course been considered previously.However, if, on one hand, the RMT formulation of (5) (see Methods) was thoroughly studied in [74], the scaling limit (7) remained insofar unaddressed.On the other hand, the vast majority of previous works relied on analytic continuation to the real-time echo from its thermal version [70,75,76].As a byproduct of our main result 1, this approach is only valid for τ < τ cr .Note also that a breakdown of the analytic continuation was already demonstrated in [77,78].
For these reasons, the phase transition eluded prior scrutiny.

Differences with standard DPQTs
The new dynamical phase transition presented here has many exotic characteristics.
The first one is that the DQPT emerges from the scaling limit (7), in analogy with the planar limit of Quantum Field Theory (QFT) and RMT [79].The transition reflects a non-analytic change in the limit f (τ ) of the dynamical free energy (6) as a function of τ (Fig. 1).
The second one is that the phase transition is third order.This is in principle unusual from the perspective of condensed matter physics, where most phase transitions are first or second order [80], and this seems to be the setting of previously reported DQPTs [48,81,82].
Third order phase transitions have been found in the study of integrable hyperbolic Richardson-Gaudin models, relevant for both superfluidity and superconductivity [83,84].See [85] for earlier and [86,87] (and references therein) for recent discussions of third-order transitions.A common thread is the lack of experimental realization so far, contrarily to the outcome of our work.
In RMT, third-order phase transitions are widespread.The first and best known example is the GWW transition [63][64][65], and a number of third-order transitions are indeed reduced to it [88].The GWW model was originally introduced in lattice Yang-Mills theory in two dimensions [89] and has since then found overarching multidisciplinary applications.It also has played a pivotal role in solving the long-standing problem of the longest increasing subsequence in combinatorics [90,91].
The Loschmidt amplitude (5) is intimately related to the GWW model: the latter may be seen as the imaginary-time version of (5).Notwithstanding this parallel, the analytic form of f (τ ) as well as the mechanism beyond the DQPT will make clear that the third order transition we unveil is novel and radically distinct from the GWW one.
The third main difference is about how the DQPT we present here is related to the zeros of the Loschmidt echo.DQPTs [48] are usually associated to zeros in the Loschmidt echo which are attained in the thermodynamic limit.Here, as mentioned above, the Loschmidt echo is strictly positive and analytic for every even N in the thermodynamic limit.The non-analyticity arises only when taking N to infinity.
The case of odd N We will also analyze the case of odd N (meaning that we take N → ∞ considering only odd values of N ).This case has an extreme behavior regarding the relationship to zeros of the Loschmidt amplitude in stark contrast to the case of even N : the Loschmidt amplitude becomes zero at finite L, N and t -a behaviour already observed in [74].Such zeros, which translate into divergences of the dynamical free energy, are responsible for the phase transition in the odd case (Fig. 2).This implies a difference in how the phase transition is approached if one takes N even or odd, see Fig. 3.In the even case, L N (t) decays with t and the dynamical free energy f (τ ), after the phase transition, is smaller than the expected curve if there were no phase transition.In the odd case however, f (τ ) is larger and L N (t) oscillates and has zeros.It is remarkable that, despite this mathematical fact, we can prove the following Main result 4. τ cr is the same for the even and odd cases.

Experimental accessibility of the phase transition
Thanks to our main result 3, one can observe clear signatures of the phase transition at very low values of N and L. For instance, as one can see in Fig. 3, there is no difference in taking L = 9 or the thermodynamic limit for N = 3.In that case, L 3 (t = 1.29) ≈ 0.034 and one already observes a 4.5% deviation from the expected curve if there were no phase transition.Note that L 3 (t) is the probability value of returning to the state |↓↓↓↑↑ . . .↑ at time t.This can be estimated by measuring in the basis {|↓ , |↑ } each one of the spins.
By making a few million samples, the Hoeffding bound [92] guarantees the detection of such deviation with very large probability.This is the number of samples used for each data point for instance in the celebrated Google's quantum supremacy experiment [93] using a quantum processor based on superconducting circuits [94].
For N = 4 and L = 11 (see Fig. 3), there is also no difference with the thermodynamic limit and the deviation for the expected curve increases to 18%, but the value L 4 (t = 1.80) ≈ 0.001 is smaller.In this case the required number of samples is one order of magnitude larger (∼ 10 7 ), which is still within experimental reach.
Since both the XY interaction and the required initialization and measurement are very naturally implementable in superconducting circuits [94,95], these platforms are precisely the most suitable candidate for an experimental observation of this new DQPT.Indeed, other DQPTs have been already observed experimentally in such settings by measuring the Loschmidt echo [96].
Recent improvements in optical lattice microscopes have allowed to study very accurately time evolution in fermionic chains after a quench in cold atoms experiments [97,98].Since the XY interaction is also naturally implementable [99] with these microscopes, one may also expect to observe the DQPT we uncover here in optical lattices.Furthermore, the implementation of a DQPT via cold atoms in a optical cavity has been reported in [100].

METHODS
The map from the XY Loschmidt amplitude (5) with PBC to a matrix model is central in our discussion [68,101] It follows from an old and well-known result in RMT [102] that the integral (8) equals the determinant of a N × N Toeplitz matrix: where J ν are the Bessel functions of the first kind.The matrix model corresponding to a finite chain is presented in Ref. 103, and amounts to replace the integral (8) with a Riemann sum.It can be equivalently written as a Toeplitz determinant, as well.
Main result 1: Dynamical phase transition in the even N case We first consider a chain with N even and L → ∞.Corrections due to finite system size are under analytic control and are discussed below.
To set up the large N limit, express (8) in the form It is customary to introduce the density of eigenvalues ρ(θ) = 1 N N a=1 δ(θ − θ a ) and rewrite S(θ 1 , . . ., θ N ) in the functional form where − is the principal value integral.For the two pieces to be same order in N , we have replaced t/N with τ according to (7).This leads to the so-called planar limit of the matrix model [79].
A novel feature of ( 12) is the presence of an imaginary coefficient.As a consequence, (14) does not admit a standard solution with z ∈ U (1).To solve the complexified minimization problem ( 14) is one of the main technical achievements of this work.
We relax the assumption on supp( ), allowing it to be a curve Γ ⊂ C subject to the requirement lim τ →0 Γ = U (1).The details are in [103], and the upshot is that Γ is the connected component of with appropriate initial condition.Letting the system evolve in time, the dynamical free energy is non-analytic when a complex zero of (z) hits Γ.Hence, there exists a critical time τ cr (see [103]), which is the unique positive solution to i.e. τ cr ≈ 0.33137171.Solving ( 14) for (z) and using (15), we evaluate Moreover, for τ > τ cr we can extract f (τ ) it in the two regimes τ → τ cr and τ → ∞.In the former limit, implying that the DQPT is third order.Conversely, for τ → ∞ we get f (τ ) ∝ ln τ (in agreement with Fig. 1).The late time asymptotic behaviour of f (τ ), as well as the critical value τ cr , guarantee that the DQPT we uncover is not merely a Wick rotation of the GWW phase transition.

Main results 2 & 3: Finite chain effects
A finite chain length L < ∞ replaces (8) with a discrete ensemble, in which the N eigenvalues are distributed among L sites.We refine the large N analysis for the discrete ensemble and consider the scaled thermodynamic limit L, N, t → ∞ with fixed rates τ = t/N and We observe that a finite ratio 1 < < ∞ induces a DQPT at a certain time τ ( ) [103], which is an explicitly known function of .Requiring that this transition takes place at a later time compared to the DQPT discussed insofar, and thus does not shield it, we impose τ ( ) > τ cr .The inequality is satisfied for 1.2.Furthermore, we quantify the validity of the L = ∞ approximation.We prove that, in the scaled thermodynamic limit (19) and in the phase τ < τ cr , the error in neglecting the dependence on < ∞ is exponentially small in L − N if ≥ 2. The derivation (detailed in [103]) builds on the mathematical work [104] and on the theory of orthogonal polynomials.
We also complement the analytic estimate with quantitative numerical evidence.Define the error The closer E N is to zero, the better is the infinite chain approximation.Fig. 4 shows that (20) converges quickly to 0 as is increased, and that the L = ∞ approximation is more accurate away from the critical point.1. × 10 -7 1.5 × 10 -7 2. × 10 -7 2.5 × 10 -7 ℰ6(ℓ, τ) We can moreover analyze the combined effect of finite L and N versus our asymptotic analysis.A quantitative probe of such net effect is where f (resp.f N,L ) is the dynamical free energy (6) in the large N limit at L = ∞ computed analytically (resp.the dynamical free energy at finite N and 2N < L < ∞).
The quantity (21) is shown in Fig. 5-6.We infer from the plots that f N,L (τ ) approaches f (τ ) monotonically from below, and that the discrepancy between the finite chain results and the asymptotic analysis are almost entirely due to finite N effects, compared to which the < ∞ effects are negligible if > 2. These conclusions agree with the fully analytic main result 3. Main result 4: Odd N and quantum speed limit Quantum speed limits (QSLs) pinpoint underlying time scales of physical processes by producing lower bounds on the expectation value of an observable or on the rate of change of a quantum state [105][106][107] (for a review [108]).We have seen above that, for every odd N , L N (t) oscillates with damped amplitude and has zeros, as shown in Fig. 2. Thus, we focus here on the earliest time at which an evolved state |ψ(t) = e itH |ψ 0 is orthogonal to the initial state |ψ 0 [109].In the current model, let us denote the corresponding value of τ as τ QSL

N
, that is where we recall that t and τ are related through (7).By construction, τ QSL N is intimately related to the zeros of the Loschmidt echo [110].
We show the existence of a finite limiting value τ QSL of the QSL, namely The proof relies on the existence, for τ < τ QSL N , of a system of orthogonal polynomials on the unit circle, that are naturally associated to the Toeplitz determinant (9).The lower bound τ QSL > 0 can be shown applying Szegő's strong limit theorem [111], which is the classical result for the asymptotics of a wide class of Toeplitz determinants, including (9), complemented with analytic arguments.
Moreover, we obtain that the limit value is precisely which is the content of our main result 4. The proof of ( 24) is technical and we defer it to Ref. 103.It uses random matrix theory, exploiting also a connection with the Toda lattice [112].
Evidence for the identity (24) is given in Fig. 7, which shows the first 27 values of τ QSL 2k+1 computed numerically.

DISCUSSION
The physics of coherent nonequilibrium real-time evolution is very rich yet still vastly unexplored.In this work we have uncovered a novel dynamical quantum phase transition, with many distinctive features.Such DQPT is detected by the real-time Loschmidt echo of the XY model but, differently from previous works, we have analyzed it in the planar limit, guided by QFT and RMT.
The DQPT takes place as a function of time in the scaling limit (7) and its characterization has required the study of a random matrix ensemble with complex weight, for which we developed tailored analytical methods [103].In particular, the DQPT is not merely a Wick rotation of the well-established GWW transition.
The DQPT arises in the thermodynamic limit, but we have managed to fully identify the finite chain effects and gain analytic control on the exponentially small discrepancies.
While the above points hold true for N large and even, we have also discussed a phase transitions induced by the zeros of the Loschmidt echo when N is odd and have established a relation with the quantum speed limit.Every analytical result has been supported by numerical evaluations, shown in the Figures.
Third order DQPTs have never been measured experimentally.We have shown here how not only it is within experimental reach, but it is predicted for one of the most fundamental strongly-correlated systems, as is the isotropic XY Heisenberg chain.Crucially, the spin interactions of this model have already been engineered with quantum simulators [55].
The experimental verification of the DQPT discovered here is certainly the most appealing potential follow-up of our wok.Another problem open for future study is to establish the universality class of the transition.Based on heuristic arguments, we expect that it is governed by Painlevé I equation.This would mean the universality class of 2d quantum gravity [113][114][115], as opposed to that of GWW, characterized by Painlevé II [116].It would be interesting to rigorously confirm or refute this conjecture.
From a broader perspective, in this paper we have opened a new route to analyze rigorously DQPTs and to uncover their universality classes, based on the power of RMT.As the natural next step, inspired by the nonanalyticity result of [78], it is conceivable, although technically challenging, that an exhaustively characterization of DQPT as done presently can be achieved in more general 1d models.Consider a one-dimensional spin chain with L qubits, L 1.We label the qubits k = 1, . . ., L and use the notation {|↑ k , |↓ k } for the natural basis in the k th copy of the single-qubit Hilbert space H .The Hilbert space of the spin chain is Span {|↑ j , |↓ j } , the L th tensor product of the single-qubit Hilbert space, with dim H ⊗L = 2 L .
The spin chain is prepared in the initial state It is possible to assemble N ≤ L spin-flip operators acting on N distinct copies of the single-qubit Hilbert space into string operators k i = k j if i = j (and likewise for σ + k1•••kN ).Thanks to the commutation relations, we can restrict without loss of generality to configurations with Acting with the operators (S1.1) on |⇑ produces the states These states are labelled by the integer N ∈ N and a partition of length at most N .(The dependence on L is left implicit).The partition κ is related to the labels {k j } j=1,...,N of the qubits acted upon through Notice that it is important to specify N , not only κ, because those partitions with length(κ) < N give distinct operators (S1.1) for distinct values of N .The states |N, κ are multi-domain wall states, consisting of strings of ↑ alternating with strings of ↓.In particular, the empty partition κ = ∅ corresponds to the state which has a single domain wall, for a chain with open boundary conditions (OBC), or a two-sided domain wall for a chain with periodic boundary conditions (PBC).

B. Quantum quench
We consider a quantum quench protocol in which we begin with a state |N, κ and a trivial Hamiltonian and, at t = 0, we suddenly turn on the spin- 1  2 Heisenberg XX Hamiltonian We are interested in the return probability, that goes under the name of Loschmidt echo, or fidelity: 3) The corresponding probability amplitude, the Loschmidt amplitude, is Let us recall that, in general, the Loschmidt echo measures, by definition, the overlap of a state e itH0 |ψ , evolved forward in time with an unperturbed Hamiltonian H 0 , with a state e itH |ψ evolved with a perturbed Hamiltonian H. See [S1] for a review.Here we focus on one of the simplest possible setups: the initial Hamiltonian is trivial, H 0 = 0, the perturbed Hamiltonian is H = H XY and the state is (S1.2).
Assuming an infinitely long chain, L → ∞, the Loschmidt amplitude is given by [S2] where χ κ is the character of the representation κ of G, e.g. a Schur polynomial if G = SU (N ), and dU is the normalized Haar measure on G, such that The derivation of (S1.5) is based on the observation that the Loschmidt amplitude (S1.4) satisfies the differential equation In the latter expression, κ ± a means the partition obtained by appending or removing a box to the a th row of the Young tableaux representing κ.Equivalently, For instance, A solution to (S1.6) at imaginary time t = −iβ was given in [S3, S4] in the form of a determinant.The derivation is algebraic and does not rely on reality of the parameters, thus it applies straightforwardly to the case at hand.The specific form of the determinant solving (S1.6) is sensitive to the boundary conditions, so let us pause to discuss these first.We will come back and complete the proof of (S1.5) afterwards.
In this work, we consider two types of boundary conditions: i) Periodic boundary conditions (PBC for short); ii) Absorbing boundary conditions on the left and Open boundary conditions on the right (henceforth ABC for short).
By absorbing boundary conditions, we mean that we set to zero the probability amplitude of a hop of the 0 th qubit to its left [S3].We may also allow for open boundary conditions (OBC) on both sides of the chain.As we will explain below, neglecting this choice results in no loss of generality for what concerns our results.
It follows from [S2] that G N,L (t; κ) solving (S1.6) with PBC is a Toeplitz determinant.For ABC, one gets instead a Toeplitz+Hankel determinant.We do not write down their form explicitly because we will not use it in the rest of the work.More details on the derivation will be given for the κ = ∅ case below.Formula (S1.5) follows from applying the celebrated Andréief's identity [S5, S6] to these determinants.
Lemma S1.1 (Andréief's identity).Let dµ(z) be a Borel measure on D and {Φ a−1 } N a=1 and {Ψ a−1 } N a=1 be two collections of µ-integrable functions on D. Then Equivalently, one could use the generalization of the Heine-Szegő identity in [S7].In each case, we arrive at (S1.5) with Furthermore, it is possible to realize a chain with OBC at both ends by breaking the U (1) translation invariance of the PBC chain, cutting it open.It has the effect of removing the center symmetry (which is a 1-form symmetry, in the gauge theory parlance), thus the integration in (S1.5) is over P SU (N ).

Remark S1.2 (OBC and global form of G).
It is oftentimes quoted in the literature that, for OBC, the pertinent group is SU (N ).From the analysis of the symmetries of the model we get that the resulting group should be centerless, hence its global form is P SU (N ).This subtlety is not relevant for the ensuing discussion.

C. Preparing the quench in the simplest domain wall state
We focus on the simplest initial state with domain wall, To reduce clutter, we denote The determinants expressing the Loschmidt amplitude with, respectively, PBC and ABC are [S2, S8]: where J ν are the Bessel functions of the first kind.For trivial κ = ∅, the integral (S1.5) becomes simply It is customary in random matrix theory to reduce integrals of the form (S1.8) to integrals over eigenvalues, thanks to Weyl's integration formula.
Lemma S1.3 (Weyl's integration formula).Let G be a compact, connected Lie group, T G a maximal torus and W G the Weyl group of G determined by T G .Denote by dU the normalized Haar measure on G and dt the Lebesgue measure on the torus.Besides, let Φ be a class function on G.
where, denoting Roots + (G) the set of positive roots of G, Weyl's integration formula allows to rewrite (S1.8) as the eigenvalue integrals where the Vandermonde factors are [S9] As claimed in the previous discussion, the equality between (S1.9) and (S1.7) follows from identifying in the Vandermonde determinants a form suitable for applying Lemma S1.1 to (S1.9), and computing the integrals An alternative derivation of the matrix model representation (S1.8) is based on a Jordan-Wigner transformation, see for instance [S10].
Remark S1.4 (Loschmidt echo and gauge theory).The matrix model (S1.9a) for the spin chain with PBC is closely related to the Gross-Witten-Wadia (GWW) model [S11-S13], originally introduced in twodimensional lattice gauge theory.More precisely, the GWW model is the imaginary-time version of (S1.9a), i.e. replacing t → −iβ.The inverse temperature β plays the role of the inverse Yang-Mills coupling in the gauge theory interpretation [S11].The imaginary-time version of (S1.9b) is the analogue of the GWW model with U Sp(2N ) gauge group.This gauge theory analogy will partly inspire the methods to study (S1.9) in later sections.
From the expressions (S1.8) we easily recover the universal, short-time behaviour of the Loschmidt echo.
Theorem S1.5.At leading order in the short time expansion, the Loschmidt echo has the universal behaviour with the average over G = U (N ) for PBC and over G = U sp(2N ) for ABC.The superscript indicates the connected correlator Proof.Expanding the integrand in (S1.8) at small t one gets Recall that we are using the normalized Haar measure on the Lie groups, so G N (0) = 1.The above expansion then implies Observing that concludes the proof.

D. Loschmidt echo on a finite chain
The expressions (S1.9) hold for an infinitely long chain to which we apply the operators (S1.1).More precisely, (S1.9) must be understood as the L → ∞ limit of a chain of fixed length (or circumference) in which the number L of qubits is progressively increased.Note that the actual length of the chain is irrelevant, because it drops out from any probability amplitude by normalization.
The finite-L matrix models, that give rise to (S1.9) upon taking L → ∞, are discrete analogues of (S1.9), in which the N eigenvalues occupy L discrete sites.The angular variables {θ a } a=1,...,N are discretized into {θ(s a )} a=1,...,N , where In a finite chain of L qubits with PBC, the circular topology U (1) is replaced with the discrete topology where we describe Z L as a the multiplicative group e i2π(s−1)/L , s = 1 . . ., L .
ii) Likewise, in a finite chain of L qubits with ABC, the semicircle topology U (1)/Z P 2 is replaced by the multiplicative group e iπ(s−1)/L , s = 1 . . ., L .Here Z P 2 is the parity symmetry z → −z.
Therefore, the matrix models for L < ∞ are the Riemann sums: The discrete ensembles (S1.12) admit a determinant presentation, as well.
Proposition S1.6.For every k ∈ Z and ξ ∈ C, let e −w cos θ(s) e ikθ(s) . Then Proof.It is a direct application of Andreief's identity (Lemma S1.1) with measure For practical purposes, instead of working with (S1.9b)-(S1.12b),when considering the ABC chain we use the change of variables x a = cos θ a , which produces where the finite-L model has domain The function in the definition of S L comes from the logarithmic form of the arccosine function, and then writing z = e iArcCos(x) .
Remark S1.7 (Finite echo in the thermodynamic limit).We stress that a very important consequence of the simple quench setup we consider is the finitude of the Loschmidt echo in the thermodynamic limit L → ∞.This is in sharp contrast with most quench protocols, in which one typically has ln L(t) ≈ −L f (t) at L 1, for some L-independent function f [S14].

E. Determinants of Bessel functions
Looking back at the expressions (S1.7) and plotting them as functions of t for several values of N ∈ N, we uncover the following pattern: • For N ∈ 2N + 1, L N,L (t) decays as e −t 2 near t 0, then vanishes at a given t = t QSL N and shows damped oscillations from that point on; • For N ∈ 2N, L N,L (t) decays as e −t 2 near t 0, then reduces its slope and has an asymptotic decay to zero with small oscillations on top of that trend.
(See also [S15], where such behaviour was first noticed).We stress that the behaviour of L N,L=∞ (t), including the oscillations for N odd, are a consequence of the analytic properties of the Bessel functions J ν (2t): • Exponentially decreasing near t 0; • Oscillating at t 1.
For the most part of the body of this note, we consider N ∈ 2N, that is, an even number of spins |↓ .The discussion of N ∈ 2N + 1 is deferred to Section S5.
As a side remark, we notice that the factors of i a−b (where i = √ −1) in (S1.7a) appear in such a way that ∀ N ∈ N, that is, they drop out of the determinant, cf. Figure S1.

S2. DYNAMICAL QUANTUM PHASE TRANSITION
In this section we discuss the large N planar limit of the models (S1.9a)-(S1.14),and prove the existence of a third order dynamical quantum phase transition.
For concreteness, we begin with the analysis of the L = ∞ chain, and come back to the experimentally accessible L < ∞ case in the next section.
The content of this section is organized as follows.
• In Subsection S2 A we introduce the planar limit and set up the problem at large N .
• Then, we take a detour to discuss the large N third order phase transition taking place at imaginary time.In Subsection S2 B we review the large N solution of G N,∞ (−iβ), which is the GWW model, and the associated third order phase transition.This is a review of known results, which we report here for completeness and to introduce the tools we will use in later subsections.After that, in Subsection S2 C we extend the result to the ABC chain in imaginary time, that is, to the large N limit of G (ABC) N,∞ (−iβ).The latter model was recently studied in [S16].
• We then proceed with the analysis of interest, with real time.We discuss the large N planar limit of the chain with ABC in Subsection S2 D and with PBC in Subsection S2 E. We will prove therein that the chain undergoes a third order DQPT.These two subsections form the core of the work and contain many of the main technical achievements.
Because of their intrinsically more technical nature, each of the Subsections S2 B to S2 E is segmented into several steps, to orient the reader through the computations.

Statement of the results
We will uncover a third order dynamical quantum phase transition in the planar large N limit, N, t → ∞ with τ = t/N fixed.The existence of this transition is insensitive to the choice of boundary conditions.Moreover, it will be manifest in the derivation that the phase transition we find is not merely a rotation of the GWW phase transition from imaginary to real time.This is confirmed both by the critical value τ cr = 0.33137 . . .and the explicit form of the dynamical free energy.
Very much in the paradigm of dynamical quantum phase transitions [S14], we observe a change in the analytic properties of the dynamical free energy in the complex plane.

A. Planar limit of matrix models
Let us begin by setting up large N planar limit of the models (S1.9a)-(S1.14)as well as their imaginarytime version.Standard reviews on the large N limit of matrix models include [S17-S19].For unitary matrix models see e.g.[S20-S23].
To treat real and imaginary time at once, throughout the present subsection we define the quantities Ĝ where {c k } k≥1 is a finite collection of real coefficients and {T k } k are the Chebyshev polynomials of the first kind.For later reference, {U k } k will be the Chebyshev polynomials of the second kind.Here w ∈ C is a complex parameter, and the real-and imaginary-time dynamics are recovered substituting real time : Moreover, the original models of interest correspond to the coefficients We consider the planar limit [S24] of the matrix models (S2.17).The planar limit is a large N limit originally introduced by 't Hooft for the study of Quantum Chromodynamics [S25] and has since then proved crucial in a wide variety of problems in high-energy physics and beyond.It owes its name to the fact that only planar Feynman diagrams are retained in this large N limit.
We start by rewriting Ĝ(PBC) where the "effective actions" are They are obtained by passing the Vandermonde determinants in the exponential.Observe that, in the limit N → ∞, every sum over the indices a, b = 1, . . ., N grows linearly in N , thus the terms 1 N N a=1 have a well-defined large N limit.For the two pieces in the effective action to compete at the same order in N , for N 1, one is led to consider the 't Hooft scaling of the parameter w, with For the imaginary-time dynamics, w → β, the planar limit corresponds to send the inverse temperature β → ∞ linearly with N .In other words, the 't Hooft scaling intertwines the large N and low temperature limits, instead of taking them independently.
We conclude that, with the scaling (S2.19), the integrands in (S2.18) are suppressed at large N as e −N 2 Seff , with the effective actions having a finite large N limit.Therefore, in the large N limit, the integrals (S2.18) are dominated by the stationary points of the integrands, which are the saddle points of the effective action.
This leads to the consideration of the system of N saddle point equations: PBC : ∂S These are the Euler-Lagrange equations for (S2.18).Explicitly, we have the saddle point equations PBC : 1 Note the factor of 2 coming out of differentiating the double sum.In the PBC case, this factor of 2 cancels a factor 1 2 from the derivative of the Vandermonde term.Remark S2.2.Note that the square-root term in the ABC effective action is sub-leading at large N al will not play a role in determining the saddle points.As a matter of fact, our results in the planar limit apply to any matrix model where ϕ : [−1, 1] → C (i) has no branch cuts on (−1, 1), and (ii) has coefficients that do not depend on N .
In the large N limit, the discrete index a = 1, . . ., N can be effectively replaced by a continuous index This substitution implements the replacement for arbitrary function F , where the N eigenvalues are grouped into a function θ(a) defined according to The procedure is exactly analogous for the substitution x a → x(a) = x aN in the ABC case.
One then arrives at the saddle point equations ∀ 0 < a ≤ 1: ABC : up to O(1/N ) corrections, that are negligible in the large N limit.The symbol − means the Cauchy principal value integral.
Clearly, it is not handy to deal with a system of N equations in the large N limit.For this reason, it is customary to introduce the eigenvalue densities (we will omit the superscript when no confusion can arise) which are normalized by definition and have compact support at large N .Besides, it follows from the definition together with the discussion above that lim Thus, in practice, we read the above equation as ρ (PBC) (θ)dθ = da and enforce the replacement and likewise for ρ (ABC) (x).
Putting all the pieces together, we are finally led to the saddle point equations These are integral equations to be solved for the eigenvalue densities ρ (PBC) and ρ (ABC) .For the XY Heisenberg spin chain we are interested with, these equations reduce to PBC: The substitutions w = β and w = it give the imaginary-time and real-time dynamics, respectively.Because we are dealing with the planar limit, we introduce the scaled quantities imaginary time : To conclude the computaton, we recall from (S2.18) that, at leading order in the large for either choice of boundary condition.In the latter expression, S eff [ρ] means the effective action evaluated on the eigenvalue density that solves the pertinent equation in (S2.23).Therefore, it follows by construction of the planar limit that lim with the right-hand side finite.
Remark S2.3 (Normalization of the free energy).It follows from this general argument that, in the planar limit, the matrix models we consider behave as ln G(t) ∝ N 2 f (τ ), with f (τ ) a continuous function of τ = t/N , independent of N .Therefore, we normalize the dynamical free energy of the system according to in agreement with the standard normalization in the random matrix theory literature.
With the aim of reducing confusion among imaginary-time and real-time quantities, we define for imaginary tme Therefore which, replacing a = aN and rearranging, yields Sending N → ∞ and b → a using the characterization (S2.22), we conclude that ρ(z) ≥ 0.
Because both matrix models we consider have compact support, the eigenvalue densities ρ (PBC) (θ) and ρ (ABC) (x) must be non-negative definite.Solutions to (S2.24) that do not satisfy such constraint are to be discarded.

B. Review of the GWW third order phase transition
The celebrated Gross-Witten-Wadia (GWW) model [S11-S13] is the imaginary-time version of (S1.9a), equivalently (S2.17a) with w = β and c k = 2δ k,1 .It was originally devised in the study of two-dimensional gauge theory on the lattice.The term exp −βTr(U + U −1 ) was derived from the Wilson lattice action for a one-plaquette model of lattice two-dimensional gauge theory [S11].
In this subsection we review the derivation of its third order phase transition, as a warm up to present the saddle point techniques.The reader familiar with the subject may move directly to the next subsection.
To begin with, we make a change of variables θ a → θ a −π for later convenience, which shifts the integration domain to [−π, π] N and reflects β → −β.According to the general discussion of Subsection S2 A, the planar limit of the model is governed by the saddle point equation where we recall that γ = β/N plays the role of the inverse 't Hooft coupling.

Statement of the result for the eigenvalue density
The solution to (S2.28) is given by [S11, S13] where 1 1 is the indicator function and θ 0 in the phase γ > 1 2 is the positive angle with sin θ0 2 = 1/(2γ).

Solution in the first phase
In order to solve (S2.Using the ansatz (S2.30) we can perform the integrals on the left-hand side and get a generating function for the Fourier coefficients of ρ(θ).Equating with the right-hand side one finds The 0 th Fourier coefficient is not fixed by the saddle point equation, but rather is fixed by normalization together with the ansatz (S2.30).
It is straightforward to extend the solution thus obtained to the saddle point equation (S2.23a), finding The solution (S2.31) is subject to the constraint which is violated around θ = ±π for γ > 1 2 , and around θ = 0 for γ < − 1 2 .Because γ is a scaled inverse temperature, we only consider the physically meaningful region γ > 0. Therefore, at γ > 1  2 , the solution (S2.31)is invalid, indicating a phase transition.We must drop the assumption (S2.30) and look for a different eigenvalue density.

Solution in the second phase
Giving up on the assumption (S2.30), we then make the ansatz that ρ(θ) is supported on (one or more) arc(s) on the unit circle.
To this aim, it is standard procedure to introduce the planar resolvent where we are adopting the standard shorthand notation ω ± (z) := lim ε→0 + ω(e iθ ± iε).
Then, expressing (S2.28) in terms of the exponentiated variable z = e iθ ∈ Γ, yields This is a jump equation for ω(z) along Γ, supplemented with the condition lim |z|→∞ ω(z) = 1, which follows from the normalization of the eigenvalue density together with the definition (S2.32).We have thus set up a Riemann-Hilbert problem, whose solution governs the planar limit of the GWW model in the phase γ > 1 2 .The solution for ω(z) is derived in a standard way.The jump equation (S2.34) implies a contour integral representation for the planar resolvent, with integration contour encircling Γ.The final result then follows from standard contour manipulations, see e.g.[S21-S23] for the details.

Free energy
Equipped with the solution (S2.29) for the eigenvalue density, we can compute F(γ) as defined in (S2.27).More precisely, it is convenient to compute The result is that F (PBC) (γ) is continuous and twice differentiable, but it fails to be smooth at γ = 1 2 .This is the third order GWW phase transition.

C. A GWW transition for the ABC chain in imaginary time
We now extend the above result to the imaginary-time version of the ABC.

Statement of the result
The outcome of the analysis is: the model has a third order transition at the same critical value γ cr = 1 2 as the imaginary-time PBC chain (i.e. the original GWW model).Moreover The result has been recently obtained in [S16] using slightly different methods.

Solution in the first phase
Taking the large N limit after the change of variables x a = cos θ a yields the imaginary-time version of (S2.24b): In fact, we can easily solve the more general problem where U k−1 are the Chebyshev polynomials of second kind.From the Stieltjes transform of the Chebyshev polynomials we infer Specializing to the case of interest, we have Again, the solution holds for 0 ≤ γ ≤ 1 2 .Beyond the critical value, the solution ceases to be non-negative close to the edge x = −1.Therefore, the GWW transition becomes a hard-edge to soft-edge transition in this setup.

Solution in the second phase
We look for a so-called "soft-edge" density, in which the inverse square-root singularity near x = −1 is replaced by a square-root behaviour near a boundary A > −1.
The solution to (S2.36) under these circumstances is known (see e.g.[S26]) and goes under the name of Tricomi's formula [S27].The final answer in the case of interest reads where A = γ−1 γ is fixed by normalization.
In the first phase, the relation (S2.35) is a corollary of Johansson's extension of Szegő strong limit theorem [S28].Our result is to prove it beyond the phase transition.

Theorem S2.6 ([S28, S29]
).Let {c k } k≥1 be a collection of coefficients satisfying and define It holds that lim Proof.The U (N ) part of the theorem dates back to Szegő [S29].For its extension to other classical Lie groups, in particular U Sp(2N ), see [S28, S30], but be aware of the difference in normalization.Szegő's theorem and its extension consider the strict N → ∞ limit, with all other parameters fixed.The agreement with the value in the first phase follows by real-analyticity arguments, see [S22] for an explicit proof.

Folding trick
The spin chain model we are considering endows us with a visual, albeit non-rigorous, argument for the factor of 2 in formula (S2.35).
Imagine to begin with a chain with PBC and 2N + 1 flipped spins, yielding a U (2N + 1) matrix integral.Because of the N 2 -scaling, we have at leading order in N .Then, we remove the middle one among the 2N + 1 qubits and replace it with an absorbing wall.At the same time, we cut open the chain at the antipodal point, thus obtaining two copies of a chain with N flipped spins and ABC on the left, OBC on the right.Because the two effects of introducing boundaries/interfaces are expected to be sub-leading in N , we predict Combining the last two equations we obtain (S2.35), which is of course only valid at leading order in the planar limit.It is indeed known [S28], and in agreement with our argument, that the ABC model has 1/N corrections, as opposed to the 1/N 2 corrections in the PBC model.
We emphasize that the spin chain argument predicting (S2.35) is general and does not rely on the specific couplings of the model, as long as it admits a description in terms of matrix integrals.As a matter of fact, the relation holds in the real-time dynamics as well.
To conclude the digression into imaginary-time dynamics and the associated third order phase transition, we compare in Figure S2 the finite N results with the asymptotic formula.The plot shows the perfect agreement of our asymptotic expressions, derived analytically, with the numerical exact result at large but finite N .

D. Dynamical quantum phase transition: ABC
At this stage we can go back to the main subject: solving the real-time dynamics at leading order in N 1.We begin from the chain with ABC and solve (S2.24b) with ξ = iτ first, that is The chain with PBC is discussed in the next subsection.
Because the right-hand side of (S2.39) is imaginary, the standard techniques do not apply, as it is obvious that they would produce a real-valued left-hand side.We therefore need to extend the methods utilized above to account for the imaginary values of the parameters in the effective action.
We begin with a bird-eye discussion of complex actions and their saddle points, and then go on and solve the planar limit of the Loschmidt echo for the chain with ABC.

Remark S2.7 (A minus sign).
To avoid dragging cumbersome factors of (−1) in the ensuing computations, we make a change of variables xa = −x a in the matrix model (S1.14) (and drop the tilde henceforth).In practice, it is tantamount to solve (S2.39) with replacement τ → −τ .

Complex saddles
We have noted that the saddle point equation (S2.39) requires a complex solution for ρ(x).It follows from the effective action lim ) having an imaginary coefficient.Therefore, in the planar large N limit, G N,L=∞ (t) ≈ e −N 2 S (ABC) eff [ρ] is dominated by complex saddle point configurations.It is known that, in general, effective actions in Quantum Field Theory and in matrix models may admit complex saddle points.We refer to [S31-S33] for extensive discussion, and [S34, S35] for a treatment of the complex saddle points in matrix models.An important feature of the models (S1.9a)-(S1.14) is that the complex saddles of the effective actions are the dominant ones in the planar limit.
Remark S2.8.The fact that the saddle point configuration may be complex does not automatically imply that the action evaluated on that saddle point is complex.In fact, it turns out not to be the case for the models we consider.
In conclusion, to solve (S2.39) we gain insight from the Quantum Field Theory approach and relax the condition on the integration domain, so that we let it pass through the complex saddle points.Remark S2.9 (Replica wormholes).The inclusion of complex saddles in the gravitational path integral is playing a prominent role in black hole physics.The recent efforts to solve the Hawking paradox and the factorization problem build on the refined analysis of complex saddle points, known as replica wormholes [S36, S37].It has been proposed that these wormholes can become the dominant saddle point, against the Hawking saddle point configuration.This change of dominant saddle is accompanied with a first order phase transition, consistent with the Page curve.See [S38] for an overview.

Contour deformation and reality condition
In the light of the above discussion, in practice we allow suppρ = Γ ⊂ C, i.e. the saddle point configuration of the matrix model (S1.14) places the eigenvalues on some arc or union of arcs Γ in the complex plane.See e.g.[S39] for an early discussion and [S40] for a recent application in the context of black holes.
Of course, consistency requires and we will henceforth only focus on solutions Γ that, for small enough values of τ (in a sense made precise momentarily) are homotopic to [−1 , 1].
Having relaxed the assumption on suppρ, we write (S2.39) as .40)where now z, y ∈ C need not be real.Notice the minus sign in the right-hand side compared to (S2.39), according to Remark S2.7.We define a coordinate σ ∈ Γ and parametrize the embedding Γ → C through σ → z(σ), which we normalize according to | ż(σ)| 2 = 1.The integrals below are understood in such a parametrization, which we leave implicit.
After having obtained a distribution ρ(x) that solves (S2.40), we determine the shape of Γ requiring the probabilities So, in particular, Γ must be a level set along which the integral is real-valued.

Solution in the first phase
Gaining insight from the imaginary-time solution in Subsection S2 C, we expect that the solution in the first phase is the analytic continuation of the imaginary-time solution.
Indeed, under the ansatz that Γ is homotopic to [−1, 1], the function where ArcSin is the arcsine function, i.e. the extension of the inverse sine function sin the reality constraint becomes which is (S2.43)upon multiplication by √ −1.Equation (S2.43) admits disconnected solutions and, by construction, we ought to retain the one component that is smoothly deformed into

Cross-check: consistency with known results
Notice that, taking τ → 0, the condition (S2.41) reduces to where we have chosen as initial point z 0 = −1.The latter equation is solved by z ∈ [−1, 1], thus giving back the contour as it should.
Furthermore, replacing iτ → γ in (S2.41), that is, applying the argument to the imaginary-time dynamics of the system, the solution is still given by z ∈ [−1, 1] (cf. Figure S3), consistently with the standard procedure and the result of Subsection S2 C.

Criticality
At every τ > 0, the contour (S2.43) intersects the imaginary axis in the upper half plane.The eigenvalue density ρ(x) in (S2.42) has a complex zero at  that is, placed in the positive imaginary semi-axis.Increasing τ > 0, z 0 (τ ) descends from i∞ and eventually intersects Γ at the critical value τ cr , which is the (lowest positive) solution to The solution is At τ > τ cr , Γ breaks up into a two-cut solution, with the gap in the support opening around z = z 0 (τ cr ) ≈ i1.51.
See Figure S4.From the random matrix theory perspective, this is a one-cut to two-cut phase transition, with the peculiarity that the cut Γ is deformed away from the real axis.This is distinct from the hard-edge to soft-edge transition we have seen for the imaginary-time (finite temperature) model in Subsection S2 C.

Solution in the second phase
In the new phase, the two disconnected arcs forming Γ join −1 to A and B to +1, respectively, for a pair of points (A, B) in the upper-half C-plane.The reflection symmetry of the problem imposes A = B and follows from the continuity of ρ(x) at the critical point.
In the new phase, we can set up a Riemann-Hilbert problem for the planar resolvent ω(z), similar to what was done in the second phase in imaginary-time, and solve it adapting the standard techniques to the present case.
We omit the technical details, which are involved but standard.One can follow [S17] step by step for the two-cut phase, noting that the procedure is solely based on complex contour deformations, which go through identically when the endpoints A, B lie in C and not in R. The upshot is that we find the two-cut solution and A = B = i 2τ , while A = − B is implicitly fixed by the normalization condition and the knowledge of A.
Let us mention that, a posteriori, one may have leveraged the insight from two-cut matrix models with hard edges supported on subsets of R, as for example in [S41], to derive (S2.45).Indeed, the naïve ansatz for the resolvent would be precisely the planar resolvent ω(z) associated to (S2.45).Then, running the argument backwards, by selecting a closed contour of integration encircling Γ one shows that such ω(z) solves the Riemann-Hilbert problem descending from (S2.39).

Dynamical free energy
Having obtained the eigenvalue density, we can use it to compute the dynamical free energy (S2.26).In the first phase we have df (ABC)  dτ where we use the form (S2.42) of the eigenvalue density.To pass to the second line we have used that Γ is homotopic to [−1, 1].
We have not been able to obtain a closed form expression for f (τ ) in the phase τ > τ cr .Nevertheless, we can leverage the continuity of df (ABC) dτ at the critical point, which follows from the continuity of ρ(x) and Γ.We then compute the expansion With these expressions at hand, we can compute dτ 2 τ >τcr close to the critical point, which is enough to obtain the order of the phase transition.
Proposition S2.11.The dynamical free energy f (τ ) is twice differentiable and its third derivative is discontinuous at τ = τ cr .
Proof.We start by computing the second derivative of f (τ ) close to the critical value τ cr : We have used the two-cut eigenvalue density (S2.45) replacing A and B with their approximate form (S2.46).Thus, we first differentiate and then evaluate the resulting integrals.The contributions involving dA dτ and dB dτ vanish because they multiply the integrand evaluated at the boundary, that vanishes.We finally get lim The same argument can be improved to compute directly the first derivative of the dynamical free energy in the second phase.We obtain the result the symbol "≈" indicating that we only retain the lowest non-trivial order in (τ cr − τ ).We conclude that the phase transition is at least third order.By applying the same ideas to d 3 f (ABC) dτ 3 τ >τcr we have explicitly checked that it is non-zero at τ cr .
We conclude that the dynamical free energy signals a third order DQPT.

Late time behaviour of the dynamical free energy
We can study the late time asymptotics τ → ∞ of the dynamical free energy f (τ ).In this regime, the symmetry arguments are enough to constraint the form of A and B at leading order in 1/τ .One finds and a direct computation akin to the ones above yields which means the dynamical free energy behaves logarithmically at late times.

Relation with earlier works
Before concluding the present analysis, it is worthwhile to mention that the transition we find is essentially the one first discovered in [S42], and further elaborated upon in [S43, S44].Indeed, even though the model we are interested in is (S1.9b), after changing variables into (S1.14)and applying Andréief's identity, we land on a Hankel determinant which only differs by the one considered in [S42-S44] by terms which are sub-leading in N , cf. also Remark S2.2.Not surprisingly, then, we find agreement with the results therein, although using different mathematical techniques, namely a saddle point analysis in the complex plane instead of the asymptotics of a system of orthogonal polynomials.

Statement of the result
We now solve the saddle point equation for the chain with PBC and prove the existence of a third order phase transition at the critical point (S2.44).Moreover, in agreement with the general argument presented at the end of Subsection S2 C, we will find by direct computation the relation We consider (S2.24a) and rewrite it in exponential variables z = e iθ , u = e iϕ : Recall that the distribution is related to ρ through (e iθ ) = ρ(θ).As in the real-time ABC chain, we ought to relax the assumption on the integration contour, in order to catch the contribution from the complex saddles of the model (S1.9a). 3

Solution in the first phase
We begin with the ansatz supp = Γ where, as opposed to the imaginary-time model of Subsection S2 B (i.e. the GWW model), we let Γ be a 1-cycle in C homotopic to the unit circle S 1 .Topologically, [Γ] ∈ H 1 (C * ) but the geometry (i.e. the shape) of the curve Γ depends on τ and is subject to the constraint lim τ →0 + Γ = S 1 .With this ansatz, and gaining insight from the imaginary-time result in Subsection S2 B, we obtain that solves the saddle point equation.We also observe that where Int(Γ) (resp.Ext(Γ)) is the interior (resp.exterior) of the closed Jordan curve Γ.
Proposition S2.12.The 1-cycle Γ is the unique connected component homotopic to S 1 of the level set Proof.It follows from the analogue of the reality condition (S2.41) by direct computation.

Criticality
The solution (S2.48) for (z) has zeros in C located at and a phase transition takes place at the lowest positive value of τ at which either of these points hits Γ.The criticality condition is then which, using elementary manipulations, becomes This is exactly the same equation for τ cr found in the ABC chain, cf.(S2.44).The phase transition then takes place at exactly the same critical value, as expected.
A direct study shows that, as τ is increased, Γ is progressively deformed away from S 1 and eventually pinches at τ = τ cr .See Figure S5.

Dynamical free energy
As in the ABC chain, we are unable to find a closed form expression for f (τ ) in the phase τ > τ cr .Nevertheless, we can mimic the computations above and evaluate df (PBC)  dτ 0<τ ≤τcr = 2τ using (S2.48).Besides, approximating near τ cr in the second phase, we are able to prove that the transition is third order.The derivation is almost identical to the previous subsection.We conclude by showing f (τ ) in Figure S6.For comparison with the imaginary-time phase transition, we plot f (PBC) (τ ) and F (PBC) (γ) together in Figure S7.The plots make manifest that the dynamical free energy has a much more marked discrepancy from the parabola τ 2 than its imaginary-time analogue, in agreement with the logarithmic behaviour at late times derived analytically.

Relation with earlier works
Expression (S1.9a) for the Loschmidt echo was thoroughly analyzed in [S15].However, that work did not consider the planar limit.Our result for f (PBC) (τ ) precisely interpolates between the large N (at fixed t) and large t (at fixed N ) regimes of [S15].

Double-scaling limit and universality class
It is an intriguing open mathematical problem to rigorously establish the universality class to which the DQPT belongs.
Setting up a Riemann-Hilbert problem, following e.g. in [S46, S47], one may take a double-scaling limit in which τ → τ cr ) as N goes to ∞, with (τ cr − τ ) ∝ N −η for some positive rational exponent η.However, it seems not possible to find a scaling such that the Riemann-Hilbert problem reduces to that of Painlevé II, the latter being the case for the double-scaling limit of the GWW model.
Therefore, we expect that the DQPT we uncover belongs to a different universality class than GWW.From the analysis of the random matrix mechanism behind the transition, we may conjecture that the corresponding equation governing the double-scaling limit is Painlevé I.That is to say, the universality class would be the same as two-dimensional quantum gravity [S48-S50], rather than GWW.It would be interesting to prove (or disprove) this conjecture.

S3. ROBUSTNESS AGAINST IMPURITIES IN THE INITIAL STATE
It may happen that, in the preparation of the initial state |ψ 0 , one or more of the spins are not aligned, thus producing a state that slightly differs by the single domain wall state |↓↓ • • • ↓↑↑ • • • ↑ .A desirable feature of a quantum quench, that makes the DQPT more easily accessible experimentally, is the robustness of the phase transition against the introduction of such impurities, or defects, in the initial state.In this section we show that it is indeed the case for the quench we are considering.Recent theoretical studies of spin chains (or lattice models) with defects include [S51-S55].where the notation D (p) stands for defect.As usual, we assume for simplicity L N and set L = ∞ throughout the analysis.Besides, we consider the chain with PBC, being the analysis for ABC completely analogous.

Statement of the result
A careful analysis of the amplitude in presence of an impurity, defined in (S3.51), using random matrix theory techniques, leads to lim .52)for every t, meaning that the presence of an impurity does not spoil the DQPT.

Derivation of the leading order behaviour
From the derivation in Subsection S1 B, in particular from (S1.5), we have that where A p means the p th antisymmetric representation of U (N − 1).We are interested in large N behaviour of the ratio where the average is taken in the random matrix ensemble (S1.9a).To reduce clutter, here and in the rest of the section we omit the dependence on N, L and t when there is no risk of confusion.
To begin with, we introduce the generating function of the defects:4 Recall that, when written in terms of the eigenvalues e iθa a=1,...,N of the random matrix U , the character in the p th antisymmetric representation is the p th elementary symmetric polynomial.The generating function then reads It is explicit in the latter form that the determinant insertion in the matrix model (S1.9a) shifts the integrand by with the second summand in the right-hand side being O(N ).Hence, its effect is negligible at large N when compared to the term N 2 S eff , and it does not alter the saddle point solution.We conclude that whence the advertised result (S3.52).

Derivation of the defect contribution
We can go on and evaluate the first non-trivial contribution to the ratio (S3.53), which is the contribution of the defect in the planar limit.As we have said, the insertion of the impurity has a sub-leading effect that preserves the saddle point solution.We can thus evaluate D(w) in the planar limit using the eigenvalue density computed in Subsection S2 E. Rewriting (S3.54b) in the planar limit, we obtain As with the dynamical free energy, it is easier to evaluate the derivative In the first phase τ < τ cr the integral evaluates to i2τ .Exponentiating and expanding, we finally arrive at τ <τcr ≈ (it) p p! valid in the large N limit.In particular, the closer the impurity to the "outer region" with all spins up, corresponding to a smaller the value of p, the smaller is its net effect.Once more, we lack an explicit expression in the phase τ > τ cr , but the formulas, and especially the scaling with N , apply equally well.

Other impurities in the initial state
The previous derivation is easily adapted to the case in which the impurity consists of a spin |↓ N +p , p positions away from the domain wall.The difference in this case would be to replace N + 1 with N − 1, which of course has no effect in the large N limit, and to replace the p th antisymmetric representation A p with the p th symmetric representation S p .The defect generating function (S3.54b) is replaced by while the rest of the argument goes through essentially unchanged.

S4. DYNAMICAL QUANTUM PHASE TRANSITION AT FINITE SYSTEM SIZE
We now move back to the problem of analyzing the effect of a finite number of qubits on the DQPT.To this aim, we consider the finite-L, discrete matrix models (S1.12) and take the scaling limit together with the planar limit with τ = t N fixed.

Statement of the results
The aim of the present section is to argue that the DQPT discussed above persists in the limit (S4.55) that accounts for the finite size effects.Two take-home messages emerge from our analysis: • A finite triggers a phase transition, which nevertheless does not spoil the DQPT previously uncovered if 1.2; • In the phase τ < τ cr and with ≥ 2, the dependence on the system size is exponentially small in (L − 2N ).

Discrete matrix models and phase transitions
It has been known since the work of Douglas and Kazakov [S60] that the discreteness of a matrix ensemble may induce a phase transition.The argument is simple: the discreteness imposes a lower bound on the distance between any two eigenvalues of the random matrix ensemble, namely the lattice spacing.This lower bound becomes a constraint on the eigenvalue density at large N .Whenever a given solution to the saddle point equation ceases to satisfy such a constraint in a region of parameter space, the model undergoes a phase transition.
For the models we consider, which are discretized versions of U (N ) and U Sp(2N ) ensembles, we have the following result, see [S61].where ρ D unif. is the normalized, uniform probability density on the domain D. The discrete domains in (S1.12), arising from the consideration of the physical spin chain, are discrete subsets of U (1) subject to the following additional constraints.Lemma S4.2.The eigenvalue densities associated to the discrete matrix ensembles (S1.12) satisfy the condition PBC: It is crucial that, by construction, a finite number of qubits induces a discretization of the matrix ensembles (S1.9).Then, in the ABC chain, the change of variables x a = cos θ a produces a non-uniform discretization of the interval [−1, 1], inherited by the uniform discretization of the semicircle.This is encapsulated in the expression (S1.15).
Proof of Lemma S4.2.Lemma S4.2 is a corollary of the more general Lemma S4.1.Since we will only need the content of Lemma S4.2, we prove it explicitly, building on [S21, S22].
The derivation of (S4.56) is as follows.Consider the discrete ensembles in (S1.12) and use the invariance of the summand to restrict the discrete angles {θ(s j )} j=1,...,N , as defined in (S1.11), to the principal Weyl chamber which, having the planar limit in mind, is conveniently rewritten in the form Here C is the length of the support onto which the eigenvalues live, which is the circumference C = 2π for PBC and the semicircle C = π for ABC.Obviously, the constraint is ineffective if one first takes L N .Instead, expressing everything in terms of the discrete complex variables z(s) and sending N → ∞ with the scaling (S4.55), one gets C| (z)| ≤ , which yields (S4.56).In imaginary time, the absolute value can be omitted, because (z) = ρ(θ) is positive and of real argument.However, to generalize to real time dynamics, we ought to be careful with the absolute values.
The rest of the present section is organized as follows.
• Before delving in the analysis of the fate of the DQPT when considering a realistic chain with L < ∞, we mention related results for the imaginary-time model, i.e. the discrete GWW model.These are collected at the beginning of Subsection S4 A.
• Then, in the body of Subsection S4 A we show explicitly how to derive a finite-size induced phase transition in the imaginary-time PBC chain, and extend the result to the imaginary-time ABC chain.
• Having gained insight in the (arguably simpler) imaginary-time model, we face the effects of a finite system size in the real-time models we are interested in Subsection S4 B.
• The final Subsections S4 C-S4 D are devoted to show that the infinite system Loschmidt echo L N,∞ (t) approximates the finite size Loschmidt echo L N,L<∞ (t) with exponential accuracy in the first phase.

A. Finite size effects in imaginary time
Exponential accuracy of the infinite system approximation in imaginary time Consider the imaginary-time version of the discrete ensemble (S1.12a), that is, the discrete GWW model.Rigorous estimates for the rate of convergence of such discrete ensemble (and its generalizations) to its continuous analogue (S1.9a), from the first phase, have been given in [S61].See also [S8] for a discussion in the context of spin chains.It turns out that the error in considering the continuous matrix model instead of the finite-L, discrete model is exponentially small in (L − N ) [S61].
Moreover, the scaling limit of the discrete GWW model has been extensively analyzed in [S21] (and later in [S22] a generalization of the model), where it is shown that the discreteness of the ensemble induces a Douglas-Kazakov-type [S60] phase transition as a function of the parameter .
It it obvious that, for finite 1, the system size has no effect on the GWW phase transition.More specifically, the finite value of produces a Douglas-Kazakov phase transition at which does not affect the large N GWW phase transition if ≥ 2 + ε for arbitrary but fixed ε > 0. We refer to [S21] for an exhaustive analysis.
It is an easy task to repeat the argument to the imaginary-time chain with ABC.From the eigenvalue density (S2.37) and imposing (S4.56b), we get the same critical value γ = −1 2 for the phase transition induced by the discreteness.
The special value = 2 has a neat physical meaning.Recall that = L/N , thus a chain with PBC and = 2 in the state |ψ 0 has exactly half of the qubits |↓ and half |↑ .For < 2, it would be convenient to rephrase the problem as starting with the state

B. Finite size effects on the DQPT
We start by recalling that the map x = cos θ produces a non-uniform discretizaion of the interval [−1, 1], and hence of Γ.The constraint (S4.56b) then reads In particular, define so that a finite-size induced DQPT takes place at the critical value τ which solves Both eigenvalue densities have local maxima at z = ±1, and attain the same extremal value √ 1 + 4τ 2 .The other local maxima over C do not satisfy the initial condition, namely do not belong to the undeformed contour in the limit τ → 0. We then have for both PBC and ABC.Requiring that this phase transition takes place at τ > τ cr , not to invalidate the derivation of the DQPT discussed above, we are led to the condition > where C. Finite vs infinite system: a numerical study To give a quantitative estimate of the error in approximating a realistic chain, consisting of a number L of qubits, with L = ∞, we define From [S61] we have for a function c(γ) > 0 which does not depend on and a suitable δ > 0. We have derived above δ = 1 at leading order in the planar limit.Thus, before the GWW phase transition, the error in approximating the imaginary-time system by an infinite number of qubits is exponentially small in the system size.
One of the outcomes of Section S2 is that, in the first phase, the real-time and imaginary-time dynamics are related by analytic continuation, a fact that does not hold beyond the respective phase transitions.From that, we are tempted to expect that the finite-L dependence in the first phase is exponential also in the real-time system.

D. Exponentially accurate thermodynamic limit
As already mentioned above and elucidated in [S8], borrowing rigorous estimates from the work [S61] one proves that, in imaginary time, the ratio for a c(γ) > 0, 0 ≤ γ < 1 2 and a suitable δ > 0. Therefore, for a sufficiently large system, the Loschmidt echo can be approximated by its infinite-size limit with error exponentially small in (L − N ).
Although the procedure of [S61] does not straightforwardly extend to the real time dynamics, Figure S9 suggests that similar estimates hold.In the rest of this subsection, we argue that this is indeed the case.
We work with the discrete model (S1.15), i.e. we henceforth focus on the chain with ABC.Following The integration contour, shown in Figure S10, is , where for arbitrarily small ε, ε > 0. Notice that we are cutting an arbitrarily small neighbourhood of the edges x = ±1 of size ε.This choice excises the accumulation points x = ±1 of the zeros of the orthogonal polynomials we will introduce momentarily.Because of the square root behaviour of the measure near the edges, the integrand vanishes linearly in ε.
The contour integral on the right-hand side of (S4.58) is solved by residues and, thanks to the suitable choice of function ν(z), yields the left-hand side.
Along the way, observe that, writing z = e iArcCos(x) , with ArcCos(x) the arccosine function, Here ν(x) is a shorthand notation.For x ∈ Σ − (resp.Σ + ), z = e iArcCos(x) is mapped to a contour slightly inside (resp.outside) the semicircle e iθ : 0 ≤ θ ≤ π , see the plot in Figure S11.where {P j } j∈N is a system of monic polynomials.In practice we use the orthogonal polynomials with respect to the measure e iN λDx dx studied in [S42], with identification λ D = −4τ .Note that, in (S4.59), we have used the invariance of the determinant under addition of rows and columns to replace the factor x a+b−2 , that arises in taking the Hankel determinant, with P a−1 (x)P b−1 (x) = x a−1 x b−1 + lower terms.
From the analytic properties of the integrand in each entry of (S4.59), it follows that where h a = h a (t), valid for τ distinct from every point at which any of the {h a } a=0,...,N −1 vanishes.Importantly, a sufficient condition at leading order in N is τ < τ cr , in perfect agreement with the rest of our discussion so far.
With the planar limit in mind, we restrict to τ < τ cr and write [S61]: The auxiliary function v(x) accounts for both the change of orientation and the separation of the diagonal term h a−1 δ ab .In terms of the variable z = e iArcCos −1 (x) , it reads which stems from the basics of random matrix theory.We emphasize again that we are using polynomials that are orthogonal with respect to the continuous measure on [−1, 1], not with respect to the discrete one.
We finally arrive at So far we have mostly adapted and extended [S61] to the non-uniformly discretized Hankel determinant (S1.15).It remains to show that the matrix K is exponentially damped in N, L 1.
We have already observed that v(x) is exponentially small as a function of L: x) .

Summary: finite versus infinite system
Let us summarize our finding for L < ∞, its thermodynamic limit with scaling (S4.55), and its relation with the system in which we set L = ∞ from the onset.
• The Loschmidt echo in the infinite system approximation undergoes a third order DQPT.This is shown in Subsections S2 D-S2 E.
• A finite system of L qubits experiences the same phase transition if the system size satisfies > ≈ 1.2.The bound is derived in Subsection S4 B.
• Furthermore, in the phase τ < τ cr , the approximation of the finite system Loschmidt echo L N (t)| L<∞ by L N (t)| L=∞ in the limit (S4.55) is exponentially accurate if > 2. This bound is proven in Subsection S4 D.
Let us stress the distinction between > , which guarantees that the DQPT is not overturned by finite-L effects, and the requirement > 2, which guarantees that the results in the L = ∞ approximation reliably describe a finite chain. in the strict large N limit, without scaling t.By analyticity of the Loschmidt amplitude, we have that the estimate extends to the limit with scaling, as long as τ remains away from τ QSL .That is, Besides, we know that the left-hand side is analytic for τ < τ cr .Imposing the equality of the two limiting values we get the constraint τ QSL − ε ≥ τ cr ∀ε > 0.
Step 3. Identify the QSL We have accomplished step 2 of the proof exploiting the relation (S5.64), which bridges between even and odd N .Let us recall that it is an approximate identity, which holds up to terms that are sub-leading in N , which is sufficient for our purposes.
Let us rewrite (S5.64) by sifting N + 1 → N : We also observe that, near the first zero, G N (t) ∝ (t − N τ QSL N ), which again is a consequence of the analysis in [S43] and we neglect sub-leading terms.The linearity near the zero implies that lim Taking one further derivative of (S5.65), we have on the left-hand side and on the right-hand side We now set N = 2n + 1, take the logarithm of both expressions, divide by (2n + 1) 2 , and eventually take the limit n → ∞.Equating the two sides, as dictated by having differentiated (S5.65), we have that

Figure 2 .
Figure 2. Different behavior of the Loschmidt echo LN (t) at even and odd N , exemplified in N = 3, 4.

Figure 3 .
Figure 3. Dynamical free energy.Left: N = 3, L = 9.Right: N = 4, L = 11.The bold blue line is the exact value for the given L, the thin orange line is the exact value at L = ∞ and the black dashed line is the extrapolation if there were no phase transition.

Lemma S4. 1 .
Let ρ(z) be the eigenvalue density associated to a discrete matrix model with discrete domain D, with |D|.Then ρ(z) is upper bounded according to |ρ(z)| ≤ ρ D unif.
L k=1 |↓ k = |↓↓ . . .↓ and act with the spin-flip operators σ + 1,...,L−N on the last (L − N ) qubits.The value = 2 is the self-dual point of the Z 2 spin flip duality.At such point, the duality of the model becomes a symmetry.

6 Figure S9 .
Figure S9.Plot of EN as a function of .

Figure S10 .
Figure S10.Choice of contour Σ for the application of Cauchy's theorem.
thus is exponentially suppressed in L in both cases.Next, we bring out {h a−1 } a=1,...,N from the corresponding columns in (S4.60), noticing that G
.19) Remark S2.1.In the original gauge theory framework, w is inversely proportional to the squared Yang-Mills coupling g 2 YM , so that 1 ξ ∝ g 2 YM N is what is customarily called 't Hooft coupling in the gauge theory literature.