Full Counting Statistics in the Transverse Field Ising Chain

We consider the full probability distribution for the transverse magnetization of a finite subsystem in the transverse field Ising chain. We derive a determinant representation of the corresponding characteristic function for general Gaussian states. We consider applications to the full counting statistics in the ground state, finite temperature equilibrium states, non-equilibrium steady states and time evolution after global quantum quenches. We derive an analytical expression for the time and subsystem size dependence of the characteristic function at sufficiently late times after a quantum quench. This expression features an interesting multiple light-cone structure.


I. INTRODUCTION
The statistical nature of measurements of observables is a fundamental principle of quantum mechanics. Measuring the same observable in identically prepared systems leads to different measurement outcomes that are described by a probability distribution that depends on both the state |Ψ and on the observable O considered. The full probability distribution P (O, |Ψ ) encodes detailed information about quantum fluctuations in the system. It is of particular interest in situations where the first few moments do not provide a good description of the distribution. Quantum mechanical probability distributions in the guise of Full Counting Statistics (FCS) have been studied for some time in mesoscopic devices [1,2]. More recently it has become possible to analyze them in systems of ultra-cold atomic gases [3][4][5][6][7][8]. This has broken new ground in the sense that one is dealing with (strongly) interacting many-particle systems and a variety of observables, typically defined on subsystems, can be accessed. This has motivated a number of theoretical works of FCS in equilibrium states [9][10][11][12][13][14][15][16][17][18][19][20], and after quantum quenches [7,[20][21][22][23][24]. A second motivation for studying FCS has been the observation that in non-interacting fermionic systems with particle number conservation the FCS of particle number within a subsystem is directly related to the entanglement entropy [25][26][27][28][29][30][31][32][33][34] and provides indirect information about the latter.
From a theoretical point of view calculating the FCS for a given observable on a sizeable subsystem poses a formidable problem and as a result only very few exact results are available even in simple equilibrium situations. Even less is known about FCS after quantum quenches. This motivates reconsidering FCS in the transverse field Ising chain (TFIC). The TFIC is a key paradigm for quantum phase transitions [35] and a simple, but non-trivial, many-body system without particle number conservation and therefore provides an ideal playground for studying FCS both in and out of equilibrium. Indeed, thanks to the mapping of the TFIC to a model of non-interacting spinless fermions with pairing term it is possible to analytically determine ground state and thermal properties, see e.g. [35][36][37], as well as describe the non-equilibrium dynamics of local observables [38][39][40][41][42][43][44] and of the reduced density matrix of a block of adjacent sites [43,[45][46][47] after a global quantum quench. A summary of these developments is given in the recent reviews Refs [48,49].
In this work we focus on the FCS of the simplest observable, the transverse magnetization within a block of adjacent spins. In the ground state this problem has been previously analysed in Refs [9,11] and generic Gaussian states have been considered as well [14]. We note that the ground state FCS of the longitudinal magnetization at the critical point has been determined in Ref. [10] and the ground state FCS of the subsystem energy was considered in Ref. [18].
This manuscript is organised as follows. In section II we first introduce the TFIC and briefly summarize the important steps for diagonalizing the Hamiltonian. We then define the FCS and the associated generating functions considered in this work. In section III we provide a novel derivation of an efficient determinant representation for the FCS in general Z 2 invariant Gaussian states. The result is equivalent to that of Ref. [14]. This result is applied in section IV to the determination of the FCS in equilibrium states. In the ground state we recover the results of Ref. [9]. Our results for the FCS in finite temperature equilibrium states are to the best of our knowledge new. In section V we turn to the main point of interest: the time evolution of the FCS after a global quantum quench. We consider the situation where the system is prepared in a pure state at a finite finite energy density and then time evolved with a Hamiltonian H that does not commute with the initial state density matrix, which leads to non-trivial dynamics. We present explicit results for general "transverse field" quenches as well as evolution starting in a classical Néel state. The main result of this work is presented in section VI: an analytic expression for the time evolution of the FCS after a transverse field quench. In section VIII we summarize our results and comment on a number of issues that deserve further investigation.

A. Transverse Field Ising chain
In the following we consider the spin-1/2 transverse field Ising model on an infinite chain The ground state phase diagram features ferromagnetic (h < 1) and paramagnetic (h > 1) phases that are separated by a quantum critical point in the universality class of the two-dimensional Ising model [35]. The order parameter that characterizes the transition is the longitudinal magnetisation GS|σ x j |GS . At finite temperature spontaneous breaking of the Z 2 symmetry of H(h) is forbidden and hence the order present in the ground state at h < 1 melts. In order for this paper to be self-contained we now briefly summarize the relevant steps for diagonalizing the Hamiltonian (1). A more detailed discussion can be found in e.g. the Appendix in [42]. The TFIC is mapped to a model of spinless fermions by a Jordan-Wigner transformation where c j are fermion operators obeying canonical anticommutation relations {c † j , c k } = δ j,k . Setting aside the issue of boundary conditions the Hamiltonian takes the form This is diagonalized by a Bogoliubov transformation where {α k , α † p } = δ p,k and the Bogoliuobov angle is The Hamiltonian takes the form where the dispersion relation is given by The ground state of H(h) is equal to the Bogoliubov vacuum state defined by

B. Full Counting Statistics and Generating Function
We are interested in the properties of the smooth and staggered components of the transverse magnetization of a chain segment of length . These are defined as Given a density matrix ρ that specifies the quantum mechanical state of our system, the probability distributions for the transverse subsystem magnetizations are given by In the following we will focus on the characteristic functions of these probability distributions, defined as By construction, the expansion of χ (u,s) (λ, ) around λ = 0 generates the moments of the associated probability distribution. The following relations are readily inferred from the definition of χ (u,s) (λ, ) These properties imply where we have defined the weights

III. GENERATING FUNCTION FOR A GENERAL GAUSSIAN STATE
In this section we show how to obtain the generating function (11) for a general Gaussian state with a novel method that is however equivalent to the one used in [14].
Our starting point is the realization that (11) depends only on the reduced density matrix of the block A of adjacent spins where we have introduced the auxiliary "density matrices" Here the "partition function" Z (a) ensures the normalisation Tr ρ (a) = 1. A fundamental property that we will exploit in the following is that both ρ A and ρ (a) are Gaussian operators in the fermionic representation of our problem, cf. section II A. Hence they are univocally determined by the correlation matrices of the fundamental fermionic operators [50][51][52]. Moreover, the trace of the product of Gaussian operators such as (15) can be expressed in terms of the associated correlation matrices [53]. This is a very useful property, see e.g. Ref. [47] for a related application, that forms the basis of our analysis.
In order to proceed we need to specify a convenient basis of operators. This is provided by Majorana fermions related to the lattice spin operators by The Majorana fermions satisfy the algebra They are related to the Jordan-Wigner fermions (2) by a 2l−1 = c † l + c l and a 2l = −i(c † l − c l ). As we are dealing with Gaussian density matrices we can follow Refs. [50][51][52] and Wick's theorem to express ρ A in terms of the subsystem correlation matrix Γ A nm Γ A nm = Tr [ρ a m a n ] − δ nm , 1 ≤ m, n ≤ 2 .
As the Pauli matrices form a basis in the space of operators over C 2 the reduced density matrix of a subsystem A that consists of neighbouring spins at sites i = 1, . . . , can be expressed in the form where α i = 0, x, y, z. We now restrict our discussion tor density matrices that are invariant under the Z 2 transformation P σ z l P = σ z l , P σ x,y l P = −σ x,y l .
In this case the Jordan-Wigner strings cancel and the reduced density matrix (RDM) is mapped to an operator expressed in terms of Majorana fermions acting on the same spatial domain We note that the case where P ρP = ρ can be dealt with by the method set out in Ref. [47]. The RDM (22) can be written in an explicit Gaussian form as where W is a skew symmetric 2 × 2 hermitian matrix. Using Wick's theorem the matrix W can be related to the correlation matrix (19) tanh The auxiliary density matrices ρ (u,s) (16) can be expressed in the Majorana basis in a completely analogous way. The corresponding 2 × 2 correlation matrices Γ (u,s) are given by The only non-vanishing matrix elements are This implies that Γ (u,s) are block-diagonal, e.g.
where σ y is the 2 × 2 Pauli matrix.
We are now in a position to write down a convenient determinant representation for the generating functions χ (u,s) (λ, ). To do so we employ a relation derived in Ref. [53]: given two Gaussian density matrices ρ 1,2 with correlation matrices Γ 1,2 the trace of their product is given by Applying this relation to our case we arrive at the following determinant representations where Γ A and Γ (u,s) are given in (19) and (26), (27) respectively.

A. Simplification in special cases
Equation (29) has been derived for a general Z 2 -invariant Gaussian state with density matrix ρ. If the state is also invariant under translations and reflections with respect to a site the generating function χ (u) (λ, ) can be simplified further. Indeed, under these conditions, the correlation matrix assumes a block Toeplitz form [45,51] where g l = Tr a 2n a 2n+2l−1 = −Tr a 2n−1 a 2n−2l , Taking advantage of the block diagonal form of the correlation matrix of the auxiliary density matrix in (27) we can cast the generating function in the form where Γ is a block Toeplitz matrix

B. Expressions for the first few cumulants
The determinant representation (29) of the generating function provides an efficient way for determining the cumulants of the probability distribution, which is the main purpose of the function itself. The cumulants are obtained in the usual way from the series expansions of ln χ (u,s) (λ, ) The first few terms of the series expansion are where we have definedΓ The first three cumulants are Specifying to the case of density matrices ρ that are invariant under translations and reflections around a site we have Tr where F and G are the × Toeplitz matrices For the first three cumulants we obtain It is straightforward to generalise these considerations to higher cumulants because Tr Γ n can always be written as the sum of the traces of products of F and G.

IV. FULL COUNTING STATISTICS IN EQUILIBRIUM
In this section we analyze the generating function χ (u) (λ, ) obtained from (29) and the associated probability distribution in equilibrium configurations. We first consider the ground state FCS, which has been previously studied by Cherng and Demler in [9]. We then turn to the FCS in finite temperature equilibrium states, which to the best of our knowledge has not been considered in the literature.

A. Full counting statistics in the ground state
In the ground state the generating function is of the form (32), (33) with entries where θ k is the Bogoliubov angle (5). By rearranging rows and columns,Γ can be brought to a block diagonal form with × matrices G and G T (41) on the diagonal and zero otherwise. This allows us to express the generating function as This is precisely the result previously obtained by Cherng and Demler [9] by a different technique. They considered the generating function The Toeplitz determinant (45) can be analyzed by standard methods [9]. The symbol τ (e ik ) of a block Toeplitz T with elements (T ) ln = t l−n is defined through the equation The symbol of the block Toeplitz matrix (45) is given by As long as the symbol has zero winding number a straightforward application of Szegő's Lemma gives, For h < 1 and λ > λ c (h) the winding number of the symbol is 1 and the above result gets modified accordingly [11].
For a detailed analysis we refer to Ref. [11]. The full counting statistics for the transverse magnetization in the entire system was studied in [54] and the result is identical to (49). Thus considering the subsystem instead of the entire system only makes a difference for h < 1 and λ > λ c (h), as discussed in [9]. We note that all cumulants can be obtained from (49) since they are defined by the expansion close to λ = 0. Consequently, the first cumulants are given by The non-zero values of C 3 and C 4 show that the probability distribution is non-gaussian.

B. Full counting statistics at finite temperature
We now turn to the FCS in finite temperature equilibrium states, for which we are not aware of any results in the literature. In this case, the correlation matrix has the same structure as for the ground state, but now where k is the dispersion relation (7). Since f l = 0, the same simplifications as in the ground state case apply and the generating function can be expressed as In Fig. 1 we show P (u) w (m) for subsystem size = 20 and several different temperatures. We employ a log-linear plot in order to make the deviations of the probability distributions from a Gaussian form (which would correspond to a parabolic form) more apparent. We can see from Fig. 1 (a) that the temperature dependence for h < 1, corresponding to the ferromagnetically ordered phase at zero temperature, is not very pronounced. In contrast we see a much stronger temperature dependence in the paramagnetic phase, cf. Fig. 1 (b). At low temperatures the probability distribution is as expected asymmetric as a result of the applied field and is seen to display an even/odd structure. The latter disappears quickly as temperature is increased, whereas the asymmetry remains until the temperature exceeds the scale set by the magnetic field. In Figs 2 and 3 we show the skewness and excess kurtosis of the probability distribution as a function of subsystem size for a range of temperatures. These are defined as the thermal expectation values Both skewness and excess kurtosis are non-vanishing for finite β and , which establishes that the distribution is not Gaussian. A very peculiar feature is that at fixed skewness and excess kurtosis are non-monotonic functions of the temperature. Furthermore, we observe that at a fixed temperature they both tend to zero as the subsystem size is increased. This signals that the corresponding probability distribution approaches a Gaussian. This is expected as for large subsystem sizes the laws of thermodynamics apply and the probability distribution is then approximately Gaussian with a standard deviation that scales as √ .
• We initialize the system in the Néel state |↑↓↑↓ . . . ↑↓ , thus breaking translational symmetry by one site. This symmetry is restored at late times after the quench and it is an interesting question how this is reflected in the probability distributions of observables.
A. Transverse field quench h0 −→ h In this quench protocol both the Hamiltonian and the initial state are translationally invariant. The characteristic function has the determinant representation (32), (33) with [45] where Using Szegő's Lemma it is straightforward to obtain the large-asymptotics in the initial (t = 0) and stationary (t = ∞) states. The t = 0 result corresponds to a ground state at field h 0 and has been discussed earlier.

Behaviour in the stationary state
The late time asymptotics of the generating function can be determined from Szegő's Lemma. For quenches into the paramagnetic phase h > 1 it takes the form The O( −1 ) corrections also follow from Szegő's Lemma. The real and imaginary parts of corrections included) are shown for a transverse field quench from h 0 = 5 to h = 2 and subsystem size = 100 in Fig. 4. For quenches into the ferromagnetic phase and λ < λ c (h 0 , h), Eq. (58) continues to hold. However, for λ > λ c (h 0 , h) the symbol exhibits non-zero winding number and the analysis needs to be modified, cf. Appendix A. The probability distribution in the stationary state is obtained by Fourier transforming χ (u) (λ, , t). Examples for several transverse field quenches are shown in Fig. 5. We again employ a logarithmic scale to make the deviations from a Gaussian form more apparent. In Figs 6 we plot the skewness and the excess kurtosis of the steady state probability distributions for a number of transverse field quenches. We observe that in all cases both skewness and excess kurtosis tend to zero for large subsystem sizes. This signals that the probability distributions approach Gaussians in the large-limit. While the steady states are non-thermal now, they still exhibit finite correlation lengths. Employing the same arguments as for finite temperature ensembles then implies that the cumulants of S z u ( ) are proportional to in the large-limit. This in turn suggests that skewness and excess kurtosis should scale as −1/2 and −1 respectively, while the standard deviation scales as 1/2 . These expectations are in perfect agreement with our findings.

Scaling collapse
At finite times the FCS and the probability distribution can be computed efficiently from the determinant representation (32). Importantly we observe that for sufficiently large values of and t there is scaling collapse 1.
The property (59)  For quenches towards the ferromagnetic regime the scaling collapse for general values of λ can be significantly worse, and then really only emerges at rather large subsystem sizes , cf. Fig. 9. Like in the case of the stationary state discussed above, c.f. section V A 1, there exists a critical valueλ c (h 0 , h) of the counting parameter such that for λ <λ c (h 0 , h) the scaling collapse is excellent, while for λ >λ c (h 0 , h) no collapse is observed at the times and subsystem sizes of interest here. For the cases we have consideredλ c (h 0 , h) coincides with λ c (h 0 , h), which is the value of the counting parameter above which the symbol has non-zero winding number. We note however, that in cases like the one shown in Fig. 9 the generating function itself is extremely small and will not give a significant contribution to the corresponding probability distribution.

Time dependence of the probability distribution
There are basically four different kinds of transverse field quenches and we now consider them in turn. 1. Quenches within the ferromagnetic phase. For such quenches the probability distribution remains very narrow and approximately Gaussian throughout, cf. Fig. 10. For the parameters considered the average relaxes quickly towards its stationary value.
2. Quenches within the paramagnetic phase. Here the initial probability distribution exhibits an even/odd structure. This can be understood by doing perturbation theory around the large h 0 limit, cf. Appendix B. After the quench the mean of the probability distribution broadens and shifts towards smaller values of m. The alternating structure is initially preserved but then gets smoothed out. At late times P (u) w (m, t) is well described by a Gaussian.
3. Quenches from the paramagnetic to the ferromagnetic phase. Here the probability distribution is initially peaked at a large value of m and displays an even/odd structure. At later times it broadens and becomes smooth, while relaxing towards its stationary profile in an strongly oscillatory manner. 4. Quenches from the ferromagnetic to the paramagnetic phase. In this case the probability distribution shows very little variation in time and remains narrow and approximately Gaussian throughout the evolution. It was pointed out in Ref. [62] that the return amplitude exhibits a nonanalyticity at some finite time t * after the quantum quench. This phenomenon was termed a "dynamical phase transition". Local operators are known to be insensitive to this phenomenon [41][42][43]. We have investigated the behaviour of P (u) w (m, t) in the vicinity of t * but have not observed any unusual effects. We conclude that the probability distribution for the smooth subsystem magnetization in the transverse field direction is also insensitive to the "dynamical phase transition".

B. Quench from the Néel state
We now turn to the time evolution of P (u,s) w (m, t) when the system is initialized in the Néel state |ψ 0 = |↑↓↑↓ . . . ↑↓ . This explicitly breaks translational invariance by one site, but retains invariance under translation by two sites. As a result the subsystem correlation matrix is now a 4 × 4 block-Toeplitz matrix where we have assumed the subsystem size to be even and Here the various two point functions are given by In the following we will determine the characteristic functions where again we have defined According to our general discussion in section III they have determinant representations of the form

Behaviour in the stationary state
We first consider the probability distributions for a finite subsystem of even size in the late time limit. As we will now show, the stationary state for the Néel quench is locally equivalent to an infinite temperature state. To see this we first note that the energy of the Néel state is It is easy to see using their explicit representation in terms of spins [49] that the expectation values of all higher conservation laws also vanish ψ 0 |I (n,±) |ψ 0 = 0.
This in turn implies that the conserved Bogoliubov mode occupation numbers are given by These characterize an infinite temperature equilibrium state. We conclude that the system will relax locally [49] to an infinite temperature steady state at late times after the quench. Using this observation it is then straightforward to work out the probability distributions P (u) (m, t = ∞) of S z u ( ) = j=1 σ z j and P (s) (m, ∞) of S z s ( ) = j=1 (−1) j σ z j . As shown in the introduction we have As we are dealing with an infinite temperature state, we may calculate P w (r) by using a grand canonical ensemble and working in the simultaneous eigenbasis of the σ z j 's. This reduces the calculation of P w (r) to the combinatorial problem of how many eigenstates there are for a given eigenvalue of S z u ( ) or S z s ( ). This is easily solved in terms of the binomial distribution The result (70)   this peak broadens and relaxes towards the Gaussian profile (70). When quenching to the ferromagnetic phase, cf. Fig. 14, an additional feature emerges: an even/odd structure evolves at short times after the quench.
The probability distribution of S z s ( ) is useful for investigating the restoration of the translational symmetry. In the initial state P (s) w (m, t = 0) features a single peak at m = − /2, which is a characteristic fingerprint of the classical Néel state (in z-direction). We first discuss quenches into the ferromagnetic phase. Here at short times after the quench

P (s)
w (m, t) develops an even/odd structure and broadens significantly. The average of the probability distribution oscillates strongly in time and decays very slowly to its stationary value, which is a Gaussian distribution centred around m = 0. This shows that translational symmetry is restored very slowly.
The behaviour for quenches into the paramagnetic phase is broadly similar. An even/odd structure develops at early times, but is less pronounced that for quenches to the ferromagnetic phase. The average of P (s) w (m, t) again oscillates strongly around m = 0, but is seen to relax much more quickly than for quenches to the ferromagnetic phase. Approximate translational symmetry gets restored more rapidly.

VI. ANALYTIC RESULTS FOR THE PROBABILITY DISTRIBUTION
We now restrict our discussion to the particular case of transverse field quenches. As we have seen above, in this case the characteristic functions χ (u) (λ, , t) exhibit a scaling collapse at late times, cf. (59). This suggests that it might be possible to obtain analytic results for the late time asymptotics by a suitable generalization of the multi-dimensional stationary state approximation method previously used to determine the asymptotics of the order parameter two-point function [41] and the entanglement entropy [46]. As we will see, such a generalization is indeed possible, even though the case at hand is significantly more complicated.
Our starting point is the following expression ln χ (u) (λ, , t) = ln (cos λ) + 1 2 Tr (ln(1 − tan λ Γ )) , which is derived from (32) by using the identity ln (det (A)) = Tr (ln (A)). The second term in (72) can be expanded This then leads us to examine integer powers (Γ ) n of the correlation matrix. Unlike in the case of the order parameter two-point function analyzed in [42] odd powers do not vanish because Γ is not a real anti-symmetric matrix. The symbol t (k) corresponding to the correlation matrix Γ is defined by Its explicit expression for a magnetic field quench from h 0 to h iŝ where θ k and ∆ k have been previous defined in (57). Following Ref. [42] we can represent the trace of powers of the correlation matrix as multiple integrals where we have defined k 0 ≡ k n and We now change variables The integration ranges in the ζ variables is determined by the constraints The integral over ζ 0 can now be carried out as the integrand does not depend on it. This gives where µ({ζ}) is the size of the range of ζ 0 under the constraints (79) µ( ζ) = max 0, min A. Multi-dimensional stationary phase approximation For large values of the integrals can be carried out using a multi-dimensional stationary phase approximation. As the symbol is independent of ζ j the stationarity conditions for the ζ j 's implies that the leading contribution to (80) derives from the region We may thus replace k j with k 0 everywhere except in rapidly oscillating terms in the symbol such as e 2iε(kj )t . In [42] this procedure was referred to as localisation rule. As in [41] application of this rule gives Obtaining a closed form expression for F ( k) is however much more involved than for the order-parameter two point function studied in Ref. [42]. We conjecture that application of the localization rule to F ( k) results in Here the sum is over all partitions of the set of integers {0, 1, . . . , n − 1} into three sets A 1 , A 2 and A 3 , where the number of elements in A 1 is constrained to be even. The size of the set B = {b 1 , b 2 , . . . } is denoted by |B| and we have defined Finally, sign(A) is the sign of the permutation required to bring the (integer) elements of the set A into ascending order. We have explicitly checked (84) for 1 ≤ n ≤ 15 but have not been able to find a rigorous proof for it. We now use the identity (for even k) to rewrite the time-dependent factors in (84). This gives where (A) r is the r'th element of the set A and Application of the localization rule to (80) hence results in an expression of the form In the next step we carry out a multi-dimensional stationary phase approximation for the 2n − 2 integrals over ζ 1 , . . . , ζ n−1 and k 1 , . . . , k n−1 . We will assume that there is a single saddle point and use where σ A the signature of the matrix A (i.e. the difference between the numbers of positive and negative eigenvalues), which is the Hessian of the function q evaluated at the saddle point In our case the saddle point conditions are where γ A,k = 4t (−1) The Hessian A is a matrix of the form and hence we have det (A) = −4 1−n and σ A = 0. The value of µ( ζ) at the saddle point for a given sequence The saddle point approximation thus gives The leading contribution to the final integral can then also be determined by a stationary phase approximation. This shows that all terms with |A1∪A2| r=1 (−1) pr = 0 are suppressed at late times by a factor of 1/ √ t. Conversely, the leading contribution to χ (u) (λ, , t) at late times arises from terms with |A1∪A2| r=1 (−1) pr = 0, which requires |A 1 | + |A 2 | to be even.

Structure of µ( ζ (0) )
At this point it is useful to investigate the structure of µ( ζ (0) ) for a given term in the multiple sum over p 1 , . . . , p |A1|+|A2| in more detail. For simplicity we focus on a particular example In this case we have where v k0 = ε (k 0 ) is the group velocity of Bogoliubov fermions at momentum k 0 and Θ(x) is the Heaviside step function. The step function in (97) is reminiscent of the light-cone structure found for two point correlation functions of local operators [66][67][68][69] and entanglement entropies [45,70,71]. Repeating the above exercise for leads to the result All other cases can be worked out analogously and lead to Heaviside step functions Θ( − 2m |v k0 | t) with m ∈ N 0 .
In principle one could determine further terms f n,0 but their contribution turns out to be negligible for all cases we have considered. The contributions f n,m>0 (λ, k 0 , t) are more difficult to simplify. While the term f 0,1 can still be obtained without further approximations, in order to obtain closed form expressions for m > 1 we have resorted to an expansion in powers of sin(∆ k0 ). This is expected to give very accurate results for small quenches, which are defined as producing a small density of elementary excitations through the quench [41,42]. The leading terms are then conjectured to be of the form f 0,1 = −i tan ∆ k0 ln 1 + ie iθ k 0 cos ∆ k0 tan λ 1 + ie −iθ k 0 cos ∆ k0 tan λ , As we will see below, the contributions described by (101) and (102) are sufficient to obtain an extremely accurate description of χ (u) (λ, , t). The constant C can be fixed by comparing the t → ∞ limit of (100) to the result obtained previously for the behaviour in the stationary state. For later convenience we define two approximations as where a = 1, 2 and where we set f 2,1 = 0. The structure of the integrand in our result (100) is reminiscent of that found in connected two-point correlation functions [41] and entanglement entropies [46]. In the latter quantities it gives rise to a "light-cone" behaviour in developing connected correlations and the spreading of entanglement respectively. In contrast to these cases the expression (100) involves an infinite number of "light-cone structures" with velocities that are integer multiples of the maximum group velocity. Since the generating function involves complicated sums over multi-point correlation functions on the interval [1, ] this is not in contradiction with the celebrated Lieb-Robinson bound [73]. The light-cone structure in connected two-point correlators and entanglement entropies can be understood in terms of simple semi-classical quasi-particle pictures [45,66]. It would be interesting to develop an analogous understanding for the novel structure observed in the generating function, but this is beyond the scope of the present paper.

VII. ACCURACY OF THE ASYMPTOTIC RESULT
Our analytic result (100), (101), (102) gives the leading contributions in the space-time scaling limit [42] , t → ∞, /t fixed. An important question is how good this asymptotic result describes the behaviour of χ (u) (λ, , t) at small and intermediate times and subsystem sizes. In order to answer this question we now turn to a comparison between our analytical results (103) and a direct numerical evaluation of the determinant representation (32), (33). The numerical errors in the latter are negligible.

A. Small-λ regime
A representative comparison between the analytical results χ (u) 1,2 (λ, , t) for small values of λ and numerics is shown in Fig 18 and 19. We see that χ (u) 1 (λ, , t) reproduces the numerics very well at late times after the quench. In contrast, the oscillatory behaviour at short times is clearly not captured. The improved approximation χ (u) 2 (λ, , t) (103) is seen to be in excellent agreement with the numerics.
By construction the oscillatory part of the analytic result is most accurate over the entire range of the "counting parameter" λ when sin ∆ k0 is small, i.e. for small quenches. For quenches where sin ∆ k0 is no longer small we still find excellent agreement between the analytic and numerical results as long as tan(λ) is small. This can be understood by noting that for such values of λ the infinite sum in (73) is dominated by the first few terms, i.e. small values of n. On the other hand, higher orders of sin ∆ k0 only emerge for larger values of n. Therefore the leading order result in sin ∆ k0 already provides a very good approximation in the small-tan(λ) regime even when sin ∆ k0 is not small. This observation is of significant practical importance: As shown in Fig. 20 in a particular example |Reχ (u) (λ, , t)| and |Imχ (u) (λ, , t)| are largest in the vicinity of λ = 0 (except at short times). This implies that the corresponding probability distribution, which is the object we are ultimately interested in, will be dominated by the small-λ regime. As a consequence (100), (101), (102) provide a good approximation for the calculation of P  to h = 1.5. We observe that the characteristic function is small unless λ is small. The behaviour for quenches within the ferromagnetic phase and quenches between the phases is similar.

B. Large-λ regime
In the large-λ regime we have to distinguish between the cases where the symbol in the stationary state has zero or non-zero winding number, c.f. section V A 1. The first case covers quenches to the paramagnetic phase. Here we find that our analytic result is again in good agreement with numerics. The second scenario applies to quenches to the ferromagnetic phase and λ > λ c (h 0 , h). We have shown in section V A 2 that there is no good scaling collapse in this regime of counting parameters for the moderate subsystem sizes and times of interest here. It should therefore not come as a surprise that the asymptotic result does not provide a good approximation in this regime. Presumably (100), (101), (102) no longer hold in this regime because the analytic continuation of the power series expansion of the logarithm (73) becomes non-trivial in this case. In practice the failure of the analytic approach to give a good account of the generating function in this parameter regime is irrelevant as χ (u) (λ, , t) itself is extremely small and makes a negligible contribution to the probability distribution. As shown in Fig. 21, the main contribution to the latter, which after all is our object of interest, arises from the small-λ regime of the generating function, which is well approximated by our analytic expressions.

C. Relative errors
In order to provide a more quantitative discussion of the quality of the approximate results (103) we consider the relative errors where χ (u) 1,2/num (λ, , t) are respectively the analytic approximations (103) and the result of the numerical computation of the determinant representation (32), (33). In Fig. 22 we plot the time dependence of the relative errors for a  maximal value of sin(∆ k0 ) within the domain of integration approximately 0.54, which means that higher orders in f 1,1 can be important. As we have argued above, this will be the case if tan(λ) is not small. In Fig. 22 (a) λ = 0.1 is taken to be small, and the quality of both approximations χ (u) 1,2 (λ, , t) is seen to be excellent. In Fig. 22 (b) the counting parameter λ = 1.4 is taken to be large. This leads to a significantly larger error, which is however still fairly small and also decays in time. We see that the analytic results provide a good approximation for all values of λ.
We now turn to a parameter regime, in which our analytic results no longer provide a uniformly good approximation for all values of the counting parameter λ. Fig. 23 shows results for a quench from h 0 = 0.2 to h = 0.8. The maximal value of sin ∆ k0 in the integration range is now 0.71 so that higher orders in f 1,1 can again be important. For small values of λ the relative errors of both analytical approximations are small and decreasing in time. On the other hand χ (u) 1,2 (λ, , t) cease to provide accurate approximations for large values of λ with λ > λ c (h 0 , h) as can be seen in Fig. 23 (b). However, we want to stress once more that χ (u) (λ, 200, t) itself is extremely small in this parameter regime and makes only a negligible contribution to the probability distribution.

D. Probability distributions
An asymptotic expansion for the probability distribution P (u) w (m, t) can be obtained by Fourier transforming the generating function, cf. Eq. (14). As expected on the basis of the discussion above, we find that the analytic result becomes very accurate at sufficiently late times for all quenches. At intermediate and short times we still find excellent agreement between the analytical and numerical results for quenches originating the ferromagnetic, see e.g. Fig. 24 (a). For quenches from the paramagnetic phase the analytic result is an excellent agreement with numerics at short and intermediate times as long as the quench is "small". In practice this covers all quenches within the paramagnetic phase as long as h is not very close to 1. For other quenches the corrections to the f 1,1 term in (102) will become significant at short and intermediate times.

VIII. CONCLUSIONS
We have analysed the full counting statistics of the transverse and staggered magnetization of a subsystem in the thermodynamic limit of the transverse field Ising chain. We derived a convenient determinant representation for the corresponding generating functions χ (u,s) (λ, , t). We first considered the FCS in equilibrium states and showed that the probability distributions are always non-Gaussian except in the limit of infinite subsystem size at finite temperature. We determined the temperature and field dependence of the generating function as well as the first few cumulants. We then moved on to the main focus of our work, the calculation of the FCS after quantum quenches. We considered two quench protocols: transverse field quenches and evolution starting from a classical Néel state. We first determined the FCS in the stationary states reached at late times. The probability distributions are again non Gaussian, except in the limit of infinite subsystem size. We analyzed the time evolution of the probability distributions P (u,s) (m, t) for a variety of quenches by numerically evaluating the exact determinant representation for the generating function (the numerical errors incurred are negligible). For transverse field quenches originating in the paramagnetic phase P (u,s) (m, t) showed interesting smoothing and broadening behaviour in time. In contrast, P (u,s) (m, t) displayed a simpler behaviour for quenches originating in the ferromagnetic phase. In the case of a Néel quench P (s) (m, t) encoded detailed information on the restoration of translational invariance. The numerical approach provided us with evidence for the existence of a scaling regime for the generating function in which we observed data collapse according to the scaling form (59). This is turn allowed us to proceed with the derivation of the main result of our work: the analytic expression (100) for the FCS after transverse field quenches in the space-time scaling limit t, → ∞, t/ fixed. This was achieved by a substantial generalization of the multi-dimensional stationary phase approximation method of Refs [42,46]. We performed a careful comparison of our analytic results to numerics (that has negligible errors) and found excellent agreement on the level of the probability distributions for all cases considered. We observed that the expression for the generating function exhibits an interesting multiple light-cone structure that has no analog in either correlation functions of local observables [41] or in the entanglement entropy [46]. An interesting open question is whether this structure can be understood in terms of the kind of semiclassical quasi-particle picture that has been successfully employed to explain the main features observed in the dynamics of both entanglement [45] and correlations [66].
Our work provides the first analytic results for FCS after quantum quenches and hopefully will pave the way for further studies. Here we have focussed on the FCS for the transverse magnetisation. It would be very interesting to determine the FCS for the longitudinal magnetisation, which is the order parameter characterising the Ising quantum phase transition. A more straightforward but interesting extension would be to study certain observables in free fermion models with long-range hopping and/or pairing [72,74,75]. Similarly, the probability distribution of the (smooth) subsystem magnetisation in the spin-1/2 Heisenberg XXZ chain should be calculable both at finite temperatures [76] and in the stationary states after certain quantum quenches [77][78][79][80][81][82][83][84][85]. For quantum quenches in the regime where bosonization provides a good approximation [86] the full time evolution of the probability distribution for certain observables can be obtained in a straightforward way. Finally, the case of integrable chains of higher spin could be studied both in equilibrium [87,88] and after a quench [84,89].

(A2)
Here T (τ ) denotes an infinite Toeplitz matrix with symbol τ . In the case where the block-size is 1, this reduces to the Szegő limit theorem ln det (T ) = The large asymptotics of Toeplitz determinants in cases where the symbol τ has winding number ±1 is given by [42,91] ln det (T ) = We have seen that the probability distributions P (u,s) (m, t) exhibit an even/odd structure in m for short times after quenches starting in the paramagnetic phase. In this appendix we show that this structure can be understood in perturbation theory around the h → ∞ limit. For simplicity we consider the probability distribution P (u) (m) in the ground state at h 1. In the limit h → ∞ the ground state is the saturated ferromagnetic state along the transverse field direction The corresponding probability distribution is a delta function at m = − /2. The other eigenstates of j σ z j are denoted by |n (0) . The leading correction to the generating function arises at second order in perturbation theory in H 1 = j σ x j σ x j+1 . The relevant corrections to the ground state are Substituting this into the expression for the generating function gives In order for (0) n|H 1 |0 (0) to be non-zero the product state |n (0) must have precisely two overturned spins. Let us denote their positions by j and j + 1. For ≥ 2 we then have This is seen to exhibit an even/odd effect as the corrections for m = /2 mod 2 are proportional to the subsystem size.