Dynamics of longitudinal magnetization in transverse-field quantum Ising model: from symmetry-breaking gap to Kibble-Zurek mechanism
Michał Białończyk Bogdan Damski

















eq

⟨
⟩

i
c.c.
i_1 i_2…i_N
j_1 j_2…j_N

#1(#1)
#1[#1]
#1|#1|
#1{#1}


    p_GS0.65fitGSEX for cross#1#2⟨ #1  #2⟩@nameJagiellonian University, Marian Smoluchowski Institute of Physics, Łojasiewicza 11, 30-348 Kraków, Poland
            We show that   the symmetry-breaking gap of the quantum Ising model in the transverse field can be extracted
            from  free evolution  of the longitudinal magnetization taking
            place after a gradual  quench of the magnetic field. 
            We perform  for this purpose  numerical simulations  of   the  Ising chains
            with either periodic or open boundaries.
            We  also study  the condition for adiabaticity of evolution of the longitudinal magnetization
             finding excellent agreement between 
             our simulations and  the prediction based on the Kibble-Zurek theory of non-equilibrium
            phase transitions. Our results should be relevant for ongoing cold atom and
            ion experiments targeting either  equilibrium or dynamical aspects of  
            quantum phase transitions. They could be also useful for benchmarking D-Wave
            machines.
            Dynamics of  longitudinal magnetization in   transverse-field quantum
            Ising model: from  symmetry-breaking gap to   Kibble-Zurek mechanism
                Michał Białończyk and Bogdan Damski
                December 30, 2023
            ==========================================================================================================================================
            
            
            
            § INTRODUCTION
            
            
            Quantum phase transitions take place when tiny changes of the external
            parameter, such as the magnetic field in spin models or the lattice height 
            in cold atom setups, can induce radical changes in  ground state properties of the
            system <cit.>. This can happen when the system is near the 
            critical point separating its   phases. 
            
            
            One of the typical features associated with quantum phase transitions is the
            symmetry-breaking phenomenon, where one of the phases of the
            thermodynamically-large system has degenerate 
            ground states, which do not respect the symmetry of
            the Hamiltonian. 
            Such a degeneracy is typically lifted  in
            finite systems, where a small energy gap between  lowest-energy eigenstates of
            the Hamiltonian is present. Rapid disappearance of this gap with increase of the system size signals the
            onset of the symmetry-breaking phenomenon that is fundamentally important for
            understanding of both equilibrium and non-equilibrium  phase
            transitions.
            It is thus  very  interesting  to address the question how one can
            access the symmetry-breaking gap in real systems to observe emergence of such a compelling  phenomenon.
            
            
            We make a step forward in this direction by providing an explicit two-stage scheme for
            studies of the symmetry-breaking gap in 
            the paradigmatic system undergoing a quantum  phase transition: the quantum Ising model in the
            transverse magnetic field <cit.>.
            In the first stage, the system is gradually driven to the desired value of
            the transverse magnetic field 
            (see <cit.> for   studies of  different aspects of dynamics of the
            quantum Ising model under similar driving).  In the second stage, the system undergoes free
             evolution  and oscillations of its  longitudinal magnetization are recorded.
            They  encode the symmetry-breaking gap.
            
            
            While our studies target a specific model, we expect that they 
            can be generalized to other, not necessarily 
            exactly solvable systems, which can be experimentally approached in cold atom and
            ion simulators of various  condensed matter models
            <cit.>.
            In fact, while pursuing  this work, we have come across a recent  paper
             discussing experimental studies of the symmetry-breaking gap in a
            cold atom cloud <cit.>. We will first present 
             our results and then compare the  two approaches. 
            
            The outline of this paper is the following. We explain the idea behind our
            work in Sec. <ref>. 
            We  illustrate in
            Sec. <ref> the
            concepts related to the symmetry-breaking gap and the 
            equilibrium longitudinal magnetization in the simplest version of the Ising model composed of just two
            spins.  We then discuss in Sec. <ref> dynamical extraction of those quantities from after-quench 
            free evolution of the longitudinal magnetization in the periodic
            Ising chain composed of several spins. Next, we extend these studies in Sec.
            <ref> to systems composed of up to several hundreds of spins by
            considering Ising chains with open boundaries, where  computations of the longitudinal
            magnetization can be more efficiently done. The overall
            discussion of our results is presented in Sec. <ref>. 
            Finally, technicalities related to the  studies of the  Ising
            chains with periodic (open) boundaries  are presented in Appendix <ref> (Appendix <ref>).
            
            
            
            § IDEA
            
            
            To explain the idea behind our proposal, we introduce the 
            quantum Ising model in the transverse magnetic field g, whose Hamiltonian 
            in the periodic chain can be written as <cit.>    H(g) = -∑_i=1^Nσ^x_iσ^x_i+1 + g σ^z_i,
               σ^x_N+1=σ^x_1,
            
            
            where σ_j^x,y,z is the Pauli matrix acting on the j-th spin,  N
            is the number of spins, and   g≥0 is assumed in this work. 
            Such a model has two phases: the ferromagnetic phase for  0≤ g<1 and 
            the paramagnetic phase for g>1.
            The breaking of the ℤ_2 symmetry of this model, associated with
            the σ^x_j→-σ^x_j symmetry of the Hamiltonian, 
            can be easily explained  at g=0, when the two ground states have all
            spins aligned in the ± x direction 
            
            |→→→⋯,  |←←←⋯.
            
            Then an arbitrarily  small perturbation  along the x direction, say
            -h∑_iσ^x_i with h→0^+, makes 
            |→→→⋯ the ground state of the system breaking the spin-flip symmetry
            supported by the Hamiltonian. Importantly, degeneracy of the ground state, in thermodynamically-large 
            systems, persists in
            the whole ferromagnetic phase.
            
            
            The idea that we have is that one can initially set g=0 and  prepare the system in one of the states 
            given in (<ref>), say  
            
            |ψ= |→→→⋯.
            
            We propose then to evolve the system according to the following  protocol 
            
            g(t)=
            {[               g_f+t/τ_Q for  -g_fτ_Q  ≤ t  ≤  0;                     g_f                for  t>0 ].
            
            with τ_Q controlling the rate of  driving   the system to the
            magnetic field g_f, where  free evolution begins. 
            
            
            The natural observable for  studies of the symmetry-breaking gap is 
            the longitudinal magnetization
            
            M_x=σ^x_i,
            
            where the choice of the site index i is arbitrary in a periodic chain.
            Such an observable is sensitive to the way in which the symmetry is broken.
             Indeed, it is well-known that 
            
            -M_x^≤ M_x≤ M_x^,
            
            where in the thermodynamically-large system  <cit.> M^_x={[ (1-g^2)^1/8 for  0≤ g≤1;           0    for  g>1 ]..
            M_x^ will be called the equilibrium longitudinal magnetization. It was
            obtained in the N→∞ limit from  the asymptotic form of the two-point correlation functions.
            We will study it in finite systems, where a different approach has to be employed (see below).
            Typical results that we obtain are illustrated in Fig. <ref>, where we  see that 
            M_x^ is analytic when N<∞. In such systems,
            the  drop of the longitudinal magnetization is not as
            steep as in (<ref>) on the ferromagnetic side.
            On the paramagnetic side,
            non-zero value of M_x^ appears–see   (<ref>) and discussion around
            it. 
            
            
            
            
            
            
            The question now is how the symmetry-breaking gap can be extracted from
            dynamics of the longitudinal magnetization. To explain that, we first note 
             that  the expectation value of the σ^x_j operator vanishes 
             in 
             all eigenstates of  Hamiltonian 
            (<ref>) for g>0. This  can be proved with Wick's theorem <cit.> in the fermionic representation 
            of the Ising model, which we briefly introduce  in Appendix <ref>.
            Non-zero value of the longitudinal magnetization  
            can be  obtained by evaluation of  (<ref>) in a proper  superposition of  two  eigenstates 
            of the Hamiltonian.
            In fact,  deceptively simple-looking initial state (<ref>) is a macroscopic
            superposition of two  eigenstates belonging to  different sectors of the spectrum of
            Hamiltonian (<ref>).
            
            To understand this statement, we note that   Hamiltonian
            (<ref>) commutes with the parity operator,
            
            P=∏_i=1^N σ^z_i,     [H,P]=0.
            
            This leads to splitting of the Hilbert space into two subspaces, where the eigenstates have either
            +1 or -1 parity <cit.>. The  symmetry-breaking gap δ is the  difference
            between the energies of the lowest-energy eigenstates    in the negative and positive parity subspaces
            in the ferromagnetic phase, where the symmetry-breaking phenomenon takes place.
            Such a quantity, however, will be of key importance in 
            our studies  for any value of
            the magnetic field g. We will thus call it the symmetry-breaking gap
            even when it will be computed at the critical point or in the paramagnetic
            phase, which  will simplify our discussion.
            
            
            These features of the spectrum are illustrated in Fig. <ref>.
            A closed-form expression for δ  was given in  <cit.>
                δ  =  
                g^N ∫_0^1 dt 4 N/πt^N-3/2√((1-t)(1-g^2 t))/1-(gt)^2N  for  0≤ g<1,
                
                   δ = 2g-2 +   g^-N∫_0^1 dt 4 N/πt^N-3/2√((1-t)(g^2-t))/1- t^2N/g^2N  for g>1,
                
                   δ = 2tanπ/4N  at  g=1.
            
            These expressions can be used for showing that
            δ vanishes exponentially (algebraically) with  the system size deeply in the ferromagnetic
            phase (near the critical point) <cit.>. Deeply in the paramagnetic phase, δ is
            well-approximated by an expression that is obtained from (<ref>) after 
            neglecting the term containing  the integral. 
            
            Besides the symmetry-breaking gap, there are  also  dynamical gaps
            Δ_+ and  Δ_-. They are  defined as energy gaps 
            in  subspaces of positive and negative parity 
            between the ground state and the first
            excited state that can be
            populated during evolution starting from (<ref>)–see Fig.
            <ref>.
            Their importance comes from the fact that commutation relation (<ref>) 
            prohibits dynamical transitions between the two parity  subspaces. As a
            result, system's excitation due to  driving (<ref>) 
            depends on Δ_±, the quench time τ_Q, 
            and the initial state  (<ref>) for time evolution. It  
            has nothing to do with the symmetry-breaking gap δ as long as there is
            no symmetry-breaking perturbation in the system, which is the case in our
            studies.
            
            
            
            If we now note that σ^z|→=|←, which implies 
            σ^z|←=|→, we see that the ground states 
            in the positive and negative parity subspaces at g=0 can be written as 
            
            |_±(g=0)=|→→→⋯± |←←←⋯/√(2),
            
            which in turn allows us to cast initial state (<ref>)  to the following form
            
            |ψ=|_+(g=0)+|_-(g=0)/√(2).
            
            
            Next, we use (<ref>) and employ time-dependent Schrödinger equation to arrive at    
               |ψ(t)=|ψ_+(t) + |ψ_-(t),
               P|ψ_±(t)=±|ψ_±(t).
            
            Finally, knowing that σ^x_i has vanishing matrix elements in  the positive and
            negative parity subspaces, which can be easily proven with the help of Wick's theorem, 
            we get 
            
            M_x(t)=ψ(t)|σ^x_i|ψ(t)=
            ψ_-(t)|σ^x_i|ψ_+(t) + c.c.
            
            This  expression can be very  complicated for fast quenches, i.e.,
            for small τ_Q in (<ref>).
            However, if we  assume that the driving  is slow enough to be nearly 
            adiabatic, then the following expression will properly approximate the exact
            result
            
            M_x(t)=cos∫_-g_fτ_Q^t dt δ[g(t)]_-[g(t)]|σ^x_i|_+[g(t)].
            
            Such an expression is  still  non-trivial because we are unaware of a closed-form expression for
            the above matrix element in an arbitrarily-sized system.
            If we now consider t>0, i.e., the free evolution stage in our problem, we
            will get from (<ref>) that 
            
            M_x(t)=cosδ(g_f) t+const M^_x(g_f),
            
            which can be used for extracting the symmetry-breaking gap and the equilibrium
            longitudinal magnetization out of either  numerical or experimental data.
            
            The question now is what is the condition for adiabaticity in our system so
            that approximation (<ref>) can be used.
            Clearly, the system is most prone to being excited near the critical point,
            where the dynamical gap Δ_± 
            is smallest,   and the quantum version <cit.> of the
            Kibble-Zurek (KZ) theory of non-equilibrium phase transitions
            <cit.> can be used. 
            A simple criteria then exists and is based on comparision between the system
            size and the non-equilibrium  length-scale ξ̂ characterizing excitations resulting from the
            quench <cit.>. The latter, according to  the KZ theory, scales with the quench time as
            ξ̂∼τ_Q^ν/(1+zν),
            
            where z and ν are the critical exponents.
            These exponents  describe disappearance
            of the dynamical gap and the correlation length, which are 
            assumed to be proportional to |g-g_c|^zν and  |g-g_c|^-ν,
            respectively, near the critical point g_c.
            Comparing the two length scales, one gets the following  condition for a crossover between
            adiabatic and non-adiabatic evolutions approaching or crossing the critical
            point
            τ_Q∼ N^(1+zν)/ν∼ N^2,
            
            where we substituted   z=ν=1 for the quantum Ising model in the transverse
            field. In other words, for quench times of the order of N^2 or larger we
            expect the evolution to be nearly adiabatic.
            
            
            
            
            
            
            
            § TWO SPINS
            
            We illustrate here the quantities  introduced in the previous section by showing 
            how the symmetry-breaking gap and the longitudinal magnetization can be calculated  in the
            simplest version of the Ising chain, the one with just two spins. 
            
            To start, we need eigenenergies and normalized  eigenstates of (<ref>),
            which for N=2 read 
            
                E__+=-2√(1+g^2), |_+=A_-|↑↑+B_-|↓↓, 
                    E__-=-2,  |_-=|↑↓+|↓↑/√(2), 
                    E__-=2,   |_-=|↑↓-|↓↑/√(2),
                    E__+=2√(1+g^2),  |_+=A_+|↑↑+B_+|↓↓,
            
            where  σ^z|↑=|↑,
            σ^z|↓=-|↓,
            
            A_±=1/√(2)√(1+g^2± g√(1+g^2)),  
            B_±=-g±√(1+g^2)/√(2)√(1+g^2± g√(1+g^2)),
            
            and   _± refers to the  excited states in the
            corresponding subspaces. 
            
            Using these results, one immediately gets
            δ= E__-- E__+=2√(1+g^2)-1,
            
            which agrees with (<ref>)–(<ref>) for N=2.
            
            
            
            It is then a standard exercise to show that the extremal values of 
            the longitudinal magnetization in the state, which is an arbitrary
            superposition of |_+ and |_-, are 
            
            M_x   =±_-|σ^x_i|_+=± M_x^
               =±1/2√(1+g/√(1+g^2))+
            √(1-g/√(1+g^2)).
            
            
            
            
            This formula  is depicted in Fig. <ref>, where one easily
            notices its departures from thermodynamic limit expression (<ref>).
            These differences do not vanish in the whole paramagnetic phase.
            In fact, for a periodic system composed of N spins, it is easy to show that the longitudinal magnetization 
             does not vanish  in the limit of g→∞
            M^_x(g→∞)=1/√(N).
            
            Such a simple result follows from elementary  observation that in the limit of g→∞ the
            positive parity ground state  approaches 
            |↑↑↑↑⋯, while the negative parity one
            approaches 
            
            |↓↑↑↑⋯+
            |↑↓↑↑⋯+
            |↑↑↓↑⋯+⋯/√(N).
            
            
            § DYNAMICS OF  PERIODIC CHAINS
            
            We will discuss here dynamics of the longitudinal magnetization in  a periodic Ising  chain. 
            There are different ways how one can approach this problem. 
            
            
            The most direct one is to work  
            in the Fock space spanned by all combinations of up/down states of all spins  such as
            |↑↓↓↑⋯. Such a  space can be
            then easily cut into positive and negative parity subspaces, where 
            time evolutions  can be independently performed.
            The advantage of such an approach is that it allows for easy computation of
            the longitudinal magnetization. The disadvantage is rather obvious: one has to
            keep track of 2×2^N-1 amplitudes, which highly limits the available
            system sizes.
            
            
            Another approach relies on mapping the system onto non-interacting fermions
            via the Jordan-Wigner transformation (<ref>). This reduces the problem to
            finding time evolution of N-1 two-level uncoupled systems (<ref>), which can be
            straightforwardly   numerically achieved. This is the most
            efficient way of getting wave-function (<ref>). The problem with this approach, however, 
            lies in the complexity of the computation of the longitudinal magnetization,
            which we  explain in Appendix <ref>.
            
            
            
            
            We will thus use the first above-mentioned approach, i.e., the direct numerical
            solution, to characterize dynamics of the N=12  periodic chain. In the
            next section, we will sacrifice the translational invariance by doing 
            calculations in much larger  chains with open boundaries,
            where such  problems with computation
            of 
            the longitudinal magnetization are absent.
            
            
            
            The results of our numerical simulations, starting from initial state
            (<ref>) and employing quench protocol (<ref>),  are presented in Figs.
            <ref>–<ref>, where we have chosen quench time τ_Q
            large-enough to ensure that  evolutions will be nearly adiabatic. We
            see there that after stopping the driving, say for 
            
            0≤ t≤ L,
            
            there are
            periodic oscillations, which will be studied in  various ways. 
            
            First, we will  extract  from numerical
            data  the difference between  positions of the first two maxima of M_x(t)
            and use it to compute   the frequency ω_12 of oscillations. We
            will also get  from such data  half of the difference between  M_x(t) in
            the first maximum and the second minimum to obtain the oscillation amplitude  A_12. 
            These quantities are illustrated in Fig. <ref>. 
            The frequency ω_12 and amplitude A_12
            estimate δ(g_f) and M_x^(g_f)–see (<ref>)
            for  justification of this statement[Note that we work with ħ=1, where  the frequency and the
            energy  can be directly compared.].
            As far as numerical simulations are concerned, the uncertainty  of getting ω_12
            and A_12 is negligible.
            This is perhaps the easiest way of
            characterization of  after-quench oscillations as it requires observation of no more than two oscillation periods.
            
            
            Second, we will fit the periodic function
            
            M_x(t)=A_cos(ω_ t+const)
            
            expecting that ω_ and A_ will estimate the same
            physical quantities as ω_12 and A_12.
            Such fitting will be done with NonlinearModelFit function from <cit.>. This 
            function also provides  uncertainty of the fitted parameters, which is  negligible in 
            our studies of  ω_ and A_.
            A few oscillation periods  provide  enough data for reaching
            the full potential of this technique.
            
            
            
            
            
            
            Third, we  will compute the discrete Fourier transform
               M̂_x(ω_j)=1/S∑_s=0^S-1M_x(sΔ t)e^ω_j sΔ t, 
               Δ t=L/S,  ω_j=2π j/L,   j=0,1,…,S-1,
            
            where S is the number of data points that we generate during free evolution.
            We will work with S being   of the order of a few tens of thousands
            (its exact value is of marginal use in the following discussion).
            We will look for the global maximum of |M̂_x| studying  its position 
             ω_max and  value |M̂_x(ω_max)|. 
             The former will estimate δ. The latter,
             after multiplication  by a factor of two, will be of the order of  M_x^. We will
             thus introduce A_max=2|M̂(ω_max)|  in analogy to the notation that we have 
             used  above.
            Both remarks  follow from (<ref>) and   (<ref>).
            
            
            
            We have applied all these techniques to numerics from Figs. 
            <ref>a–<ref>a getting results that are 
             presented in Table <ref>. We see there   very good agreement between the oscillation
            frequencies ω_12 and  ω_ and the 
            values of the symmetry-breaking gap δ from Table <ref>.
            Moreover, similar  agreement is found between the oscillation amplitudes
            A_12 and  A_fit listed in Table <ref> and M_x^ from Table <ref>.
            
            If we, however, use data from these  tables to compare  ω_max and
            A_max to the symmetry-breaking gap and the equilibrium longitudinal
            magnetization, we will find much larger relative discrepancies.
            Estimation of the symmetry-breaking gap through ω_max is mostly affected by 
             the spectral resolution of the Fourier transform  
            
            2π/L.
            
            This is seen 
            in  Figs.  <ref>b–<ref>b, where the horizontal 
            spacing between the data points ranges  between
            0.006  and 0.02, which is non-negligible  relative to δ.
            Estimation of the equilibrium longitudinal
            magnetization through A_max is additionally influenced by population of different Fourier
            modes decreasing |M̂(ω_max)| with respect to its adiabatic value.
            This can be seen through Plancherel
            theorem. Finally,  mismatch between the symmetry-breaking gap and the
            frequency grid {ω_j} additionally affects  |M̂(ω_max)|.
            
            A good agreement that we have found using the first two approaches is the
            result of choosing so large  quench time τ_Q, that  
            evolutions are nearly adiabatic.
            If we now  take a look at  
            Figs. <ref>c–<ref>c, we will be able to notice
            small deviations from perfectly adiabatic solution  (<ref>).
            Similar features will be seen in the simulations of chains with open
            boundaries, which
            brings us to the next section of this work, where they will be discussed.
            
            
            
            
            § DYNAMICS OF   CHAINS WITH OPEN BOUNDARIES
            
            We will study in this section dynamics of   Ising chains with open
            boundaries. Such systems are described  by the following Hamiltonian
            H̃(g) = -∑_i=1^N-1σ^x_iσ^x_i+1 - ∑_i=1^N g σ^z_i.
            
            
            There are at least two reasons for their consideration.
            
            First, they are  more experimentally-relevant than the 
            periodic chains  as it is a very complicated task to engineer a periodic
            coupling between the spins  (see e.g. <cit.> for an
            elaborate proposal how this might be achieved in a cold ion simulator of spin systems).
            
            Second, computations of the longitudinal magnetization can be  efficiently done
            in the  chains with open boundaries.
            It is so because the parity operator (<ref>)
            does not appear during the diagonalization 
            of Hamiltonian  (<ref>) and so both parity subspaces are diagonalized 
            with one and the same set of transformations 
            (see Appendix <ref> for explanation why this is not the case in periodic
            chains and 
            note that H̃ commutes with the parity operator).
            This allows  for efficient evaluation of the longitudinal magnetization in
            systems composed of up to a few hundreds of spins, which is a major step
            forward with respect to our studies of periodic chains.
            This flexibility with respect to the system size allows us for systematic
            studies of the transition from the non-equilibrium regime, where the
            Kibble-Zurek theory describes the system's excitation, to the nearly  adiabatic regime, which we
            have already begun to investigate in Sec. <ref>.
            
            
            
            Moving on to the actual calculations, we mention that technical details
            related to 
            diagonalization of Hamiltonian (<ref>) and time evolution that it generates are
            presented in Appendix <ref>. We only summarize  here  some basic
            formulae that are necessary for understanding of the following discussion.
            
            After Jordan-Wigner mapping (<ref>), the Hamiltonian is
            diagonalized by the Bogolubov transformation so that it finally reads
               H̃ = ∑_i=1^N E_i γ_i^†γ_i - 1/2,
               {γ_i, γ_j^†} = δ_ij,  {γ_i, γ_j } = 0,
            
            where the energies of single-particle excitations, defined in our work  with respect to
            the ground state energy, are sorted   in ascending order,
            E_1 ≤ E_2 ≤…≤ E_N,  so that 
            the symmetry-breaking
            gap is
            
                δ=E_1.
            
            
            They are given by 
            
                E_i = 2√(g^2-2gcosθ_i + 1),
            
            where  θ_i's are obtained from <cit.>
                gsin(N+1) θ = sin(N θ),    0<θ < π.
            
            
            For g>N/(N+1), θ_1,…,θ_N are the  real roots of
            (<ref>).
            For 0≤ g ≤ N/(N+1), most interestingly, there is one purely imaginary solution of (<ref>),
            Re(θ_1)=0. Besides that,
            there are N-1 real roots of (<ref>): θ_2,…,θ_N.
            Equation (<ref>) cannot be solved analytically for an arbitrary value of
            the magnetic field g. Its numerical solutions, relevant for the subsequent discussion, are
            collected in Table <ref>. 
            
            
            
            
            Besides the symmetry-breaking gap, we are also interested in the longitudinal
            magnetization, which is now position dependent. Its  equilibrium value in the  Ising chain with open 
            boundaries  was 
            recently discussed in <cit.>, where it was analyzed how 
            the “ends” of the chain affect its value. To minimize
            their influence on our results, we will focus our attention on the center of the system
            by calculating 
            
            M_x=σ^x_N/2.
            
            The  equilibrium values of such defined  longitudinal magnetization, for the parameters relevant for the subsequent
            discussion, are listed in Table <ref>.
            The technical details of
            computation of (<ref>) are discussed in Appendix <ref>.
            
            An important thing now is to note that if we write the  Schrödinger-picture
            wave-function at the time the quench stops as 
               |ψ(t=0)=∑_ a_|,
               ∑_ |a_|^2=1,
               |=γ_1^†^i_1γ_2^†^i_2⋯γ_N^†^i_N
            |GS,
               i_n=0,1  and γ_n|GS=0  1≤ n≤ N,
            
            then at times t>0
            M_x(t)=∑_a_
            a_  e^ t∑_n=1^N(j_n - i_n)E_n|σ^x_N/2|
            with
            ∑_n=1^N (i_n+j_n)  being  odd.
            
            If this condition  is not satisfied, then the matrix element in (<ref>) identically  vanishes,
            which can be shown with Wick's theorem. This feature impacts Fourier spectra
            of the longitudinal magnetization.
            
            
            
             §.§ Quenches stopping in  ferromagnetic or paramagnetic phase
            
            We start discussion of our numerical simulations  from    
            the non-equilibrium quench that has been stopped  in   the 
            ferromagnetic  phase (Fig. <ref>).
            The first thing that catches our attention in this figure is “roughness” of free evolution of M_x.
            To understand it, we need to  look at modulus of its  discrete Fourier transform. 
            
            By doing so, we first notice   a series of peaks enumerated by A, B, etc. whose
            maxima are  placed at ω_A, ω_B, etc. listed in Table
            <ref>. If we now use data from  Table <ref> and equations (<ref>) and (<ref>),
            we can note that within spectral resolution of the Fourier transform  ω_A, ω_B, ω_C, and ω_D
            can be identified with E_1, E_3-E_1-E_2, E_2, and E_3, respectively. 
            
            Occupation of the Fourier modes around ω_A is responsible for the 
            oscillation period  2π/δ(g_f)  marked in the upper panel
            of Fig. <ref>.
            Next, we note that (i)
            ω_B, ω_C, etc. are larger  by at least a factor of ten   than ω_A
            and (ii) |M̂_x(ω_A)| is larger by at least a factor of fifteen 
            than |M̂_x(ω_B)|,|M̂_x(ω_C)|, etc.
            The (i) observation means that there will be high frequency oscillations on
            top of the base  oscillation, whose frequency is approximated by ω_A. 
            The (ii) remark implies that they will
            have small amplitude relative to the amplitude of the base  oscillation. 
            Both features are nicely seen in the upper panel of Fig.
            <ref>.
            They explain small fluctuations of the data presented there.
            
            
            
            
            A quite different situation is encountered when the quench stops in the
            paramagnetic phase. This is illustrated in Fig. <ref>, where we see beats.
            A simple explanation of this observation comes again from the discrete Fourier
            transform, where we see  two leading peaks centered around ω_A and ω_B.
            The beats result from the fact that |M̂_x(ω_A)| is comparable
            to |M̂_x(ω_B)|. Moreover, we note that ω_A, ω_B, and
            ω_C are of the same order of magnitude, which eliminates high frequency oscillations
            seen in Fig. <ref>. Looking more quantitatively at the Fourier
            transform from Fig. <ref>, we notice that
            ω_A, ω_B, and ω_C can be identified with
            E_1, E_2, and E_3 within the spectral resolution of the Fourier
            transform (Tables <ref> and <ref>). 
            
            
            Finally, at the risk of stating the obvious, we mention that 
             we recover adiabatic results, akin to 
            those presented in Figs. <ref> and <ref>, by 
            increasing  the quench times from Figs. <ref> and
            <ref>.
            In the opposite limit of fast transitions, free evolution of the longitudinal
            magnetization becomes noisy and so less  interesting in the context of our  studies.
            
            
            
            
             §.§ Quenches to critical point
            
            We will discuss now the transition to the adiabatic regime for evolutions ending
             at the critical point. Such evolutions are depicted in Fig. <ref>.
             For small τ_Q, we see a train of narrow peaks, whose amplitude
              decreases as time goes by. The  magnetization in
             between the peaks is nearly zero. 
            As  evolutions slow down, the decay of the peaks' amplitude slows and the peaks' width increases shrinking  time
            intervals, where the system is  unmagnetized in the longitudinal direction. By the time those intervals
            disappear, the curve describing dynamics of the longitudinal magnetization has a triangular-like shape.
            Further increase of the quench time brings the expected single-frequency
            dynamics characteristic of the adiabatic evolution (<ref>).
            
            
            The rather unusual shape of  oscillations of the longitudinal magnetization for the fastest
            quenches from Fig. <ref> comes from substantial population of several Fourier
            modes. This is illustrated in Fig. <ref>, where the subsequent Fourier
            peaks are centered at the energies of consecutive single-particle excitations (Tables
            <ref> and <ref>).
            
            
            
            The question now is why the period of these 
            oscillations is approximately given by 2π/δ(g=1).
            This would be an expected result for  adiabatic  evolutions, where the
            oscillation pattern would be cosinusoidal (<ref>). It may  thus
            be a bit surprising that multi-frequency  oscillations from Fig. <ref>
            exhibit the same period. This can be understood 
            by noting that 
            
            E_i≈ (2i-1)δ(g=1)≈(2i-1)π/N
            for the lowest-energy modes–see the second row of Table <ref>, 
            (<ref>), and the inset of Fig. <ref>. 
            If such a relation would hold for
            all E_i's, then the oscillations of the longitudinal magnetization during
            free evolution would be perfectly adiabatic with the 2N period regardless 
            of the quench time τ_Q. 
            Relation (<ref>), however, is approximate, which we also illustrate 
            in the inset of Fig. <ref>. This  explains quasi-periodicity
            of the oscillation pattern for the fastest quench in Fig. <ref>.
            For a bit slower, but still non-adiabatic quenches depicted in this figure,
            (<ref>)  properly captures these single-particle
            excitation energies that give the main contribution to the free dynamics of
            M_x(t). This  is sufficient for explanation of the oscillation period from
            Fig. <ref>.
            
            
            
            This  oscillation period 
            is twice smaller than the oscillation  period at the critical point
            of the periodic system, which can be trivially   shown with
            (<ref>).
            The linear dependence of the oscillation period 
            on the system size, albeit with a different prefactor,  was also observed in our earlier studies 
             of the quantum Ising model <cit.>, where we investigated  free
            dynamics of the transverse magnetization  after  quenches 
            moving  the system from the paramagnetic phase to the critical point. 
            Finally, we mention  that approximation  (<ref>) works  only at  the critical
            point, which explains why different dynamics have been observed in Sec. <ref>.
            
            
            To study quantitatively the crossover  from the non-equilibrium to the adiabatic regime, 
            which we depict in Fig. <ref>,
            we
            need some measure of the deviation of non-equilibrium evolution from the
            adiabatic limit. A good measure should be easily numerically and
            experimentally accessible. It should be also stable against fluctuations of
            the data for M_x(t). Several options seem to be  available.
            
            
            
            
            
            First, one may  analyze  modulus of the discrete Fourier transform. For example, one can study 
            how the global maximum around the symmetry-breaking gap grows with increasing τ_Q.
            Alternatively, one may  research  how  the other
            extrema disappear in such a limit. This choice, however, is problematic for
            the reasons explained by the end of Sec. <ref>. 
            For example, there are limitations imposed by the  spectral resolution of the discrete Fourier transform 
            (<ref>). To overcome them, either long free evolution times are needed or
            some fitting procedure allowing for precise interpolation of the properties of
            extrema of |M̂_x| from sparse data. This is a
            complication affecting both numerical and experimental studies. The latter
            would be also affected by the fact that the Fourier transform is not directly
            measured and so its extraction out of M_x(t)
            will necessarily
            bring some inaccuracies that may play a role in the Kibble-Zurek scaling analysis.
            
            
            Second, one may use  a more straightforward approach by studying the amplitude  and
            spacing  of the first two peaks of M_x(t), just as in Sec. <ref>.
             Such a method, however, is
            susceptible to fluctuations of the data. This can be improved by averaging results 
            collected for several peaks, but this would again require long free evolution
            times, which is problematic.
            
            
            
            
            
            
            Third, one may fit  (<ref>) to  free evolution of the longitudinal magnetization and study 
            such obtained amplitude of oscillations A_ (Fig. <ref>).  
            Such a procedure uses all information contained in M_x(t)–not only the one stored 
            in the extrema of either M_x or |M̂_x|–and so long evolution times are not needed.
            Moreover, it   should work well with irregular
            data averaging out the fluctuations, which is of interest in the context of high-precision 
            numerical and experimental research. This is the approach that we will employ.
            
            Before moving on, we again mention that the fitted amplitude A_  converges to 
            M^_x(g_f) in the adiabatic limit.
            For faster quenches, however, it  underestimates the real amplitude of oscillations,
            which is seen  in Fig. <ref>. This has no effect on
            our studies, which is perhaps best illustrated by the excellent agreement 
            between the scaling exponent extracted out of the fitted amplitude and the
            predictions of the Kibble-Zurek theory (see below). 
            
            
            To proceed, we define the crossover quench time τ_Q^ by the condition
            |A_(τ_Q)-M^_x(g_f)/M^_x(g_f)|< η  τ_Q>τ_Q^,
            
            where η is the threshold set on  the relative  
            difference between the fitted amplitude of oscillations and its asymptotic in τ_Q value.
            We will use in this formula the amplitude obtained by fitting (<ref>) to M_x(t) for 0≤ t≤12N,
            which corresponds to roughly 6 oscillation periods in the  chain with open
            boundaries. 
            Moreover, we will set η=10%, which should be large-enough to be
            experimentally-relevant and small-enough to describe the crossover to the adiabatic limit.
            
            
            Our results for A_, in the experimentally-relevant system composed of N=50 spins 
            <cit.>, are presented in 
            Fig. <ref>a, where we see that the fitting procedure produces a
            perfectly smooth curve monotonically approaching M^_x(g_f). The threshold
            η  is illustrated in Fig. <ref>b. 
            
            
            Repeating such analysis for system sizes 20≤ N≤ 300, we have obtained detailed
            results for the crossover quench time τ_Q^, which we present in Fig.
            <ref>a. As we anticipate   from (<ref>) that 
            τ_Q^(N)∼ N^a,
            
            where a>0 is the scaling exponent, we  display results for τ_Q^ on a double logarithmic plot in Fig.
            <ref>b.
            Instead of  a straight line, we find in this figure a curve slightly bending upwards as the system size grows.
            This means that the  exponent a  increases with N. To quantify this observation, we fit
            lnτ_Q^ = aln N + b
            
            to numerical data from four  different ranges of the system sizes. The results are collected in Table <ref>,
            where we see that a approaches the value of 2 for the largest system sizes that we consider. 
            This is in a very good agreement with the Kibble-Zurek
            scaling argument (<ref>), which is supposed to work best in 
            the large-system limit. Finally, we notice that the increase of a with N is monotonic, leaving no
            doubts about  stability of the procedure of extraction of τ_Q^ from
            the free evolution data for the longitudinal magnetization.
            
            
            
            
            
            
            
            
            
            § DISCUSSION
            
            
            
            The goal of this work was  to investigate how the symmetry-breaking gap,
            which is of crucial importance in the discussion of quantum phase transitions, 
            can be studied with the help of quantum quenches.
            We have chosen for this purpose an exactly solvable model, the quantum Ising model in the transverse field,
            and analyzed   its dynamics  after  quenches
            induced by the gradual change of the magnetic field. These  quenches
            start from the easy-to-prepare broken-symmetry ground state at the vanishing magnetic field.
            They  bring the
            system to the desired value of this field, where the symmetry-breaking gap
            can be  read from the  subsequent   free-evolution dynamics of the longitudinal
            magnetization. In this way a small symmetry-breaking gap can be seen through 
            large-amplitude oscillations of the longitudinal magnetization.
            
            We have discussed different ways of analyzing the oscillatory dynamics of such 
            magnetization showing that one can also extract the 
            equilibrium longitudinal magnetization out of them. All this can be accurately
            done if the quench is slow enough, which we have  studied
            in the context of the Kibble-Zurek theory
            of non-equilibrium phase transitions. An excellent agreement between
            predictions of this theory and the dynamics of the longitudinal magnetization
            has been found. 
            
            
            
            
            
            Although our studies have been done in the Ising model, they can be
            extended to other systems exhibiting the symmetry-breaking phenomenon. 
            For example, the Ising-like ones  with long-range interactions 
            that are typically found in cold ion and atom  emulators of spin systems
            (see e.g. <cit.>). 
            These systems provide a promising platform for experimental realization of the
            studies discussed in our work for two reasons.
            First, their size is  finite  rather than thermodynamic making  their 
            symmetry-breaking gap large-enough to be experimentally measurable. 
            This should not be taken for granted because  it is not the case in traditional
            condensed matter setups discussed in the context of phase transitions. 
            Second,  there has been 
            substantial   progress in the experimental studies of the dynamics  of 
            such systems (see e.g. <cit.>). 
            
            Another promising platform for
            implementation of our ideas is provided by D-Wave machines, which can be also
            used for simulations of   spin models (see e.g. <cit.>). 
            Quite interestingly, 
            D-Wave-based investigations of non-equilibrium Kibble-Zurek dynamics of the one-dimensional quantum Ising
            model in the transverse field were  recently reported  in <cit.>.
            It should be thus possible to use  our predictions  for critical assessment 
            of the performance of such devices.
            
            
            Talking about experimental realizations, the symmetry-breaking phenomenon
            was recently experimentally investigated in <cit.>.
            These studies were done in a cold atom cloud, where each atom was simulating 
            the 16-spin Lipkin-Meshkov-Glick model. This is the Ising-like model
            with the  nearest-neighbor spin-spin interactions replaced with identical couplings 
            between all the spins. Besides exploration of a different Hamiltonian,
            these interesting studies differ from our work  in 
            the following aspects. First, the quenches   
            start in the paramagnetic phase. Second, the initial state for
            them occupies one of the parity subspaces  and so the symmetry-breaking 
            perturbation is used to populate the other one as well. Third, perhaps most importantly, 
            they are limited to one, rather small system size, and they do not explore 
            the non-adiabatic Kibble-Zurek dynamics of the longitudinal magnetization,
            which is of considerable  importance in our work.
            
            Finally,   we mention that we  hope that this work will trigger interest in 
            the experimental studies of the symmetry-breaking phenomenon, which could lead
            to  quantitative insights into the very nature of quantum phase transitions.
            This would be most interesting 
            in systems  that can be neither analytically nor numerically studied  in the foreseeable
            future.
            
            
            
            
            § ACKNOWLEDGMENTS
            
            We thank Marek Rams for  useful discussions, sharing with us the results of
            his numerical simulations,  reading the manuscript, and providing us with most
            useful feedback.
            MB and BD were supported by the Polish National Science Centre (NCN) grant
            DEC-2016/23/B/ST3/01152. 
            
            
            
            § LONGITUDINAL MAGNETIZATION IN JORDAN-WIGNER-TRANSFORMED   PERIODIC ISING CHAINS
            
            We will explain here  the source of difficulties in computation of the longitudinal magnetization 
            in the free fermion  representation of the periodic Ising model.
            While doing so, we will also derive expressions for dynamical gaps
            Δ_±, which have been introduced in Sec. <ref>.
            
            Such a representation  comes from  employment of the Jordan-Wigner transformation
               σ^z_i=1-2 c_i^† c_i,  σ^x_i=(c_i+c_i^†)∏_j<i1-2 c_j^† c_j,
               {c_i, c_j^†}=δ_ij,  {c_i, c_j}=0,
            
            after which  Hamiltonian (<ref>) reads
            
             H(g) = -∑_i=1^N-1  f_i,i+1 - g ∑_i=1^N (c_i  c_i^†  -  c_i^†  c_i)+ f_N,1 P,
              f_i,j =  c_i^† c_j -  c_i  c_j^† -  c_i c_j +   c_i^†  c_j^†.
            
            The complication one encounters now is that the parity operator  is not quadratic in 
            fermionic c_i and c_i^† operators. This  can be found by combining
            (<ref>) with (<ref>). As a result, Hamiltonian (<ref>)
            is non-quadratic. Its exact analytical diagonalization  is still possible, but 
            the price one has to pay is that one has  to split the Hilbert space into
            positive and negative parity subspaces imposing different boundary conditions
            on the fermionic  operators in those subspaces <cit.>. 
            This is realized by going to the momentum space 
            
            c_j =e^-π/4/√(N)∑_k c_k e^ kj,
            
            and choosing different quantization 
            schemes for momenta in the positive and negative parity subspaces, 
            
            k_+=±π/N,±3π/N,…,±N-1/Nπ, 
            
            and 
            
            k_-=0,±2π/N,±4π/N,…,±N-2/Nπ,π,
            
            respectively. These expressions are valid for even N 
            (see <cit.> for comprehensive discussion of the Ising
            diagonalization intricacies). 
            Therefore, different transformations are used to diagonalize 
            the two parity subspaces in periodic chains. No such complications 
            appear in the chains with open boundaries, where the 
            parity operator is absent  in  Jordan-Winger-transformed Hamiltonian
            (<ref>).
            
            Combining the results of 
             <cit.>
            and
            <cit.>, time-dependent wave-function (<ref>) can be  obtained from
            
                |ψ_+(t)⟩ =∏_k_+>0u_k_+(t) - v_k_+(t)
                c_k_+^† c_-k_+^†| vac⟩,
                   |ψ_-(t)⟩ =e^2 tc_0^†∏_0<k_-<πu_k_-(t) - v_k_-(t)
                c_k_-^† c_-k_-^†| vac⟩,
            
            where the state |vac  is annihilated by all c_k
            operators and  time evolution of the Bogolubov modes is governed by
            d/dt(
            [ v_k; u_k ])
            =2(
            [ g-cos(k)  -sin(k);  -sin(k) cos(k)-g ])
            (
            [ v_k; u_k ]).
            
            
            Computation of the time-dependent longitudinal magnetization in our periodic system is now reduced to evaluation 
            of 
            
            M_x(t)=ψ_-(t)|c_1+c_1^†|ψ_+(t)+
            If we now put (<ref>) and  (<ref>) into (<ref>) and then 
            transform  c_k_± operators to the position space inverting (<ref>)–so that all
            operators are defined  on the same Hilbert space–we will quickly realize how
            complicated the resulting  expression is. This obstacle discourages us from
            using the free fermion representation in our studies of periodic chains.
            
            Finally, having (<ref>) at hand, one can find by diagonalization of the 2×2
            Hamiltonian that the energy gap for excitation of the pair of ± k fermionic modes is
            
            4√(g^2-2gcos(k)+1).
            
            The smallest value of that gap, in each of the parity subspaces, is 
            the dynamical gap Δ_±–see also Fig. <ref>. Thereby,
            Δ_+=4√(g^2-2gcos(π/N)+1),  Δ_-=4√(g^2-2gcos(2π/N)+1).
            
            
            § DIAGONALIZATION AND TIME EVOLUTION OF ISING CHAINS WITH OPEN BOUNDARIES
            
            We will briefly summarize here technicalities related to diagonalization
            of the  Ising chain with open boundaries, computation of its longitudinal magnetization, and
            its time evolution.
            
             Diagonalization. We follow here <cit.> providing  an early take on this subject. 
            The diagonalization begins  with Jordan-Wigner spin-to-fermion mapping (<ref>)
            transforming   Hamiltonian (<ref>) to the following quadratic form
               H̃ = Ψ^† HΨ,
               Ψ^† = (c_1^†… c_N^†  c_1 … c_N),
                H =  ( [  A  B; -B -A;    ]),
            
            where A and B are N × N tridiagonal matrices
               A_ij = g  δ_i,j - 1/2  δ_i,j+1 -1/2  δ_i+1,j,
               B_ij = 1/2  δ_i,j+1 -1/2  δ_i+1,j.
            
            We mention in passing that we have corrected a misprint from  <cit.> in the expression for
            B_ij. 
            
            Next,  for every value of the magnetic field g, one  can perform the Bogolubov transformation
               Ψ = βΓ,
               Γ^† = (γ_1^†…γ_N^† γ_1 …γ_N)
            
            choosing real orthogonal matrix β in  such a way that (<ref>) is diagonal. This leads
            to (<ref>) and the related equations (<ref>)–(<ref>). Two remarks are in order now.
            
            First, as a self-consistency check of our calculations, we have 
            verified  that energies of single-particle excitations, which  we have
            obtained
            from numerical 
            diagonalization of H, very well  agree  with the results coming from 
            (<ref>) combined with (<ref>).
            
            Second, the symmetry-breaking gap can be analytically  calculated 
            from (<ref>) only at the critical point, where 
            
                δ(g=1)  = 4 sinπ/4N+2≈π/N.
            
            This  is about twice larger than the symmetry-breaking gap in the periodic
            chain  (<ref>).
            
             Equilibrium longitudinal magnetization. 
            To compute the equilibrium longitudinal magnetization M_x^, we evaluate 
            (<ref>) in the state 
            | GS ⟩ + γ_1^†| GS ⟩/√(2),
            
            where | GS⟩ is the ground state annihilated by all γ_i operators. 
            This leads to the following expression for the longitudinal magnetization 
            after employment of  Wick's theorem  
            
                M_x^ = ⟨ GS |γ_1  σ^x_N/2| GS ⟩.
            
            The operator γ_1 σ^x_N/2 can be conveniently written as
            
                γ_1 σ^x_N/2 = γ_1   a_N/2  b_N/2-1 a_N/2-1… b_1   a_1,  a_i = c_i + c_i^†,
                 b_i = c_i - c_i^†.
            
            Using Wick's theorem again, one can show that 
               ⟨ GS |γ_1  σ^x_N/2| GS ⟩ = Pf(G),
               G=(
            [              0       γ_1a_N/2     γ_1b_N/2-1     γ_1a_N/2-1              ⋯         γ_1a_1;                             0   a_N/2b_N/2-1   a_N/2a_N/2-1              ⋯       a_N/2a_1;                                            0 b_N/2-1a_N/2-1              ⋯     b_N/2-1a_1;                                                                                         ⋮;                                                                                         0 ]),
            
            where Pf  stands for Pfaffian, the lower triangle of the G matrix can be obtained by the relation G=-G^T, and 
            the expectation values  are calculated in the ground state |GS. 
            Pfaffians of  skew-symmetric matrices can be efficiently  computed using Hausholder transformation <cit.>. 
            
             Time evolution. We work in the Heisenberg picture. Our evolutions start
            at  time t_0 from the equal superposition of the two  lowest-energy eigenstates of
            H̃[g(t_0)]| GS[g(t_0)]⟩ + γ_1^†| GS[g(t_0)] ⟩/√(2).
            
            Two remarks are in order now.
            
            First,  initial state (<ref>)  for time evolution is constructed in the same way as 
            for  evolutions in periodic chains.
            In particular, the two states in  (<ref>) have different
            parities. In fact, it is perhaps worth to say again  that  Hamiltonian H̃ 
            for the  Ising chain with open boundaries  commutes with the parity operator.
            Therefore, its
            eigenstates can be labeled with the ±1  parities.
             Moreover,  expectation values of the σ^x_i operators in all
            eigenstates of H̃ are zero. The very same properties are found  in  periodic
            chains, which have been  discussed in Secs. <ref>–<ref>. 
            
            Second, using  quench protocol (<ref>), one gets g(t_0=-g_fτ_Q)=0 and the initial state (<ref>) is 
            given by (<ref>).  It is numerically convenient 
            for us, however,  to begin  evolutions from the  slightly non-zero g, which
            we do by choosing t_0 such that g(t_0)=0.001. 
            
            Time-dependent longitudinal magnetization  (<ref>)  is
            then  expressed as
            
            M_x(t) = Re⟨ GS[g(t_0)] |γ_1  σ^x_N/2(t) | GS[g(t_0)] ⟩,
            
            where the operator γ_1 is defined at time t_0.
            The matrix element in this equation can be  computed just as (<ref>) expect a_i and b_i
            operators are now time dependent. Thus,  we need to know their time evolution,
            which can be extracted from 
            
                Ψ(t) = U(t)Ψ(t_0),
            
            where the N× N unitary matrix U(t) can be obtained by solving 
            d/dt U(t) = -2 H  U(t)
            
            with the initial condition U(t_0)=1. Equation (<ref>) can be derived from 
            the Heisenberg equations for the c_i(t) and c_i^†(t) operators.
            We solve it numerically with the  Suzuki-Trotter method of order two with the 
            time step smaller or equal to 0.01<cit.>. 
            We have checked that such a procedure produces well-converged results.
            
            Having U(t) and the Bogolubov matrix β at the  time t_0,
            we can relate c_i(t) and c_i^†(t) operators to γ_i and γ_i^† appearing 
            in the diagonal form of H̃[g(t_0)]. 
            Namely,
            
                Ψ(t) = W(t)  Γ,   W(t) = U(t)  β,
            
            where the  matrix W has the following structure
            
                W = 
                ([ C D; D C ])
            
            with C and D being N× N complex matrices.
            Transformation (<ref>)  can be used  to compute all  correlation functions, 
            from time-dependent version of (<ref>), 
            needed for getting  (<ref>).
            For example, after straightforward manipulations one can show that 
            
                a_ia_j = ∑_k=1^N 
                C_ik  D_jk + 
                C_ikC_jk + 
                D_ik  D_jk +  
                D_ik C_jk .
            48 natexlab#1#1 bibnamefont#1#1 bibfnamefont#1#1 citenamefont#1#1 url<#>1 urlprefixURL [Sac(a)]SachdevnoteS. Sachdev Quantum Phase Transitions (Cambridge
              University Press, 2011).[Con()]ContinentinoBooknoteM. Continentino  Quantum Scaling in Many-Body Systems: An
              Approach to Quantum Phase Transitions (Cambridge University Press, 2nd
              edition, 2017).[Sac(b)]SachdevTodaynoteS. Sachdev and B. Keimer, Phys. Today  64, 29 (2011).[Lew()]LewensteinAdvnoteM. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen De,
              and U. Sen, Adv. Phys.  56, 243 (2007).[Jac(a)]JacekPRL2005noteJ. Dziarmaga, Phys. Rev. Lett.  95, 245701 (2005).[Dor()]DornerPRL2005noteW. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. 
              95, 105701 (2005).[Pol(a)]PolkovnikovPRB2005noteA. Polkovnikov, Phys. Rev. B  72, 161201(R) (2005).[Ral()]RalfPRA2007noteS. Mostame, G. Schaller, and R. Schützhold, Phys. Rev. A 
              76, 030304(R) (2007).[Pol(b)]PolkovnikovPRL2008noteR. Barankov and A. Polkovnikov, Phys. Rev. Lett.  101,
              076801 (2008).[Sen(a)]SenPRB2009noteS. Mondal, K. Sengupta, and D. Sen, Phys. Rev. B  79,
              045128 (2009).[San(a)]SantoroPRB2009noteD. Patanè, L. Amico, A. Silva, R. Fazio, and G. E. Santoro,
              Phys. Rev. B  80, 024302 (2009).[Jac(b)]JacekPRA2007noteL. Cincio, J. Dziarmaga, M. M. Rams, and W. H. Zurek, Phys. Rev.
              A  75, 052321 (2007).[Sen(b)]SenPRA2009noteK. Sengupta and D. Sen, Phys. Rev. A  80, 032304 (2009).[Arn()]ArnabPRB2010noteA. Das, Phys. Rev. B  82, 172402 (2010).[Kol()]KolodrubetzPRL2012noteM. Kolodrubetz, B. K. Clark, and D. A. Huse, Phys. Rev. Lett.
               109, 015701 (2012).[Fra()]FrancuzPRB2016noteA. Francuz, J. Dziarmaga, B. Gardas, and W. H. Zurek, Phys. Rev.
              B  93, 075134 (2016).[San(b)]SantoroJstat2015noteA. Russomanno, S. Sharma, A. Dutta, and G. E. Santoro, J. Stat.
              Mech. (2015) P08030.[Pus()]PuskarovSciP2016noteT. Puskarov and D. Schuricht, SciPost Phys.  1, 003
              (2016).[Apo()]ApollaroSciRep2017noteS. Lorenzo, J. Marino, F. Plastina, G. M. Palma, and T. J. G.
              Apollaro, Sci. Rep.  7, 5672 (2017).[Mic()]MichalOnenoteM. Białończyk and B. Damski, J. Stat. Mech. (2018)
              073105.[Ado()]AdolfoPRL2018noteA. del Campo, Phys. Rev. Lett.  121, 200601 (2018).[Mar()]MarekArxivnoteM. M. Rams, J. Dziarmaga, and W. H. Zurek, Phys. Rev. Lett.
               123, 130603 (2019).[Blo()]BlochSci2017noteC. Gross and I. Bloch, Science  357, 995 (2017).[LMG()]LMGexpnoteV. Makhalov, T. Satoor, A. Evrard, T. Chalopin, R. Lopes, and S.
              Nascimbene, Phys. Rev. Lett.  123, 120601 (2019).[Lie()]Lieb1961noteE. Lieb, T. Schultz, and D. Mattis, Ann. Phys. (N.Y.)  16,
              407 (1961).[Pfe()]PfeutynoteP. Pfeuty, Ann. Phys.  57, 79 (1970).[Wic()]WicknoteG. F. Giuliani and G. Vignale,  Quantum Theory of the
              Electron Liquid (Cambridge University Press, 2005).[BDJ()]BDJPA2014noteB. Damski and M. M. Rams, J. Phys. A  47, 025303 (2014).[BDP()]BDPRL2005noteB. Damski, Phys. Rev. Lett.  95, 035701 (2005).[Kib()]KibbleRevnoteT. W. B. Kibble, Phys. Rep.  67,183 (1980).[Zur()]ZurekRevnoteW. H. Zurek, Phys. Rep.  276, 177 (1996).[del(a)]delCampoRevanoteA. del Campo, T. W. B. Kibble, and W. H. Zurek, J. Phys.:
              Condens. Matter  25, 404210 (2013).[del(b)]delCampoRevbnoteA. del Campo and W. H. Zurek, Int. J. Mod. Phys. A  29,
              1430018 (2014).[Jac(c)]JacekAdv2010noteJ. Dziarmaga, Adv. Phys.  59, 1063 (2010).[Mat()]MathematicanoteWolfram Research, Inc., Mathematica, Version 11.0, Champaign, IL
              (2016).[Kor()]KorenblitNJP2012noteS. Korenblit  et al., New J. Phys.  14 095024 (2012).[Cab()]CabreraPRB1987noteG. G. Cabrera and R. Jullien, Phys. Rev. B  35, 7062
              (1987).[Osk()]OskarnoteO. A. Prośniak, Phys. Scr.  94, 085201 (2019).[Luk(a)]LukinNature2017noteH. Bernien  et al., Nature  551, 579 (2017).[Mon()]MonroeNature2017noteJ. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H.
              Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, Nature  551, 601
              (2017).[Luk(b)]LukinNature2019noteA. Keesling  et al., Nature  568, 207 (2019).[Por()]PorrasPRL2004noteD. Porras and J. I. Cirac, Phys. Rev. Lett.  92, 207901
              (2004).[Sch()]SchaussQST2018noteP. Schauss, Quantum Sci. Technol.  3, 023001 (2018).[Har()]HarrisSci2018noteR. Harris  et al., Science  361, 162 (2018).[Jac(d)]JacekSciRep2018noteB. Gardas, J. Dziarmaga, W. H. Zurek, and M. Zwolak, Sci. Rep.
               8, 4539 (2018).[You()]YoungPRB1997noteA. P. Young, Phys. Rev. B  56, 11691 (1997).[wim()]wimmernoteM. Wimmer, ACM Trans. Math. Software  38, 30 (2012).[tro()]trotternoteN. Hatano and M. Suzuki, Finding Exponential Product Formulas of
              Higher Orders, in  Quantum Annealing and Other Optimization Methods,
              edited by A. Das and B. K. Chakrabarti (Springer, 2005), pp. 37-68.
