Diffusion from Convection

We introduce non-trivial contributions to the diffusion constant in generic many-body systems arising from the quadratic fluctuations of their ballistically propagating, i.e. convective, modes. Our result is obtained by expanding the current operator in the vicinity of equilibrium states in terms of powers of local and quasi-local conserved quantities. We show that the only finite contribution to diffusion constant arises from the second-order terms in this expansion. Our formalism implies that whenever there are at least two coupled modes with degenerate velocities, the system behaves super-diffusively in accordance with the non-linear fluctuating hydrodynamics theory. Finally, we show that our expression saturates the exact diffusion constants in quantum and classical interacting integrable systems.

Introduction: From its inception statistical physics has steamed to derive the universal laws of hydrodynamics and thermodynamics from the microscopic origins. Its strength lies in the universality, since the uncovered mechanisms have an unprecedented reach. One of the outstanding problems is to understand the transport behavior of strongly interacting many-body systems and identify the relevant degrees of freedom regulating the dynamics at large hydrodynamical scales. While it is well understood that ideal transport arises as a consequence of local conservation laws which prevent the current from decaying [1][2][3], much less is known about diffusive transport. Obtaining the diffusion constant of a quantum interacting system is, in general, a very nontrivial task, even with the powerful numerical methods available for one dimensional systems [4][5][6]. Recently many non-trivial results were obtained in quantum integrable systems [7][8][9][10][11][12], holographic matter [13][14][15][16][17], random models [18][19][20][21] and some general chaotic systems [5,22]. In this letter, we use the insights from integrable theories to derive a closed-form expression corresponding to a part of the diffusion coefficients, which arises due to the presence of non-trivial interactions. While we require that the system exhibits convection, i.e. ballistic transport, for some degrees of freedom, we make no assumption about the integrability structure, which makes our theory applicable to systems such as non-integrable anharmonic chains [23], non-integrable cellular automata [8,24], and general many-body systems with at least one conserved quantity [25], as for example quantum fluids with translational invariance symmetry [26,27].
Our result is based on the power series expansion of the current operator within the hydrodynamical cell in terms of conservation laws. We show that the diffusion constant is related to the projection of the current onto the second-order term in this expansion. This contribution to diffusion should be understood as a consequence of the dispersion of convective modes in thermal or nonthermal ensembles. On the level of the second-order term in this expansion, the Euler scale equations for the eigenmodes of the system can be used to derive a closed-form expression for the diffusion constant. Our result can be connected with the recent advances on the computation of diffusion constant in quantum integrable models [9,10,28], and demystifies the "magic formula" from [29], which connects the diffusion constant to the curvature of the self-Drude weight. Furthermore, we show that our expansion saturates the diffusion constants in integrable theories, which were previously derived using different methods, relying on their integrability structure. A natural connection can also be made with a lower bound on the diffusion constant in terms of quadratically extensive quantities [30], and in terms of the curvature of Drude weights [31].
Our approach is reminiscent of the non-linear fluctuating hydrodynamics (NLFHD) theory [23,27,32,33], which also employs the expansion of the currents in terms of conservation laws. One of the main achievement of NLFHD is to show that Kardar-Parisi-Zhang [34] or Lévy super-diffusive universal transport arises whenever a quadratic coupling of the modes in the second-order expansion is present. Generalization of this result manifests itself within the framework of our theory as a divergence of the diffusion constant in the presence of non-zero coupling of the modes with degenerate velocities. Hydrodynamics on the super-lattice: The starting point of our hydrodynamical approach is to replace microscopic degrees of freedom with macroscopic modes describing the behavior on large space-time scales. This will enable us to identify the relevant degrees of freedom which contribute to diffusive transport. We are interested in systems on infinite chains i ∈ Z; however, our approach can also be understood as a lattice regularization of an emerging field theory. The lattice is decomposed into the fluid cells of length ∆x = . The operators on the hydrodynamical lattice are given as sums of (quasi-)local operators o i within the fluid cell o(x) = Assuming that in the large time limit many-body systems locally equilibrate [35], the only relevant information is contained in the dynamics of conservation laws. Therefore we can define the set of maximum entropy en-arXiv:1911.01995v1 [cond-mat.stat-mech] 5 Nov 2019 sembles, where the equilibrium has been reached within every hydrodynamical cell where β = (β 1 , β 2 , · · · ). Here q i (x) are the modes, i.e. the densities of local conserved quantities Q i , with i = 1, . . . , n c , within the hydrodynamical cell. We also use the Einstein's repeated indices summation convention . Note that the summation x;∆x over the super-lattice is preformed with steps of size ∆x = . Note that in the expressions without derivatives the homogeneous background is assumed β i (x) = β i . The N -point connected correlation function is defined by [36] where the coefficients a i (x) are set to 0 after taking the derivatives. Expectation values of operators between distinct fluid cells vanish with respect to any exponentially clustering density matrix [37] where we introduced the normalized N -point con- For a given density matrix ρ(β) it will also be important for the following to define a set of normal modes n i [33]. Normal modes n i are linear combinations of the densities of conserved charges, as expressed by the rotation matrix R j i = R j i (β), q i = (R −1 ) j i n j . They form localized wave-packets on top of the density matrix ρ(β) which move with distinct velocities v k on the super lattice. At the Euler scale, the Fourier normal modeŝ n i (k, t) = x;∆x e −ikx n(x, t) satisfy the continuity equation with ω i (k) = v i k. Its explicit representation in the real space then reads n i (x, t) = ∆x 2π π/∆x −π/∆x dke ikxn i (k, t). The normalization of R is set by requiring orthonormalizatoin of normal modes (RCR T ) ij = δ ij , where we introduced the susceptibility matrix C ij = q i (0)q j (0) c n . In equation (4) there can be corrections of the order O( −1 ), due to the spreading and diffusion of the normal modes wave-packets, which are absent at Euler scales. We shall show in the following that in order to compute the diffusion constants only the Euler scale dynamics of the normal modes is needed (see also [38]). Transport coefficients by current expansion: The current operator is determined through the continuity equation ∂ t q l,i (t)+j l,i+1 (t)−j l,i (t) = 0. The main idea of our derivation is to perform a hydrodynamical expansion of the current operator δj k = j k − j k on the super-lattice in terms of the local and quasi-local conserved charges, which holds on the level of correlation functions with respect to some stationary density matrix ρ(β). Note that the subtracted current average j k in the definition of the variation of the current δj k is taken with respect to the homogeneous background β k (x) = β k , while the homogeneous limit of δj k is taken only after taking the derivatives with respect to the expectation values of the charges q i (y) = q i (y) . Importantly, the velocities v i and the rotation R of normal modes are obtained by diagonalizing the matrix (∂ qi(x) j k (x) ). The remainder R can contain two types of additional terms which contribute to diffusion. The first one arises due to the presence of non-local conserved quantities whose norms scale quadratically Q 2 c ∼ L 2 with the system size L, and will be discussed bellow in relation to the lower bound on diffusion [30]. The second contribution corresponds to the gradient of (quasi-)local charges R = M j i (q j (x + ∆x) − q j (x))/∆x, and can be set to 0 by imposing P T (parity-and time-reversal) symmetry to the currents and charges [28]. Importantly, stochastic or driven systems may in principle allow for such terms due to the breaking of PT symmetry, see [39].
To demonstrate the utility of our expansion (5) we employ it to obtain the Drude weights, namely the coefficients parametrising the ballistic direct conductivities [40,41], which are defined as D kl = lim t→∞ (2t) −1 t −t ds x j k,0 (s)j l,x (0) c , for some mode k and l. In terms of the hydrodynamical currents the matrix of Drude weights reads and can be represented using the hydrodynamical expansion as The only contribution to the Drude weight (6) that remains finite after the hydrodynamical expansion (5) in the hydrodynamical limit → ∞ corresponds to the linear order, which reproduces a well-known result [42,43] In the normal mode basis the Drude weights read where we introduced extensive normal modes N i = R j i Q j . The physical interpretation of this result is the following. The current operator at the origin j l (0) creates excitations on top of the density matrix ρ(β), see also Fig. 1. The only stable excitations, i.e. the excitations that survive the hydrodynamical limit → ∞, and which are able to reach the current operator at point x, j l (x) are the densities of conserved charges. In the normal mode basis excitations have a well defined velocities v k and gives a contribution to the Drude weight when it reaches the current operator at j l (x). We now focus on the object of primary interest: the Onsager (diffusion) matrix corresponding to the late time sub-ballistic contribution to the current j k (0) generated in the system after the equilibrium state is perturbed by a constant gradient of the charge q l . The Onsager coefficient is related to the diffusion constant via Einstein's relation [28,44] The second order term in the expansion (5) produces a finite convective contribution L c kl to the Onsager matrix. The contribution is once again easily understood in terms of the dynamics of normal modes (reminiscent of a kinetic-like picture [10]). At time t the nonzero contribution is produced by excitations created by the currents at the origin, that reach the operator j(x). The weight of contribution is given by the overlap of two normal modes and the current, see FIG 1. We obtain (see [37] for details) a remarkably simple expression for the convective Onsager matrix where we introduced the renormalized coupling coeffi-cientG expressed in terms of the quadratic matrix G jk is rotated by R's when changing basis from charge densities q i to normal modes n i , while the renormalization |v i − v j | corresponds to the rescaled time t/∆x for which the two normal modes n i and n j overlap with the current at position x. Similarly as for the Drude weight the convective Onsager matrix can be nicely represented in a closed form in terms of normal modes (see [37] and also [38]) Lower bounds on diffusion: The diagonal elements of Onsager matrix L c kk are generally expected to be a lower bound for the exact ones. There have been several proposals in the past years to lower bound the diffusion coefficients. Here we establish a direct connection between our expression (10) and two recent results.
The first bound on diagonal diffusion coefficients (1) was derived in [30].
Here v LR is the Lieb-Robinson velocity and Q is a conserved quantity which scales quadratically Q 2 c ∼ L 2 with the system size L. A set of such conservation laws can be obtained by multiplying two local integrals of motion Q = i≥j α ij N i N j , where we choose traceless N i . Provided that all of the quadratically extensive quantities are of this form, we can show that our expression supersedes this lower bound L c kk ≥ L s kk . Plugging a generic Q expressed in terms of squares of normal modes and maximizing it with respect to the free parameters α ij , we obtain a lower bound [37] for the full discussion). On the other hand the diagonal elements of convective Onsager matrix L c kk in (12) can be lower bounded by choosing |v i − v j | ≤ 2v LR , which demonstrates that L c kk ≥ L s kk . However, additional quadratically extensive quantities, not given by product of local and quasi-local charges, can in principle exist. These terms, included in the extra contribution R in (5), would give an extra contribution to the full Onsager coefficients L kk .
The second lower bound on the spin/charge diffusion constant in the zero magnetization sector/half-filling has been proposed in [31] and further studied in [12,45]. It corresponds to the curvature of the Drude weight with re-spect to the magnetization S z filling ν(T, h) = 4T S z as L s,s ≥ ∂ 2 ν D(h)/v LR . Taking the second derivative of the normal-mode representation of the Drude weight (8), we obtain the expression which is reminiscent of (12). Noticing that magnetization is the eigenmode N i = m with zero velocity v i = 0, and that the only non-zero matrix elements in L c m,m include exactly one instance of the magnetization normal mode, we reproduce the lower bound expression by replacing the renormalization |v i | of the matrix elements (12) by the Lieb-Robinson velocity v LR . Diffusion in integrable systems: Following recent developments in the hydrodynamic description of integrable systems [10,28,[46][47][48][49][50][51], the transformation R as well as other transport matrices such as A, B, and G can be computed exactly. In integrable systems normal modes are the stable quasi-particle excitations that fully describe the thermodynamics of the system [52,53]. Their properties in stationary states are described by their density ρ tot i and occupation function A simple computation yields a neat expression for the G-matrix in the quasi-particle basis, labelled by the set of indices i = (θ, a), where θ denotes the quasi-particle momentum and a the quasi-particle type (see [37] for its derivation): The dressed scattering phase shift T dr = (1 − T n) −1 T can be thought of as a length of the jump of the quasi-particle upon scattering with another quasi-particle, if both of them are immersed in a thermal background [10,50]. The convective coefficients (11) in the normal mode basis read and fully reproduces the exact Onsager matrix in integrable models L c ij ≡ L ij obtained in [9,28] (see also [54]) via the exact form-factors expansion, and is also equivalent to the expression obtained for the classical hard-rod gases [55,56].
The operatorial expansion of the current (5) in terms of conserved charge densities q i can be thought of as a generalization of the thermodynamic form factor expansions in terms of particle-hole excitations in integrable systems [28]. In this representation, a single particle-hole pair excitation on top of the equilibrium state yields the Drude weights, and the two particle-hole pairs give the full diffusion constant [28]. Analogously to the absence of the higher order contributions in the hydrodynamical expansion, more than two particle-hole pairs give no contribution to diffusion. Since L c ij ≡ L ij we can conclude that diffusive transport in integrable models is given purely by the dispersion of their ballistic modes.
Our result also gives a simple explanation of the curious "magic formula", which was established rigorously at infinite temperatures [29] and numerically at finite temperatures. The "magic formula" relates the curvature of self-Drude weight D self s,s = dt j s (0, t)j s (0, 0) c where j s is the spin (or charge) current, with the spin diffusion constant D s,s = L s,s /C s,s = ∂ 2 ν D self s,s (ν). This is now readily understood by noticing that the curvature of the self-Drude weight is equivalent to equation (10) with q j = S z and j k = j s , and that in integrable models L c kl ≡ L kl . NLFHD and super-diffusion: Much like our approach, NLFHD is based on expanding the expectation values of the extensive currents in terms of the local conservation laws up to the second-order and adding phenomenological diffusion constants and white noise terms related by the fluctuation-dissipation relation [23,27,32,33]. Consequently, when there is a single conserved charge, we obtain a noisy Burgers' equation that is known to belong to the KPZ universality class with dynamical exponent z = 3 2 [34,57]. This implies that the exact diffusion constant is divergent L = ∞. Similar behavior can be identified in the presence of multiple modes, and detailed analysis reveals a plethora of new super-diffusive universality classes arising in the presence of appropriate multi-mode couplings [33]. This behavior can also be understood within our framework, since the self-coupling term corresponds to the presence of nonzero matrix elements with degenerate velocities (11), resulting in the divergence of the convective Onsager coefficients (10). More specifically, the diffusion constant of a mode diverges whenever the mode coupling matrix G jk i has non-zero elements for at least two modes with degenerate velocities. Such situations were typically exempted from past NLFHD applications, where it was usually assumed non-degeneracy of velocities, but they are expected to be important in integrable chains, where super-diffusion can also be observed for the spin or charge degrees of freedom close to half-filling [58][59][60][61]. Conclusion: We have introduced an operatorial expansion of the currents in a many-body system in terms of hydrodynamical densities of conservation laws in generic stationary states. We have shown that the second-order terms of this expansion give rise to a finite contribution to the Onsager matrix and therefore to the diffusion constants of the system. Our framework unifies previous results on diffusive transport in one dimension and shows that in integrable systems convective contributions saturate the exact diffusion constants, demonstrating that other mechanisms are absent. Nevertheless, convective contributions account only for a part of the full diffusion constant in classical probabilistic dynamical systems [8,24], and are completely absent in the spin chain with strong dephasing [39]. These partial results call for clarification of which contributions are non-vanishing in certain systems and whether generic Hamiltonian systems support non-trivial quadratic charges that are not included in the convective part of diffusion.
For the case of integrable models, we have provided another non-trivial verification of the expression for the diffusion constants and proved the so-called magic formula. Moreover, we have shown how the diffusion constant of a certain charge can diverge whenever there is degener-acy in the velocity of the modes. While this mechanism is valid for generic systems, it might prove to be useful to explain the emerging KPZ super-diffusive dynamics of spin and charge in quantum and classical isotropic chains [29,62,63]. In these models the spin/charge mode at half-filling has zero velocity and it constitutes an accumulation point for all other modes with finite velocities [64] that approach zero as v i ≡ v θ,a ∼ a −1 [45,62] for any θ. An interesting direction for future research is to develop these ideas on a more rigorous footing.
While we here explored the effects of convective modes on the Onsager matrix, they can also contribute to the level of higher-order cumulants of the currents. Our construction indicates that the contributions are hierarchically ordered with respect to the order of the cumulant, analogously to the BBGKY hierarchy in the Boltzmann equation [65,66]. Finally, our expansion can also be applied in higher dimensional systems, where normal diffusion is typically found. It would be interesting to compare our convective contribution with known results [67].

S1. CORRELATIONS ON SUPER-LATTICE
In order to identify the effects of multi-point correlation functions on transport coefficients we consider an infinitely large spin chain, which we divide into the parts of equal length ∆x = , i.e. the hydrodynamical cells. We will denote the macroscopic coordinate which determines the position on this supper-lattice by x.
The hydrodynamical density o(x) is given by in terms of the quasilocal operator Quasilocality means that the spectral norm is upper bounded by for some α, ζ > 0, and the C * norm • . The operators o [i,i+k] , are supported on the sublattice [i, i + k], and act as an identity operator everywhere else. The multi-point connected correlation function is defined as with Z(a 1 (x), · · · ) = tr (exp(log(ρ) + a i (x)o i (x)) in the limit a i (x) = 0. We are considering a set of states ρ which satisfy exponential clustering where we assumed the ordering i ≤ j ≤ k ≤ l. Note that for any two operators a and b a trivial inequality holds true Using the properties above we can show that two hydrodynamical operators decay exponentially with the distance if |x − y| > 2 , i.e. the operators are not located in the neighboring cells. In order to demonstrate the property (S7) we can divide the connected correlator into two contributions In order to upper bound the second term s 2 we can use the trivial bound (S6) and quasilocality (S3) In order to evaluate the first term, we take into account that |x − y| > 2 , implying that the minimal distance between the operators in o 2 and densities o 1 is l 2 . This enables us to show that Now we proceed to show that the contribution from the neighboring cell is sub-polynomial in hydrodynamical scale for arbitrary κ. We divide the sum into three parts with 0 < κ 2 < κ 1 . Similarly as before we can lower bound the second sum by using a trivial lower bound (S6) and quasilocality Using similar arguments as before, we get Using a trivial bound (S6), we can upper bound s 3 by Since κ = 2κ 1 + κ 2 can be arbitrarily small, we arrive at the result (S11).
We are now in position to prove extensivity and orthogonality of arbitrary multipoint connected correlation function on the super-lattice The N -point connected correlation function is at most extensive, since the expectation value o is proportional to the volume of the cell ∝ . And the N -point correlation function corresponds to the N − 1 point derivative of the expectation value. In the absence of phase transitions divergences are absent, implying extensivity of N -point correlation function.
We proceed to show that if not all of the operators neighboring positions, connected correlation function vanishes exponentially. This can be easily demonstrated by noticing that the clustering property is not altered by adding local quantities to log(ρ), i.e. the equations (S11) and (S7) hold even if a i (x) = 0. Let's assume that o i (x) and o j (y) do not occupy the same or the neighboring cells. Using (S7) we have that a 1 (x), ..., a n (x)) ), (S16) and taking the derivatives produces at most polynomial factor N −2 . Similarly the derivative of the upper bound (S11) produces the correction of the same order.

S2. DERIVATION OF DRUDE WEIGHTS
In the derivation of Drude weight and diffusion kernel, we will use the dynamics of normal modes on Euler scale, which reads up to 1 corrections. In order to establish how the current expansion works in actual computations, we will apply it to the computation of Drude weight. The Drude weights D i,j are defined by For computing this object, we need the first term in the hydrodynamic expansion (5) only. Insert this term into (6) gives where we used the relation ∂ β i (x) = C ij ∂ qj (x) , which follows from the clustering property (S15). Alternatively, changing to the normal mode basis with the convention RCR T = 1, one can write the Drude weight as where N i = R j i Q j is the total charge in the normal mode basis. Note that in the computation involving the current expansion (5), we always take the homogeneous limit of the averages only at the end of computations.

S3. DERIVATION OF THE ONSAGER MATRIX
In this section, we present a derivation of the convective contribution to the Onsager matrix L c u,v , and hence diffusion constant. Recall that the Onsager matrix is given by the following expression where D u,v corresponds to the Drude weight. In order to better understand how each term scale with , let us rewrite it in terms of the hydrodynamic current j We will take the hydrodynamic limit → ∞ only in the end. In order to study corrections to the Euler scale hydrodynamics, which corresponds to the nonzero Drude weight, and can be interpreted as a consequence of the first term in the expression for the current (5). It is useful to consider a part of the current that characterizes the sub-Euler contribution The Onsager matrix now reads Inserting the expression for the current into the equation (S22) we obtain where, in the fifth and sixth lines, we used and defined for an arbitrary hydrodynamic operator o. To proceed, let us first deal with a building block M o ij and rewrite it in terms of normal modes. Using the solution of n(x, t) (S17), we have Focusing on the leading contribution when → ∞, we rescale the variables k → ∆xk and k → ∆xk , in terms of which M o ij becomes where we also used the fact that we are considering the almost-homogeneous limit 1, hence we can assume translation invariance in correlators. We further perform the rescalings x → x ∆x , t → t ∆x , which gives The integration over t can be done as follows which allows us to obtain a compact expression of M o Notice that in the final line, the hydrodynamic observable o(x, t) is replaced by the ordinary local observable o(x, t).
We further observe that the curvature term ∂ 2 o(0,0) ∂β i ∂β j can be written as according to which (S25) becomes where H matrix corresponds to H ij v = ∂ qj (0,0) ∂ qi(0,0) j v (0, 0) . The G-matrix is defined by We are finally in the position to derive the exact convective contribution to the Onsager matrix. Putting everything together, we have Note that the convective Onsager matrix can also be written as where The H-matrix can be expressed explicitly in terms of correlation functions using and reads where

S4. HIGHER ORDER CORRECTIONS
Here we elaborate on possible higher order contributions to the Onsager matrix from the hydrodynamical expansion of the current. For simplicity we will consider only the third order correction steaming from a single charge It will be clear that higher order contributions vanish in the hydrodynamical limit → ∞. Following the same steps as in the derivation of the Onsager matrix from the second order, we obtain the contribution with Going to the normal mode basis, we obtain After rescaling k → ∆xk, k → ∆xk , k → ∆xk , and x → x ∆x , t → t ∆x , we obtain Summation over x and integration over time produces Let's assume that v i ≥ v j ≥ v r . Integrating over k , we get Integration over k and k finally yields Obviously the contribution vanish in the limit → ∞, provided that the degeneracies are absent.

S5. DERIVATION OF THE G-TENSOR IN INTEGRABLE SYSTEMS
Here we compute the G-tensor for integrable systems directly from known generalised hydrodynamics expressions [43] . As in the main text we use the Latin indices i = (θ, a) to denote the pairs of quasi-momentum θ and the particle type a, wherefore the calculations boil down to simple matrix-like manipulations. To reiterate, the G-tensor reads We first normalize R according to RCR T = 1. Let us define the normalizationN so that R =N (1 − nT ). Then, recalling that C = (1 − nT ) −1 ρ(1 − n)(1 − T n) −1 in the quasi-particle basis in GHD, where ρ i = ρ tot i n i , we realize that if we chooseN asN where χ i = ρ i (1 − n i ) is the quasi-particle susceptibilities, then R is properly normalized as well as diagonalizing A. Then the G-matrix becomes

S6. QUADRATIC LOWER BOUND IN NORMAL MODE BASIS
Here we relate the lower bound in terms of quadratic charges derived in [30] to the convective Onsager matrix L c kk . For simplicity, we restrict the discussion to the infinite temperature state, however the generalization to finite temperatures should be possible by invoking exponential clustering property.
The first step in the derivation is to generalize the lower bound to multiple charges, by considering the norm of the operator O = 1 where we take the normal modes N i and the extensive current operator J = x j x to be traceless and Hermitian. This gives us the lower bound Note that the three point function in the above expression is indeed equivalent to the three point connected correlation function due to the tracelesness of the operators, however the connected correlation involving four copies of normal modes corresponds to the two point connected correlation function of the terms in the brackets. In order to maximize the contribution on the right hand side we take a derivative with respect to α * ij , which produces the set of equations Now we proceed to prove that where δ ijkl is 1 only if all of the four indices are the same and 0 otherwise. We also assumed that i ≥ j and k ≥ l. This equality follows from the general consideration where the quasilocal extensive operators are assumed to be traceless. Note that the limit L → ∞ is considered before sending the support of operators to ∞. In order to evaluate the above sum, we divide it into two parts. The first part corresponds to the case in which at least three of the densities overlap at least at one point. This immediately implies that all four densities overlap at some point. If the largest support of the density in the above sum is r β = max(r 1 , r 2 , r 3 , r 4 ), then the absolute value of the sum over α k , k = β of such terms can be upper bounded by tr (1) ≤ r 3 β exp(−η(r 1 + r 2 + r 3 + r 4 )), η > 0, where we introduced the maximal inverse quasilocality length η = max(η 1 , η 2 , η 3 , η 4 ) of the operator densities appearing in the above equation. Furthermore we can upper bound r β < (r 1 + r 2 + r 3 + r 4 ), implying that the summing over r 1 , r 2 , r 3 , r 4 yields a finite contribution. The only summation that remains is the one over β. This produces a factor which is proportional to the system size L. Secondly, let's consider a product of two norms of operators