The impact of different unravelings in a monitored system of free fermions

We consider a free-fermion chain undergoing dephasing, described by two different random-measurement protocols (unravelings): a quantum-state-diffusion and a quantum-jump one. Both protocols keep the state in a Slater-determinant form, allowing to address quite large system sizes. We find a bifurcation in the distribution of the measurement operators along the quantum trajectories, that's to say, there is a point where the shape of this distribution changes from unimodal to bimodal. The value of the measurement strength where this phenomenon occurs is similar for the two unravelings, but the distributions and the transition have different properties reflecting the symmetries of the two measurement protocols. We also consider the scaling with the system size of the inverse participation ratio of the Slater-determinant components and find a power-law scaling that marks a multifractal behaviour, in both unravelings and for any nonvanishing measurement strength.


Introduction
Entanglement [1,2] plays an important role in the unitary dynamics of many-body quantum systems in many different contexts: It marks quantum phase transitions [3]; It behaves differently in the dynamics of thermalizing, integrable and many-body localized quantum systems [4][5][6]; It even allows to detect the existence of topological boundary modes [7][8][9][10].Recently, attention has been moved also to the behaviour of entanglement in situations beyond the unitary dynamics, where the evolution of monitored systems is considered.The interplay between the intrinsic dynamics of the system and that induced by the quantum measurement process can lead to a variety of scaling regimes for the asymptotic entanglement entropy, giving rise to the so called entanglement transitions.
In this paper we consider a free-fermion chain undergoing a dephasing Lindbladian, and describe it as an average over random quantum trajectories in two different ways (unravelings).Specifically, we employ either a quantumstate-diffusion (QSD) unraveling [31] or a quantum-jump (QJ) unraveling.Both unravelings are chosen in such a way that the state is kept in a Gaussian form (a Slater determinant), thus a numerical analysis of systems with up to L = O(10 2 ) sites is easily affordable.
For this model, the behaviour of the asymptotic entanglement entropy averaged over the quantum trajectories has been carefully scrutinized and is still under debate.After the pioneering work in Ref. [31], where a saturation of the entanglement entropy with the system size (arealaw) was predicted, Refs.[39,73] claimed the presence of a transition from a logarithmic increase to an area-law behaviour.More recently, Ref. [36] challenged this result, suggesting that there is only area law and the transition is actually a crossover, due to the exponential growth with the inverse measurement strength of the size where the entanglement entropy saturates.
Here we focus instead on the distribution of the expectations of the measured operators along the quantum trajectories.Properties of similar distributions have already been studied in [45,74,75].In general, the properties of this distribution are not related to the entanglement properties, although with an exception [74].Comparing the behaviour of this distribution in the two unravelings, we observe differences as well as similarities.In the QSD unraveling this distribution is symmetric around n = 1/2, while in the QJ one this symmetry is absent.This reflects the fact that, while QSD is invariant under particle-hole symmetry, QJ is not.In both unravelings the distribution displays a bifurcation, moving from unimodal to bimodal shape.One can see this phenomenon already at finite size, being related to the frequent local measurements forcing the state to be locally similar to an eigenstate of the measurement operators, in analogy with the quantum Zeno effect [76][77][78][79].In the case of QSD, where the bifurcation occurs for a measurement strength γ ⋆ QSD ≈ 0.2, the two maxima in the bimodal phase are symmetric and stem continuously at the transition point, in a way formally reminiscent of the Landau picture of second-order phase transitions.At the bifurcation point, the single maximum bifurcates into two maxima, with a discontinuity in the derivative analogous to the pitchfork bifurcation in classical dynamical systems [80][81][82].In the QJ case, the bifurcation occurs at γ ⋆ QSD ≈ 0.23.In contrast with the QSD case, here the maxima are not symmetric, and the smaller one appears discontinuously at the bifurcation point, in a way formally reminiscent of the mean-field analysis of first-order phase transitions.
The bifurcation originates from the interplay between the unitary dynamics and the measurement operations.In particular, when the measurement strength is large, the second prevail and the state is expected to be similar to a product state.Consistently with that, in both cases, for large measurement strength, the distributions show strong maxima near n = 0 and n = 1, marking that the state is near to a separable one.This also agrees with what is known for the asymptotic entanglement entropy [31,36,39,73], that for QSD tends to 0 for large measurement strengths.
Finally we study the localization of the state, to probe whether the bifurcation corresponds to a delocalizationlocalization transition or not 1 .In order to do that, we first note that, for both unravelings, the state is a Slater determinant.We consider the components of this determinant and evaluate the inverse participation ratio (IPR) of these components in the space basis.This is a standard probe for localization, used for instance in studies of Anderson localization [85].We find no sharp transition, but we see that the IPR scales with the system size L as a power law ∼ L −α with 0 < α < 1.For both unravelings, α depends smoothly on γ similarly to what happens in the two-dimensional case [47].This is a mark of multifractal behaviour of the components of the Slater determinant [86,87]: The system is anomalously delocalized and never achieves perfect localization, as confirmed by results on the conductivity [83].Such multifractal behaviour also occurs at the transition between extended and Anderson-localized phases, so we can say that our model is always in a critical Anderson regime.
The paper is organized as follows.In Sec. 2 we present the model Hamiltonian, its diagonalization and its symmetries.In Sec. 3 we describe the two unraveling protocols that, on average, give rise to the same Lindblad equation.In Sec. 4 we show that for both unravelings the state along a trajectory is always in a Slater determinant form.In Sec. 5 we compare the two unravelings and show that, although the asymptotic entanglement entropy behaves similarly, the distributions of the expectations of the measurement operators have different symmetries, but show both a bifurcation (Sec.5.1); we also discuss the scaling of the IPR averaged over disorder and its multifractal behaviour (Sec.5.2).Finally, in Sec.6 we draw our conclusions.

Model Hamiltonian
We consider a model of spinless free fermions on a one-dimensional lattice with L sites, described by the Hamiltonian where ĉ( †) j are anticommuting fermionic operators acting on the j-th site2 .After defining the vector Ψ ≡ (ĉ 1 , . . ., ĉL ) T , the Hamiltonian (1) can be written as where we introduced the L × L matrix H such that A generic fermionic quadratic model as the one in Eq. ( 2) can be cast in a more transparent form by defining the fermionic operators such that Φ ≡ γ1 , . . ., γL T = U Ψ, (5) where U denotes the unitary transformation which diagonalizes this matrix H: Thus, we finally get: which represents a model of free γk -fermionic particles with dispersion relation ω k ≥ 0. The eigenstates of Ĥ can be put in a simple Gaussian form (Slater determinants) and are fully characterized by the two-point correlation matrix of ĉ-fermions, ).This makes such kind of quadratic fermionic systems particularly suited to numerical computations up to large sizes L ∼ O( 103 ).It is also worth noticing that Gaussianity of states is preserved by the application of any operator that can be written as the exponential of an other operator that is quadratic in the fermionic operators (see, e.g., Refs.[31,[88][89][90]).
We finally remark that the Hamiltonian in Eq. ( 1) conserves the total fermion number N = j ĉ † j ĉj and is invariant under the particle-hole transformation, the one that exchanges every ĉ † j with the corresponding ĉj , and thus exchanges every nj with 1 − nj .

Monitored dynamics
To characterize the effects of measurements on the unitary quantum dynamics generated by the Hamiltonian in Eq. ( 1), we construct a suitable stochastic Schrödinger equation, whose exact form depends on the kind of measurement protocol (unraveling) one wants to study [91,92].The measurements provide noise and on each stochastic realization of the noise the state evolves along a so-called quantum trajectory.In what follows, we focus on two distinct unravelings that give rise to the same master equation for the mean density matrix of the system, once averaged over the stochastic realizations.

Quantum state diffusion
In the first measurement protocol, we assume the number ⟨n j ⟩ t of fermions on the j-th site (with nj = ĉ † j ĉj ) to be continuously measured with a rate γ.This generates a QSD dynamics, that is, a collection of Wiener processes.For small time evolution steps δt, it can be cast in the following form [31]: neglecting normalization constants and assuming ⟨•⟩ t ≡ ⟨ψ(t)| • |ψ(t)⟩.The δW i t variables are normally distributed with zero mean and variance γ δt.
We emphasize that this unraveling is invariant under the particle-hole transformation.Indeed, by exchanging nj ↔ 1 − nj in Eq. ( 8), one gets back the same equation with an immaterial multiplying factor in front of the state and the sign of δW j t changed.This different sign is immaterial, since the δW j t are Gaussian variables with zero mean, and all the moments of any distribution over the randomness are left unchanged by the transformation.

Quantum jumps
Secondly, we consider an unraveling based on occasional yet abrupt measurements of the local operators mℓ = 1 + nℓ , where 1 denotes the identity operator.We assume that these measurements can occur on any site ℓ with a probability p ℓ ∈ [0, 1] that we define below.This kind of measurements are modeled by a QJ dynamics that describes the evolution of the system for discrete time steps δt.At each step one of the following operations on the state can randomly occur: , where γ denotes the measurement rate, we project the system state with one of the L measurement operators 3 : 2. With probability p = 1 − j p j we evolve with a non-Hermitian effective Hamiltonian Ĥeff where Ĥeff = Ĥ − 3 2 iγ j nj .Since for spinless fermions 1 + nj = e nj ln 2 (at most one fermion per site is allowed), it is immediate to see that this kind of dynamics can be implemented by preserving Gaussianity [43,44].Moreover, as one can easily check, it is not invariant under particle-hole transformation.

Average over trajectories
When averaged over the stochastic realizations, in the limit δt → 0, both QSD and QJ unravelings lead to the same Lindblad master equation for the mean density matrix ρ t = |ψ(t)⟩ ⟨ψ(t)|, where the overline indicates the ensemble average, [93] In this equation, mj denote the (Hermitian) measurement operators and γ quantifies the strength of the coupling between the system and the measurement apparatus.Note that in the QSD protocol (Sec, 3.1) we set mj ≡ nj , while in the QJ protocol (Sec.3.2) we set mj = 1 + nj .In fact, substituting mj → 1 + nj , one recovers the same Lindblad equation for the measurement operator nj .
In other words, the two stochastic dynamical protocols above are unravelings of the same Lindblad master equation (11) [43,44].As a consequence, any linear function of the density matrix ρ t must be independent of the chosen unraveling.We mention, for example, the expectation value of a generic observable Ô, Finally we point out that, independently of γ, the dynamics described above preserves the total number of fermions N = j nj .As a consequence, one can work in a subspace of the full Hilbert space with a fixed filling.In what follows we work in the half-filling subspace and fix

Methods
To reduce the numerical complexity of the problem, we exploit the conservation of the number N of fermions, as stated above.We choose as initial condition the Néel state with N = L/2 particles: where |Ω⟩ is the vacuum of ĉ-fermions, such that ĉj |Ω⟩ = 0, ∀j = 1, . . ., L. Then, we let it evolve according to one of the two measurement protocols described before, up to time t.The time evolved state |ψ(t)⟩ can be cast in the form of a generic Gaussian state.The full information of such state is contained in a L × N matrix U t , defined by such that U † t U t = I N ×N and U t U † t = C(t), being C(t) the two-point correlation matrix with elements C ij (t) = ⟨ĉ † i ĉj ⟩ t .The state Eq. ( 14) is in a Slater-determinant form [94].
As explained in Ref. [31], for the QSD protocol, the dynamics can be split in two steps accounting for the unitary and the dissipative one.Note that, in general, the operators generating these steps do not commute, therefore this approach is based on a Trotterization of the full dynamics in the limit of small time steps δt.The dynamical step thus reads with Since the whole dynamics does not preserve the norm of the state, to restore the unitarity, at the end of the step in Eq. ( 15) we have to perform a QR-decomposition where For the QJ protocol, we exploit a similar procedure.The matrix after the projective step Eq. ( 9) on the ℓ-th site is obtained from V t+δt = T U t with while the state after the non-Hermitian step Eq. ( 10) is obtained by applying V t+δt = e −i H eff δt U t , Even in this case, at the end of each step, we have to restore the unitarity of V t+δt by performing a QRdecomposition [see the above Eq.( 17) and the related discussion].

Comparison of the two unravelings
In this section we discuss the differences in the dynamical behaviour of the system, under the action of the two unravelings described in Sec. 3. In particular we focus on: (i) the distribution of the expectation values of the local number operator, once this is measured in time over the various trajectories (Sec.5.1); (ii) the localization properties of the states averaged along the trajectories (Sec.5.2).In all the following we are going to perform evolutions up to a time t f = 10 3 , and choose δt = 0.05 for QSD and δt = 0.16/L for QJ, in order to ensure convergence.

Distribution of ⟨n j ⟩ t
We first consider the statistics of the expectation values ⟨n j ⟩ t ≡ ⟨ψ t |n j |ψ t ⟩ detected along the quantum trajectories in the different unravelings.Because of translation invariance of the average dynamics, in the long-time limit the statistics of this expectation value does not depend on the site j, therefore we can drop the explicit dependence.
To further simplify the notation, hereafter we also drop the time dependence and pose n ≡ ⟨n j ⟩ t .As discussed in Sec.3.3, because of the linearity of the trace operation, the ensemble average over the two stochastic processes (QSD and QJ) should coincide [see Eq. ( 12)].We checked this Fig. 1 Temporal behaviour of n ≡ ⟨n j ⟩ t for the two unravelings (orange and blue curves for QSD and QJ, respectively) for γ = 0.1 (a) and γ = 1.0 (b).In both cases, curves oscillate around the mean value n = 0.5.However, for small measurement rates they remain concentrated around the mean value, while for large γ they shift close to the extremes.The various colors refer to different values of γ.These plots have been obtained for L = 128, but the curves are almost independent of the system size (see Fig. 3).
numerically and indeed we found that n = 1/2, independently of the unraveling: This agrees with the fact that we work in the half-filling subspace.
Let us start showing some examples of ⟨n j ⟩ t versus t along single trajectories, for a single noise realization.Examples of such trajectories are plotted in Fig. 1 for (a) γ = 0.1 and (b) γ = 1.0.The colors refer to the two stochastic processes, QSD (orange) and QJ (blue).For both unravelings, the trajectories fluctuate around the ensemble average value n = 1/2, however the nature of these oscillations depends on γ.In fact, while for small measurement rates they stay close around 0.5 (panel a), 3 Distributions of n along a time trajectory, averaged over Nr = 80 noise realizations for QSD (top panels) and QJ (bottom panels) protocols, for γ = 0.1 (left column) and γ = 0.4 (right column).The various colors refer to different system sizes.
for larger γ values they shift and spend more time near the extremal values 0 and 1 (panel b).
This can be better appreciated in Fig. 2, where we show the probability distribution P (n) along a time trajectory for different values of γ and for both QSD (panel a) and QJ (panel b) unraveling protocols.We evaluate this distribution as the histogram of ⟨ψ t |n j |ψ t ⟩ along the trajectory, averaging over sites j and noise realizations, in order to get smoother curves.Furthermore, we evaluate each distribution over a set of trajectories lasting a time t f , large enough that all transients have vanished.The main emerging feature is that, for increasing γ, the distribution crosses from a unimodal to a bimodal character, reflecting the differences in the trajectories discussed in Fig. 1.We refer to this behaviour as a "quantum bifurcation" [75].
Looking at the shape of P (n), we can see some differences between the two unravelings.In fact, because of the particle-hole symmetry (see Sec. 3.1), in the QSD dynamics the distribution is symmetric around the ensemble average value n = 1/2.This behaviour is reminiscent of the minima of the free energy in the Landau model.In contrast, for the QJ protocol we do not observe the same symmetry, consistently with the fact that such unraveling breaks the particle-hole symmetry (see Sec. 3.2), but we can still identify a crossover from a unimodal asymmetric distribution to a bimodal asymmetric one, where a second local maximum appears.In any case, we carefully checked that, even in this latter case, the average value n = 1/2 is consistent with the half-filling assumption and the conservation of particle number.
We should also emphasize that the numerical data reported in Fig. 2 are for a system with L = 128 sites.However they are almost independent of the system size, as explicitly shown in Fig. 3, in analogy with the quantum Zeno effect [76][77][78][79].As expected, a more evident dependence on L emerges for smaller values of γ (left panels, γ = 0.1), while for larger measurement strengths (right panels, γ = 0.4) differences in the curves from L = 16 to L = 128 become hardly visible on the scale of the figure.
To gain further insight on the properties of the bifurcation, in Fig. 4 we show the position of the two local maxima of the distribution, n max − and n max + , defined such that P (n max + ) ≥ P (n max − ).In Fig. 4(a) we consider the case of the QSD protocol: We see that, at the bifurcation point γ ⋆ QSD ≈ 0.2, the two maxima stem symmetrically (i.e., n max + − n = |n max − − n|) and continuously from the single maximum, and the derivative with respect to γ is discontinuous at the bifurcation.This behaviour is formally similar to the mean-field analysis of a second-order phase transition (the maxima here correspond to minima of the free energy there), or a pitchfork bifurcation in classical dynamical systems [80][81][82].In Fig. 4(b) we consider the case of the QJ protocol: The maxima are not symmetric, and we call the one corresponding to the smaller value as "secondary maximum".Here the bifurcation occurs in a different way compared to the QSD case, because at the bifurcation point γ ⋆ QJ ≈ 0.23 the secondary maximum appears discontinuously.For γ ≈ 0.35 the global and the secondary maxima swap with each other.This behaviour is formally similar to what happens in a first-order phase transition, with the maxima of our distribution corresponding to the minima of the free energy in that case.
We notice that, for both unravelings, at γ ≈ 0.4 the position of the global maximum approaches the value n max + = 1.This is consistent with the fact that, for large values of γ, the measurement prevails and the state of the system is close to a product one, as the results on the asymptotic entanglement entropy for the QSD unraveling confirm [31,36,39,73].We have checked that the entanglement entropy behaves similarly also for the QJ case (such kind of similarity in the entanglement behaviour for different unravelings of the same Lindblad dynamics contrasts with what has been observed with the additional presence of nearest-neighbour coherent pairing terms, which modify the symmetry properties of the Hamiltonian [43,44]).Nevertheless we stress that the entanglement transition and the bifurcation are different and unrelated phenomena [75], thus the above analysis of the probability distributions is not expected to provide information on the entanglement behaviour.

Inverse participation
A related question is whether the change in the behaviour of the distributions of ⟨n j ⟩ t could be related with some change in the space-delocalization properties of the state.To this purpose, we study the time averaged IPR of the components of the Slater determinant, defined as with In general, the IPR is such that 1/L ≤ µ IPR µ ≤ 1, where the two bounds refer to a perfectly delocalized and a localized state, respectively.Fig. 5 The factors α QSD (orange dots) and α QJ (blue squares) vs γ, as obtained by a numerical fit of the IPR data with a power-law as a function of L, cf.Eq. ( 22) (see the inset in log-log scale, for the QSD protocol).
Let us first focus on the QSD unraveling.To characterize the localization properties when varying the measurement strength, we numerically evaluate the IPR as a function of the system size and for different values of γ, as shown in the inset of Fig. 5. Different lines refer to various values of γ, from γ = 0.04 to γ = 1, with steps of δγ = 0.04.Qualitatively similar results have been obtained for the QJ protocol (for this reason, we decided to show only data for the QSD protocol).We then fitted the data with a power law and found that both unravelings exhibit a similar scaling with L, but with a power-law exponent α(γ) that depends on the considered measurement protocol (main frame of Fig. 5, where the orange dots refers to QSD, while the blue squares to QJ).In both cases, the single-particle states are perfectly delocalized (α = 1) only in the γ → 0 limit.For γ > 0, the scaling exponents deviate from 1 and differ from each other.This behaviour is typical of anomalous delocalization of states with multifractal properties [86,87].The fact that we do not observe discontinuities for α as a function of γ, seems to indicate that delocalization properties are independent of the possible presence of an entanglement transition.In particular, α continuously decreases with increasing γ, marking that the system becomes less and less delocalized, but never attains localization.This result fits with the findings of [83], where the authors see that the conductivity decreases as a power law with γ but never vanishes.

Conclusions
In conclusion we studied the dynamics of a free-fermion chain under dephasing, considering two protocols of random measurements that, once averaged over the quantum trajectories, provide the same Lindbladian: From one side we considered the QSD unraveling that preserves the particle-hole symmetry, and from the other side a QJ unraveling that breaks such symmetry.In both protocols the state of the system can be always cast as a Slater determinant, allowing us to analyze quite large sizes.We focused on the expectations of the measured operators along the trajectories and studied the properties of their distributions.When the measurement strength γ lies below a certain threshold, the distributions are unimodal, while above this threshold they become bimodal, giving rise to a bifurcation transition.This bifurcation threshold is located at γ ⋆ QSD ≈ 0.2 for QSD and at γ ⋆ QJ ≈ 0.23 for QJ, and QSD preserves particle-hole symmetry providing distributions symmetric around n = 1/2, while QJ breaks this symmetry and gives rise to asymmetric distributions.For QSD, two symmetric maxima appear at the bifurcation, in formal analogy with second-order phase transitions.The two maxima stem from the single one continuously, with a discontinuity in the derivative in γ, as occurring for the pitchfork bifurcation in classical dynamical systems.On the other hand, in the QJ case the two maxima are asymmetric and the secondary one appears discontinuously, as occurring in first-order phase transitions.Quite curiously, the distributions display a very weak dependence on the system size.
Finally we analyzed the behaviour of the IPR of the components of the Slater determinant, averaged over time and randomness.We considered the scaling of this quantity with the system size L, to understand the space localization properties of the state, in analogy with what is done in Anderson-localization problems.We always find a power-law scaling of the form ∼ L −α , where the exponent α equals 1 for γ → 0, marking perfect delocalization.For γ > 0, we find that α depends continuously on γ and for both unravelings 0 < α < 1, marking a multifractal behaviour of the Slater-determinant components.This is what happens at the localization-delocalization transition in Anderson-localization problems, so we can say that our model is always in an Anderson critical phase.
Perspectives of future work focus on the study of properties of the measurement-operator distributions and of the IPR in nonintegrable monitored models [51,95,96], in order to understand the possible existence of quantum bifurcations or localization transitions in these cases.

Fig. 2
Fig. 2 Distributions of n along a time trajectory, averaged over Nr = 80 noise realizations for QSD (panel a) and QJ (panel b) protocols.The various colors refer to different values of γ.These plots have been obtained for L = 128, but the curves are almost independent of the system size (see Fig.3).

Fig. 4
Fig. 4 Position of the local maxima of P n vs γ for QSD (panel a) and QJ (panel b) protocols.For QSD the two maximua stem continuously and symmetrically at the bifurcation point γ ⋆ QSD ≈ 0.2, with a discontinuity in the derivative with respect to γ.For QJ the maxima are asymmetric: the secondary one appears discontinuously at γ ⋆ QJ ≈ 0.23, while the global and the secondary one swap each other at γ ≈ 0.35.Numerical parameters: L = 128, Nr ≥ 48.