Measurement induced transitions in non-Markovian free fermion ladders

The Referee writes: The current work tries to understand the interplay between non-Markovianity of the bath and the steady-state entanglement of the system. To facilitate the discussion, the authors consider a one-dimensional complex fermion ladder with particle-number conservation, where one of the two chains is measured and serves as a non-Markovian bath. The dynamics is designed to preserve the Gaussainity of the state so that the problem is numerically tractable. The authors examine two different cases, the persistent measurement where the measurement probability is 1, and the sporadic measurement where the measurement probability is less than 1. In both regimes, non-area law behaviors are observed when the bath dynamics is fast (t2 being large). The authors also invoke a quantitative measure of the non-Markovianity.

I found both the question and the setup in this manuscript interesting. Attempt to provide a quantitative characterization of the non-Markovianity in a many-body system is also very nice. However, I found that some of the results are not clear to me and not conclusive enough. The quality of the numerics needs be improved.

I would like to hear the author's responses to the following questions first before making any decision.

Our Response: We thank the Referee for carefully reading our manuscript. Below, we address in detail their comments.

The Referee writes: 1. Is it possible to draw some phase diagrams for the two cases of persistent and sporadic measurement?

Our Response: In Fig. 2b and 7 we drew phase diagrams in the $t_{12}$--$t_2$ space, based on the relative difference in entanglement entropy/negativity. As specified in the text, this quantity can identify the area law phase, but does not distinguish between volume and logarithmic scaling for the non-area-law phase. For that a scaling with $L$ is necessary, which requires more extensive numerical simulations and limits the parameter space we are able to explore within a reasonable time frame.

On the other hand, we can draw more precise phase diagrams in the $t_2$ space for different values of $p$.

The Referee writes: 2. For the persistent measurement, since the bath has no entanglement after the measurement, it seems to me that the entanglement of the system fermion is the same as the entanglement of the system plus the bath. Namely, it is effectively probing a single-chain entanglement transition. And then there is no volume-law phase in this case. Is it possible to make some connection between the current setup and that of Ref 53 to get a more clear conclusion? E.g. whether the non-area law exists in the thermodynamic limit, and is the transition also of KT type?

Our Response: We thank the Referee for raising this point. In the regime $t_1\gtrsim t_{12}\gg t_2$, we agree with the Referee that for $p=1$ we obtain a setup equivalent to that of Ref. [53]. In fact, for large $t_{12}/t_2$ every site on the inner chain is coupled to an essentially static site on the outer chain. In this regime, measuring the outer chain site is equivalent to performing a measurement -- which may not be exactly a measurement of the particle number occupation -- on the inner chain site, without any memory effect (given that the dynamics of the outer chain is very slow). For all other regimes, we do not expect to find such a similarity with the single chain of Ref. [53] since memory effects play an important role.

The Referee writes: 3. From appendix A, the half-chain entanglement (both the von Neumann entropy and the negativity) exhibits a significant fluctuation, e.g. Fig 11 and 14. The label of vertical axis is the same as eqn7, which implies it is averaged over trajectory although the text says it is for a single trajectory. In any case, the results in the main text, e.g. Fig 2 and 6, have no visible error bar. This is surprising to me. I am wondering how the error bar can be invisible given the fact that the single-trajectory result has a strong fluctuation.

Our Response: The label has a typo and we thank the Referee for spotting this. All plots in Fig. 11 to Fig. 14 refer to the case of a single trajectory. While a single trajectory may have large fluctuations in time, the average over trajectories shows small fluctuations and thus a small error bar. The cases for which the error is large compared to the values of the observable (see Fig. 16a for example) are the ones where the value of the observable is very small, so that the error bar becomes hardly visible due to a low absolute error.

The Referee writes: 4. The result in Fig 7 is obtained for L=16, which is much smaller than the system sizes used in the previous plots. I am curious what happens at larger system sizes?

Our Response: We have results of the behavior of the negativity at large system sizes only for $t_{12}=\pi/2$ and some selected values of $t_{2}$. Obtaining a phase diagram like in Fig. 7 for large system sizes would demand extensive numerical simulations.

However, we can infer the qualitative behavior at larger $L$ from the results in Fig. 8, which indicate a transition from area law to CFT phase as $t_2$ and $p$ are increased for fixed $t_{12}=\pi/2$. For $t_{12}=\pi$, we expect the area law to not be present at all.

Moreover, we have added simulations at larger system sizes and a scaling analysis with $L$ along a cut line in $t_2$ for $t_{12}=\pi/2$

The Referee writes: 5. In Fig 8, is it possible to do some finite-size scaling for subfig b,d,f, as the authors have done in appendix C, to infer what happens in the thermodynamic limit?

Our Response: We thank the Referee for the good suggestions. We performed this analysis, which is now reported in Fig. 8 b), d), f). The trend we observe is that increasing the systems sizes where we perform the fit, the contribution from the logarithmic part increases, suggesting that the linear contribution is a finite size effect (which would disappear in the thermodynamic limit). However, this is not very clear for large $t_2$ and small $p$, where it seems that these finite size corrections are more important.

We also added an Appendix D reporting an analysis of the residuals for logarithmic and linear fit.

The Referee writes: 6. The dynamics has a periodicity in $t_{12}$, which can be seen through most of the plot. However, the plots of the measure of non-Markovianity $N(\phi)$, Fig 9 and 10, do not have such a structure. What is the reason for that?

1. When $t_{12} = \pi$, the two chains are decoupled and the system chain undergoes unitary dynamics. Thus, I would expect $N(\phi)$ to have a dip near $t_{12} = \pi$? However, this is not the case from the plot. Do I miss something and how to properly understand the result?

Our Response: The periodicity (and the decoupling of the chain at $t_{12}=\pi$) stems from the periodic boundary conditions of the ladder in the large $L$ limit. We suppose that the deviation in $\mathcal{N}(\Phi)$ originates from the small system size $L=4$ that we use in the calculation of $\mathcal{N}(\Phi)$.

The decoupling of the chain is exactly true for $t_{12}=\pi$ and $t_1=t_2$ and for large $L$. Thus, in the small sizes limit we do not expect to see a sharp drop in $\mathcal{N}(\Phi)$

The Referee writes: 8. The authors mention two difficulties of calculating $N(\phi)$, and thus only present the result for L=4. I am curious whether one can still utilize the Gaussianity to compute it for a slightly larger system size, e.g. L=10. Specifically, the Gaussianity of a density matrix $\rho$ implies that it maps Slater determinant states to a linear superposition of Slater determinant states, i.e. $e^{-c^\dagger M c} |n\rangle= \sum_m A_{n,m}|m\rangle$, where $|m\rangle$ are Slater determinant states that form a complete orthonormal basis, M is determined via the fermion two-point function and $A_{n,m}$ can also be calculated by utilizing Gaussianity. This way we can write down $\rho_1 - \rho_2$ as a matrix directly, e.g. in the basis of Fock states, without calculating the state of the two chains. I am wondering whether this is helpful in increasing the accessible system size and curious of the authors’ thoughts.

Our Response: The Referee raises a good point, which was also raised by Referee 1.

Gaussianity is in general a multiplicative but not additive property. The state of the system is Gaussian only along each trajectory $\alpha$, but in order to calculate $\mathbf{N}(\Phi)$ we need access to the average density matrix $\rho=\sum_{\alpha}\rho_{\alpha}/N_{tr}$ with $N_{tr}$ the number of trajectories. This density matrix is not Gaussian, and it means in particular that we cannot write $\rho_1 - \rho_2$ as a Gaussian state, so that we cannot calculate the $\mathbb{L}_1$ trace distance defined in Ref. [90] (Breuer et al.) using Gaussian states.

The suggestion of Referee 1 also cannot be applied to Gaussian states, as the relative entropy requires, e.g., the computation of $\log(\rho_1)=\log(\sum_{\alpha}\rho_{\alpha}/N_{tr})$ which again cannot be written as a Gaussian state -- unless one expands the logarithm, in which case it is possible in principle but requires computing infinitely many terms (or, alternatively, establish rigorous bounds on such truncated sums, that we are presently not aware of).

A solution is to use a $L_2$ trace distance $\text{Tr}(\rho_1-\rho_2)^2=\text{Tr}(\rho_1^2-2\rho_1\rho_2+\rho_2^2)$. In fact $\rho_1\rho_2=\sum_{\alpha,\beta}\rho_{\alpha}\rho_{\beta}/N_{tr}^2$, where $\alpha$ labels the trajectories with $\rho_1$ as initial condition and $\beta$ labels the trajectories with $\rho_2$ as initial condition. Each term in the sum is a Gaussian density matrix, so that $\text{Tr}(\rho_1-\rho_2)^2$ can be written as a double average (over $\alpha$ and $\beta$) of Gaussian density matrices, meaning it can be simulated using Gaussian states.

We performed this calculation and included some of the results in the paper, see the new Fig. 11. We remark that this technique does not allow to increase much the size of the system. In fact, while everything scales polynomially with the system size, there is an additional overhead due to the simulation of the quantum trajectories, absent in the exact diagonalization numerics. This overhead is roughly a factor $\sim10^2$. Moreover, computing $\sum_{\alpha,\beta}\rho_{\alpha}\rho_{\beta}/N_{tr}^2$ requires inverting the correlation matrix for each value of $\alpha$ and $\beta$, i.e. roughly $\sim10^4$ times. These additional computationl costs add up, making numerics with Gaussian states impracticable for system sizes larger than $L\sim12\div16$.

In agreement with the ED analysis, we observe that the dynamics is always non-Markovian and that for $t_{2}=5$ the degree of non-Markovianity is enhanced for smaller values of $p$. We do observe different trends (discussed in detail in the manuscript) of the non-Markovianity as function of $t_2$ compared to ED simulations, but we cannot make more definitive statements, since the two analysis are based on different trace distances.

The Referee writes: 9. $N(\phi)$is computed for the system. The introduction emphasizes that previous studies focus mainly on Markovian bath. Therefore, I am wondering whether we should instead measure the non-Markovianity of the bath, which addresses the effect of the bath memory on the system dynamics more directly. In this case, one may also need to modify the protocol a bit to actually control the memory of the bath. Maybe I misunderstand the main question the introduction posts, but could the authors clarify my confusion?

Our Response: The Referee raises a very good point. Information back-flow from system to bath also means there is an information back-flow from bath to system, so non-Markovianity of the bath means non-Markovianity of the system. However, calculating the non-Markovianity of the bath dynamics using $N(\phi)$ has the same (exponential) computational cost as the non-Markovianity of the system.

We also tried computing the correlation functions of the bath, to show that they do not decay like a Markov delta function but rather show memory effects and finite correlation times. However, this numerical computation as well is exponentially hard to perform as it needs to be exact diagonalization techniques and cannot be calculated using Gaussian states.

The referee writes: I think the idea of the paper is interesting and has potential to lead to insightful results or new phenomena. The results as presented in this version seem, to me, a bit preliminary: there is numerical evidence for different behaviors, but not enough resolution to locate a transition or to distinguish a logarithmically-entangled phase from a volume-law.

The $I_{L/8}$ diagnostic seems to rule out the volume-law but that is a bit indirect, as the "volume-law" ansatz ($I_{L/8}->0$) is based on the random circuit model, from what I can tell; whereas this particular free-fermion ladder could in principle realize a volume-law entangled phase with different universality (say, one based on "rainbow" states). Absent any hints of saturation in the entropy or negativity it is hard to rule out a volume-law phase.

I think the difficulty in interpreting the numerics comes in part from the fact that the system sizes studied are a bit small ($L\leq256$, counting both legs of the ladder); in this field I think simulating ~1000 fermions should be fairly straighforward, and that could help sort out some of the finite-size issues. Other diagnostics, such as tripartite mutual information and purification of a probe qubit, could also be helpful.

For the non-Markovianity diagnostic, the size limit (L=4) is particularly problematic. I would suggest a modified measure based on the relative entropy instead of the trace distance. The relative entropy is also monotonic under channels but has the advantage of being efficiently computable for fermionic Gaussian states (since log(rho) is a quadratic Hamiltonian). I don't know if any desirable properties get lost but it seems like this should capture the desired property and greatly alleviate the size limit.

Our response: We thank the Referee for the positive overview of our work and for the useful comments.

Regarding the diagnostic based on $I_{L/8}$, we also employ a diagnostic based on $I_{L/4}$ at the same time, making our analysis more accurate than with just one of this two tools.

We have performed additional simulations, extending system sizes up to $L=256$ for each chain (so 512 in total), which is a system size comparable with other free fermions simulations. We remark that an additional difficulty of performing numerical simulations on this system arises from the measurements part of the evolution, whose complexity scales as $\sim L^3$ and is typically larger than the complexity of the unitary dynamics. Moreover, we are analyzing a 2D phase diagram in a system that depends on many parameters ($t_{12}$, $t_2$, $p$) - a task which requires greater resources.

The idea of the Referee to look at the purification of a reference qubit is interesting, but there are no benchmarks in the literature on how to apply this method to a mixed system. Indeed coupling the probe qubit to the system alone, without also coupling it indirectly to the bath, is complicated and could lead to spurious effects in the purification of the qubit.

Studying the tripartite mutual information is very useful and interesting for a precision study of the transition, but this is beyond the scope of our work, since our main goal is to distinguish between the phases, for which a combination of entanglement entropy and bipartite mutual information is already sufficient. We have commented on this in the conclusion.

The suggestion of the referee is interesting and we comment more in detail on non-Markovianity later in the reply.

The referee writes:

In addition, I have some questions. 1. It's not clear to me exactly how the diagnostic of non-Markovianity Eq.(22) is computed here. Starting from different initial states, there are different probability distributions for outcomes, etc; are these drawn independently for the two states? Is the measure computed between pairs individual realizations, and then averaged? Or are the states averaged over realizations first? It would be good to clarify this better. Depending on how this is done, I think it may return non-zero values even in purely Markovian cases, since the dynamics under consideration is not a channel. A good sanity check would be to apply this to the p=1 model but also reset all bath modes to some fixed configuration after measurement (this eliminates all memory from the bath and makes the dynamics Markovian).

Our response: We realize that the text was not very clear on this point. We do not use a quantum trajectories approach, but rather perform an evolution of the density matrix of the full system (system chain + bath chain) according to the Lindblad master equation that describes the dynamics. This is done using exact diagonalization techniques, hence the small system sizes: $L=4$ corresponds to a density matrix of a system with 8 sites, which requires manipulating matrices of size $2^{16}$.

We do not perform any averaging, as the Lindblad dynamics already describes the average over quantum trajectories. The Lindblad evolution is computed for many random pairs of initial density matrices $(\rho_1,\rho_2)$. At each time the degrees of freedom of the bath chain are traced out, and the distance between the two reduced density matrices $d_{\rho} (\rho_1, \rho_2)$ is computed. The integral over the regions of positive $\partial_td_{\rho}(t)$ is then maximized over $(\rho_1,\rho_2)$.

The referee writes: 2. I think the volume-law ansatz for $I_{L/4} ~ L^{1/3}$ [Table I] is based on the case of two contiguous intervals (i.e. touching at one endpoint). I think for the case of antipodal intervals the scaling should be O(1), based on the configuration of minimal membranes.

If the authors can address these questions and improve the paper accordingly then I am in principle in favor of publication.

Our response: The scaling of the bipartite mutual information $I_{L/4} ~ L^{1/3}$ reported in Ref. [84] (Zabalo et al.) refers to antipodal intervals, see Fig. 1b of the main text of Ref. [84], pag. 23 of the Supplemental Material and Fig. 17 in the Supplemental Material. Actually, we are not sure that there is an independent method (e.g., beyond minimal cut) to derive such scaling. Indeed, in our manuscript we calculate the bipartite mutual information for antipodal intervals; even though the work in Ref. [84] refers to stabilizer circuits, we would expect a volume law arising in our model to exhibit the same behavior.

We remark that the possibility raised by the Referee of considering contiguous intervals is also interesting: we are unsure if this would be able to distinguish the two phases at $p=1$, but such configuration might be useful for $p<1$ case (at least, in the context of negativity).

The referee writes: the results are inconclusive in terms of establishing a phase diagram, phase transition, or different universal properties of the phases considering the fermionic Gaussian nature of the model, the sizes studied seem small relative to standards in the field, which may contribute to weakness (1) some lack of clarity in the meaning of "(non)-Markovianity" in the present context Extend size of simulations and/or try other diagnostics (e.g. tripartite mutual information) to better control finite-size effects Consider alternative non-Markovianity diagnostic based on relative entropy (see Report) Check comments 1 and 2 from report and correct as needed

Our response:

The behavior of the mutual information is a good diagnostic tool for $p=1$, and together with the scaling of the entanglement entropy, allows to distinguish between area-law, CFT and volume-law phases. We focus on distinguishing between the phases rather than performing a precision study of the transition, so that these diagnostic tools combined with the numerical simulations for the entropy should suffice to achieve that.

For what concerns the $p<1$ case, we performed simulations of the negativity for larger system sizes, up to $L=256$. These are helpful in distinguishing volume-law from CFT phases. It would also be interesting to investigate tripartite correlations for mixed states, but we are not aware of any result in this scenario that we could use as a benchmark, so that this tool would be of limited utility.

The suggestion of the Referee to use a non-Markovianity diagnostic based on relative entropy is interesting and is also suggested by Referee 2.

Gaussianity is in general a multiplicative but not additive property. The state of the system is Gaussian only along each trajectory $\alpha$, but in order to calculate $\mathbf{N}(\Phi)$ we need access to the average density matrix $\rho=\sum_{\alpha}\rho_{\alpha}/N_{tr}$ with $N_{tr}$ the number of trajectories. This density matrix is not Gaussian, and it means in particular that we cannot write $\rho_1 - \rho_2$ as a Gaussian state, so that we cannot calculate the $\mathbb{L}_1$ trace distance defined in Ref. [90] (Breuer et al.) using Gaussian states.

The suggestion of Referee 1 also cannot be applied to Gaussian states, as the relative entropy requires, e.g., the computation of $\log(\rho_1)=\log(\sum_{\alpha}\rho_{\alpha}/N_{tr})$ which again cannot be written as a Gaussian state -- unless one expands the logarithm, in which case it is possible in principle but requires computing infinitely many terms (or, alternatively, establish rigorous bounds on such truncated sums, that we are presently not aware of).

A solution is to use a $L_2$ trace distance $\text{Tr}(\rho_1-\rho_2)^2=\text{Tr}(\rho_1^2-2\rho_1\rho_2+\rho_2^2)$. In fact $\rho_1\rho_2=\sum_{\alpha,\beta}\rho_{\alpha}\rho_{\beta}/N_{tr}^2$, where $\alpha$ labels the trajectories with $\rho_1$ as initial condition and $\beta$ labels the trajectories with $\rho_2$ as initial condition. Each term in the sum is a Gaussian density matrix, so that $\text{Tr}(\rho_1-\rho_2)^2$ can be written as a double average (over $\alpha$ and $\beta$) of Gaussian density matrices, meaning it can be simulated using Gaussian states.

We performed this calculation and included some of the results in the paper, see the new Fig. 11. We remark that this technique does not allow to increase much the size of the system. In fact, while everything scales polynomially with the system size, there is an additional overhead due to the simulation of the quantum trajectories, absent in the exact diagonalization numerics. This overhead is roughly a factor $\sim10^2$. Moreover, computing $\sum_{\alpha,\beta}\rho_{\alpha}\rho_{\beta}/N_{tr}^2$ requires inverting the correlation matrix for each value of $\alpha$ and $\beta$, i.e. roughly $\sim10^4$ times. These additional computationl costs add up, making numerics with Gaussian states impracticable for system sizes larger than $L\sim12\div16$.

In agreement with the ED analysis, we observe that the dynamics is always non-Markovian and that for $t_{2}=5$ the degree of non-Markovianity is enhanced for smaller values of $p$. We do observe different trends (discussed in detail in the manuscript) of the non-Markovianity as function of $t_2$ compared to ED simulations, but we cannot make more definitive statements, since the two analysis are based on different trace distances.