Full Counting Statistics of Charge in Quenched Quantum Gases

Unless constrained by symmetry, measurement of an observable in a quantum system returns a distribution of values which are encoded in the full counting statistics. While the mean value of this distribution is important for determining certain properties of a system, the full distribution can also exhibit universal behavior. In this paper we study the full counting statistics of particle number in one dimensional interacting Bose and Fermi gases which have been quenched far from equilibrium. In particular we consider the time evolution of the Lieb-Liniger and Gaudin-Yang models quenched from a Bose-Einstein condensate initial state and calculate the full counting statistics of the particle number within a subsystem. We show that the scaled cumulants of the charge in the initial state and at long times are simply related and in particular the latter are independent of the model parameters. Using the quasi-particle picture we obtain the full time evolution of the cumulants and find that although their endpoints are fixed, the finite time dynamics depends strongly on the model parameters. We go on to construct the scaled cumulant generating functions and from this determine the limiting charge probability distributions at long time which are shown to exhibit distinct non-trivial and non-Gaussian fluctuations and large deviations.


I. INTRODUCTION
Symmetry and universality are two of the most powerful concepts in theoretical physics.The former can dramatically reduce the complexity of systems, places stringent constraints on allowed physical processes and gives rise to conservation laws via Noether's theorem.The latter instead explains why vastly different systems can display near identical features and how simple physical principles can underpin many seemingly complex phenomena.These concepts have been extensively studied in the context of closed quantum systems which are close to equilibrium leading to the discovery of many ubiquitous properties and the development of numerous powerful and widely applicable techniques of analysis.In recent years however questions of universality and its emergence in far from equilibrium systems have come to the fore and in this context one-dimensional integrable models have been widely studied [1][2][3][4][5][6][7][8].These models possess special symmetry properties endowing them with an infinite number of mutually commuting conserved charges thereby placing strong constraints on their dynamics [9].At the same time these properties facilitate exact analytic solutions of the models through Bethe ansatz techniques allowing in-depth analysis of their thermodynamic properties [10].
Despite this analytic control however, uncovering the non-equilibrium properties of integrable models remains challenging and is still a highly active area of research.Nevertheless many universal features have been established, in particular concerning the dynamics within a subsystem where it has been shown that the system relaxes locally to a stationary state described by a generalized Gibbs ensemble (GGE) [1][2][3][4][5][6][7].Building upon this a number of exact techniques have been developed including the quench action method [11,12] which allows one to determine this GGE explicitly and generalized hydrodynamics (GHD) which describes the long time and large scale dynamics of an inhomogeneous state.[13][14][15].
A significant driver of interest in the non-equilibrium dynamics of integrable models has been the advent of numerous experimental platforms which allow for the simulation of isolated many-body quantum systems with a high degree of accuracy and control [16][17][18].Chief amongst these are ultra cold atomic gas setups which have the ability to faithfully simulate integrable systems.The Lieb-Liniger model [19] of interacting bosons, the Gaudin-Yang model [20,21] of interacting fermions and the sine-Gordon field theory [22] are all well known integrable models which can be simulated within cold atom experiments [17,18,[23][24][25][26][27][28][29][30][31][32][33][34].The relaxation of such models to GGEs [35] as well as the vailidty of GHD [36][37][38] has been observed and tested extensively in such experiments ensuring that despite their apparent fine tuned nature integrable models are the appropriate description of these systems.
In this work we shall study the interplay between nonequilibrium dynamics of integrable models and the fluctuations of their charge within a subsystem in the context of cold atom gas experiments.We do this by calculating the full counting statistics (FCS) of the particle number in the Lieb-Liniger (LL) and Gaudin-Yang (GY) models.The systems shall be taken out of equilibrium by quenching them from an initial Bose-Einstein condensate (BEC) state which is not only experimentally relevant and conceptually simple but also allows us to characterize the dynamics exactly.We shall show that there exists a simple relationship between the cumulants of charge in the initial state and the GGE to which these systems relax at long times.This relationship is in fact universal, relying only on the presence of integrability and not on details of the model.However, while the cumulants in the final state are fixed by the initial state their decomposition in terms of quasi-particle excitations differs from model to model meaning that the finite time dynamics will differ.Specifically, we find that the connected m-point function The quasi-particle picture of charge fluctuations within a region A. The quench causes pairs of correlated quasi-particles with equal and opposite momenta, k, −k to be emitted from the initial state.The pairs propagate throughout the system, carrying with them the charge.A quasiparticle pair's contribution to charge fluctuation differs depending on whether both members of the pair (solid lines) or a single member (dotted lines) are inside A. At short times the fluctuations inside A are given by pairs while times which are long compare to the subsystem size only a single member of a pair can contribute.
of the charge, N in a subsystem, A at the initial and final times are related by a factor of 1/2 m−1 , that is lim We provide a general derivation of this result below along with an explanation in terms of the qusaiparticle picture [39,40](see Fig. 1) and explicit checks in our models of choice.Moreover, since the leading order (linear in subsystem size) behavior of the cumulants of charge fluctuations in the GGEs are the same in all cases they define the same limiting coarse-grained or continuous probability distributions (PDs) which can be explicitly determined.Nevertheless, the microscopic PDs (retaining information about sub-leading corrections) can be different.This fact arises from the transmutation of particle statistics in the different parameter regimes of each model and manifests in a simple and intuitive way.For example, in the strongly attractive Fermi gas the fermions are bound into tight pairs forming hard core bosons with double the charge, a fact which is reflected in vanishing probabilities for observing an odd number of particles in a subsystem.The rest of the paper is arranged as follows: In Sec.II we review our object of study, the full counting statistics and some basic properties of cumulant generating functions and their associated probability distributions.We also briefly recall some properties of integrable models and prove a general relationship between the FCS calculated in a GGE and in the diagonal ensemble.In Sec.III we introduce the models we shall study as well as the particular initial states of interest.In Sec.IV we determine the cumulants of charge in Lieb-Liniger model for both repulsive and attractive regimes.In Sec.V we carry out an analogous calculation in the repulsive and attractive regimes of the Gaudin-Yang model.In the penultimate section we study the analytic properties of the scaled cumulant generating functions for both our models.We use this knowledge to determine the limiting charge probability distributions in both the initial and final states.In the last section we summarize our work and discuss open questions and potential future directions.

A. Full counting statistics and charge probability distributions
We wish to characterize the fluctuations of the particle number in the LL and GY models which have been quenched from BEC states.While these charges are conserved quantities for the full system, when we restrict to a subsystem they exhibit nontrivial fluctuations and dynamics .To study this we calculate the full counting statistics of the charge within a region A of length ℓ.We take ℓ to be large but much smaller than the full system, 1 ≪ ℓ ≪ L where L is the full system size and in the end take the thermodynamic limit ℓ, L → ∞.For a charge operator N whose restriction to the subsystem is denoted NA the FCS are defined as where ρ A (t) = tr Ā[ρ(t)] is the reduced density matrix of A at time t and ρ(t) = |Ψ 0 (t)⟩⟨Ψ 0 (t)| is the density matrix of the full system which has been quenched from the intial state |Ψ 0 ⟩ → |Ψ 0 (t)⟩ = e −iHt |Ψ 0 ⟩.While the FCS of charge and related quantities like the work distribution have been studied previously in certain interacting integrable models [52,[62][63][64][65][66][67], recent developments have allowed to study these quantities in far from equilibrium systems as well as their associated current statistics [54,[68][69][70][71][72][73][74][75][76].
From the FCS we can determine many properties of the dynamics of the particle number fluctuations within the subsystem and its interplay with the relaxation of the system to its long time steady state.Specifically, we can calculate the connected correlation functions of the charge via the cumulants of the FCS, for m ∈ N and where ⟨•⟩ c means the connected part of the correlation function e.g.
Moreover, by continuing β → iβ we can obtain the charge probability distribution P (n, t), which gives the probability that a measurement of NA at time t returns the value n.The latter is a natural quantity to study from the experimental point of view and particularly so in the case of cold atom experiments.While the expectation value of an observable such as an order parameter is a central quantity to understand the nature of a system, the full distribution of measurement outcomes is also enlightening and can unveil universal behavior.In cold atom experiments large numbers of measurements are required to be performed and so the full probability distribution is naturally obtained [77][78][79][80].
By our analytical methods which we shall introduce shortly, the computation of the cumulants or the cumulant generating function log Z β (A, t) is only feasible if the length of the subsystem ℓ is also infinitely large.In particular, the moment generating function for an extensive quantity in a large subsystem can be generally written as where C s (β, t) is the scaled cumulant generating function (SCGF) defined as and whose derivatives w.r.t.β are the scaled cumulants denoted by κ s m .Therefore it is the scaled cumulants and their generating function, i.e., quantities with extensive scaling w.r.t. the subsystem size, which can be computed.From the SCGF, we can obtain the large ℓ limiting probability distribution (PD) for the charge fluctuations via the rate function I(z) through where ≍ is understood as and we assumed the same scaling for the time variable t.Importantly the rate function I(z) can be obtained as the Legendre-Fenchel transform of the SCGF according to the Gärtner-Ellis theorem.The rate function also governs the large deviations of the limiting coarse-grained probability distribution, which is a continuous distribution even if the microscopic PD is discrete.Whereas in a strict sense our methods allow for the computation of I(z) via the SCGFs, we shall also use Eq. ( 4) and approximate Z iβ (A, t = 0, ∞) at the initial and the steadystates by exp (ℓC s (iβ, t = 0, ∞)).Doing so retains the discrete nature of charge fluctuations and allows for the onset of some interesting microscopic effects, which vanish in the ℓ → ∞ limit but are expected to be observed in large but finite subsystems.

B. Review of integrable models and the Thermodynamic Bethe Ansatz
In the subsequent sections we shall explicitly calculate the scaled cumulants κ s m (t) as well as SCGF C s (β, t = 0, ∞) in the LL and GY models however before this we shall derive some generic properties which will allow us to relate the scaled cumulants in the initial state, denoted by κ s m (0), and in the GGE, denoted by κ s m (∞).To do this we keep things general and recall briefly some basic properties of integable models: Integrable models posses an extensive number of conserved quantities charges whose associated operators Q(k) , k = 1, 2, . . .commute with the Hamiltonian.These imbue the model with a stable set of quasi-particle excitations which are indexed by a discrete species index n = 1, . . ., N s and parameterized by a continuous rapidity λ These states are also simultaneous eigenstates of all the conserved charges of the model, namely, The first three conserved charges are chosen to coincide with the particle number N = Q(0) , momentum P = Q(1) and Hamiltonian Ĥ = Q(2) .For simplicity we shall denote the charge, momentum and energy of a quasi-particle of species n and rapidity λ by q n = q (0) n (λ) and ϵ n (λ) = q (2) n (λ).In the thermodynamic limit and at finite density the model can be treated using the methods of the thermodynamics Bethe Ansatz (TBA) [10].Within this approach, the quasiparticle content of a stationary state of the model can be encoded in the distributions ρ n (λ), ρ h n (λ), ρ t n (λ) and ϑ n (λ) which are respectively the distribution of occupied quasiparticles of species n, the distribution unoccupied quasiparticles, the total density of states ρ t n (λ) = ρ n (λ)+ ρ h n (λ) and the occupation function ϑ n (λ) = ρ n (λ)/ρ t n (λ).All of these functions are related to each other by the Bethe-Takahashi equations which take the form where (•) ′ denotes differentiation w.r.t.λ and * is the convolution f * g(y) = dyf (x − y)g(y).Here T nl (λ, µ) is the scattering kernel which characterizes the scattering between quasi-particles of type n and l with rapidities λ and µ.In the cases of interest this is symmetric in the species index T nl (x) = T ln (x).
When the state contains many excitations the bare quasi-particle properties become dressed due to the interactions.These dressed quantities are denoted by (•) dr and satisfy the integral equations Explicitly, the dressed charge is obtained from the solution of The quasi-particle velocity on the other hand is given by the ratio of dressed quantities, v n (λ) = ϵ ′ dr (λ)/p ′ dr (λ).This can be recast into the equation where we have used ρ t (λ) = p ′ dr (λ).These sets of integral equations typically need to be solved numerically but in certain cases can be solved analytically as is the case for the LL model discussed below [81][82][83].
C. Full counting statistics using TBA Having reviewed these basic properties and set up our notation we can now move on to determine the relationship between the FCS at t = 0 and in the GGE.We begin with the latter and define where we have neglected the o(ℓ) contributions and have introduced the generalized Gibbs ensemble Herein it is necessary to include both local and semilocal conserved charges in the GGE in order to obtain the correct description of the long time state [84][85][86][87][88].The Lagrange multipliers above, β (k) , are determined by matching the expectation values of Q(k) in ρ GGE to those in the initial state.To compute C s GGE (β) with β ∈ R we utilize the result of Ref. [68] (see also [59,89]) claiming that where f GGE denotes the free energy density of a GGE characterized by the chemical potentials β (k) and β shifts one of the chemical potentials that corresponds to the conserved quantity under investigation.We can rewrite the above equations in the typical manner of the thermodynamic Bethe ansatz by exchanging the trace for a path integral over the distributions ρ n (λ), ρ h n (λ) [10], where and S Y Y is the Yang-Yang entropy which is proportional to ℓ as well and which counts how many microstates |λ⟩ correspond to the same macrostate given by the distributions ρ n (λ).Evaluating this functional integral via a saddle point approximation we obtain with η n,β being given by while η n,0 is evaluated in the same way at β = 0 and can be related to the occupation function of the GGE via ϑ n (λ) = 1/(1 + η n,0 (λ)).In order to evaluate C s GGE (β) explicitly one must know the functional form of g n (λ) which can be done by making use of the quench action method.Using this one can show that it is obtained from the extensive part the squared overlap between an eigenstate |λ⟩ and |Ψ 0 ⟩ in the thermodynamic limit [90].Namely, For the states that we consider the left hand side is known explicitly and the function g n (λ) is simply extracted from this.
We now turn to the calculation of the FCS in the initial state.In this instance, as discussed further in Appendix A, the initial value can be obtained by considering the diagonal ensemble.After expressing this as a path integral we find where the factor of 1 2 in front of S Y Y arises from the fact that the overlap with our initial state is non-zero only for parity invariant eigenstates which have half the entropy contribution.Here we have denoted C s (β, 0) by C s 0 (β), which we express more explicitly as with ηn,β being given by log ηn, Either comparing Eqs. ( 19) and ( 23) or the TBA systems ( 20) and ( 24), we then find that where C s GGE is given by (20).Therefore comparing with (16) we also find that the cumulants are related as The above relations are one of the main results of this paper.We note that while ( 27) is a technical result whose proof relies upon the details of the TBA and quench action formalism, it also admits a very intuitive interpretation based upon the quasi-particle picture of entanglement dynamics [39,40], characterized in Fig. 1, which proceeds as follows.The initial state of the system can be expressed as a collection of correlated pairs of quasiparticles, of opposite momenta, which are excited by the quench.A pair of quasi-particles can then be viewed semi-classically as emerging from a single point in space.
Its constituents propagate ballistically in opposite directions throughout the system thereby spreading correlations over a wider region.The fluctuations of charge within the subsystem at t = 0 are therefore governed by the distribution of pairs of quasiparticles which necessarily carry charge 2 and are indexed by the absolute value of the momenta.At long time, the quasiparticles have spread throughout the system and it is not possible for both members of an initially correlated pair to be inside the subsystem, thus the fluctuations of charge in the subsystem are governed by single quasi-particles of charge 1 and indexed by the momentum rather than its absolute value.Nevertheless, since the quasiparticle occupation function is a constant of the dynamics the distributions of pairs and unpaired quasi-particles are the same.Thus we arrive at the expression (26) where we can recognize the 2β as originating from the charge of a pair of quasiparticles while the factor of 1 2 comes from the fact that one sums over only the absolute value of the momenta.
A notable consequence of the relations is that since the right hand sides of ( 26) and ( 27) are a property solely of the initial state then the left hand side is independent of the model specifics.In the examples below we shall argue and then show explicitly that for BEC initial states κ s m (0) = d ∀m where d is the average density i.e., the charge is Poisson distributed and accordingly which is exactly the same expression as the one obtained from the steady-state GGE using (18) as we shall demonstrate shortly.Evidently (27) does not tell us how the charge distribution evolves between its initial value and the GGE value which necessarily depends on the model.However we can also characterize this in integrable models using recently obtained results on the charged moments and full counting statistics using the method of space-time duality [59,60].In particular, as discussed above the scaled cumulants obey a quasi-particle picture [39,40] meaning that charge fluctuations are spread throughout the system via the ballistic propagation of pairs of quasiparticles which are shared between A and Ā resulting in where we have introduced K m n,x (λ) with x = 0, ∞ being the rapidity and species resolved scaled cumulant in the initial state or GGE, i.e.
Additionally, min[2t|v n (λ)|, ℓ] is the characteristic function which counts the number of quasi-particle pairs shared between A and Ā, while ℓ−min[2t|v n (λ)|, ℓ] can be interpreted as the number quasiparticle from pairs which are contained solely within the subsystem.Thus, within the subsystem the charge fluctuations are governed either by pairs of quasi-particles both of whom are inside A and whose contribution is given by K m n,0 (λ) or single quasi-particles which originated from outside A or whose partner has exited A and which contribute K m n,∞ (λ).The first cumulant returns the expectation value of the charge and is conserved which can be seen using ( 28) and ( 29) at m = 1.The second cumulant gives the charge susceptibility which has a compact expression Higher cumulants become more complicated and involve dressing of lower cumulants.For example, the third cumulant is given by where we see that it depends also on the dressed second cumulant.We omit higher expressions for cumulants which can nevertheless be systematically computed.
Before moving on, we make two brief remarks on the preceding discussion.First, the result (27) relies upon the validity of the diagonal ensemble for calculating the initial state value which is discussed further in Appendix A.
To check the validity of ( 27) we shall compute the scaled cumulants in the initial states by independent means and show that the results of the diagonal ensemble, involving interaction-dependent functions, yield the same result.This is not always the case however and the diagonal ensemble cannot be applied to initial states such that the subsystem is an eigenstate of the charge and hence has no fluctuations.This is the case for the magnetization of the GY model quenched from the BEC state.However for the models and initial states which we consider the diagonal ensemble correctly reproduces the initial FCS of the particle number.Second, although some analytic results exist for the distributions functions which determine (20) and (23), their explicit evaluation is feasible only in certain special limits of our models.In spite of this, at generic interactions, it is straightforward to numerically compute the (first few) scaled cumulants and check the validity of (27).The numerical computation of the SCGF is also possible but instead of comparing the numerically obtained SCGFs against (26) on a strictly finite interval, we shall rather construct the SCGFs at arbitrary interactions as a Taylor series defined by the cumulants, which is justified by the exponentially decaying behavior of the scaled cumulants.

III. MODELS AND SETUP A. Lieb-Liniger model
The standard model for describing one-dimensional interacting bosons in the context of cold atom experiments is the Lieb-Liniger model.The Hamiltonian is given by (32) Here b † (x), b(x) are canonical bosonic operators satisfying [b(x), b † (y)] = δ(x − y) they describe bosons of mass m which interact via a density-density interaction of strength c on a system of length L. We shall consider both the repulsive c > 0 and attractive c < 0 cases.From here on we take m = 1/2 for simplicity and assume periodic boundary conditions.The model has a single U (1) charge, the particle number N = dx b † (x)b(x).
In the repulsive case there is only one species of quasiparticle with bare charge, momentum and energy given by As outlined above, the properties of a stationary state of H are encoded in the distributions ρ(λ), ρ h (λ), ρ t (λ) and ϑ(λ) which are connected via the Bethe equations where the scattering kernel is given by a n (x) = In the attractive case the spectrum of the model is entirely different and there are an infinite number of quasi-particle species corresponding to bound states of n bosons.These quasi-particles have the following bare charge, momentum and energy We denote the distributions of the n-boson bound states by ρ n (λ), ρ h n (λ), ρ t n (λ) and ϑ n (λ).The Bethe equations in this case then take the form where for n > m the two particle scattering kernel is We shall study the dynamics of the system quenched from the BEC state of N particles given by which is the ground state of the model at c = 0 and is also an eigenstate of N .When considering just the subsystem however ρ A (0) contains states in all particle number sectors less than N and the charge probability distribution can be determined by means of a simple argument.Since the wave function is constant the probability of measuring a charge in A for the system with N = 1 is given by ℓ/L.For higher particle number since the state has no spatial correlations the detection of a boson is an independent event with a constant rate ℓN/L and therefore has a binomial distribution.In the thermodynamic limit we take N, L → ∞ with d = N/L the density held fixed and thus we end up with a Poisson distribution with rate ℓd and ⟨N m A ⟩ c 0 = ℓd or κ s m (0) = d, ∀m.From this observation it is possible to construct the reduced density matrix ρ A (0). Introducing the boson on the restricted space b † where |0 A ⟩ is the vacuum inside A. This has the required Poisson distribution of charge and also has the same correlation within the subsystem as |Ψ 0 ⟩⟨Ψ 0 |.The entanglement entropy between A and Ā can also be determined from this expression and originates solely from the charge fluctuations in the subsystem.The von Neumann entanglement entropy is given by Shannon entropy of the charge probability distribution and so scales as S A ≃ 1 2 + 1 2 log 2πℓd.

B. Gaudin-Yang model
For the case of interacting spinful fermions in one dimension we shall study the Gaudin-Yang model which is also widely used to describe cold atom gas experiments.The Hamiltonian is ) where ψ † σ (x), ψ σ (x) are two species of fermionic operators with σ =↑, ↓ which obey {ψ † (x) σ , ψ σ ′ (y)} = δ σσ ′ δ(x−y).These describe spin 1/2 fermions of mass m interacting via a local inter-species density-density interaction of strength c.As before we shall take m = 1/2 and assume periodic boundary conditions but we shall consider both repulsive c > 0 and attractive c < 0 interactions.In this case there are two U (1) conservation laws corresponding to particle number , however we shall only concentrate on the former.
The quasi-particle content of the model depends on whether one is in the repulsive regime or attractive.In the repulsive regime there are an infinite number of species and we denote their distributions by ρ(λ), ρ h (λ), ρ t (λ), ϑ(λ) and σ n (λ), σ h n (λ), σ t n (λ), ϑ n (λ) where n ∈ N. The former types of quasi-particles are associated to spin up fermions, they have unit charge and magnetization and momentum and energy p(λ) = λ, ϵ n (λ) = λ 2 .The latter quasi-particle types, also called strings are associated to the spin degrees of freedom, they carry zero charge, magnetization −2n and have zero energy and momentum.The integral equations for these distributions are where for T nm (λ) was introduced above In terms of these distributions the U (1) charges of a state are given by N = dλρ(λ) and M = N − 2 ∞ n=1 n dλ σ n (λ).In the attractive regime there is an additional type of quasi-particle which is a bound state of two fermions forming a spin singlet.We denote its distributions by ρ(λ)ρ h (λ), ρt (λ), θ(λ).They carry charge 2 and magnetization 0 while their momentum and energy are p(λ) = 2λ, ε(λ) = 2λ 2 .The corresponding integral equations are For the fermionic model we shall again study the dynamics emerging from a BEC state wherein the bosons are formed from pairs of fermions in a singlet at the same point in space, where N L is a normalization factor.
In contrast to the previous case this is not an eigenstate of the model at any c but is an eigenstate of particle number ⟨Φ 0,N | N |Φ 0,N ⟩ = 2N and magnetization ⟨Φ 0,N | M |Φ 0,N ⟩ = 0. Once again upon tracing out Ā, particle number is no longer conserved however the magnetization is still 0. In this instance the wavefunction is not flat due to the Pauli exclusion of the fermions however for large enough subsystem size and at finite density we can apply the same arguments as before and determine that the charge distribution is also Poisson.In particular, it is easy to check analytically the first few cumulants, which, in the thermodynamic limit, in fact yield κ m (0) = 2d if d denotes the density of singlet pairs.Furthermore the initial reduced density matrix has the form In other words, we have again obtained constant cumulants, which can be attributed to a Poisson distribution.It is, nevertheless, important to recall the fact that for any subsystem, the total magnetization is zero, which means that the distribution is only Poissonian for the effective 'quasi-bosons'.Therefore the fermionic distribution is microscopically different, namely it has vanishing probabilities for odd fermion numbers.In the thermodynamic limit, this distribution can have the same cumulants and scales to the same limiting continuous PD, determined by the rate function of the Poisson distribution with density 2d.We shall revisit this feature in Sec.VI.The von Neumann entanglement entropy can be obtained here also and similarly coincides with the Shannon entropy of the Poisson distribution.

IV. FCS IN THE LIEB-LINIGER MODEL
We now turn to the explicit calculation of the SCGF in the Lieb-Liniger model treating the repulsive and attractive regimes separately.The quench dynamics emerging from the state (38) has been studied previously in both the repulsive [81] and attractive cases [82,83] with many properties already determined analytically and with thorough discussion on the steady-state GGE in the repulsive case [91].Following, we shall make use of these exact results to obtain the FCS.

A. Repulsive interactions
The one remaining ingredient that is required to determine the FCS is the overlap function (22).For the particular initial state (38) this has been obtained in [92].In the repulsive case where there is only a single quasiparticle species this is given by As a result the function C s GGE (β) is given by where η β is given by the solution to and as stated above ϑ = 1/(1 + η 0 ).This type of integral equation admits an analytic solution [81] which is given where I m (x) is the modified Bessel function of the 1st kind.From this we can determine exactly both the FCS in the initial state and in the GGE.It can then be checked that in the initial state we recover ⟨N m A ⟩ c 0 = ℓd + o(ℓ).This statement holds true for arbitrary interactions c > 0 as can be confirmed numerically, but the result can be particularly easily seen in the Tonks-Girardeau limit of c → ∞ in which case η β = e −β λ 2 /4d 2 .Accordingly we have that for the initial state in this limit for β ∈ R, i.e., the SCGF of the Poisson distribution which has cumulants all equal.That is, we have recovered the know SCGF of the Poisson distribution at this limit.Switching now to the long time limit we have that the but it is again easy to check numerically that the cumulants are in fact predicted by (28) at arbitrary c > 0. This implies that Eq. ( 20) actually gives the same function 2d(e β/2 − 1) for any c > 0 and for β ∈ R, which can also be confirmed by numerical comparisons.Note that this way we confirmed the predictions of ( 27) and (26) in the repulsive regime of the LL model.From these analytic results we can also write down the full time dynamics of the cumulants in the Tonks-Girardeau limit and find where we used the fact that in this limit v(λ) = 2λ.At finite c since the function K m ∞ (λ) and the quasiparticle velocity depend on the interaction strength the same relations do not hold.Instead, to analyze the finite c and t behavior we plot the first four cumulants using as a function of the scaled time ζ = t/ℓ in Fig. 2 for c = 1 and d = 10.From this we see that while the first two cumulants are monotonic the higher cumulants are not.To understand the non-monotonicity we also plot the rapidity resolved cumulants K m 0 (λ) for d = 1, c = 1.From this we see that K m≤2 0 (λ) are positive functions while K m>2 0 (λ) are negative for some rapidities.By inspecting (29) and using (28) one sees that the positiveness implies that the cumulant is monotonic in time while one which is not can allow for non-monotonic behavior.In addition one can note that higher cumulants take longer to relax to their GGE values (dashed lines).This also can be understood by examining the behavior of the rapidity resolved cumulants.For m > 1 K m 0 (λ) all vanish at the origin but have at least one set of extrema close to it.These extrema are closer to the origin at higher m indicating that the slower quasi-particles contribute more to the higher cumulants resulting in a slower approach to their asymptote.Lastly we note that for m > 1 the number of pairs of extrema present in K m 0 (λ), m − 1, is related to the number of extrema in κ s m (t) as a function of time, namely m − 2.
To investigate the effect of changing the interaction we plot in Fig. 3 the second scaled cumulant κ s 2 (t) as a function of time for c = 1/2, 1, 10. From this we can also see from these plots that for lower interaction strength it takes longer to reach the asymptotic value.This could be anticipated from the fact that at c = 0 the initial state becomes the ground state of the model.To understand this microscopically however we also plot the corresponding K 2 0 (λ) as a function of λ and see that they are peaked closer to the origin for smaller c.Thus the second cumulant is dominated by slower modes for lower interaction strength c resulting in a slower approach to its asymptote.In Fig 4 we plot also the fourth cumulant and see again that it approaches its asymptote slower for lower interaction strength which is a consequence of the fourth cumulant being dominated by slower modes.An analogous effect occurs if one holds c constant and instead changes d.This is related to the anomalous relaxation phenomenon called the quantum Mpemba effect [93][94][95]

B. Attractive interactions
We turn now to the case of attractive interactions.Here there are an infinite number of quasiparticle species and the overlap function for these can be obtained from (51) via The resulting cumulant generating function is where η n,β are determined by These latter equations can be rewritten using some standard TBA identities to a more convenient form [82,83] log where s(λ) = sech πλ 2|c| .As with the repulsive case these admit an analytic solution.It can be shown that with the remaining functions obtained through for n > 1 and where η 0 (λ) = 0.These expressions become cumbersome to treat for n > 1 but simplify considerably, as does the SCGF in the infinite interaction limit |c| → ∞.Taking this limit in the analytic solution one can we find that η 1,β = e −β λ 2 /4d 2 while all other η n>1,β diverge, indicating that they do not contribute to (59).The resulting expression for the SCGF is the same as in the Tonks-Girardeau limit (55) and so the cumulants also agree.This remains true also at arbitrary time meaning that κ s m (t The lack of bound-state contribution in the infinite interaction limit can be understood on energetic grounds, and the fixed energy of the initial state cannot support the formation of higher bound states.Interestingly however, in the same limit the local density-density correlation function only receives a contribution from the two-particle bound states [82,83]. As the interaction strength is decreased bound states of increasing length contribute to the SCGF.In Fig. 5 we plot the second cumulant κ s 2 (t) as a function of time for different interaction strengths.We see that in contrast to the repulsive case the attractive system takes much longer to relax to its asymptotic value.This is a result of the contribution of the bound states to the cumulants being dominated by the slowest quasi-particles.This is seen also in Fig 5 where we plot K 2 n,0 (λ) for n = 1, 2, 3 at |c| = 2 and see that as n increases the mode resolved cumulant becomes more peaked about λ = 0. Similar behavior is also seen for the higher cumulants.

V. FCS IN THE GAUDIN-YANG MODEL
We now consider the FCS in the Gaudin-Yang model, again splitting our analysis into two parts dealing with the repulsive and attractive regimes separately.In both cases the key component of our analysis, the overlap functions have been obtained previously by considering appropriate limits of the integrable quenches in the Hubbard model [96,97].

A. Attractive interactions
We start by considering the attractive regime and taking into account the quasi-particle content of the model reviewed in Sec.III B we find where the functions η n,β , ζ β , ζβ satisfy the integral equations Unlike in the LL model these equations do not admit an analytic solution that we are aware of and so must be treated numerically.There are however two limiting cases of interest.The first is the non-interacting limit c → 0 − in which case the dependence on the η n,β functions drop out and we find [10] for β ∈ R, which is the SCGF of the Possion distribution with rate 2d and is therefore twice the result we found in the Tonks-Girardeau limit of the LL model (55).Accordingly we immediately find that the cumulants are κ s m (0) = 2d where we recall that the total fermionic density is 2d.Additionally, the model has a second analytically tractable limit of interest, |c| → −∞.In this case not only do η n,β not contribute but now neither does ζ β and the FCS are determined solely by the twofermion bound states.This results from the fact that in the strongly attractive limit the fermion pairs in the initial state which make up c † 0 form bound states which cannot be broken apart.Thus the long time steady state consists only of these quasi-particles and no others.For this we find ζβ = e −2β λ 4 /4d 4 and the corresponding initial state FCS are for β ∈ R, i.e., the same result as in the free fermion limit, namely the SCGF of the Poisson distribution.In these two limits the SCGF can similarly be computed at the steady-states, namely we have and for β ∈ R, which equal the SCGF of the LL model at infinite repulsion upon the substitution d → 2d.
To investigate the cumulants at finite interaction we numerically integrate Eqs. ( 66)-( 68) which requires that we truncate the number of strings which are included i.e. η n,β (λ), σ n (λ) = 0, ∀n > N string and also that we impose a cutoff on the allowed rapidities, |λ| ≤ Λ. Doing so we can confirm that the cumulants of the steady state in fact do not depend on the interaction strength.Similarly to the case of the LL model, this implies that (65) yields 4d(e β/2 − 1) with d denoting the density of fermions, which can be confirmed numerically.
Finally, the full time evolution of the cumulants in the non-interacting limit can then be found using the previously result of (57) and similarly for c → −∞ whereas at finite c we again resort to numerical integration of the TBA system.In Fig. 6 we plot the second cumulant as a function of time for different interaction strength.Here we note note that in contrast to the LL model the initial state is never an eigenstate of the model and so regardless of interaction strength the cumulant relaxes at finite scaled time.Nevertheless, the approach to the GGE value still depends upon the interaction strength with the different dynamics being attributed to the changing role of the quasi-particles.Specifically for large negative interaction the two particle bound states are dominant and as shown in Fig. 6 their contribution is peaked nearer the slower quasi-particles leading to an almost linear initial decrease.As the interaction strength is lowered the unbound particles also start to contribute and their mode resolved cumulant is more spread in rapidity space leading to a faster initial decay.The string contribution is negligible in all cases.

B. Repulsive interactions
As a final example we look at the repulsive GY model.In this regime the spectrum of the model does not include the two particle bound states and we have that where the functions η n,β and ζ β satisfy a set of integral equations of the form (25).After using some standard TBA identities we find [10] log log where µ c is a Lagrange multiplier, introduced to fix the average density to be 2d and which behaves as lim c→0 µ c /c 2 = 2d 2 .As in the attractive case we cannot find an analytic solution to these equations however we can again consider the non-interacting limit c → 0 + .In this limit the dependence on the η n,β functions drop out and 1 [10].From this we then find that which agrees with the expression obtained from the (69).Thus the non-interacting limits of the FCS approached from both the repulsive and attractive sides agree.Unlike the attractive case the opposite limit of c → ∞ presents some problems.Indeed, exactly in this limit the initial state does not appear in the Hilbert space of the Hamiltonian since it precludes the possibility of two fermions being at the same position.By considering a lattice regularization of the model and taking the limit of infinite repulsion prior to the continuum limit it can be shown that the distribution functions become constants [96].Consequently, in the repulsive GY model, as one increases the interaction strength the distribution ρ(λ) becomes more and more spread in rapidity space.In the absence of an analytic solution to (74) and ( 75) must be treated numerically by first imposing some cutoff on both the string number and the rapidity space as was done in the attractive regime.However, the aforementioned spreading in rapidity space of the distributions results in a strong dependence on these truncation schemes leading to unreliable data except at small c where the results are qualitatively the same as the non-interacting limit.Despite this however, we expect that the relationship between the cumulants remains valid.The rapidity resolved second cumulant for the two particle bound states (solid lines, increasing interaction values corresponding to higher peaks), K2 0 (λ) as a function of λ for different interaction strengths.Also shown is the rapidity resolved second cumulant for the unbound states (dashed lines), K 2 0 (λ).In the c = 0 − case only unbound states are present (sole dashed curve distinguishable from the λ axis) while in the opposite limit |c| → ∞ only bound states are present.For any finite interaction the unbound states are strongly suppressed and not discernible on the same scale.

VI. CHARGE PROBABILITY DISTRIBUTIONS
In the previous sections we demonstrated that the scaled cumulants of charge fluctuations in the steadystate are universal for BEC quenches in both the LL and GY models.That is, they do not depend on the interaction strength of these quantum gases and they are determined by scaled cumulants in the initial states, which are essentially identical in both models.We now address the characterization of the full charge probability distribution in the steady state of the LL model and in the initial and steady-states of the GY model.The standard way of achieving this task is using the Gärtner-Ellis theorem, i.e., computing the Legendre-Fenchel transform of the SCGF and using Eq.(8).During the previous analysis, nevertheless, we also argued that the SCGF in the steady-states of the LL and GY models can be written by the analytic function independently of the interaction strengths, and where the density d denotes the density of bosons or the 'quasibosons' in the initial states of the LL and GY models, respectively.In fact, C s an equals Eq. ( 56), that is, the C s GGE (β)| c→∞ expression explicitly computed from the QA equations in the LL model.Similarly, C s an (β) equals (71) as well as (72) [that is, the C s GGE (β)| c→0 − and the C s GGE (β)| c→−∞ limits in the GY model] apart from a factor of 2, which further underpins the validity of C s an (β) and Eqs. ( 77) at arbitrary interactions.Finally we note that the numerical solutions for the SCGF on a finite interval, originating from the QA equations also confirm (77).
Given the explicit analytic expression for the SCGF of the steady-states, and focusing first on the LL model, the rate function I(z) can be easily obtained, namely i.e., twice the rate function of the Poisson distribution and thus P (ℓz, ∞) has non-trivial and non-Gaussian fluctuations and large deviations.The rate function above unambiguously characterizes the limiting coarse-grained PD P (ℓz, ∞) in the asymptotic sense (cf.( 9)).
Nevertheless, as we have anticipated, it is also instructive to approximate the steady-state PD by using that is by neglecting the o(ℓ) terms in the 2nd cumulant generating function.To use (79) the knowledge of the 2nd SCGF is required, but given the analytic nature of (77), we can invoke that In Fig. 7 we display the rescaled limiting probability distribution P (z, ∞) ≍ const × e −ℓI(z/ℓ) obtained from the rate function and the discrete probabilities resulting from (79) to demonstrate their agreement.It is important to stress, that although the PDs look Gaussian, this is not the case, since the cumulants and the asymptotics of the PDs are markedly different.One main feature that can be observed is that the variance of the PD in the steady-state is 1/ √ 2 times that of the initial state, or in other words, the distribution gets squeezed during the time evolution.
A remark can be made in connection with approximating the steady-state PD of a large but finite system via (79).Notably, the Fourier integral can be analytically computed resulting in lengthy expressions of trigonometric functions and exponential integrals.However, due to these expressions it is easy to show that C s an (iβ) does not define a PD in the strict sense, as for small ℓ negative probabilities can occur.Nevertheless increasing ℓ, as expected, the probabilities become positive numbers sum up to one and the discrete probabilities also reproduce the predicted results for the cumulants, which we have checked explicitly focusing on the first four.We now turn to the discussion of the PDs in the GY model and first consider the PD in the initial state.Similarly to the LL model, the scaled cumulants in the initial states are constant, that is, κ s m (0) = 2d, where d in this case is the density of fermion pairs.In the infinite subsystem size limit, the limiting continuous PD is therefore described by the rate function of the Poissonian, i.e., I 0 (z) = 2d − z + z log(z/(2d)).Nevertheless, the microscopic, discrete PD is anticipated to be different from a Poissonian, due the vanishing probabilities of odd particle numbers.As already stressed, in a strict mathematical sense our TBA/QA based methods allow only for the characterisation of the continuous, coarse-grained limiting probability distributions via the rate function.However, we can provide an approximate discrete probability distribution through (79) that takes into account vanishing probabilities and scales to the coarse-grained PDF as well, in the following way.Namely, when directly computing the 2nd SCGF at the free fermion limit using explicitly the integral in (69)  the Fourier integrals (79) vanish for every odd number.In addition, the integrals can again be explicitly evaluated in terms of analytic expressions and hence it can be checked that for large enough subsystems, the probabilities are positive, sum up to one and the scaled cumulants match the expected value 2d.That is, real physical behavior is correctly predicted by C s 0,GY and (79).We again stress that such a microscopic effect, i.e., the vanishing of probabilities for odd numbers, shall not affect the limiting and continuous PD P (ℓz, 0), however in large but finite subsystems, this feature is anticipated to be present and the corresponding discrete PD is assumed to be well approximated by C s 0,GY (iβ) via (79).Finally, we discuss the charge probabilities in steadystate of the GY model after the BEC quench.The coarsegrained PD P (ℓz, ∞) ≍ e −ℓI(z) can be simply obtained by specifying that I(z) equals (78) upon the d → 2d substitution.In the following, we again attempt to focus on some microscopic features of the discrete charge distribution in a finite subsystem.Specifically we are interested in whether vanishing probabilities for odd particle numbers can occur in the steady-state PD as it occurred for the PD of the initial state where it was expected based on physical considerations, but was indicated by the numerically computed 2nd SCGFs as well.
In particular, the feature of vanishing odd probabilities does take place in the c → −∞ limit, as in this case, only the fermion, spin-singlet bound states have non-vanishing densities in the steady-state.That is, the unit charge is two, hence probabilities for odd numbers must vanish [98].Indeed, when taking a look at the 2nd SCGFs originating from the QA equations and plotted in Fig. 8, it can be seen that in the c → 0 − , i.e., the free fermion limit, the analytically computed function C s (iβ) obtained from the integral in Eq. ( 71) performing the β → iβ substitution equals C s an .Nevertheless, for non-zero interaction strengths including the case of infinite attraction, we obtain collapsing non analytic functions.This function can be computed using the integral in (72) after substituting β with iβ.The function C s ∞,GY one obtains this way equals 4d(e iβ/2 − 1) on the [−π/2, π/2] interval, and outside this interval it is given by the properties Although the non-analytic behavior may again be attributed to the improper treatment of the branches of the logarithm, interestingly, together with the symmetry properties of the function they capture the aforementioned microscopic physical effect, and provide an approximate discrete PD via (79).That is, similarly to the fermionic BEC state, the non-analytic 2nd SCGF with the particular symmetry properties imply vanishing probabilities for odd fermion numbers, when the charge fluctuations in large but finite subsystem are considered via (79), which is the physically anticipated behavior at c = −∞.Moreover, just like in the case if the free fermion BEC initial state, the Fourier integral (79) involving C s ∞,GY can be computed analytically and when the subsystem size ℓ is large, the PD consists of positive probabilities summing up to one, and increasing ℓ the scaled cumulants converge to the prescribed value κ s m (∞) = 2d/2 m−1 as we have explicitly checked for the first four cumulants.
The situation is less understood at finite interaction strengths in the attractive regime.Although the numerically obtained 2nd SCGFs collapse to C s ∞,GY (iβ), in this case it is not clear whether or not the symmetry properties of the 2nd SCGF can be associated with a microscopic effect.Namely, it is not obvious why the probabilities for observing an odd number of particles in a subsystem should be zero, since at finite c, the root densities for the species associated with charge 1 are non-vanishing.Their contribution to the cumulants is, however, strongly suppressed for any non-zero interaction (see Fig. 6) and so it is possible that in the thermodynamic limit the odd probabilities are vanishing.Unfortunately this question cannot be definitively answered within the framework of the QA techniques we are using, therefore we leave it open for further investigations.In the repulsive regime of the model, however, we do not anticipate ambiguities and the presence of such microscopic effects due the lack of fermion bound states.
To conclude this section, we again wish to stress that by our methods we are able to give precise predictions for the continuous limiting PD P (ℓz) via the rate function I(z) for all the non-trivial cases, that is for the steadystate of the LL and the GY model.Going beyond the continuous description of the limiting PD, or, eventually the ℓ → ∞ limit, is nevertheless out of the range of the methods we have.Therefore, we find it worthwhile to note that despite this, in many instances we could still approximate discrete and microscopic PDs describing the charge fluctuations in large but finite subsystems via (79).In particular, by enhancing the periodicity of 4d(e iβ/2 − 1) as explained above, physical microscopic structures can be recovered in PDs of the fermionic BEC initial states as well as of the steady-state of the GY model at infinite attraction which manifest in vanishing probabilities for odd charges.This periodicity property of the 2nd SCGF can emerge naturally when the function is numerically computed from the QA equations and is a result of evaluating the logarithm at complex values and hence it may be purely accidental.In fact when the derivative of the 2nd SCGF w.r.t.β is computed at infinite attraction in the GY model, an analytic function with period 4π is obtained and hence the 2nd SCGF is assumed to inherit the same properties.Therefore, such indications of the numerically computed 2nd SCGFs are always to be reconciled with additional physical considerations, which nevertheless are lacking at the moment for the case of the steady-state fluctuations of the GY model at finite attractive interactions.

VII. CONCLUSIONS
In this paper we analyzed out-of-equilibrium charge fluctuations in two paradigmatic models of one dimensional quantum gases.In particular, we investigated the Lieb-Liniger model which describes interacting bosonic particles, and the Gaudin-Yang model in which the interacting particles are fermions and explored their entire parameter space by considering generic repulsive and attractive interactions.We focused on two particular initial states, the Bose-Einstein condensate state and its fermionic analog for the Lieb-Liniger and Gaudin-Yang models, respectively.These choices allowed for a rather complete characterization of non-trivial charge fluctuations over the course of the entire time evolution.This achievement is due the applicability of powerful methods relying on the integrability of the physical systems and the peculiar structure of the initial states.
In particular, using novel analytical techniques [59,60,68] and inspired by the quench action method [11], we determined all the scaled cumulants of charge fluctuations in the steady state as well their time evolution.These quantities characterize the fluctuations in a very large subsystem by keeping the leading order (linear in subsystem size) behavior.Whereas the time evolution of the scaled cumulants generically depends on the interaction strength of the models, surprisingly, their value in the steady-state was found to be independent of it.Moreover, it was also revealed that the scaled cumulants in the steady-state are uniquely characterised by the corresponding scaled cumulants in the initial state in a remarkably simple fashion.In particular, the scaled cumulants κ s m (0) in both the conventional and the fermionic BEC states are κ s m (0) = d and κ s m (0) = 2d, respectively, where d is the density of bosons and fermion pairs in the two states.In the steady-states the scaled cumulants κ s m (∞) are instead given by the universal relation κ s m (∞) = κ s m (0)/2 m−1 .This relation was established based on the explicit determination of the scaled cumulants in the initial and steady-states.A formal derivation of this relationship was also provided based on comparing the full counting statistics in the diagonal ensemble and in the generalized Gibbs ensemble and by showing that the former can capture the fluctuations in the initial state.
Thanks to the exponential decay of the scaled cumulants in the steady-states we could naturally invoke the scaled cumulant generating function whose determination based on only the knowledge of the cumulants, in general, might not be straightforward.Nevertheless, the analytic function obtained this way also agrees the generating functions obtained directly from the quench action method at specific interaction strength, when analytic computations are feasible or at generic interaction strengths, when the agreement can be checked numerically.Computing the Legendre-Fenchel transform of this function we obtained an explicit expression for the rate function which characterizes the limiting continuous probability distribution in an infinitely large subsystem.In accordance with the scaled cumulants in the steady-state, this function predicts non-trivial and non-Gaussian fluctuations and large deviations for the charge.Additionally, we also considered the probability distributions in large but finite subsystems, where their discrete nature is still visible.An interesting feature occurs in the fermionic BEC initial state as well as in the steadystate of the Gaudin-Yang model at infinite attraction.In these cases the probabilities of observing an odd number of charge in a subsystem identically vanish.Nevertheless, resolving such discrete and microscopic features in strictly finite subsystems is beyond the scope of our current methods and hence was accomplished by relying on additional physical input.
Our work admits several pathways for further investigations.A surprising finding is the universal relationship between the scaled cumulants of the initial and steadystates, which is a consequence of the integrability of both the models and the initial states [99].It would be important to identify the precise conditions for the onset of this phenomenon and to give a more rigorous explanation than what is presented in this work.A noteworthy remark is that, for initial states with vanishing scaled cumulants, such predictions clearly cannot hold.However, universality was manifest in our examples in another way as well, namely by the lack of dependence on the interaction strength of the models.It is an interesting question whether such interaction-independence can be observed in other, possibly non-integrable models or in other quench protocols, or it again requires the integrability of both the post-quench Hamiltonian and the initial state.By applying the semi-classical arguments of the quasi-particle picture it may be possible to go beyond initial states which have a pure pair structure at least for free models as has been carried out already for entanglement dynamics [100,101].
Finally we would like to highlight that similarly to the fluctuations of the charge, other conserved quantities can be investigated as well.This may require some care in the particular case of the Lieb-Liniger model, as some charges can have divergent expectation values and suitable linear combinations have to be considered [91].Nevertheless, the methods applied in this work can be applied for other charges, at least for the ones with extensive initial cumulants like the energy in a subsystem.Last but not least it would be worthwhile to accomplish the challenging task of developing analytic or semi-analytic methods that can capture the sub-leading behavior of fluctuations.This would be important to characterize certain microscopic effects, such as the fluctuations of conserved quantities in a subsystem for which the cumulants admit sub-extensive scaling e.g.those associated with the KPZ universality class [102].
rewriting of ⟨Ψ 0 |e β N |Ψ 0 ⟩ in the diagonal ensemble.Using the saddle point approximation and taking the logarithm of (A5) we can write down the SCGF as a difference of two effective free energy densities, i.e., where we adapted the notation of the logarithmic overlaps from Sec. II, took into account the 1/2 factors due to the pair-structure of the initial state, and ρsp and ρ sp are the two saddle-point root distributions of the nominator and the denominator of (A5).

j
, j = 1, . . ., M n .An eigenstate of the model is specified by the rapidities of the quasi-particles which are present, and we denote it by |λ⟩ = |λ

FIG. 6 .
FIG. 6.(a) The scaled second cumulant, κ s 2 (t) as a function of the scaled time ζ = t/ℓ for different values of the interactions strengths in the attractive regime.The increasing interaction strengths correspond to curves with increasing values at the ζ = 0.2 intersection.(b)The rapidity resolved second cumulant for the two particle bound states (solid lines, increasing interaction values corresponding to higher peaks), K2 0 (λ) as a function of λ for different interaction strengths.Also shown is the rapidity resolved second cumulant for the unbound states (dashed lines), K 2 0 (λ).In the c = 0 − case only unbound states are present (sole dashed curve distinguishable from the λ axis) while in the opposite limit |c| → ∞ only bound states are present.For any finite interaction the unbound states are strongly suppressed and not discernible on the same scale.

FIG. 8 .
FIG. 8.The real and imaginary (inset) parts of the 2nd scaled cumulant generating functions C s (iβ) in the steadystate of the Gaudin-Yang model with a fixed density d = 2.4 and various interaction strengths.The functions are obtained from the solution of the QA equations and performing numerical integration (c = −1.8,red dashed line) or using exact expression for the integral (20) (c = 0 and c → −∞ indicated by black ticks and blue continuous line, respectively).The function for c = 0 coincides with C s an (iβ), i.e., the analytic function defined by a Taylor series via the cumulants.
log⟨Ψ 0 |e β N |Ψ 0 ⟩ = ℓC s 0 (β) + o(ℓ) (78)7.The probability distributions (large red dots and black line) of charge fluctuations in the steady state of the LL model after a BEC quench in a large but finite subsystem.The red dots were obtained using Eq.(79) whereas the black continuous line corresponds to the rescaled limiting PD obtained from the rate function(78).The small blue dots show the Poisson PD of the initial state.The subsystem size ℓ equals 60 and the particle density is unity d = 1.