Rigorous bounds on dynamical response functions and time-translation symmetry breaking

Dynamical response functions are standard tools for probing local physics near the equilibrium. They provide information about relaxation properties after the equilibrium state is weakly perturbed. In this paper we focus on systems which break the assumption of thermalization by exhibiting persistent temporal oscillations. We provide rigorous bounds on the Fourier components of dynamical response functions in terms of extensive or local dynamical symmetries, i.e. extensive or local operators with periodic time dependence. Additionally, we discuss the effects of spatially inhomogeneous dynamical symmetries. The bounds are explicitly implemented on the example of an interacting Floquet system, specifically in the integrable Trotterization of the Heisenberg XXZ model.


Introduction
Response functions, or susceptibilities, can be used to probe the symmetry breaking phenomena. In analogy with static susceptibilities, which probe transition from the ordered to disordered phase and, for example, characterize the spontaneous breaking of the spacetranslation symmetry in crystals [1], the dynamical susceptibilities carry information about the dynamical phases of matter. For instance, ideal conductivity, which is a particular manifestation of ergodicity breaking, can be related to the Drude weight, which corresponds to the zero-frequency behavior of the dynamical response function [2]. More generally, non-vanishing dynamical response functions at finite frequencies imply breaking of time-translation symmetry of equilibrium states.
Obtaining rigorous or explicit results in quantum strongly interacting many-body systems is a formidable task even in the presence of "exact" solvability. Typically one has to rely on numerical simulations, which, however, are bound to fail, either due to finite size effects or because of the entanglement growth. One of the most important rigorous results explaining the origins of non-ergodic behavior is Mazur's lower bound on asymptotics or time-averages of dynamical correlation functions [3]. Through the scope of this bound, non-ergodicity can be understood as a consequence of underlying extensive or local symmetries of the system. In particular, local symmetries can be related to localization phenomena in many-body systems [4][5][6], while extensive symmetries lead to ideal transport at arbitrary temperature and lack of thermalization in integrable models [2,[7][8][9].
The emergence of time-translation symmetry breaking in strongly interacting systems has recently been related to extensive dynamical symmetries [28], calling for the development of a rigorous framework. In this paper we provide such a framework by deriving strict lower bounds on AC dynamical response functions in terms of dynamical symmetries. The lower bounds are derived for autonomous Hamiltonian dynamics as well as for periodically driven (Floquet) systems. We apply our results to a nontrivial example of a many-body Floquet system that breaks the discrete time-translation symmetry, specifically to the integrable Trotterizaiton of the spin-1/2 XXZ model [29,30].

Breaking of time-translation symmetry
Breaking of the time-translation symmetry is associated with a failure of a perturbed stationary state to return to stationarity, even in the infinite-time limit t → ∞. Analogously, we can consider the space-translation symmetry breaking from a dynamical point-of-view: it occurs, when the information about spatial inhomogeneities induced by the perturbation is retained indefinitely after the perturbation has been switched off. In general we will consider stationary states ρ(µ), where µ denotes the set of chemical potentials. Typically, only few chemical potentials describe the complete steady-state manifold pertaining to the system. In Hamiltonian systems, one of the chemical potentials is the inverse temperature µ 1 = β. Assuming that the thermal state of the system ρ(β) = exp(−βH) is slightly perturbed at t = 0, say by a Hermitian operator B, i.e. H → H − (δ/β)B, we probe the dynamics of the local observable A by considering the first order in the expansion of its expectation value. Here, we have introduced the averages with respect to the perturbed, A δ = tr (Aρ)/ tr (ρ), and thermal state, A ≡ A 0 . The response of expectation values to a perturbation can be interpreted in terms of the canonical Kubo-Mori-Bogoliubov inner product Physically sensible perturbations of extensive Hamiltonians H are either extensive, e.g. a sum of local operators acting on a spin lattice, or themselves local. In order to probe the frequency dependence of the response, it is useful to introduce the dynamical response function At zero frequency it represents the time-averaged perturbative correction (1) to A 0 . Now, the system is called non-thermal, provided that f AB (0) = 0 and f HB (0) = 0. The condition f HB (0) = 0 ensures that the energy is conserved in the first order of the perturbative expansion, and in turn implies that the thermal ensemble with associated expectation values should remain the same under the assumption of thermalization. However, due to f AB (0) = 0, the time average of the expectation value does not coincide with the old thermal average A 0 . If, in addition, the finite-frequency response function does not vanish, i.e. f AB (ω) = 0, for some ω = 0, the system does not relax to any stationary distribution, for arbitrarily weak perturbation of the equilibrium ensemble. This, then, characterizes the breaking of time-translation symmetry.
Similarly, the breaking of space-translation symmetry on the lattice can be probed by spatially modulated extensive observables with local densities a n , via the frequency and wavevector-dependent response function As before, the nonvanishing susceptibility f AB (k, ω) = 0 for a nonzero value of the wavevector, k = 0, implies that the system will retain memory related to the spatial modulation of the perturbation. Note that our definitions (3) and (5) differ from the standard definitions of dynamical response functions by an extra factor of 1/T in the Fourier transformation. This factor ensures that f AB (ω) is finite, and not a Dirac delta singularity, for a perfectly harmonic response at frequency ω.

Local and extensive dynamical symmetries
Analogously to how the non-thermal behavior can be understood as a consequence of local and extensive conservation laws [8], the time-translation symmetry breaking relates to the existence of local and extensive dynamical symmetries [28].
For a concise presentation let us revisit the notions of locality and extensivity, as they play a prominent role in our story. We consider extended quantum systems defined on a regular lattice of N sites with finite-dimensional local Hilbert space. Prominent examples of such systems are spin chains or spin lattices. The notions of (effective) locality and extensivity are, in general, state-dependent [31]. Operator a is local with respect to the thermal state ρ (or generic clustering state that is invariant under the time evolution), if its Kubo-Mori-Bogoliubov norm is finite in the thermodynamic limit Operator A is extensive with respect to the state ρ, if its norm is proportional to the volume of the system (i.e. the number of local physical sites, N ) in the thermodynamic limit, and has a nonvanishing overlap with at least one local operator b Extensivity, as defined here, is sometimes also referred to as pseudo-locality [32]. While extensive dynamical symmetries are of central importance when considering extensive perturbations of stationary states, local dynamical symmetries are crucial when dealing with local perturbations. Such symmetries should be employed to study discrete time crystals that can arise in many-body localized systems [15]. In systems with multiple extensive conserved quantities Q j , i.e. [H, Q j ] = 0, the array of equilibrium ensembles is naturally enlarged and local observables are expected to relax to their equilibrium values, described by the associated set of generalized inverse temperatures (or chemical potentials) β j [33,34]. The set of stationary states can be rigorously established using extensive charges as flows on the space of operators [31]. Such systems behave non-thermally, and are well-described by non-thermal maximum-entropy states, the so-called generalized Gibbs states ρ GGE = exp(− j β j Q j ), provided that the complete set of extensive integrals of motion Q j has been identified. Upon a slight perturbation of the state ρ GGE by an extensive operator, the expectation values of local observables will relax back to a slightly perturbed stationary state ρ GGE . Similarly, if the system admits local (non-extensive) integrals of motion q j , [H, q j ] = 0, it fails to relax to the thermal state even if the perturbation of the initial stationary ensemble is local.
In order to be considered as extensive (or local) dynamical symmetries, operators Q j should satisfy the definition of extensivity (or locality), as well as the eigenoperator condition Clearly, extensive dynamical symmetries admit a simple periodic time-dependence Q j (t) = exp(iω j t)Q j . The notion of (dynamical) symmetries can be trivially extended to non-autonomous timeperiodic (Floquet) systems. In this case the dynamics of local observable a is generated by a periodic, time-dependent Hamiltonian with a period τ . In the case of time dependent Hamiltonians we will focus on the stroboscopic time evolution of observableŝ where ←− exp is the time-ordered exponential, the arrow denoting the direction of increasing time in the time-ordering. In the case of non-autonomous systems, local or extensive dynamical symmetries can again be defined by locality or extensivity and quasi-periodicity of their time-dependenceM Quasi-periodicity arises since ω j τ /(2π) is in general not a rational number, meaning that for any integer n > 0, (M τ ) n [Q j ] = Q j . Thermalization and equilibration in quantum systems can be understood through the scope of eigenstate thermalization hypothesis (ETH). It is thus instructive to clarify the role of dynamical symmetries also in this regard. Clustering property of the initial state ensures that the dynamics is restricted to states within a narrow window of energy and expectation values of other extensive conservation laws. After a short time, the off-diagonal elements of the time-evolved observables vanish due to dephasing and suppression in the thermodynamic limit, and we observe the onset of stationarity, which can be described either by a representative eigenstate, or equivalently, by a statistical ensemble. Systems with dynamical symmetries avoid equilibration by circumventing two assumptions of the ETH. As a consequence of the eigenoperator condition (9), there exists an eigenstate |ψ corresponding to the energy E +ω j , for any eigenstate |ψ with energy E, provided that Q j |ψ = 0 (then, in fact, Q j ∝ |ψ ).This implies that (i) dephasing does not occur for eigenstates with the energy difference of the order O(1) and (ii) that the corresponding off-diagonal elements are not exponentially suppressed in the thermodynamic limit, resulting in perpetual oscillations of some local (or extensive) observables.

Bounds on susceptibilities
The purpose of this section is twofold. We first generalize the bound provided by Mazur [3] to AC response functions and secondly, provide general bounds on the off-diagonal components of the susceptibility matrix. We will first focus on time-independent Hamiltonian systems. Let us start by considering the following operator where Q j are either extensive or local dynamical symmetries. To obtain the lower bound on the Fourier components we consider the Kubo-Mori-Bogoliubov norm (1) of the operator (13), with respect to some stationary state ρ(β), which is clearly non-negative Writing out all of the components explicitly, we obtain the bound To study the susceptibilities in the thermodynamic limit N → ∞, the latter has to be taken after dividing each term by the system size N . Equivalently, we can, in this case, substitute the Kubo-Mori-Bogoliubov bracket A, B by lim N →∞ A, B /N . Only then we take the infinitetime limit T → ∞ of each term in the inequality (15). The optimal set of coefficients α j can be obtained by maximizing the quadratic form on the right-hand side of (15), yielding the set of equations where the Kronecker delta on the left-hand side emerges from time-averaging Remember that in the standard definition of response functions the prefactor 1/T is absent, whence the integral over time instead results in the Dirac delta peak. We remark that the overlaps Q j , Q l in Eq. (16) are nonzero only if ω j = ω l , which follows from the time-translation invariance of the thermal state. It thus suffices to consider only the dynamical symmetries Q j that correspond to the same frequency ω j ≡ ω l in the eigenoperator condition (9). In order to represent the results compactly, we now introduce a hermitian matrix of kernels Q = Q † , with elements Q j,l = Q j , Q l . Additionally, we require the knowledge of overlaps between the observable A and dynamical symmetries Q j , namely A j = δ ω,ω j Q j , A , gathered into the vector A. Plugging the solution of the set of equations (16) into (15) and taking the limit T → ∞, yields the lower bound The left-hand side of the above equation can be further simplified to if a weak requirement that the Fourier components of the dynamical response function asymptotically approach their averaged value is satisfied. The lower bound thus reads The result (21) straightforwardly generalizes to the off-diagonal elements of the susceptibility tensor. Specifically, considering a pair of different observables A, B and defining the overlaps B j = δ ω,ω j Q j , B , we obtain the following bound where the phase φ is chosen appropriately, to render the right-hand-side of the first inequality (22) real and nonnegative. In order to arrive at the bound (22), we substituted A → αA + βe iφ B in the inequality (21) and took the derivatives w.r.t. α andβ at α = β = 0. Note that the above result can be further generalized to systems with spatially modulated response functions where the overlaps A k , Q k , and B k contain only the contributions of spatially modulated extensive symmetries Q j , associated with the wave vector k, i.e.

S[Q
Here,Ŝ denotes conjugation by a single-site lattice shift (one-site translation). In the discrete time (Floquet) case the dynamical susceptibility can be defined as where l ∈ Z counts the steps in the stroboscopic time evolution. In this case the lower bounds can be obtained by replacing the continuous-time variable with stroboscopic steps, t → lτ , and time averages with sums, lim T →∞

2L
L l=−L , in the above derivations.

Driven Heisenberg chain at root-of-unity anisotropies
To illustrate how persistent oscillations of extensive many-body observables occur, we consider a Floquet driven spin-1/2 chain of even size N ∈ 2N, introduced in Ref. [29,30] as an integrable Trotterization of the anisotropic Heisenberg (XXZ) model. For simplicity, we assume that the thermal ensemble corresponds to the infinite-temperature state ρ(0) = 2 −N 1, although the results can be generalized to arbitrary time-translation invariant state generated by extensive conserved quantities of the model. We fix a unit period τ = 1 here, and consider a discrete (stroboscopic) evolution where time t is an integer. Specifically, for an arbitrary observable A we define dynamical map where Each of them comprises local unitary quantum gates acting on pairs of neighboring sites. Operators s α , for α ∈ {x, y, z}, denote the spin-1/2 matrices. The circuit representation of a single time step is shown in Fig. 1 Here we focus on an interesting regime of the model, characterized by (i) extensive symmetries that break the spin-flip (Z 2 ) symmetry, and more importantly, (ii) extensive dynamical symmetries that break its U(1) invariance, corresponding to the conservation of magnetization. While the first set of symmetries prevents decay of spin-current fluctuations and leads to ideal spin transport [30], the second set, which is the focus of this paper, prevents relaxation of local observables that couple different magnetization sectors, provided that the magnetic field is nonzero, i.e. h = 0. Since the dynamical symmetries are related to the integrability structure of the model, it proves useful to parameterize the Floquet operator in terms of anisotropy parameter γ and staggering δ Through this mapping the local unitary gate (28) is related to the trigonometric R-matrix, which constitutes the local conservation laws of the XXZ model; see Appendix A for the review of the integrability structure. In order to reproduce the continuous-time Heisenberg evolution, the propagator (26) should be expanded in δ to the leading order. Dynamical symmetries Y (λ) of the model are continuously parametrized: the discrete index j in Eq. (12) is substituted by a complex spectral parameter λ ∈ C. They stem from the so-called semicyclic complex-spin representation of the XXZ model's symmetry group [28,35,36] that lacks the lowest-weight state and exists only for γ ∈ π /m m + 1, ∈ 2N, < m , Appendix B for elaborate details on their construction. In short, the semicyclic representation of the symmetry group is spanned by tridiagonal m × m matrices L α , where α ∈ {0, +, −, z}, that encode coefficients in the operator-basis expansion of Y (λ) according to 1 Here, s 0 = 1 is a 2 × 2 identity matrix and s ± = s x ± is y , while |0 is the highest-weight state of the semicyclic representation, used to project out the auxiliary degree of freedom, upon which operators L α act. The absence of the lowest-weight state in the auxiliary space results in periodic action of the operator L + , sketched in Fig. 2, and lies at the origin of the U(1) symmetry breaking. 2 Due to this periodic action, each term of the continuously-parametrized extensive symmetries Y (λ) possesses a surplus of m spin-raising operators s + . In the presence of magnetic fields this results in the oscillatory time-dependence with a frequency, proportional to the order of the anisotropy m and the magnetic field strength h. Note that the evolution equation (32) associates Y (λ) † to a negative frequency ω = −hm. Since all terms in Y (λ) † contain a surplus of m spin-lowering operators s − , the continuous family {Y (λ) † } of adjoint dynamical symmetries is independent of {Y (λ)} and should be included in the bounds on the Fourier components of dynamical susceptibilities. In Appendix B we show that Y (λ) are extensive in the size of the system, whenever |Reλ − π/2| < π/(2m). In this strip in the complex plane one can compute their overlap Y(λ, µ) = Y (λ), Y (µ) according to a conjectured and thoroughly checked formula which is essential in the computation of the bounds on the dynamical susceptibilities.
1 Strictly speaking, Eq. (31) is inaccurate, since in each term of the operator-basis expansion, one of the operators in the product L α 1 . . . L α N actually differs from the rest. Nevertheless, this representation suffices for illustrative purposes. The accurate matrix-product form of Y (λ) is presented in Appendix B. 2 The terms highest-and lowest-weight state come from the fact that L + , in terms of algebraic relations that it satisfies, corresponds to a spin-lowering and not spin-raising operator on the auxiliary space (see Appendix B). To relate to the discussion of the latter, we will consider a set of extensive observables overlapping with dynamical symmetries Y (λ) and their adjoint counterparts Y (λ) † . The linear-response susceptibility of observable A evidently possesses nonzero Fourier components associated with frequencies ω = ±hm, arising from terms m−1 r=0 s ± n+r . They can be bounded from below by means of the system of integral equations which generalize the bound (21) to the case, where dynamical symmetries are enumerated by a continuous parameter λ instead of a discrete index (following Ref. [32]). The projection A(λ) = Y (λ), A onto the dynamical symmetries reads and can be discerned from the matrix product form of Y (λ), shown in Appendix B. Assuming that the system of dynamical symmetries is complete, the lower bound, computed by solving the linear integral equation (35), saturates. We can then conjecture the asymptotic behaviour of the temporal auto-correlation function to be This result is indeed corroborated by the numerical evidence obtained via the time-dependent density-matrix renormalization group method (tDMRG) shown in Fig. 3.
The time translation symmetry breaking is not specific to the Floquet driven XXZ model, but can be observed also in its continuous-time limit; see Ref. [28] for an elaborate discussion. The corresponding dynamical symmetries are simply Y (λ) evaluated at δ = 0, while the dynamical equation then reads

Conclusion
In this article we derived rigorous bounds on AC dynamical response functions, providing insight into the spontaneous time-translation symmetry breaking of the underlying dynamics. The research outlines a rigorous approach towards the study of quantum time crystals, and probing of the spatial symmetry breaking phenomena, such as dimerization, by shedding the light on its microscopic origins. From the broader point of view our study centers on types of information that can be preserved in interacting many-body quantum systems, and might thus prove lucrative from the point of view of quantum information processing and storage. From this perspective, finding concrete examples of systems, which go beyond the toy models presented in this paper, is of paramount importance. There are two somewhat related properties that should be at the forefront of this endeavour. The first one is stability to perturbations, and the second one the viability of implementing these systems in experimental setups.
While we focused on close-to-equilibrium setup in this paper, the dynamical symmetries have a profound effect on the long-time description when the system is initially prepared in the state that is far from equilibrium as well. It was conjectured that in this case the appropriate description is provided by time dependent maximum entropy ensembles [28], however any clear theoretical account of this phenomena is yet to be provided.
It was recently proposed [37][38][39] that less-local conservation laws might have effects on the equilibration rate and normal conductivity. The effects of e.g. quadratically extensive dynamical symmetries is yet to be explored in this context. into the trigonometric R-matrix of the spin-1/2 XXZ model [30]. Explicitly where P n,n+1 is a permutation (transposition) that exchanges the one-particle states of two neighbouring spins, while the trigonometric R-matrix reads For h = 0, the propagator U = U o U e with half-steps (27) can be recovered from the integrability structure as where we have introduced a continuous family of transfer matrices in which λ ± = λ ± δ/2 denote shifts in the spectral parameter. As a consequence of the Yang-Baxter equation

B Extensive dynamical symmetries in the driven XXZ model
The extensive dynamical symmetries in the driven XXZ spin-1/2 model originate in the transfer matrixT in which R-matrices have been substituted by Lax operators L a,n (λ) = 1 cos[γS z a ] + 2 cot λ s z n sin[γS z a ] + csc λ sin γ s + n S − a + s − n S + a , that form the higher-spin transfer matrices of the XXZ model [8]. Here, 1 is the identity operator, s z n and s ± n = s x n ± is y n act on the n-th spin-1/2 degree of freedom, while operators in bold act over the auxiliary space, labeled by index a. They generate a complex spin-s algebra, also known as quantum group U q (sl 2 ), where q = exp(iγ). Focusing on root-of-unity q (i.e. q k = 1, for some k ∈ N, which yields γ = π/k, ∈ 2N), we choose a semi-cyclic of the quantum group. Its dimension is set by the order of the root of unity, defined as m = min{k ∈ N | q k = 1}. The oscillatory behaviour of the dynamical symmetries originates in the absence of the lowest-weight 3 state, i.e. S − |m − 1 = β|0 . Before examining them in detail, we state several important properties: 1. They can be computed according to Eq. (49) only for odd orders m of the root-of-unity parameter q (see, however, Ref. [35] for the details on how to treat the anisotropies parametrized by the root-of-unity q of even order).
The first property is a consequence of the fact that the auxiliary representation (48) is only irreducible for odd m. A detailed discussion of this peculiar character of the semi-cyclic representations has been presented in Ref. [35]. For simplicity, we will restrict ourselves to the set of anisotropies, given in Eq. (30). Regarding the second property note that, in the absence of magnetic fields Y (λ) are exactly conserved. This follows from the specific type of Yang-Baxter equation, the so-called RLL relation which implies [T (λ), T (µ)] = 0, provided that h = 0. We proceed to examine the second and the third property in detail.

B.1 Explicit form of the dynamical symmetries
Due to cyclicity of the trace over the auxiliary space in (46), the derivative on β can always be translated to the leftmost position in the string of Lax operators. Then, using and denoting L 0 a,n (λ) = L a,n (λ) | β,s=0 , the operator (49) can be rewritten as Here and below, symbolŜ denotes the conjugation by a one-site lattice shift, e.g.Ŝ(s α n ) = s α n+1 , for which the periodic boundary conditions implyŜ N = 1.
We observe that, in each term, the highest-weight and the lowest-weight vectors in the auxiliary space have to be coupled by the string of Lax operators. For any n ∈ {1, 2, . . . N } we have (recall that β = s = 0) where the local densities read Except on the first 2r + 1 (respectively 2r) sites, these local densities act trivially, i.e. as identities. Importantly, they contain a surplus of m spin-raising operators s + n acting on the physical degrees of freedom, as the string of Lax operators still needs to connect the states |1 and |m − 1 in the auxiliary space. This can only be achieved via the auxiliary spin operator S − [see the representation (48)], which is coupled to s + n in the Lax operator (47). This explains the lower boundaries on the sums in (56), since contributions with a surplus of less than m spin raising operators vanish. Due to the surplus of exactly m spin raising operators in each term of Y (λ) we now have where relation [s z n , s + n ] = s + n has been used. We proceed to examine the extensivity of dynamical symmetries Y (λ).

B.2 Extensivity of the dynamical symmetries
For simplicity, let us consider infinite-temperature susceptibilities, so that the state that enters the rescaled inner product (2) corresponds to the featureless identity matrix ρ = 2 −N 1, while the inner product itself becomes of the Hilbert-Schmidt type. To determine the extensivity of the dynamical symmetries, we need to consider the kernel of overlaps that will be referred to as the Hilbert-Schmidt kernel. For ease of notation, we have conjugated the spectral parameter in the first factor of the inner product. The local densities (57) Its computation is further facilitated by employing the matrix product structure of the local densities. In particular, we can define the auxiliary transfer matrix T a 1 ,a 2 (λ, µ) = 1 2 tr n [L 0 a 1 ,n (λ)] T L 0 a 2 ,n (µ) , where (•) T denotes the partial transposition of the operators over the physical degrees of freedom, namely those, indexed by n in equation (47). Two indices a 1 and a 2 denote two copies of the auxiliary space. The auxiliary transfer matrix now facilitates the computation of overlaps in the Hilbert-Schmidt kernel (60), for example q [2r,±] (λ), q [2r,±] (µ) = m − 1, m − 1| T(λ ∓ , µ ± )T(λ ± , µ ∓ ) r−1 |1, 1 4 sin(λ + ) sin(λ − ) sin(µ + ) sin(µ − ) , while the extensivity of the dynamical symmetries Y (λ) now depends on the spectrum of its projection onto the subspace W = lin{|k, k | 1 ≤ k ≤ m − 1} that is invariant under its action, namely T(λ, µ)W ⊂ W. In particular, the sums over the support-size index r in (60) converge, if the spectrum lies inside the unit circle in the complex plane. Linear extensivity of dynamical symmetries Y (λ), i.e. finiteness of Y(λ, µ), now follows from the following observation, proven in [32]: Proposition 1 For | λ − π/2| < π/(2m) the eigenvalues of the auxiliary transfer matrix T(λ, µ) projected onto its invariant subspace W = lin{|k, k | 1 ≤ k ≤ m − 1} are strictly below 1 in the absolute value.
Evidently, the imaginary shifts of the spectral parameters by ±δ/2 do not cause violation of this condition as it constraints only the real part of λ. The explicit form of the Hilbert-Schmidt kernel Y(λ, µ) can be conjectured, similarly as was done in [30], and is given in Eq. (33). It has been extensively numerically and analytically checked that it also reproduces the normalized inner product of the semicyclic symmetries in the continuous-time limit δ → 0; see Ref. [35].