Exact multistability and dissipative time crystals in interacting fermionic lattices

The existence of multistability in quantum systems beyond the mean-ﬁ eld approximation remains an intensely debated open question. Quantum ﬂ uctuations are ﬁ nite-size corrections to the mean-ﬁ eld as the full exact solution is unobtainable and they usually destroy the multistability present on the mean-ﬁ eld level. Here, by identifying and using exact modulated dynamical symmetries in a driven-dissipative fermionic chain we exactly prove multistability in the presence of quantum ﬂ uctuations. Further, unlike common cases in our model, rather than destroying multistability, the quantum ﬂ uctuations themselves exhibit multistability, which is absent on the mean-ﬁ eld level for our systems. Moreover, the studied model acquires additional thermodynamic dynamical symmetries that imply persistent periodic oscillations, constituting the ﬁ rst case of a boundary time crystal,to the best of our knowledge, a genuine extended many-body quantum system with the previous cases being only in emergent single-or few-body models. The model can be made into a dissipative time crystal in the limit of large dissipation (i.e. the persistent oscillations are stabilized by the dissipation) making it both a boundary and dissipative time crystal.


M
ultistability in driven-dissipative models usually means the presence of two possible stationary states of the system that can be distinguished by local observable measurements.Although on the level of the mean-field approximation or for classical systems it can be easily established whether or not it exists 1,2 , its actual existence, in particular in lowdimensional strongly interacting systems remains quite controversial with both theoretical and experimental works reporting differing conclusions  . The xisting approaches usually rely on sophisticated theoretical techniques for including perturbations of finite-size corrections to the mean-field or large scale efficient numerical simulations such as t-DMRG 10 or projected entangled pair states (PEPS) 21 .The general lore in the literature is that the lower dimension, the more likely it is that the full quantum fluctuations (finite-system size corrections) will destroy the multistability and restore the generic unique stationary state of the model [27][28][29][30] .However, in general, the question remains unsettled in any dimension (see e.g., ref.21 for an advanced numerical study in two dimensions).Therefore, exact results on this controversial problem are desirable.
Here, we report an exactly solvable model describing a broad class of driven-dissipative and interacting fermionic chains that shows quantum multistablity related to the strong symmetries of the Liouvillian.In general, the strong symmetries are not necessary for the existence of multistability (as degenerate stationary states may emerge in the thermodynamic limit only without a strong symmetry), nor are they sufficient since all of the degenerate stationary states implied by the strong symmetries may have the same expectation values for local observables.More generally, multistability is taken to be an effect not due to any manifest symmetry.However, in our case emergent strong symmetries do guarantee multistability in local observables.In contrast to other potentially bistable models 27 , here, rather than being detrimental, quantum fluctuations are essential for the multistability.Further, we show that for certain types of jump operators the Liouvillian has dynamical symmetries that lead to sustainable oscillations in the non-equilibrium long-time behavior of the system, a phenomenon known as dissipative time crystal.Our approach is based on identifying a novel modulated 76 spectrum generating algebra (SGA) 77 , which is non-local and fermionic.Since standard dynamical symmetries 31 are extensive and local SGA 78,79 and the SGA here is extensive and semi-local in spin, we call it a semi-local dynamical symmetry 80 .This dynamical symmetry, being fermionic, cannot be relegated to a non-Abelian symmetry, unlike previously known cases based on closed algebras settling the question of whether such operator relations are possible 81 .Later, for sake of simplicity, we specialize the general dissipativedriven fermionic model to a quadratic model and show that these models have an infinite set of emergent (thermodynamic) superextensive raising operators that we call super-extensive dynamical symmetries.We provide evidence that the total effect of all these operators is that the model displays very slow finite-size decay due to the presence of strong symmetries 82 , which guarantee degenerate stationary states (null space of Liouvillian).These have attracted lots of interest recently due to their utility for quantum information storage [83][84][85][86][87][88][89][90][91][92][93][94][95][96] .Moreover, as this behavior is robust to a wide-class of perturbations and occurs in the thermodynamic limit any for generic parameter values, it also constitutes the first example of a boundary time crystal, to the best of our knowledge, in an short-range interacting model.We also derive an analytical lower bound for the decay of such behaviors when the dynamical symmetry is not exact.Finally, we support our theoretical findings using a phase space approach to numerically calculate the Liouvillain spectrum as well as the correlation functions dynamics of the open quantum system dynamics highlighting the emergence of the dissipative time crystalline and the quantum multistabiliy as well as the scaling with the system size.

Results and discussion
The model.Consider the following interacting Kitaev chain model 97 , with c y j and c j being the (Dirac) fermionic creation and annihilation operator, w is the hopping amplitude, Δ is the p-wave pairing correlation, μ is the on-site chemical potential, and V is a novel interacting (beyond-quadratic in general) term we use to model the stronly interaction Majorana fermions, physically originated from electron-electron interactions 98 .The number operators are n j ¼ c y j c j and the chain is subject to the periodic boundary condition, i.e., c N+1 = c 1 .We note that complex values of coupling coefficients can be physically realized with e.g., constant phase gradients 99 or laser coupling as recently employed in the simulation of a bosonic ladder 100 .
This model is commonly mapped onto a spin-1/2 transverse field Ising model at V = 0 101 .Note that the model we consider here is not mappable to integrable XYZ spin chains or noninteracting models in contrast to other interesting results 102,103 .
As will be apparent shortly, it is convenient to define Majorana fermion operators as, These fulfill the following anti-commutation relations, Furthermore, here we will study a dissipative case with an incoherent Markovian driving modeled by a Lindblad master equation, where L is a quantum Liouvillian of the form and the Lindblad jump operators are local Majorana dissipation, as incoherent local single-fermion drive and dissipation of the form The quadratic version of this model (V = 0) has been previously studied for its topological properties 104 .We will also consider a special type of dissipative pairing as can be engineered with Pauli blocking 105 .
It is useful to define a vector space of operators with the standard Hilbert-Schmidt inner product, as hhAjBii ¼ trðA y BÞ, with a corresponding norm that we will use.
In order to solve the dynamics of ρ(t) it is useful to diagonalize L. Define λ k to be the eigenvalues of L and ρ k , σ k to be the corresponding right and left eigenoperators, respectively, Due to the semi-group properties of the Lindblad master equation all the eigenoperators are either stable or decaying, i.e., Re(λ k ) ≤ 0. Further, they always appear in complex conjugate pairs as fλ k ; λ Ã k g.Since the jump operators L μ are Hermitian in our model, the Lindblad equation is unital, and the identity matrix is always a non-equilibrium steady state (NESS), i.e., ρ 0 ¼ 1.
We are interested in the dynamics of observables O(t) when we initalize the system in ρ(0).Formally, the solution can be written as Purely imaginary eigenvalues are therefore necessary but not sufficient for persistent oscillations in physical observables due to possibly vanishing overlap of local observables with right eigenoperator, or the zero overlap of the initial state with left eigenoperators, or the presence of dense and incommensurate purely imaginary eigevalues λ k in the Liouvillian spectrum.
Emergent dynamical symmetries in the thermodynamic limit.
The origin of dynamical symmetries can be understood by studying the quadratic version of the model in ( 1) in the noninteracting limit, i.e., H 0 = H(V = 0).The model is then the standard Kitaev Hamiltonian that in the thermodynamic limit can be diagonalized with a Fourier transform followed by a Bogoliubov transformation.We obtain (up to an irrelevant shift), with the lowering operators of, Here, c(k) is the Fourier transform of the fermion annihilation operator at momentum k and d(k) is its Bogoliubov transformation where the (not-normalized) coefficients are defined as and the energy is, The momentum is restricted to the first Brillouin zone k coinciding with the topological phase of the Kitaev chain.Using (3) it is clear that even Majorana fermions γ 2j commute with any product of an even number of odd Majorana fermions γ 2j−1 .That means the interaction term of the Hamiltonian in (1) commutes with d κ from which it directly follows that, where q .Thus d κ is a modulated fermionic dynamical symmetry of the model.They are essentially similar to Goldstone modes, except they exist at finite frequency and momentum.The dynamical symmetry implies that observables with non-zero overlap with it can persistently oscillate at frequency E κ 81 .Note that there are special values of parameters for which the dynamical symmetries exist for finite systems, e.g., for μ = 0, κ = π/2 is always a solution for mod(N, 4) = 0.For more general μ, w (with |μ| < |w|); however, the dynamical symmetries exist only in the thermodynamic limit as for finite systems there is no solution for κ for general μ, w.In other words, these dynamical symmetries are emergent for systems that are large enough to have a continuum of momenta k in the 1st B.Z. Hence, they are thermodynamically emergent symmetries.Now consider the case Γ = 0. Since the jump operators L j,2 are as sums of products of odd Majorana fermions one has [L j,2 , d κ ] = 0. Therefore, S ≔ d κ is an emergent strong dynamical symmetry of the open quantum system with dissipators L j,2 in the thermodynamic limit 31,106,107 , which implies that operators S n ρ 1 ðS y Þ m are eigenoperators of L with purely imaginary eigenvalues i(n − m) E κ .Here, n, m = 0, ± 1 as S 2 ¼ 1.Note that S is local in the fermion basis in the sense of being a sum of local densities.
Here, the eigenstate dephasing 108 is not possible because the purely imaginary eigenspectrum is equally spaced.In order to qualify as a boundary time crystal the system must, in addition to the previous conditions, have persistent oscillations in local observables 47 .This is the case because S is local and has non-zero overlap with local observables meaning that the expectation values of a local observable (by ( 7)), will be finite and persistently time-periodic in the thermodynamic limit.This is clearest when considering the example O = γ 2j with the asymptotic state ρðtÞ ¼ 1 þ e iE κ t S þ e ÀiE κ t S y .This therefore constitutes, an exact example of a boundary time crystal in an extended many-body system with local interactions.Besides, we conjecture that the model is a dissipative time crystal as well (persistent oscillations induced by dissipation) because the closed model likely has oscillations that dephase due to the multitude of incommensurate eigenvalues 108 .We will discuss robustness of the boundary time crystal state in the next section.
Semi-local dynamical symmetries and stability of the dynamics.We will now consider the case when Γ ≠ 0 and in order to simplify the following discussion we will assume without loss of generality that the parameters of the model are such that the dynamical symmetry exists for some finite values of N, e.g., at μ = 0 as discussed in the previous subsection.
Let us define m j ¼ 1 2 1 À n j and the parity operator We now note some useful identities.The Hamiltonian is parity-symmetric, i.e., [H, P 1,N ] = 0, and the Lindblad jump operators are parity-antisymmetric, i.e., {L μ,1 , P 1,N } = 0 and furthermore satisfy {L μ,1 , Remark 1: As a side note we will show that this operator leads to a semi-local dynamics symmetry when mapped to the spin systems.
Here, for reasons that will become apparent, we consider the standard Wigner-Jordan mapping 109 , Following the mapping in (15), operator A 0 gets transformed to a semi-local dynamical symmetry in the spin basis, i.e., its densities commute only with operators on one side of the chain with the following form, Such semi-local symmetry operators have been studied recently in the context of generalized hydrodynamic corrections in quadratic and integrable models where their existence was associated with the topological nature of the models 80 .Topology likely plays a role in our model as it is intimately related to the Kitaev chain.These new kinds of dynamical symmetries should be distinguished from both local extensive 78,[110][111][112] and strictly local ones [113][114][115] .We emphasize, that in the model studied here there is no obvious transformation that would allow mapping this semi-local dynamical symmetry into a semi-local non-Abelian symmetry while preserving the spatial locality of H.
The existence of A 0 immediately implies the existence of a super-extensive (quadratic) charge Q ¼ ðA 0 Þ y A 0 for which it can be easily shown that trðQ 2 Þ=2 N / N 2 by observing that the expression has only finite trace for those term in the doubled sum when local densities are on the same site.For example, for μ = 0, in the Majorana basis it may be written as, This Hermitian operator defines a symmetry S = e iQ as [H, S] = 0. We remark that for the closed system (Γ = ϵ = 0), the existence of Q and its growth with system size immediately implies that the memory of the initial state decays as 1/N 2 in local observables O with trðQOÞ ≠ 0 as quantified by e.g., infinite temperature autocorrelation functions 〈O(t)O〉 via the Mazur bound 116 .For the open system, S is a strong symmetry 82,117 , which in turn implies that the non-equilibrium stationary state λ 0 = 0 is degenerate.There, ρ 0 0 ¼ 1 À ðA 0 ðA 0 Þ y Þ is the other eigenoperator of the Liouvillian null space for λ 0 = 0 where we define 1 :¼ 1=ð2 N tr½ðA 0 ðA 0 Þ y ÞÞ, in order for it to be Hilbert-Schmidt orthogonal to the left null space eigenmode ω 0 ¼ 1, i.e., trðρ 0 0 Þ ¼ 0. It is worth reminding that due to unitality of the Liouvillian, one NESS is always ρ Furthermore, based on the properties of A 0 , ρ 1 ¼ ðA 0 Þ y , ρ À1 ¼ A 0 are the bi-orthogonal eigenmodes corresponding to the purely imaginary eigenvalues of λ ±1 = iE κ .Due to unitality the left and right eigenmodes are each others conjugate transposes, i.e., σ Ç1 ¼ N ðρ ± 1 Þ y up to a normalization constant N .
It is obvious that the presence of purely imaginary eigenvalues is not visible in any local observable O because A 0 is non-local as it does not contain any local terms hence, ρ ±1 do not have Hilbert-Schmidt overlap with such observables.
The multistability situation is different however, as ρ 0 0 has a non-zero overlap with a local observable O (without loss of generality tr(O) = 0).In order to estimate its multistability, we first note that a general stationary state must be density matrix and hence a convex combination of the two stationary states N and the eigenvalues of ρ 0 0 sum into 0 by construction, c ∝ 1/N.However, since for local observables trðOρ 0 0 Þ / N 0 we have for the expectation values, hOðt !1Þi / hhOjcρ 0 0 ii / 1=N; ð18Þ which gives a lower bound on the |〈O(t → ∞)〉|.We note that this is only a lower bound because many eigenvalues have vanishing real part in the limit of N → ∞ and hence the actual expectation values can decay slower with N. Thus the model is bistable only for the finite-sized chains where quantum fluctuations are important.This is unlike typical mean-field multistability scenarios occurring in the thermodynamic limit (i.e., N → ∞) where quantum fluctuations do not play a role.Expanding for small ε = O(1/N) around k = κ ± ε we have ½L μ ; P 1;N d κ ± ε ¼ OðεÞ (by Taylor series expansion) and, likewise, [H, this implies that the dissipative gap (the real part of the eigenvalues of the Liouvillian) closes as ∝ 1/N into the same purely imaginary eigenvalues.They are hence metastable [118][119][120] .We will confirm this in the next section with a concrete example.The closing of the Liouvillian gap is associated with an algebraic temporal relaxation of the dynamics, but as we will see in the next section, it can also lead to larger quantum fluctuations (i.e., slower decay of multistability with N).Therefore, these additional dynamical symmetry lead to a different behavior than the standard power-law decay of multistability and oscillations that are usually studied for Liouvillian gaps closing [121][122][123] .
It is important to note that, even though A 2 κ ¼ 0, A κ+ε A κ ≠ 0 and thus these operators do have overlap with local operators (this follows from P 2 1;N ¼ 1), leading to boundary pseudo-time crystal behavior in local observables as will be shown in the numerical results subsection.
The dynamical symmetries are responsible for both the boundary (pseudo-)crystal and the multistability.As long as the dynamical symmetries we discovered are present these phenomena will be as well.The semi-local dynamical symmetries are completely nonperturbatively stable to all perturbations that are odd in the Majorana fermions i.e., D x ¼ ∑ j x ð1Þ j γ 2jÀ1 þ x ð2Þ j;k γ 2jÀ1 γ 2kÀ1 for arbitrary parameter set x. Namely, lets introduce an arbitrary such perturbation H → H + D x , L k,μ → L k,μ + D k,y .We will still have ½H; A 0 ¼ E k κA 0 and ½L k;μ ; A 0 ¼ ½L y k;μ ; A 0 ¼ 0 and A 0 remains a strong dynamical symmetry.What if the perturbation W is even in the Majorana fermions i.e., H → H + sW, L k,μ → L k,μ + sW?In that case the purely imaginary are stable at least to the second order in the small parameter s according to the results of ref. 106 .Hence, an degree of sub-leading stability remains.
Moreover, according to the results of ref. 106 if we can engineer ϵ to be large a quantum Zeno dynamics will emerge.More specifically, the purely imaginary eigenvalues will be stable up to corrections 1/ϵ 2 that we can engineer to be small.In other words, the dissipation will stabilize the oscillations to any perturbations.This renders the system a dissipative time crystal in the sense of dissipation inducing persistent oscillations in the our many-body system 31 .
Numerical results.In this section, for sake of simplicity, we will focus on the quadratic model, noting that the general conclusion, according to the discussion in the previous sections, holds for the interacting models, as well.
Besides, as the open system with L μ,2 jump operators leads to an exact time-periodic behavior as discussed here we merely focus on L μ,1 dissipators to show the finite-size scaling.
To check the existence of the pure imaginary eigenvalues λ k and find the multiplicity of the null space, we employ the third quantization method for calculating the Liouvillian spectrum for N-fermion chains where mod ðN; 4Þ ¼ 0 subject to the periodic boundary conditions, as described in the model subsection.As the conserved dynamics is quadratic and the jump operators are linear in fermionic basis the Liouville super-operator (L) can be diagonalized in terms of 2N normal master modes acting on the Fock states of density operators 124 .The eigenvalues of the super-operator (λ k ) can be obtained directly from the spectrum of the shape matrix, aka rapidities (β i ) as where v ! is a 2N-long binary string.An alternative approach, also finding a closing spectral gap in related models, can be found in 125 .
The whole Liouvillian spectrum can be exactly calculated by considering all v !within (1, 4 N ).Due to the linear growth of the shape matrix with N, in opposed to an exponential one, one can obtain detailed information about L without being limited to small N chains hence, an equal treatment of the finite-sized systems and larger one approaching the thermodynamic limit (cf. the methods section further details).Figure 1a-c shows the Liouvillian spectrum (λ k ) of noninteracting Kitaev model in ( 1) for μ = 0, w = 1, Δ = i at different chain lengths of N = 4, 12, and 100, respectively.For cases (b) and (c) the spetrum is zoomed in close to the imaginary axis to highlight the slowly decaying eigenvalues.The red dots show the pure imaginary eigenvalues at λ ± = ± 2i, and the green dot corresponds to the degenerate NESS at λ 0 = 0.
To examine the multistability and the long-time behavior of the system we study the two-point correlations and their time evolution.As the whole dynamics, including both the conservative and the dissipative part, is quadratic the state is Gaussian hence its first and second moments (two-point correlation functions) are sufficient to describe the system, fully.Using the Heisenberg picture we can derive the following equations of motion for the two-point correlation functions As can be seen, correlations make a closed set of coupled nonlinear equations that can be numerically solved knowing the initial states.
The time evolution of a local two-point correlation (hĉ 1 ĉ2 i) of such chains is presented in Fig. 2 showcasing the emergence of both non-stationary steady states, aka dissipative time crystal, and the multistability, manifested by two distinct values at long-time limit.Since the long-time solutions in the bistable region depend on the initial state (cf.( 7)), the equations of motion as in ( 20) and ( 21) are evolved for randomized initial states, two of them shown as red and blue lines in each panel.
To examine the scaling of the local observable with the system size (N) in Fig. 3 we plot the long-time value of jhĉ 1 ĉ2 ij as a function of the chain length.
The dots show the results of the numerical calculations when the initial correlations are chosen to be hĉ n ĉm i ¼ 1 þ i and hĉ n ĉy n i ¼ 0, i.e., having one particle on each site.As can be seen the correlation for smaller system sizes follows a power-law behavior decaying as N −1 (red line in Fig. 3) consistent with the lower bound prediction of the semi-local dynamical symmetries and stability of the dynamics subsection.For longer chains, i.e., larger N, more semi-local dynamical symmetries (semi-local finite-frequency Goldstone modes) start emerging and the decay with N slows down, consistent with the results of the emergent dynamical symmetric in the thermodynamic limit subsection.More specifically, with ω k 2 R and the complex decay rate goes down as Oð1=NÞ).This implies both multistability and the persistent oscillations.

Conclusions
In this paper we have shown that for large classes of driven strongly interacting models there exist spectrum generating algebras that are semi-local in the spin basis, which we therefore named semi-local dynamical symmetries.They generically are manifest in the thermodynamic limit only.Physically, they correspond to particle excitations that are invisible to the interaction.They also imply non-local (quadratic) conservation laws, which are promoted to strong symmetries when the system is subjected to pair dephasing.Being quadratic, these operators directly imply memory of the initial condition, i.e., degenerate stationary states and multistability that decays with system size.This means that, unlike previously studied cases of multistability, here the multistability is present in the quantum fluctuations and in the finitesize systems (beyond mean-field) rather than being destroyed by them.Our work implies that genuine (fully exact) multistability for quantum many-body systems requires emergent symmetry structure in the thermodynamic limit, even though a manifest one is not necessarily present for the finite-size system.The system in the thermodynamic limit obtains further emergent dynamical symmetries, which are finite-frequency and finite-momentum quasi-particles dressing the original dynamical symmetry excitations.Therefore, we call these emergent dynamical symmetries semi-local finite-frequency and finitemomentum Goldstone modes.For the dissipative system they imply strong symmetries.These Goldstone modes imply that, as we approach the thermodynamic limit, decay times of oscillations in the local observables diverge, but their amplitude goes to zero at least for initial states that have low-enough entanglement.Hence the system is a boundary time time crystal according to the thermodynamic requirements of ref. 47 .However, the oscillations are clean and periodic both for the isolated (closed L μ = 0) and dissipative system, therefore our system is not a dissipative time pseudo-time crystal in the sense of the dissipative time crystals introduced in ref. 31 , which would imply that dissipation is the one inducing periodic oscillations absent for the isolated system.
To the best of our knowledge, this is the first exact and fully non-perturbative result on the long-debated problem of multistability and multistability in driven-dissipative many-body quantum systems.Although we studied pairing fermionic models, the approach of thermodynamically emergent dynamical symmetries implying quasi-particles that are invisible to certain kinds of interactions is general and can be applied to both bosonic and spin systems.Our work provides an approach for proving presence of multistability in more general and widely studied quantum optical setups.
In future works, we plan to apply our approach to many-body spin and bosonic systems where multistability and persistent oscillations have been experimentally observed in the thermodynamic limit.Exploring the underlying connection between emergent collective behaviors and topology in such systems with quantum synchronization 106,[126][127][128] are other interesting directions of the future studies.

Methods
Shape matrix, rapidities, and Liouvillian spectrum.As described in ref. 124 the two parts of the super-operator, i.e., the conserved dynamics LH and the nonunitary parts LD has the following forms and where ĉi ; ĉy i are the super-operator (a-fermion) annihilation and creation operators in the operator Fock space, respectively.Here, we focus on the K þ , i.e., the even sup-space, only.
We define the 4N × 1-vector of a-fermionic operators as With this definition we can write the Liouville super-operator To write the non-unitary parts easier, we define a matrix with entries M jk ¼ ∑ N μ¼1 l μj l Ã μk hence, M = M † is a Hermitian matrix.Using these definitions we have Considering the antisymmetric properties of H, it becomes apparent that L 11 ¼ ÀL y 22 .Therefore, we have The shape matrix A defined in 124 is simply a rotation of this matrix hence, the eigenvalues of A and L + are the same.If η i , i ∈ {1, 2, ⋯ , 2N} are the eigenvalues of L 22 then the eigenvalues of L + appear in pairs as ðη i ; Àη Ã i Þ.Also from the form of L 22 it is clear that the eigenvalues appear in complex conjugate pairs as ðγ i ; γ Ã i Þ; i 2 f1; 2; Á Á Á ; Ng.Finally, one can conclude that the spectrum of the shape matrix A appear in quadruple of (ξ, ξ * , − ξ, − ξ * ).The rapidities are defined as the subset of the eigenvalues with positive  real parts from which the full spectrum of L þ can be obtained using (19).From this spectrum it becomes evident if there are any kernels, corresponding to multistability, or pure imaginary eigenvalues, corresponding to a non-stationary NESS.
Dispersion of a chain with periodic boundary conditions.Let's consider an Nlong chain with periodic boundary conditions (PBC), i.e., c m+N = c m .One can use the following Fourier transformation to find the spectral form It is straightforward to see that the anti-commutator relations of the Fourier series has the following form fcðkÞ; cðk 0 Þg ¼ 0 ; fcðkÞ; cy ðk 0 Þg ¼ δ mn δðk À k 0 Þ: ð30Þ Replacing each term by its Fourier transform, we get the following spectral form þΔe Àik cðkÞcðÀkÞ þ Δ Ã e Àik cy ðkÞc y ðÀkÞ ð32Þ If the jump operators are identical for all fermionic sites as L j ¼ ffiffi ffi g p c j þ δc y j , then M will read as follows Finally, we can use (27) to determine the rapidities from the eigenvalues of . This leads to the following dispersion relation for rapidies as β(k) where

Fig. 1
Fig. 1 Liouvillian gap closure and the signature of the non-stationary phase.Liouvillian spectrum of noninteracting Kitaev chain for different chain length a N = 4, b N = 12, and c N = 100 subject to periodic boundary conditions.In all cases the chemical potential, hopping amplitude, and pairing potential are μ = 0, w = 1, Δ = i, respectively.The red and green dots show the pure imaginary and zero eigenvalues, respectively.

Fig. 2
Fig. 2 Multistability and non-stationary behavior vs. the system size.Time evolution of the first and the second site correlation (jhĉ 1 ĉ2 ij) showcasing the multistability for two randomized initial states and the oscillatory behavior for various chain length of a N = 4, b N = 12, and c N = 100 subject to the periodic boundary condition.Blue and red lines correspond to two different randomized initial conditions.

Fig. 3
Fig.3Observable scaling with the system size.Local observable scaling vs. the chain length N when the initial correlation in all cases is the same as two-site correlation of hĉ n ĉm i ¼ 1 þ i and initial occupation of hĉ y n ĉn i ¼ 1.The dots are the results of the moments equations of motion integration, the dashed line is guide to the eye, and the red line is N −1 scaling for comparison.
cðkÞ cy ðÀkÞ :ð33ÞThe choice of this spinor is useful since we can readily write the Fourier transform of Majorana fermions as a direct rotation , e-superscripts refer to the odd and even Majorana fermions.Substituting this back into the spectral Hamiltonian we can re-write it in terms of the Majorana fermions as