Mixing times and cutoffs in open quadratic fermionic systems

In classical probability theory, the term {\it cutoff} describes the property of some Markov chains to jump from (close to) their initial configuration to (close to) completely mixed in a very narrow window of time. We investigate how coherent quantum evolution affects the mixing properties in two fermionic quantum models (the"gain/loss"and"topological"models), whose time evolution is governed by a Lindblad equation quadratic in fermionic operators, allowing for a straightforward exact solution. We check that the phenomenon of cutoff extends to the quantum case and examine with some care how the mixing properties depend on the initial state, drawing different regimes of our models with qualitatively different behaviour. In the topological case, we further show how the mixing properties are affected by the presence of a long-lived edge zero mode when taking open boundary conditions.

hypercubes [12], simple one-dimensional exclusion processes [13,14] or the Glauber dynamics of statistical models such as the two-dimensional Ising model in its high-temperature phase [15].
It seems natural at this stage to ask whether an equivalent phenomenon exists in the quantum context. This was in fact already adressed in [16], where an appropriate distance to equilibrium was defined in the quantum information language and the existence of cutoff established in some specific cases. However these results remain tied to some restrictions on the types of systems considered as well as on the nature of the initial states, and leave several open questions, for instance relating to the initial state dependence of the mixing properties.
Rather than generic theorems, our focus in this work will be to study mixing properties in a paradigmatic example of open quantum system, consisting in a free fermionic Hamiltonian linearly coupled to an external bath. More specifically, there are two types of system-bath coupling we shall consider : one corresponds to gain/loss of particles through interaction with the environment, and the other may be considered as a toymodel for Liouvillians with non-trivial topological properties [17]. Both these couplings have in common that they reduce in the classical limit to the master equation for the hypercube random walk, and are therefore good candidates for studying the interplay of quantum coherence and classical cutoffs. Another advantage of such models is that they can be solved exactly using free-fermion techniques [18], and can therefore be used to make analytical predictions on the mixing properties (in [19] free-fermionic chains linearly coupled to a bath were already studied in this perspective, generic bounds on the mixing times were provided, but no discussion of a cutoff phenomenon was made).
Our work is organized as follows. In Section 2, we introduce the models and review their classical limit. In Section 3, we describe the diagonalization of the Liouvillians in terms of a complete set of master modes. We also put forward the topological features of one of our models, in particular the existence of an edge zero mode when open boundary conditions are taken. In Section 4 we turn to the mixing properties, and construct in terms of the master modes a family of factorized initial states, from which the exact time evolution may be computed, and which we argue from numerics that they are good representatives of the "worst" and "best" conditions for mixing among all possible initial states, that is, those which lead at a given time t to the least or most possible mixing. From there, we conclude in Section 5 about the existence of a cutoff in all cases (defined, as customarily, from the worst mixing state at all times), and further describe the dependence of the mixing time on the choice of initial state. An interesting exception, discussed in Section 5.4, is the case of the "topological" model with open boundary conditions, where the existence of the zero mode results in a destruction of the cutoff, a phenomenon we do not know of and analog in classical problems. We also discuss in Section 5.3 the relation with other physical quantities, in particular the von Neumann entropy and local observables. Our findings are summarized in Section 6.

The models 2.1 Fermionic chains with linear dissipation
We consider the evolution of an open quantum system, that is a quantum system coupled to an external environment. In the Markovian description, where the coupling is supposed to be weak enough, and the environment's dynamics fast enough such that the latter does not keep any memory of its interaction with the system, the dynamics of the system's density matrix is known to be well described by an equation of the Lindblad form [6] dρ dt = Lρ := −i[H, ρ] + L D ρ , where the Liouvillian L is formed of a Hamiltonian part accounting for the system's unitary evolution, and a part L D describing the coupling with the environment. The Hamiltonian is taken here to be that of a free fermionic chain, where the c † j , c j are fermionic creation/annihilation operators, satsifying the canonical anticommutation rules Through a Jordan-Wigner transformation [20,21], the Hamiltonian (2) can also be rewritten as that of a XY spin-1/2 chain with transverse magnetic field, namely where the matrices σ x,y,z j act as Pauli matrices on the j th site of the chain, and as identity elsewhere. The case α = 1, in particular, corresponds to the Ising chain in a transverse magnetic field, or, expressed in terms of Majorana fermions , the Kitaev chain [22]. As for the dissipative part L D , we will consider in this work two choices of Lindblad jump operators L µ , both linear in the fermions. As a result the Lindblad equation (1) is quadratic, and can be diagonalized exactly [18]. The first case we will consider, dubbed "gain/loss" in the following, corresponds to two types of jump operators for each site of the chain, corresponding to gain and loss of fermions through interaction with the environment. Such models have been considered in many places in the past literature, see for instance [23,24]. Another case we will consider, dubbed "topological" for reasons that will be explained below (see Section 3.3), corresponds to the following choice of jump operators on each site : Contrarily to the gain/loss case, it is not clear how to realize the above operators in a realistic physical setting. Nevertheless, we will study them as a particularly simple toy model for dissipative fermionic systems exhibiting topological features, such as those considered in [17,25,26] (we note in particular that our model corresponds to a particular choice of the model considered in [26], namely ∆ = 1).

The classical part, and the cutoff phenomenon
A reason for the choices (5), (6) of Lindblad operators is that both lead, in the purely dissipative limit (that is, when the Hamiltonian part is removed from (1)), to a well-known classical Markovian problem. The latter is obtained by restricting to density matrices diagonal in the basis of fermion occupation numbers, whose diagonal entries we label as ρ n1,...n L , n i ∈ {0, 1}. As can easily be checked, both gain/loss and topological models lead to the same master equation which is the master equation for a classical nearest-neighbour random walk on the L−dimensional hypercube {0, 1} L (the latter is also equivalent to the Glauber dynamics of the classical Ising model at infinite temperature [15]). The walk, whose position is indexed by the L-uple n 1 , . . . n L , performs nearest neighbour jumps of rate Lγ, corresponding to each of its components changing value at rate γ. As t → ∞, it reaches the uniform stationary state, where all the ρ n1,...n L are equal to 1/2 L . Starting from an initial configuration ρ(0) (corresponding to a set of non-negative densities ρ n1,...n L (0) normalized to {ni} ρ n1,...n L = 1), a good way to quantify how fast equilibration occurs is through the total variation distance to equilibrium [10,11] One way to think of the total variation distance between two distributions is as the maximum difference between probabilities associated to a single event. The total variation distance to equilibrium has several well-known properties which hold for any Markov chain [11], in particular it is always comprised between 0 and 1, and is non-increasing with time. Given (8), a very important quantity is which is at a given time the maximal distance to equilibrium over all possible initial configurations. Looking back at the particular case of the random walk on {0, 1} L , the maximal distance is obtained at any time by taking ρ(0) to be any purely localized state, where one of the ρ n1,...n L (0) equals 1 and the others equal 0 [12]. As can be seen on Figure 1, and reviewed in more detail in Appendix A, the distance jumps from 1 to 0 around a time t mix (L) = ln L 2γ , and this jump occurs in a time window of width O(1), which becomes much smaller than t mix as L increases: this characterizes what has been coined a cutoff by the mathematical community [7,9,10]. More generally, a sequence of Markov chains indexed by some size or dimensionality L are said to exhibit a cutoff if there exist a sequence of mixing times t mix (L) (typically increasing with L) such that for any > 0, one has [11] 3 Diagonalization of the Liouvillians

The gain/loss model
Lindblad equations of the form (1), quadratic in fermion operators, have been presented and diagonalized in [18]. While generically the diagonalization goes through re-expressing the Liouvillian as a quadratic "Hamiltonian" acting on a superspace of operators and reduces the diagonalization of the latter to that of a 4L-dimensional matrix, in the present case translation invariance makes it natural to reduce the problem to a block-diagonal form without stepping to the super-operator formalism. We therefore introduce the momentum space creation and annihilation operators and it shall be clear depending on the context whether we are using real space or momentum space operators (we will try as much as possible to reserve the letter j for the sites of the chain and k for the momenta). A first step is the diagonalization of the Hamiltonian (2). This is achieved by introducing the Bogoliubovrotated fermions [21], which satisfy the canonical anticommutation relations

and in terms of which the Hamiltonian (2) becomes
Turning to the dissipative part L D , it is easy to check that we can rewrite it for the gain-loss model as Gathering (14) and (15), we can now compute the action of the full Liouvillian L on the rotated fermions η k , η † k . Because of their anticommuting nature, these are not simply annihilated by the components of L with momentum k = k. Therefore, it will turn convenient to introduce the modified fermions so that [η and therefore the modified fermionsη k ,η † k can be used to define a complete set of eigenmodes (master modes) of the Liouvillian. These are indexed by a sequence of 2L binary integers , and read where the bracket notation indicates that if for a given i both ν ± i are = 1, the factor [η kiη † ki ] should be understood as the commutator [η ki ,η † ki ]. The associated eigenvalues of the Liouvillian, namely LC ν = λ ν C ν , read Since all β ± k have a positive real part, the identity (or, rather, ρ ∞ := 1 2 L id), is the only mode with eigenvalue 0 and therefore corresponds to the unique steady state. Relaxation towards the steady state occurs exponentially at late times, with a rate given by the eigenvalue with smallest real part (in absolute value), the so-called spectral gapλ = γ. We accordingly define the relaxation time as the inverse of the spectral gap,

The topological model
We now turn to the topological model (6). The Hamiltonian part is the same as for the gain/loss model, so we look directly at the action of the dissipative part. In terms of the rotated fermions (13), we check which leads to the following action of the full Liouvillian Let us define new fermion operators and similarlyΓ k = Γ k (−1) Q ,Γ † k = (−1) Q Γ † k . Γ k and Γ † k satisfy the same canonical anticommutation relations as the fermions η k , η † k . The eigenmodes of the Liouvillian can be constructed from (25), (26), and these are conveniently reexpressed in terms of (27) as One indeed checks L id = 0 (31) and from there the master modes can be constructed as for the gain/loss model, with eigenvalues of the form (22). Depending on the value of γ/g, α, h, the band structure of the eigenvalues −2β ± k may draw different regimes. This is illustrated on Figure 2, where for simplicity we have specialized to α = 1, corresponding to the Ising chain in a transverse magnetic field. The expression of the spectral gap, indicated by the black arrow on the figure, depends on the regime under consideration. For α = 1, it is given bȳ In contrast to the gain-loss case, we note that the gap closes here for g → 0, as a result of there being many other steady states in this limit.

What is topological about the topological model ?
It is well known that in the regime of parameters |h| < 2, α = 0 the Hamiltonian (2) is in a topologically non-trivial phase, with a gapped bulk and gapless edge modes [22]. This is best seen in terms of the Majorana , in the extreme limit α = 1, h = 0: the Hamiltonian is then a sum of bilinears of the form w 2j w 2j+1 , which in the case of open boundary conditions leaves the modes w 1 and w 2L unpaired. These edge modes commute with the Hamiltonian, anticommute with the fermion number parity (−1) Q , and are therefore responsible for an exact degeneracy of the spectrum between sectors (−1) Q = ±1. Furthermore, these survive throughout the topological phase, where they can be expressed are a power series in h [27]. They are then exponentially located at the edges, and the degeneracy holds exactly in the L → ∞ limit.
We will now see that analogous features hold for the Liouvillian of the topological model, in the same regime of parameters. Recasting the action of the dissipative part L D in terms of spins, we see that it commutes with the operation Taking open boundary conditions for the Hamiltonian (2), σ x L further commutes with H in the limit α = 1, h = 0, so the operation (36) commutes with the action of the full Liouvillian L. Another important property of the operation (36) is that it anticommutes with the "parity of a-fermions" operator [18], which can be defined through its diagonal action on any product C of fermion operators as ZC := (−1) Q C(−1) Q (in other terms Z counts the parity of the number of fermion operators in the product C). Since Z further commutes with the action of the Liouvillian, we see that Ψ plays the role of a zero mode, and results for the Liouvillian in a doubly degenerate spectrum between sectors of parities Z = ±1.
Switching to different values of α and h, we see that these features persist throughout an extended phase, namely whenever the Hamiltonian is in a nontrivial topological phase. More precisely, the degeneracies hold for |h| < 2, α = 0, up to corrections exponentially small in the system size L. This is illustrated on Figure  3. In Section 5.4 we shall come back to these features, which will turn out to have interesting consequences for the mixing properties.

Mixing in the presence of quantum coherent evolution
Having at hand a complete set of master modes for both the gain/loss and topological models, we can now turn to the mixing properties of these models. In particular, we will be concerned with the fate of classical cutoff described in Section 2.2 when quantum coherent evolution is turned on (g = 0), and the dependence of the mixing properties on the choice of initial conditions. The notions of distance to equilibrium, mixing times and cutoffs have been extended to the quantum context in the past literature [16]. For a given quantum channel (Liouvillian) with a unique stationary state, the total variation distance (8) should be replaced by In the former case, the inset shows the first two eigenvalues in both sectors for increasing system sizes.
the trace distance which, in the case where the matrices ρ and ρ ∞ are hermitian (as they will here), is simply half the sum of their absolute eigenvalues. The trace distance shares some important properties of the total varation distance, in particular it is comprised between 0 and 1, and non-increasing with time.
In this section we will examine the time-dependence of the distance (37) as a function of the initial configuration ρ(0), with a particular interest in finding the "worst choice" initial conditions ρ(0), which maximize the distance (37) at a given time t.

Preamble: a look at product chains
Both the gain/loss and topological Liouvillians split into L momentum sectors, each carrying master modes which are eigenmodes of the Liouvillian and which commute or anticommute between different momentum sectors (for the topological model these were explicitly defined in equation (28)- (30), while for the gain/loss model they are simply the Bogoliubov fermions . Let aside anticommutation, our models therefore resemble the so-called product chains, where the time evolution can be decomposed as the tensor product of L independent channels. In the classical setup, product chains (the hypercube random walk being an example) are known to generically exhibit a cutoff under mild assumptions [10,28]. Some results were also established in the quantum context [16]: for instance, if the time evolution can be decomposed as a tensor product of identical quantum channels with a unique nondegenerate steady state, a cutoff was proven to hold when restricting to separable initial states, of the form It is however not known in general, whether such states indeed maximize the distance (37) at all times. Another question raised in [16] is whether the worst choice of initial state at a given time t is always a pure state, or whether, rather counterintuitively, states with some level of mixing might take slower to reach equilibrium (in the classical setup such a phenomenon has been presented in [32], where under the Glauber dynamics extra updates may result in delaying the mixing).
In order to illustrate these ideas, and before embarking into the study of the models defined in Section 2, let us briefly discuss a very simple example of product chain acting on L independent qubits, where the Liouvillian acts on each qubit as : For the sake of illustration, let us restrict to initial states of the product form, namely ρ(0) = ρ 1 (0) ⊗ . . . ⊗ ρ L (0). Each of the ρ k (0) can be decomposed as where 0 ≤ p k ≤ 1, and therefore evolves as The eigenvalues of ρ k (0) are p k , 1 − p k , and therefore p k is a measure of the purity of the initial state : p k = 0 or 1 for pure quantum states, and p k = 1 2 for a completely mixed state. At time t the eigenvalues of ρ k (t) take the form 1 2 ± (p k − 1 2 )f (θ k , t), and the distance to equilibrium is where of course ε k = ± should not be confused with the eigenenergies k of (14). At any time this distance is maximized by maximizing (in absolute value) the terms (p k − 1 2 )f (θ k , t) for each individual k. We observe here a peculiarity of two-level systems, namely that the purity of the initial states appears as a prefactor at all times, and, therefore, that the maximal distance is indeed obtained starting from a pure state 1 . More explicitly, it corresponds to choosing for each k p k = 1, and θ k = π 4 , so the eigenvalues of ρ k (t) are 1 2 ± 1 2 e −γt . The distance to equilibrium as a function of time then has the same form as that of the classical random walk on the hypercube discussed in Section 2.2 and Appendix A, and is indeed seen to exhibit a cutoff.

Constructing initial states
We now move back to our models, which in contrast with product chains cannot be written as a tensor product of independent channels, due to the anticommutation between the operators C ± k in different sectors To get an idea of the difficulties raised by the fermionic nature of the problem, let us look at the gain/loss model, starting with a single momentum sector. A generic initial density matrix can be written in the form (38), (17)- (20) it is easy to read off its time evolution, and eigenvalues. As in the toy-model discussed above, the slowest mixing (namely, the maximal distance) is obtained by choosing p k = 1 (or 0) and θ k = π 4 . However, putting all momentum sectors back together, it is not anymore a valid choice to simply consider a product of such density matrices: contrarily to the case of product chains these do not commute with one another and therefore their product, being non-hermitian, does not correspond to a physical initial configuration. Similar remarks can be paralleled for the topological model.
In the following we will construct several classes of physical initial density matrices in terms of the master modes C + k , C − k , C ± k , which will turn out to be relevant for both the gain/loss and topological models (we recall that for the former the correspondence is Single-sector commuting density matrices A natural way to work around the problem of anticommutation is to start from products of single sector commuting density matrices, namely linear combinations of id and C +− k in each sector : with 0 ≤ p k ≤ 1. For both the gain/loss and topological models, these evolve in time as and we write the corresponding single-sector eigenvalues Paired commuting density matrices Another possibility is to combine the anticommuting master modes of different sectors into commuting objects. There are many ways to do this, but in the following we will restrict to the so-called paired density matrices built from two sectors with momenta k 1 , k 2 . We therefore define whereΓ k ,Γ † k were defined in (27) for the topological model, while for the gain/loss model they should just be taken to beη k ,η † k . Eq. (44) is the most general choice of a commuting combination corresponding to a pure state, that is with eigenvalues (1, 0, 0, 0) (times the identity in other sectors). We will not justify here that starting from a pure state in paired sectors is indeed what maximizes the distance to equilibrium, but numerical studies below will confirm this fact. The explicit time dependence of ρ C k1,k2 (t) will be worked out separately for the gain/loss and the topological models in the next sections, and we will generically write the corresponding eigenvalues as Density matrices for the full system The single-sector and paired commuting density matrices can now be combined into density matrices for the full system, namely arbitrary products of them can be taken. Furthermore, we may still multiply such products by one noncommuting single-sector density matrix ρ k (t). We therefore define where the parentheses around k in the left-hand side, and around ρ k (t) in the right-hand side, mean that the non-commuting density matrix ρ k (t) may or may not be present in the product. We must then have 2m + n(+1) = L. In each sector (or pairs of sectors) the initial density matrices have internal parameters (p k , ϕ k , ϕ k1,k2 , etc...), which we leave unspecified for the moment. Writing the eigenvalues of ρ k (t) as p (1,2) k (t), and those of the single-sector/paired commuting matrices as (43) and (45) respectively, we obtain the eigenvalues of (46) as products of the latter, which results in the following expression for the trace-norm distance to equilibrium where once again parentheses account for the presence or absence of the non-commuting matrix ρ k (t). Let us emphasize that the the density matrices (46) are only a very small subset of all the possible initial states. In particular, these are completely (or almost completely, because of the pairing between sectors) factorized. In the following, we will however observe from numerics that such density matrices still encompass the worst-choice initial state at any time.

The gain/loss model
We are now ready to examine the mixing properties of our models, starting with the gain/loss model. Before turning to numerics, we want to find the conditions which maximize the distance at any time t for density matrices of the form (46). As can easily be checked, the distance (47) is generically maximized by separately optimizing the individual eigenvalues in each sector, that is taking the p (c) as close as possible to 0 or 1. Let us therefore see how this goes sector by sector.
Single-sector commuting density matrices The time evolution of single-sector commuting matrices has been discussed in the previous section, and the associated eigenvalues found to be given by (43). Their contribution to the distance is maximized by taking p k = 1 (or 0), that is, starting from a pure state in the corresponding sectors. Therefore, Paired commuting density matrices Paired commuting density matrices take the form (44), where in the present caseΓ † k ,Γ k =η † k ,η k . Plugging in the time evolution of the latter, we check that the corresponding eigenvalues take the form irrespectively of the value of ϕ k1,k2 , as well as of the value of k 1 and k 2 .
Single-sector non-commuting density matrices The case of single-sector non-commuting density matrices was already briefly discussed in the beginning of Section 4.2. These evolve as and the maximal contribution to distance is obtained by setting p k = 1 (or 0) and θ k = π 4 , irrespectively of the value of ϕ k . The associated eigenvalues then read : Gathering (48), (49), (51), it is easy to see that the density matrix of the form (46) which maximizes the distance to equilibrium at all times corresponds to a product of one single-sector non-commuting density matrix, and commuting density matrices in all other sectors (given the similarity between (48) and (49), it does not matter which of those are paired and which are single-sector). On Figure 4, we represent the associated distance as a function of time for a system of finite size L = 4 (blue curve), and plot in comparison the distances computed from exact diagonalization starting from randomly drawn initial conditions, corresponding to either pure or mixed states. We also represent as a red curve the distance for a product of commuting (paired or single-sector) density matrices, which, following the same lines as above, corresponds to the fastest mixing (that is, minimizes the distance to equilibrium at time t) when restricting to pure quantum initial state of the form (46). Several observations can be made from there. First, it is apparent that at any time the least mixed state is indeed of the form (46), with one single-sector fermionic density matrix, and commuting matrices in other sectors. This observation continues to hold for larger systems sizes and we conjecture it to be true for generic L, however we limit ourselves to presenting results for a small system here, as for larger systems the Hilbert space of possible initial states becomes longer to explore, and the upper bound much harder to saturate from randomly drawn samples. In Section 5, we will study the large L behaviour of the corresponding distance d(t), with particular interest in whether a cutoff develops in this limit. Another observation is that, conversely, the lower bound for the distance (once restricted to starting from pure quantum states) is not given by a density matrix of the form (46) (red curve), which means that the states with fastest mixing may be non-factorizable. The dashed light gray lines correspond to randomly drawn initial density matrices with various levels of mixing, while the darker gray lines correspond to initial pure quantum states. The colored lines are analytical predictions for initial density matrices of the type (46). We recall the notation k j = 2πj L for the momenta, but emphasize that in the present case the results are insensitive to any permutation of the latter.

The topological model
Let us now follow the exact same steps for the topological model.
Single-sector commuting density matrices As has been discussed in Section 4.2, the time evolution of the single-sector commuting density matrices has the same form as for the gain/loss model. Once again, their maximal contribution to the distance is obtained by choosing p k = 1 (or 0) in (41), with eigenvalues given by (48).
Paired commuting density matrices Paired commuting density matrices take the form (44). Their time evolution is read off by rewriting theΓ k ,Γ † k in terms of the master modes using (28)- (30), and we find that as a function of time the corresponding eigenvalues take the form where the explicit form of f ± k1,k2 (t) depends on ϕ k1,k2 . The corresponding contribution to the distance is maximized at a given time t for a certain value ϕ * k1,k2 (t) of ϕ k1,k2 , and the associated f * ± k1,k2 (t) will be given in Section 5 for the particular case of the Ising chain, α = 1, h = 0.
Single-sector non-commuting density matrices The generic form for a single-sector (non necessarily commuting) initial density matrix is The dashed light gray lines correspond to random initial density matrices with various levels of mixing, while the darker gray lines correspond to initial pure quantum states. The colored lines are analytical predictions for initial density matrices of the type (46), where we recall the notation k j = 2πj L .
As above we can compute the time evolution ρ k (t) by recastingΓ k ,Γ † k in terms of the master modes, and find the associated eigenvalues under the form where once again the function f k (t) depends on the parameters θ k , p k , ϕ k . The maximal distance at time t is obtained by maximizing |f k (t)|, which corresponds to a choice of parameters The corresponding explicit expression of f k (t), which we denote f * k (t), is too lengthy to be detailed here, but will be given in Section 5 for the particular case α = 1.
As for the gain/loss model we now compare the distances built from (48), (52), (54) with numerical results from randomly drawn initial density matrices, see Figure 5. At difference with the gain/loss model, the eigenvalues (48), (52), (54) now depend on the momentum sectors, and differ between the paired or single sector commuting cases. The distance to equilibrium therefore depends on the repartition of commuting matrices between paired and single-sector, as well as on the associated distribution of momenta, the choice of pairings between those, etc... We do not attempt at a generic discussion here (in Section 5 this discussion will be simplified by restricting to the particular case α = 1, h = 0, where the dependence in momentum vanishes), but restrict to plotting for L = 4 and different sets of parameters the distances associated with some relevant initial states, including the slowest and fastest mixing cases (blue and red curves, respectively). As can be observed from the blue curves in Figure 5, the maximal distance is attained at any time by a product of one non-commuting single sector density matrix, and commuting matrices in all other sectors, with the maximal of these being paired (which leaves out, for L even, one commuting single sector). We also note, as attested by the various crossings between blue curves, that the repartition of momenta for which this maximal distance is attained may vary over time.
Turning to the fastest mixing states (restricted to start from a pure quantum state), two regimes emerge, which seem to coincide with the regimes γ/g > 2 − |h|, γ/g < 2 − |h| emerging from the band structure on Figure 2. In the "small dissipation" regime (γ/g < 2 − |h|, left panel on Figure 5) things seem to go similarly as in the gain-loss model, namely the fastest mixing state is not of the form (46) and therefore corresponds to a non-factorizable state. In contrast, in the "large dissipation" regime (γ/g > 2 − |h|, right panel on Figure 5), fastest mixing does seem to be attained for a product of single sector commuting density matrices, and we further observe a separation of timescales between the fastest and slowest mixing states. This separation of timescales will be studied in more detail in the next section.

Study of the mixing times and cutoffs
In the previous section we have constructed a family of factorized density matrices (46), which, despite being a very restricted subset among all the possible initial states, do achieve at all times the slowest mixing (that is, they maximize at all times the distance (37) to equilibrium), as well as, in some regimes of parameters, the fastest mixing. We will now exploit the analytical formula (47) for the associated distance to equilibrium in order to investigate the mixing properties and existence of cutoffs in the different regimes of our models. Then, we will turn in Sections 5.3 and 5.4 to other related aspects, namely the behaviour of other physical observables with respect to mixing, and the effect of the edge mode discussed in Section 3.3 when taking open boundary conditions in the topological model.

The gain-loss model
As we have seen in Section 4.3, the slowest mixing in the gain-loss model is obtained at all times starting from a density matrix of the form (46), with one non-commuting single sector density matrix and commuting matrices in all other sectors. Using (47), the corresponding distance to equilibrium as a function of time reads (where once again the signs ε k = ± have nothing to do with the energies k ), and can be reinterpreted as the distance to equilibrium for a classical random walk on an anisotropic hypercube, with jump rate γ in one direction and 2γ in the other directions. The large L asymptotics of (56) can be tackled analogously as in the isotropic case, and we show in Appendix B that for L large, t large, where erf is the Gauss error function. The function (57) develops a cutoff for where the relaxation time t rel , defined in (23), is the inverse of the spectral gap. The asymptotic profile of the cutoff around this time is Several remarks are in order : first, the profile (59) is formally the same as for the classical hypercube random walk, and in particular retains a finite width when L → ∞. Second, it is independent of the coherent coupling strength g. Nevertheless, a striking difference with the classical case is that all timescales, and in particular the mixing time (58), are divided by two as a consequence of the constraint imposed on initial states by the fermionic nature of the problem, and which can be viewed as a kind of exclusion constraint (see section 4.2).

The topological model
For the topological model, we have seen in Section 4.4 that the least mixed state is obtained at any time from a factorized matrix of the form (46) with one non-commuting factor and paired commuting density matrices in other sectors, plus one residual single sector commuting matrix in the case where L is even. Furthermore, the way momenta {k 0 , . . . k L−1 } should be distributed between all these factors in order to maximize the distance (47) may generically depend on time in a complicated fashion as attested by the multiple crossings between the blue curves on Figure 5.
Here, we will specify to the particular case α = 1, h = 0, that is that of the Ising chain in the absence of an external magnetic field, which brings the simplification that the energies k in (14) and hence the Liouvillian eigenvalues β ± k do not depend on k, namely k = 2g for all k. In practice this means that all the blue curves on each panel of Fig. 5 collapse into a single one, for which we will be able to compute the large L asymptotics. The distance reads from (47) where the functions f * k (t) and f * ± k ,k (t) are those appearing in (54) and (52), taken for the worst-choice values of the initial state parameters, and which we here find to be As for the gain-loss model, we can reinterpret (47) as the distance for a classical random walk on an anisotropic hypercube, with the additional difference that the jumps are not exactly Poisson processes anymore, as the functions (61), (62) are not purely exponentially decaying, but only become so in the late time limit.
In Appendix B we show from there that in the large L limit all of these distances exhibit a cutoff phenomenon at times ∝ ln L, and that the asymptotic expressions of the cutoff profiles can be expressed in terms of the Gauss error function erf. In order to describe these profiles in more detail we now distinguish between the small disipation (γ < 2g) and large dissipation regimes (γ > 2g), recalling that the spectral gap (35), becomes for α = 1, h = 0,λ We will also be interested in the asymptotics of the "fast mixing" red curves of Figure 5, which is given in both regimes by This is exactly the distance to equilibrium for a classical random walk on the isotropic hypercube, and therefore develops a cutoff at times As has been discussed in section 4.4 these red curves do not necessarily correspond to the fastest mixing state (only in the large dissipation regime might it be the case), but looking at Fig. 5 it seems reasonable to expect that in both regimes these give a good indication of the spreading of mixing times for all possible initial states (conditioned to be pure quantum states at t = 0).  Figure 6: Distance to equilibrium in the topological model for the Ising chain in zero magnetic field (α = 1, h = 0) starting from two different classes of initial states, in the regimes γ > 2g, γ < 2g and for different sizes L. The blue curves correspond to the least mixed state at all time ("worst choice" initial conditions), while the red curves correspond to a product of single-sector commuting matrices, which for γ > 2g seems to be the fastest-mixing state.

Small dissipation regime (γ < 2g)
In the small dissipation regime, we find following Appendix B the asymptotics which develops a cutoff for The corresponding distance, as well as the "fast mixing" distances are plotted on the left panel of Figure 6 for various sizes. The situation here is quite similar to that of the gain/loss model : the slowest and fast mixing curves develop a cutoff at the same value mixing time (67), which is half that of the classical problem. Assuming, as seemed reasonable from Figure 5, that the distance for generic (pure) initial states remains "close enough" to these two cases, we may conclude that in this regime all initial states develop a cutoff at time (67).

Large dissipation regime (γ > 2g)
In the large dissipation regime, the asymptotics of the functions (61), (62) further simplifies to Following Appendix B, we check that the distance (60) develops a cutoff at and the cutoff profile is given by An interesting feature appearing here, illustrated on the right-hand panel of Figure 6, is the separation of timescales between the slowest mixing states and the "fast mixing" states (which, as concluded from Figure 5, are likely to be the fastest mixing states in this regime). In other terms, the time required to reach equilibrium strongly depends on the initial state in this regime. Let us point that the fact that both the fastest and slowest mixing states show a cutoff phenomenon does not imply a cutoff phenomenon for any initial state. Rather, we conclude from the above that any initial state (conditioned to be a pure quantum state at t = 0) should display a weaker version known as pre-cutoff [10,16] , characterized by the two sets of timescales t mix (L) (eq. (65)) and t (fast) mix (L) (eq. (70)) and the fact that for any > 0,

Other physical observables
In the previous sections we have observed that both the gain/loss and topological models exhibit a cutoff phenomenon, as defined by the trace-norm distance to equilibrium of the slowest-mixing state. More generally, it seems that the distance from an arbitrary initial state satisfies a pre-cutoff with two timescales t (fast) mix (L) and t mix (L). In the gain/loss model and the weak-dissipation regime of the topological model these two timescales coincide, so there is in fact a cutoff from any initial state.
A natural question at this stage, is whether the notion of cutoff extends to other physical observables, for instance to the growth of the von Neumann entropy S(t) = −Tr(ρ(t) ln ρ(t)) or to the equilibration of local observables. Let us start with the entropy. Starting from initial states of the form (46), it is easy to see that it can be decomposed at all times as a sum over individual (pairs of) sector contributions, which, for both the slowest-and fast-mixing states of the gain/loss model, converges for L → ∞ to the asymptotic expression We plot on Figure 7 the trace norm distance and the (rescaled) entropy for the fastest-mixing state as a function of time for various system sizes. While the former develops a cutoff at times t mix (L) ∝ ln L, the latter relaxes exponentially with a timescale 1/2γ for any system size, and therefore is insensitive to the cutoff. A similar conclusion can be made for local observables. This is in fact very easily seen in the classical case: starting from the configuration localized on site {0} L of the hypercube, a meaningful local observable is for instance the expectation value of the i th coordinate. As can easily be checked this evolves in time as 1 2 − e −γt 2 , so relaxes exponentially towards equilibrium with a timescale incommensurate with the mixing time. Now, there is little reason to expect that things should go differently in the quantum case. In order to illustrate this fact but keep the calculations at their simplest, let us consider the example of the gain/loss model with α = 0, and look ath the evolution of the number of particles at site j = L. In this case the Bogoliubov fermions η † k , η k coincide with the original momentum-space operators c † k , c k , and we have simply Starting for instance from the situation where there is one particle on site j = L and none on the others (as we have argued earlier, this initial state, like all pure quantum initial states, is expected to develop a cutoff as L → ∞), we have c † k c k (0) = 1/L ∀k, k , which we can plug into (76) in order to obtain the various occupation numbers at time t. The results are plotted on the right-hand panel of Fig. 7 for j = L and j = L/2. In all cases, the occupation numbers converge to their equilibrium values with a timescale 1/2γ, which is the same as for the entropy and is incommensurate with the mixing time.

Effect of the edge mode
We finally move on to considering the effect of topological features on the mixing properties. For this sake we take the topological model with open boundary conditions, as studied in Section 3.3, where it was shown that there exists throughout the regime |h| < 2, α = 0 a zero mode Ψ, commuting in the L → ∞ limit with the action of the Liouvillian, and resulting in a twofold degeneracy of the Liouvillian spectrum between sectors Z = ±1. While these features hold throughout the topological phase up to corrections exponentially small in L, they are exact at any size at the zero-field Ising point α = 1, h = 0, where the action of Ψ is simply given by (36). As a degeneracy of the spectrum means in particular closing of the spectral gap, we expect that taking open boundary conditions will have a drastic effect on the mixing properties. In fact, the whole discussion on mixing and cutoffs, including the definition (9) of the distance, is valid only in the case where the Liouvillian has one single non-degenerate steady state, which in the present case discards the particular point α = 1, h = 0. For that case, the various initial states decay at large time towards linear combinations of the form 1 2 L (id + f σ x L ), −1 ≤ f ≤ 1, and the distance as defined in (9) may take any value between 0 and 1/2 in the t → ∞ limit. Now, in the rest of the topological phase, the spectral gap is strictly speaking > 0 for any system size L, and any initial state will eventually decay towards the unique steady state ρ ∞ = 1 2 L id. The associated timescale is however now exponentially large in the system size, and cancels the cutoff effect observed for periodic boundary conditions. This is illustrated on Figure 8, where we compare the distance from randomly drawn initial states for open and periodic boundary conditions (we also plot in comparison the distances for the gain/loss model with open boundary conditions, which as should be expected does not show any drastic modification compared with the periodic case).
The long-lived edge mode also manifests itself at the level of "local" (in terms of spins, not of fermions) physical observables. The consequences of edge zero modes on the out-of equilibrium physics have been thoroughly studied in the past litterature, for closed Hamiltonian dynamics [29] but also in a dissipative (albeit quite different from the one considered here) setup [30,31]. Here we illustrate these effects by comparing the time evolution of the expectation value σ x L for open and periodic boundary conditions, see Figure 9. While in the periodic case σ x L = σ x 1 quickly relaxes to its equilibrium value, in the open case relaxation takes a much longer time, exponentially increasing with the system size. This is a result of σ x L being "almost conserved" by the Liouvillian evolution, and we indeed check in comparison a much quicker relaxation at the other end of the chain, for the expectation value σ x 1 .

Conclusions
In this work we have examined the mixing properties of open quantum fermionic systems whose time evolution is governed by the Lindblad equation, and in particular the quantum counterpart of the so-called cutoff phenomenon, a well-studied aspect of classical Markov chains. A general conclusion of our analysis is that the notion of cutoff generally carries over to the quantum framework (see also [16]): considering two freefermionic models which in the classical limit reduce to the well-studied hypercube random walk [10,12], we showed that the "worst-choice" distance to equilibrium, defined by maximizing over all possible initial conditions, develops a cutoff at a time t mix (L) ∝ ln L as L → ∞. Going further, we explored the initial-state dependence of the mixing properties, and provided some evidence that a cutoff, or pre-cutoff (depending on the model and regime) should exist for arbitrary initial states. With respect to the classical case, the quantum realm however reveals a number of new and potentially interesting features. A first aspect, well illustrated by the gain/loss model, is the halving of the mixing times with respect to the classical problem, as a result of some "exclusion constraints" due to the fermionic nature of our models. This is in contrast, in particular, with the case of product quantum channels considered in [16], where, under some assumptions, the mixing time was related to the relaxation time through t mix (L) = t rel 2 ln L (here, conversely, t mix (L) = t rel 4 ln L). Another aspect is the separation of timescales in some regimes between the slowest-mixing and fastest-mixing states. Finally, we have seen in the so-called "topological model" how the presence of an edge zero mode dramatically alters the mixing properties when open boundary conditions are taken, a phenomenon we know of no counterpart in the classical context (let us warn however, that the classical limit of open boundary conditions in the quantum chain does not correspond to open boundary conditions on the classical hypercube).
There are several interesting directions to go from there. The most natural next step would be to investigate the mixing properties of systems with a non-trivial steady state, for instance one with longrange order, entanglement [1], or current-carrying [3,33,34]. Another possible direction would be to study whether the present analysis, and in particular the existence of cutoffs, could be extended to non-markovian dynamics such as that governing the evolution of subsystem density matrices in isolated quantum systems after a quantum quench [35,36], or in random quantum circuits [37] (see also [38]). A second direction concerns the relation with other physical observables : even though we have observed in Section 5.3 that the most natural local observables as well as the von Neumann entropy are insensitive to the presence of a cutoff, namely they do not develop a sharp jump as the trace-norm distance to equilibrium does at the mixing times t mix (L), it remains an intriguing question whether the cutoff phenomenon might transpire into other physical quantities. On a more technical level, it would of course be interesting to move on to interacting systems, for instance using the mapping of some interacting Liouvillians onto Bethe ansatz solvable quantum Hamiltonians [39][40][41]. The study of mixing in this framework however seems to be a challenging issue. In fact, even in classical Markov chains which can be mapped onto quantum integrable systems (for instance the symmetric or assymetric exclusion processes [42][43][44]), relating mixing properties to the Bethe ansatz language seems to be a difficult task.

4
. Similarly, under the distribution π(x), the distribution of |x| is a normal distribution of mean L 2 and of variance L 4 . We may therefore write where Φ µ,σ is the cumulative distribution function of the normal law with parameters (µ, σ), that is, We check from there that the distance d(t) jumps from 1 to 0 at a time with increases logarithmically with L (see Figure 1): more precisely, expanding to second order in e −γt , we obtain d(t) = erf e −(γt− 1 2 ln L)

B Asymptotic expressions for the various distances
Our goal here is to use techniques similar to those presented in Appendix A to derive the asymptotic behaviour of the distance d(t) for "anisotropic" random walks as appearing in the main text. We consider a random walk on the hypercube {0, 1} L , where now the different components change value at different rates. In fact, these are not necessarily Poisson processes, and we define the probabilities P(x i = 1) = 1 2 (1 − f i (t)), where f i (t) are monotonously decreasing from 1 to 0 as t grows from 0 to ∞, so that the stationary distribution π is the uniform one. In practice we will consider the case of interest in Section 5.2, which is the case where L is even, and Accordingly, we decompose |x| = x 1 + x 2 + |x + | + |x − |. The probability of a given configuration at time t reads Proceeding as in Section A, computing the distance d(t) implies a restricted summation over configurations with P(x) ≥ 1 2 L , which corresponds to where we have defined By virtue of the central limit theorem, when L is large, y is described by a normal law of parameters Similarly, under the uniform distribution π(x), y is distributed under normal distribution of parameters (92) Decomposing the distance as we therefore obtain d(t) = x1,x2=0,1 p(x 1 , x 2 )Φ µ,σ (y max (x 1 , x 2 )) − 1 4 Φ µ0,σ0 (y max (x 1 , x 2 )) , where Φ µ,σ was defined in (81).