Scale invariant distribution functions in quantum systems with few degrees of freedom

Scale invariance usually occurs in extended systems where correlation functions decay algebraically in space and/or time. Here we introduce a new type of scale invariance, occurring in the distribution functions of physical observables. At equilibrium these functions decay over a typical scale set by the temperature, but they can become scale invariant in a sudden quantum quench. We exemplify this effect through the analysis of linear and non-linear quantum oscillators. We find that their distribution functions generically diverge logarithmically close to the stable points of the classical dynamics. Our study opens the possibility to address integrability and its breaking in distribution functions, with immediate applications to matter-wave interferometers.


Introduction
Scale invariance is a defining property of continuous phase transitions, which are invariant to the renormalization of the space and time coordinates. This scale invariance can be used to find the universal properties of the neighboring phases, through the renormalization group (RG) method. By construction, the RG approach does not directly apply to systems described by a small number of degrees of freedom, whose dimension cannot be rescaled continuously. A fundamental question is whether these few-body systems can show a universal behavior, and how to detect it 1 .
To address this question, we consider scaling transformations that act on physical observables, and look for the invariance of their distribution functions. A trivial example is offered by constant distribution functions, which do not change when the observables are rescaled. Interestingly, systems at thermal equilibrium generically belong to this universality class: under a scaling transformation of the variables, thermal fluctuations, and thus the temperature, effectively increase. In the asymptotic limit, the rescaled distributions tend to an infinitetemperature ensemble, where all possible values are equally probable, and the probability distribution is a constant. A natural direction to look for non-trivial scaling laws is offered by systems that do not thermalize, such as integrable models following a quantum quench. Most previous studies considered quenches in many-body systems and analyzed the scaling of the spatio-temporal coordinates 2 . Here, we study sudden quenches in few-body quantum oscillators and show that they give rise to probability distributions with a novel type of scale invariance.
Our definition of scale invariance is analogous to the common case, but involves distribution functions, rather than correlation functions. In connection to phase transitions, it is common to define the scale invariance through the two point correlation function F (x 1 −x 2 ) ≡ φ(x 1 )φ(x 2 ) , where φ is some physical property of an extended system and x 1/2 are two positions in space. The system is then said to be scale invariant for large x if F satisfies the scaling ansatz F (x) ≈ λ α F (λx) (1) where α is a critical exponent 3 . Eq. 1 is satisfied, for instance, if the correlation function decays at large distance as a power-law, F (x) ∼ x −α . In this paper, we instead consider the distribution function P of a physical observable x, and show that under appropriate condition they can be scale invariant as in Eq. 1 in the vicinity of a stable fixed point. Specifically, the above mentioned thermal case corresponds to a situation where P (x) = Ae −αx 2 , which tends to a constant for small x. As we will show, the distribution functions of quenched oscillators are generically characterized by a logarithmic divergence 4 P (x) ≈ κlog(x), which is scale invariant because P (λx) = P (x) + κlog(λ) ≈ P (x), for x → 0.
1 One possible strategy that was discussed in the literature is to use the time axis as a scaling variable [1]. The corresponding RG approaches focus on the dynamics of individual orbitals, and helps understanding the transition between regular motion and chaos [2,3]. Here, we instead consider ensembles of initial conditions, and study the statistical properties of their long-time dynamics.
2 See for example Refs. [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] for quantum quenches of integrable many-body systems. 3 For a more rigorous definition of scale invariance, one may wish to consider the finiteness of the ratio between the left and right hand sides of Eq.(1), in the limit of x → ∞. Note that in a scale invariant system, higher-order correlation functions are scale invariant as well. Their scale invariance is defined by extending Eq. 1 to multi-variable functions. 4 Note that the logarithmic divergence does not pose any problem in terms of normalizability of the distribution function because At an intuitive level, the scale invariance can be simply understood by considering the linearized equations of motion close to a stable point. Being linear, these equations are invariant under the scaling transformation x → λx, where x is the distance from the stable point and λ is a constant 5 . To obtain a scale invariant ensemble, it is then sufficient to complement these equations with a scale invariant initial state, such as a particle with a fixed momentum, whose position in real space is completely uncertain. The key result of this work is that this simple phenomenon survives non linearities and is intimately related to the model's integrability.

The harmonic oscillator
We open our discussion with the analysis of an isolated harmonic oscillator H 0 = (x 2 + p 2 )/2, where x and p are canonical conjugates. Here the simplest example of a scale invariant state is offered by |p = 0 , which satisfies x|p = 0 = const. In a semiclassical description (which is exact for an harmonic oscillator), this state corresponds to the Wigner distribution P (x, p) = P 0 δ(p), where P 0 is a normalization constant 6 . Under the effects of H 0 , this ensemble rotates in phase space: each point follows a circular trajectory around the stable point x = p = 0, with constant angular velocity. Thus, after time averaging, one obtains a distribution function that is inversely proportional to the circumference of the circle, or Here the factor 2 in the numerator accounts for the orbits starting from x and −x, which contribute to the same circumference. We can now use Eq. 2 to compute the (time-averaged) marginal probability of x where x 0 is an arbitrary cutoff, and o(1) is a constant term that does not diverge as |x| → 0. Eq. 3 shows that the distribution function of x diverges logarithmically and is therefore scale invariant (see the Introduction) 7 .
In this work, we show that the logarithmic divergence found in Eq. 3 is universal, it is not affected by non-linearities. This result is non-trivial because, for any finite x, there exists a time after which the non-linearities have a significant effect. The logarithmic divergence is nevertheless preserved, as long as the fixed point x = p = 0 is stable and the dynamics in its surroundings is characterized by invariant tori. For a scale invariant initial state, the 5 In this sense, the present scale invariant states can be associated with a Gaussian fixed point. At equilibrium, these fixed points offer the simplest example of scale invariant critical points. An interesting question for further studies is whether distribution functions can show non-Gaussian fixed points that are scale invariant as a consequence of non-linear terms. 6 See Ref. [20] for an introduction to phase-space methods for quantum mechanics. 7 This analysis can be extended to a generic harmonic oscillator with mass m, and natural frequency ω0: by working with normalized variables, it is straightforward to see that P (x) does not depend on m and ω0 (see Appendix A.1) time-averaged P (x, p) is inversely proportional to the circumference of the appropriate torus, which is in turn proportional to the distance from the stable point. The integration over one variable will then generically lead to a logarithmic divergence 8 .
For large S, the Hamiltonian in Eq. 4 is well approximated by a semiclassical description [49,50], where the spin operators are substituted by two continuous variables, n and φ, defined by S z /S = n, S ± /S = √ 1 − n 2 exp(±iφ)/2. The canonical variables n and φ respectively correspond to the number and phase differences of the two-site Bose-Hubbard model. Under this transformation, the Hamiltonian in Eq. 4 is mapped to The classical dynamics associated with this Hamiltonian has two fixed points on the line n = 0, respectively, at φ = 0 and φ = π. Their dynamical stability depends on the ratio between J and µ: for J < |µ|, the system is stable only around φ = 0, while for J > |µ| the system becomes stable around φ = π as well. This transition is associated with an equilibrium meanfield phase transition (for µ < 0), or with the disappearance of macroscopic self-trapping (for µ > 0) [51,52]. As we will see, this point determines a discontinuous change in the scaling properties of the distribution functions. To achieve a scale invariant distribution function we consider the initial states |S z = 0 . This state corresponds to the ground state of the Hamiltonian in Eq. 4 with J = 0. Thus, the present dynamics is equivalent to the experimentally-relevant situation of a quantum quench in which J is suddenly changed from 0 to a finite value [53][54][55]. In the semiclassical description of Eq. 5, this initial state is mapped to an ensemble with n = 0 and a uniformly distributed φ ∈ (−π, π), or equivalently P (n, φ) = δ(n)/2π. Fig. 1 shows the evolution of this ensemble, obtained by the numerical solution of the Hamilton-Jacobi equations derived from Eq. 5, for J = 0.2µ. The marginal distribution P (φ) is shown in the lower panel and evolves from P (φ) = P 0 = 1/2π to the universal shape P (φ) = −(1/π 2 ) log(φ), as predicted by Eq. 3. This result confirms that the nonlinear terms present in the Hamiltonian in Eq. 5 do not affect the logarithmic divergence close to the stable point.  We now compare the above-mentioned semiclassical calculations with the exact diagonalization of the quantum Hamiltonian in Eq. 4 with S = 1000. In the quantum model, the logarithmic divergence can be observed in the distribution of the operator m y ≡ S y /S = √ 1 − n 2 sin(φ), which can be approximated by m y ≈ φ, in the vicinity of the stable point n = φ = 0. The time-averaged distribution function of m y is defined quantum mechanically by where |ψ(t) = e −iĤt |S z = 0 , and H is the Hamiltonian in Eq. 4. As shown in Fig. 2(a), the resulting distribution function diverges logarithmically around m y = 0. In actual systems, this divergence is rounded at 1/S, which plays the role of the infra-red cutoff of our theory (see Appendix A.3 for details). The inset of Fig. 2(a) shows that the prefactor of the logarithm suddenly jumps at J/µ = 1: At this point, the number of stable points across the m z = 0 line jumps from 1 to 2, leading to a doubling of the prefactor of the asymptotic distribution function 9 . A similar argument can be used to determine the universal scaling of other physical observables (see Appendix A.4).
will generically tend to the time-averaged expression 9 A closer inspection of Fig. 2 shows that for J < µ, P (my) shows a cusp at finite my. This cusp is associated with two additional stable fixed points at Sz = 0, which correspond to the two ferromagnetic equilibrium states. The presence of these stable points is at the origin of the macroscopic quantum self-trapping effect. As approaching J = µ, the cusp shifts to smaller my and, for J > µ, it joins the divergence at my = 0, doubling the prefactor of the logarithm. A similar behavior can be obtained by the numerical solution of the semiclassical equations of motion associated with (5)

Breaking of integrability
The logarithmic divergence of the distribution function is due to the presence of closed orbits in the vicinity of a stable point. These orbits are protected by the integrability of Eq. 4, which involves the same number of degrees of freedom (S x , S y , S z ) as of conserved quantities (S 2 , S z , and H). To study the effects of integrability breaking terms, we now turn to two models where the number of degrees of freedom is larger than the number of conserved quantities: the Dicke model and the kicked rotor. The Dicke model [56] is a canonical model of quantum optics. It describes the interaction between a cavity mode (a) and a large ensemble of spins (or, equivalently, a single large spin S). In the thermodynamic limit of S → ∞, the Dicke model undergoes a phase transition from a normal to a super-radiant phase [57,58], at a critical value of the cavity-spin coupling, λ = λ c . This transition was throughly described both at equilibrium and out-of-equilibrium, and recently observed in cavity-QED 10 . The Dicke model has the same number of conserved quantities as the model defined in Eq. 4, but one additional degree of freedom. As a consequence, the Dicke model can give rise to a chaotic motion, whose onset occurs in the close vicinity of the phase transition [60].
We numerically simulate the Dicke model using the semiclassical equations of motion derived in Ref. [60], which are valid for S 1 (see Appendix A.5). Our initial state corresponds to a pure state where |S z = −S/2, and the photon is largely squeezed, to mimic a scale invariant state. At long times, the probability distribution of the squeezed quadrature diverges logarithmically (See Fig. 2(b)). The prefactor of the logarithm is constant for all λ < λ c , and equals that in Eq. 3. At the critical coupling λ c , the system becomes chaotic and tends to thermalize: correspondingly, the logarithmic divergence suddenly disappears (inset).
We next move to a canonical model used to describe the transition between regular and chaotic dynamics, the kicked rotor (see Appendix A.6). This model has a fixed point at x = p = 0, whose vicinity becomes chaotic at a critical value of the kick strength K c = 4. In Fig. 2(c), we show the long-time distribution obtained from an initial ensem-ble with a uniformly distributed momentum p ∼ U (0, 2π) and a constant position x = 0. We observe that the distribution function of p develops a logarithmic divergence close to p = 0. Interestingly, we find that the prefactor is not constant, but follows the empirical law κ = (−2P 0 /π) 1 − K/K c . This curve is non-analytic at K c , at the onset of chaos, where the logarithmic divergence is washed out. These findings strengthens the relation between the integrability and the logarithmic divergence of the probability distribution 11 .

Beyond Hamiltonian systems: dissipation
We now turn to study the effects of dissipation, relevant to the experimental realization with matter-wave interferometers [53,62]. We model this effect by where η is the dissipation constant 12 . In the limit of η → 0 these equations of motion are equivalent to the Hamilton-Jacobi equations associated to Eq. 4. The dissipative term is invariant under the scaling transformation φ → λφ, n → λn: As demonstrated by the numerical calculations of Fig. 3(a) (for J/µ = 0.2, η/µ = 0.1) the distribution of φ is still logarithmically divergent, although the prefactor becomes time dependent.
To understand this behavior, we go back to the phase space picture, where each point follows a spiral motion (inset of Fig. 3(b)). Close to the stable point, the motion is described by a damped harmonic oscillator, whose solution gives φ(t) = φ 0 e −ηt cos(ωt). As a consequence, the phase-space density grows as e ηt and the time-averaged distribution is given by As shown in Fig. 3(b), this expression is in quantitative agreement with the numerical solution of the full non-linear model.

Conclusion: full scaling theory
A logarithmic divergence is invariant under the scaling transformation x → λx, and this property can be used to address the effect of generic perturbations. Under the scaling transformation, all non-linear terms appearing in the equations of motion tend to zero ("irrelevant"). These terms do not affect the logarithmic divergence of the distribution functions 11 Note that the present semi-classical analysis does not take into account the dynamical localization due to quantum coherence [61]. The consequences of this effect on the logarithmic divergence requires further investigation. 12 Note that our dissipative term differs from the expression used in Ref. [62], where a force proportional to −η(dφ/dt) was considered. Our linear term has a phenomenologically similar effect, but simplifies the calculation of the correspondent fluctuating forces. ( Figs. 1 and 2(a-b)). Linear perturbations are invariant under the scaling transformation ("marginal"): These terms modify the prefactor of the logarithmic divergence, and eventually lead to its disappearance (Figs. 2(c) and 3). Finally, if a term does not depend on x, it effectively grows under the scaling transformation, and destroys the logarithmic divergence ("relevant"). A natural example is offered by the random forces associated with a coupling to a thermal bath. These forces generically drive the system towards an equilibrium distribution function, of the form P eq (φ) = P 0 exp(−E(φ)/T ), where the E(φ) is the energy. This expression is analytical around φ = 0, indicating that P (φ) does not diverge. To study this effect numerically, we consider Eq. 8 with an additional stochastic force f (t). According to the fluctuation-dissipation theorem, this force satisfies f (t) = 0, and f (t)f (t ) = 4ηT δ(t − t ), where T is the temperature of the bath. As shown in Fig. 4 (at temperature T = 0.1), the system flows towards a thermal distributions, and the logarithmic divergence is destroyed.
The logarithmic divergence of the distribution function is therefore a clear indicator of the absence of thermalization in quenched oscillators. Our scaling theory can be used to analyze the effect of generic perturbations (see Appendix A.7). This approach shows a possible way to generalize our findings to many-body systems: the Lipkin-Meshkov-Glick and Dicke models are exact mean-field solutions of interacting systems with infinite-range interactions. By considering the perturbations induced by a finite-range, it will be possible to study the crossover to extended many-body systems. Finally, by including the effects of disorder, one can attempt to describe the non-Gaussian distribution functions that were recently found in quantum quenches of many-body-localized systems [63].  Figure 4: Same as Fig. 3, in the presence of a thermal noise at temperature T = 0.1. The system thermalizes and the logarithmic divergence of the distribution function is destroyed.

A.1 Harmonic oscillator with non-unit mass and frequency
In the main text we considered an harmonic oscillator with natural frequency ω 0 = 1, and mass m = 1, whose phase space orbits are circles. Let us now consider an Harmonic oscillator of the form H = (p 2 /m + mω 2 0 x 2 )/2. Its equations of motion are given by For convenience, we now introduce the rescaled variablesx = √ mω 0 x andp = p/ √ mω 0 , whose equations of motion are If we rescale the time tot = tω 0 , we are back to the case discussed in the main text. Thus, using Eq. 3, we find that the time-averaged distribution function ofx (for smallx) is HereP 0 is determined by the initial conditions, given by P (x,p) = P ( we finally obtain Importantly, Eq. 13 does not depend on m or ω 0 , giving a first hint about the universality of this result.

A.2 Two-site Bose-Hubbard model
The two-site Hubbard model is described by the Hamiltonian where µ is the chemical potential, and J the tunneling element. Because the model commutes with the total number of particles, we restrict ourself to the subspace with a fixed N = ψ † 1 ψ 1 + ψ † 2 ψ 2 . The Hamiltonian in Eq. 14 is conveniently described in terms of N spin-1/2 variables, σ i , whose z component describes the site occupied by the i th particle [64,65]. This mapping is formally achieved through the Schwinger boson representation of spin operators S α = 1/2 i,j=1,2ψ † i σ i,j α ψ j , where α = x, y, z and σ α are Pauli matrices. By introducing the total spin operator S = N i=1 σ i , one can exactly map Eq. 14 to the Lipkin-Meshkov-Glick model, Eq. 4 of the main text, with S = N/2.

A.3 Finite size scaling
Our derivation of a scale invariant distribution functions relies on a semiclassical description of a quantum model. Specifically, the analysis of the Lipkin-Meshkov-Glick model of Eq. 4 referred to the limit S → ∞, where the quantum spin becomes a semiclassical rotor. In this appendix we consider the effects of a finite S. For this purpose, we study the steady-state distribution functions of the model for different values of S. As shown in Fig. 4, the logarithmic divergence is already evident for S = 250. Because the minimal value of m y = S y /S is 1/S, the distribution function is terminated at this value. As S increases, the cutoff becomes smaller, and the logarithmic divergence more pronounced. Thus, a finite S has a similar role to the infra-red (IR) cutoff of a scale invariant theory, which is usually determined by the finite size of the system.

A.4 Other observables
In the main text, we focused on the probability function of the variables φ and n and we showed that they diverge logarithmically around the stable fixed point n = φ = 0. The distribution function of other physical observables can be directly computed from P (n, φ). Let us consider for instance the operator m x = S x /S = √ 1 − n 2 cos(φ). Close to the stable point n = φ = 0, this quantity can be approximated as m x ≈ 1 − (n 2 + φ 2 )/2. Following the same arguments as in Sec. A.1 we obtain This result is numerically confirmed in Fig. 6.

A.5 Dicke model
The Hamiltonian of the Dicke model [56] is Here S is a spin operator (as in the main text) and a is a canonical bosonic operator satisfying [a, a † ] = 1. In the limit of S → ∞, this model undergoes a phase transition [57,58] at For large S, the Dicke model in Eq. 16 is well approximated by the semiclassical Hamiltonian (Eq. 65 of Ref. [60]) were x, p x , y, and p y are two pairs of canonical coordinates (associated with the two quadratures of of the cavity boson, and of the spin, respectively).
The correspondent equations of motion are (Eqs. 68-69 of Ref. [60]) where This model shows a transition between regular and chaotic motion at λ ≈ λ c . In our numerical calculations, we considered S = 10 6 . The initial state of the spin was chosen to represent the quantum |S z = −S (which corresponds to the ground state of the model for λ = 0). In the semiclassical pictures, this state is represented by a Wigner distribution in which y and p y are extracted from Gaussian ensembles with zero average and variances 1/(2ω 0 ) and ω 0 /2, respectively. The state of the boson was chosen to represent a vacuum squeezed state with x = p = 0, x 2 = 10 6 /4 and p 2 = 10 −6 , satisfying the minimal uncertainty relation between canonical variables. The model's parameters are chosen such that the frequency of the x and y oscillators are incommensurate: ω 0 = 1/ √ 2 and ω = √ 3. We observed empirically that the case ω = ω 0 gives rise to a distinct behavior, which requires further investigation. The equations of motion were solved using the Euler method with time-step discretization of dt = 0.01, and the distribution functions were averaged over times up to t = 100.

A.6 Kicked rotor and Chirikov standard map
The Hamiltonian of the kicked rotor is (See Ref. [66] and references therein) where δ is the Kronecker delta function. Note that in previous literature, the model is often defined with an opposite sign of K, or equivalently after the transformation x → x + π.
The stroboscopic dynamics of the model (i.e. the evolution of the system after a discrete number of time periods) is governed by the Chirikov standard map p n+1 = p n − K sin(x n ) (21) x n+1 = x n + p n+1 .
Due to the periodicity of the model, it is then common to define the dynamics on a torus, where x and p are restricted to the interval (0, 2π).
The dynamics of the model in Eq. 22 is characterized by three distinct regimes: For K < K c ≈ 0.9716 the model is localized between invariant tori (i.e. p does not grow with time); For K c < K < 4 the model has a mixed phase space, where the dynamics is diffusive for most some conditions, and localized in vicinity of the stabel point x = p = 0; For K > 4 the region around the stable point becomes chaotic.

A.7 Extended Lipkin-Meshkov-Glick model
In this section we explain how to apply the scaling analysis to predict the effect of non-linear terms on the logarithmic divergence. For this task, we consider the a generalization of Eq. 4, which includes two additional terms Within the semiclassical approach, the first term, S z = n enters into the equations of motion of dφ/dt as a constant term. This term grows under scaling and destroys the logarithmic divergence. In contrast, S 2 x = (1 − n 2 ) cos 2 (φ) is a non-linear perturbation and does not affect the logarithmic divergence. These predictions are verified numerically in Fig. 7, where we consider the initial state |S z = 0 with S=1000, evolve it in time with the Hamiltonian of Eq. 23, and compute the (time averaged) distribution probabilities of the operator m y = S y /S. As predicted by the scaling analysis, the coupling α destroys the logarithmic divergence, while β leaves it unchanged.