Reduced Dynamics of Full Counting Statistics

Dear SciPost Editors,

We hereby resubmit our manuscript to SciPost Physics. We thank both Referees for their detailed and useful reports, which helped us improve the manuscript in several meaningful ways. The two Referees made contradicting recommendations, one positive and one negative.

Referee 1 had reservations about the novelty of our work due to the impression that we were simply applying a discrete approximation to an existing method. We believe the root of this stems from a misreading of the text, and we apologize for the lack of clarity that led to it. In Sec. 2, we were not explicit enough in stating that generalizing the Nakajima--Zwanzig formalism to full counting statistics is the first main result of our paper. To our knowledge, this has never been done previously; nor has any other exact formulation of reduced dynamics in the presence of counting fields been made in the literature.

We then go on to generalize the transfer tensor formalism to counting statistics in Sec. 3, which is our second main result. This formalism is a discrete version of Nakajima--Zwanzig that describes dynamics stroboscopically. As such, applying it to many types of experimental and numerical data is more straightforward. Importantly, the connection between the Nakajima--Zwanzig formalism and transfer tensor method is nontrivial, and was not guaranteed to continue working when counting fields are thrown into the mix. Nevertheless, we show that it does work. This lets us put the new technique to work and test its efficacy in Sec. 5.

The fact that reduced dynamics can be exactly formulated in the presence of counting fields is remarkable. As Sec. 5 demonstrates, it can bypass crucial bottlenecks in existing methods for numerical simulation by extending them to longer times. Thus, our manuscript introduces a powerful new tool for researchers studying non-Markovian open dynamics in research areas spanning condensed matter physics, chemistry, and material sciences.

Referee 2 had many positive comments and some profoundly interesting questions about our novel methods. We have used these comments and questions to better our manuscript. We believe that physics is now more transparently exposed as a result of this exchange.

A detailed response to both Referees is below. We feel that the modifications, made based on their comments, have resulted in a far more approachable manuscript. We hope the Referees agree with this sentiment and recommend publication in SciPost Physics.

Sincerely,

Felix A. Pollock, Emanuel Gull, Kavan Modi, and Guy Cohen

Response to Report 1 -- Reviewer's comments quoted in italics

The paper is devoted to the development of techniques that should allow the numerical study of quantities relevant for the description of quantum transport in a regime allowing for strong coupling and memory effects.

The authors make reference to the full counting statistics approach, with the aim to determine the generating function and from it the statistics of quantities such as currents.

We thank the referee for this summary of our work, and recognizing that the results have significance in an active and practical area of research.

They introduce the evolution map for the overall state in the presence of the counting field. The trace of this quantity should provide the generating function from which observable quantities can be evaluated. To obtain this map they consider the standard Nakajima-Zwanzig approach, applied to the evolution including counting field. They further point to a discrete representation in time of the time evolution map, borrowed from what is referred to as transfer tensor method.

We must respectfully point out that, to the best of our knowledge, neither the Nakajima--Zwanzig formalism nor the transfer tensor formalism has ever been applied to counting-field modified dynamics. We had not clearly stated in our manuscript that both of these approaches are novel and first introduced here. In fact, we believe the present manuscript is the first to derive any kind of formally exact reduced dynamics (RD) technique applicable to the study of full counting statistics. We have amended some of our phrasing regarding this point, which thanks to the referee we now realize was too circumspect.

They further point to a discrete representation in time of the time evolution map, borrowed from what is referred to as transfer tensor method. If I understand correctly, in this approach the evolution is rewritten by replacing the exact evolution with the one in which the map obeys a composition law (is divisible in the language used in the paper), together with correction terms, so that keeping all contributions one still has the exact evolution.

The transfer tensor method (TTM) and the Nakajima-Zwanzig equation (NZE) both lead to a reduced description in terms of a family of completely-positive maps. Both methods are exact, but TTM is stroboscopic, while NZE is expressed at the continuum limit. In practice one often has stroboscopic dynamical data, which makes the TTM practically appealing.

We note that TTM originates in the work of Cerrillo and Cao [PRL (2014)], where it was seen as an ansatz. However, it was formally derived as an exact formulation by two of the present authors in [Quantum (2018)].

The aim of the paper is to compare the results obtained with this approach, which is a way of discretizing, with numerically exact quantum Monte Carlo techniques.

With respect, we believe this is a misunderstanding of our aim, probably due to our less-than-clear presentation in the first version of the manuscript. The primary aim of the paper is to introduce and test RD techniques able to address full counting statistics. The approach is not "a way of discretizing", but rather a way to take advantage of short memory timescales in order to obtain long-time dynamics from short-time data.

Our second main achievement is that this method allows us to obtain results of long-time dynamics from the short-time data. The Monte Carlo results are therefore used not only as a benchmark, but as the source of the short-time data. With the aid of RD, an inexpensive Monte Carlo simulation up to short times can now be augmented, granting access to long dynamics. We believe the confusion regarding our purpose of the paper is related to the previous points, and should now be resolved.

Dealing with transport properties in strong coupling regime is certainly a topic of interest, but I do not find the paper clearly written, and the real utility of the technique is also not apparent. Neither the theoretical framework is fully worked out, nor the real advantage with respect to other techniques is spelled out, the authors only consider rough compatibility with other results.

We appreciate this critical remark by the referee. It has allowed us to clarify the purpose of our paper in the current version, in which we hope the utility of the technique is now fully transparent. In terms of its advantages, as we now emphasise in the text, this is the first technique of its kind (i.e. RD for counting statistics), and there is nothing to directly compare it to.

While this goes beyond the scope of the present manuscript, the underlying technique used to obtain the short-time data (here inchworm Monte Carlo) could be compared to other numerically exact methods for simulating nonequilibrium quantum transport; and it has been in several of our previous works. However, in the context of its use in RD of full counting statistics, we are also not aware of other methods able to perform such simulations for generic fermionic transport. We nevertheless now mention more related work that has been done at certain analytically tractable limits of related models, and for the simpler spin--boson model. We also contrast it with the methods used here to avoid confusion by future readers.

On these grounds I do not support publication of the submitted paper. See below for further comments.

Given the clearer presentation in the current version of the manuscript, together with our answers above, we hope that the referee will now have a favourable view of publishing the work in SciPost Physics.

In the paper is not at all clear how the two theoretical Ansatz (both formally exact but deemed to be used as approximations) of Sections 2 and 3 are combined together. How is the Nakajima-Zwanzig expansion used in determining the transfer tensor? How are the approximations in both parts of the treatment combined and validated? The relationship between the two theoretical Ansatz is generally not clear also in the rest of the paper, e.g. in the Conclusions they are mentioned as kind of alternatives. What figure of merit do the authors consider to validate the approximations?

We apologize for the confusion and thank the referee for pointing out our poor presentation in the first version of the manuscript. The NZE and TTM are both fully equivalent and formally exact, for either regular dynamics or FCS dynamics. Approximations arise only when we truncate the memory kernel at a finite time, and then use the short-time data to obtain long-time results. This can be further complicated by noise in the short-time data, if present. We have amended the text above Eq. (13) to avoid this confusion.

Our reason for deriving the continuous time approach in addition to the transfer tensor approach is not to derive the transfer tensor approach. Rather, it is two-fold: first, it provides a more explicit connection to the (more fundamental) continuous time dynamics that will be most familiar to many readers; second, it opens the door to the application of other numerical techniques that may be better suited to the continuous time domain (as in the aforementioned Ref. [44]). We now state this explicitly in a new paragraph in the beginning of Sec. 2, in order to avoid similarly confusing future readers.

Section "Constructing the dynamical maps from data" appears a bit as a digression, what is its role in the paper?

As noted earlier, the immediate practical benefit of the RD techniques presented here is as a way to extend the range of simulations. We therefore found it prudent to describe the process of constructing dynamical maps from the output of such simulations. Without this, it may not be clear to a reader who is unfamiliar with the formalism of dynamical maps how our work might be used in practice. Furthermore, this is the actual process used to obtain the results we've presented, and its description is necessary to make our work fully reproducible.

What is the actual relevance of the norm of the transfer tensor? Is it just a way to estimate weight of different contributions? What can be learnt from its (apparently weak) $\lambda$-dependence?

The reviewer is correct that the norm is simply an estimate of the magnitude of contributions from the different terms in Eq. (19), i.e. from memory effects at different time scales. We now state this where the figure is first mentioned in the main text. As we elaborate in the third paragraph of Section 5, "the weak response in the norm suggests that the effective range of the memory kernel depends on $\lambda$ only weakly, at least in this regime." This implies that memory cutoffs used for computing regular dynamics should still be valid at finite $\lambda$. We now mention this implication in the text.

Figure 4: the lines appear to converge at variance with the steady state current reconstructed from long time Monte Carlo data. The comment to the Figure in the text does not clarify the meaning of the obtained calculations.

We have now added a sentence to the last paragraph of Section 5 explicitly stating that this convergence is what "indicates that by using the FCS-enabled TT method presented here, steady state currents can be reliably obtained from inchworm Monte Carlo simulations."

I do not see the purpose of Appendix A. It provides the standard derivation that can be found in many textbooks. Furthermore where does the divisibility condition mentioned in the first line plays a role?

The purpose of this Appendix is to demonstrate that the derivation of the NZM follows for a broader class of one-parameter families of maps than is usually considered. Note that no reference is made to unitary time evolution; instead, the only property of $X_{t:s}$ that is assumed is Eq. (24). Having derived a NZM-like equation from these minimal assumptions, we then use the fact that $\mathcal{R}_{\lambda, t}$ satisfies said assumptions to prove the validity of the expressions presented in the main text. Our goal was to present a formal underpinning for Eq. (8) and the results that follow, which are (to the best of our knowledge) novel. The reviewer's question about the divisibility condition is a good one---in fact, this property is not used in the current version of the derivation, and we no longer mention it.

It appears that the index $t_m$ in eq.35 of Appendix B is a typo.

We thank the referee for spotting this, and have corrected it to $t_l$.

What is the meaning of the $\simeq$ symbol in Eq.(13) as well as (14)? Is it not an identity?

As previously stated below Eq. (14) (and now more explicitly above Eq. (13)), the $L$ quantities represent numerical approximations to matrix elements of the density operator at different times for various initial conditions. These could be produced, for example, through Monte Carlo simulation of a particular model. As such, when combined as prescribed in the aforementioned equations they will form an object that is close, though not exactly equal, to the actual dynamical map in question.

Some sentences are not clearly formulated:

p.2 "effective non-Markovian equations of motion for quantities with the dimensionality of the small interacting region only" what is the meaning of dimension and region?

The region referred to here is $S$, introduced at the start of the paragraph. And the dimension is that of the Hilbert space acted on by the operators corresponding to the quantities in question. We have clarified this in the new version.

p.3 "the conservations of complete positivity" it sounds like completely positive is a quantity to be conserved, it is rather a possible property of a map

We have changed this to "the complete positivity of its solutions". We thank the referee for pointing it out.

Response to Report 2 -- Reviewer's comments quoted in italics

The authors introduce a new method for computing the statistics of non-local observables in open quantum systems. The method is a combination of the transfer tensor approach (ref [65]) with the counting field technique, that allows to describe the statistics of environmental observables by looking at a suitably modified system dynamics (this principle is explained, for instance, in ref [15]). The method presented in the manuscript is used to considerably reduce the time required for computing the generating function of the statistics of particle currents. The manuscript is detailed and well written, and the subject is interesting for a broad audience, including the quantum transport and quantum thermodynamics communities. The FCS combined with the transfer tensor approach represents a step forward in assessing the role of strong coupling and non markovianity in the statistics of particle and heat currents. In addition, it represents a novelty with respect to more common approaches based on non-equilibrium Green's functions and path integration.

I think the paper is suitable for publication on SciPost Physics, I have some minor comments and questions:

We thank the referee for this summary and for the recommendation to publish in SciPost Physics. The comments regarding the novelty of our derivation also pertain to the main criticism by the first referee, to whom we failed to convey this point in the previous version of the manuscript.

1) In equation (3) the observable $A$ is generic, however, the generating function can be written in the form (3) only if the operator $A_{t_0}$ commutes with the initial density matrix $\rho^{SE}_{t_0}$, am I right? This is not a great limitation, since typically the baths are supposed to be at equilibrium and their state commutes with the particle number operator. The authors could add a clarification on this point.

We thank the referee for pointing to this oversight---we have now corrected it.

2) Eqs. (3), (8), (9) and (14) can be used both for local and non-local observables. Are there some qualitative differences between the two cases? For instance does the superoperator $\mathcal{S}_{\lambda,t}^S$ maintain its explicit dependence on $\lambda$ even if $A$ is an observable living exclusively in the baths Hilbert spaces?

Yes. In general, unless the Hamiltonian trivially separates into independent system and bath components, $\mathcal{R}_{\lambda,t}$,

and hence $\mathcal{S}_{\lambda,t}^S$, will maintain a $\lambda$ dependence for most times.

This is because the exponentiated $A_t$ operators appended by the $\mathcal{A}^{-}_{\lambda, t}$ in Eq. (6) are sandwiched between non-product system--bath operators, in such a way that the subsequent projection and partial trace in Eq. (9) will not remove them.

3) Above equation (15) the authors discuss the role of positivity and trace preservation in reducing the probability that the approximation scheme gives an unphysical dynamics. About the case with non-zero $\lambda$ they comment "For finite values of $\lambda$, there are no analogous universal conditions that the dynamical maps must satisfy". I agree that, as a general statement, the sentence above is correct and there are not other conditions that can be imposed apart from the (15). However, in more specific cases, for instance when the 2 baths are initialized in different Gibbs states and $A$ is the Hamiltonian of one of the two, the generating function (and then the generator of the modified dynamics) has to satisfy a steady state fluctuation theorem. In this last framework (and in all the other cases in which integral or detailed fluctuation theorems are satisfied) do the authors think that these symmetries could be useful to stabilize the construction of the dynamical maps?

This is a beautiful question! In fact, at steady state the fluctuation theorem connects pairs of the Fourier components (or population change distributions) of the time derivative of the generating function to each other. This clearly entails constraints on the dynamics, as the referee suggests. However, Markovian master equations successfully obey these constraints. It therefore appears that rather weak conditions on the TT equivalent of the integral over the memory kernel (i.e. having the asymptotic rates obey generalized detailed balance relations) would be sufficient to guarantee fluctuation theorems, and there is no unique way to enforce them.

Another interesting route would be to consider a generating function with a complex counting field $\lambda$. There, fluctuation theorems guarantee certain symmetries between the generating function's value on the real axis; and its value on a parallel, reversed axis displaced from the real axis by a temperature and voltage dependent constant. It is not immediately clear how this property might be retrofitted onto constraints on the TT.

Either of these could be the starting point for an interesting and promising exploratory research project, and we thank the referee for the suggestion. We have noted the possibility in the manuscript.

4) The plots in fig. 3c show interesting non-markovian features of the current. I found the discussion about the physics of these plots a bit lacking. For the unoccupied initial condition, what is the physical reason for which the current suddenly increases before stabilizing to lower values?

We thanks the referee for this interesting question, and have amended the discussion to answer it. Regarding the current's rapid initial rise before reaching the plateau at a lower value, this is easily understood when you consider that our definition of the current is the time derivative of the population in the left lead. This is commonly referred to in the literature as the "left" current $I_L$, and we apologize for not clarifying this explicitly in the text (now amended). When the dot is initially empty, it is rapidly filled by the left lead at short times at timescales of order $~1/\Gamma$, until it reaches its steady-state population---at the parameters used here, this means a half-occupied state. At that point, transport is partially blocked, and the steady state current is reached.

It's worth pointing out that a complementary "right" current $I_R$ could be calculated by considering the time derivative of populations in the right lead. The right current is suppressed at short times for this initial condition, and the "average" current $(I_L + I_R)/2$ therefore does not feature the sudden rise and fall at short times.

Why do the effects of the noise seem to be stronger if we start from a magnetized initial condition?

The answer to this is related to the previous question: due to the high peak at short times, the scale of the top panel is higher. This makes it appear as if the effects of noise are different, though they are actually quite similar in practice. We now comment on this in the text.

5) At page 11 the authors write "Therefore, we expect that more robust schemes for constructing the transfer tensor will be needed to enable high-precision dynamical applications". However, from their plots it seems that the smoothing is really effective in stabilizing the results and greatly improves the precision. Do the authors expect this technique to be insufficient when considering even larger cutoffs or in different parameters regimes? I suggest to add a comment on this issue to make more clear what is the range of effectiveness of their approach.

The results we found indeed look promising, and smoothing appears to be a simple and practical way to deal with (weak) noise. However, while smoothing effectively removes the noise and stabilizes the dynamics, it also introduces a systematic bias. The behavior of this bias, and the relative merits of alternative techniques for robustly constructing transfer tensors from noisy data, is something that remains to be explored. We were simply stating that this needs to be done if one is interested in fully controlled results with rigorous data analysis.

Some typos:

• page 10, cutoof -> cutoff;
• fig. 3 the cutoff times should be 0.4 $\Gamma$, 0.8 $\Gamma$ and 1.2 $\Gamma$;

We have fixed these. Many thanks to the reviewer for pointing them out.

* Clarify role of paper in end of introduction and beginning of section 2.
* Rephrasing of eqs. (2)-(4).
* Correction to eq. (7).
* Rephrase and expand the first paragraph of subsection "Constructing the dynamical maps from data" to clarify its role.
* Mention fluctuation theorems as a possible direction for a normalizing condition in same subsection.
* Physical meaning of the norm; physical discussion of figure 3; and implications of weak dependence on counting field and agreement between TT and QMC; are all now explained in more detail in section 5.
* Several additional minor aesthetic changes and corrections.