Introduction-Conventional linear-response transport probes how the expectation values of densities and currents respond to external perturbations. In general, the point of transport experiments is to shed light on the intrinsic dynamics of a system. However, in many cases, theoretically natural questions about (say) the charge or the lifetime of a system's elementary excitations cannot be resolved at the level of linear response. To settle these questions, one must probe higher-order correlation functions or depart from equilibrium, using techniques such as nonlinear response [1][2][3][4][5][6][7][8][9][10] or noise spectroscopy [11][12][13][14].

An especially informative way of capturing these higherorder correlations is full counting statistics (FCS), i.e., the full distribution function of charge transferred across a cut. FCS was initially introduced in the context of mesoscopic systems [15][16][17][18][19][20][21], but contemporary experimental platforms built on synthetic matter, such as ultracold atomic gases or superconducting qubit arrays, allow for much more direct access to this observable. For instance, one can take simultaneous snapshots of all the particles in a quantum gas microscope experiment [22][23][24][25][26][27][28][29][30][31], which, in turn, yield the FCS. While early studies of FCS were focused on lowtemperature properties, recent work has explored how the FCS evolves after quenches [32][33][34][35][36][37][38][39]; a notable finding is that in integrable models as well as certain constrained nonintegrable models, the FCS remains strongly non-Gaussian at all times [38][39][40][41][42][43][44][45]. Many basic questions remain unanswered, however, including the conditions under which the FCS of quantum fluctuations has an effective classical description.

So far, experiments probing FCS in synthetic matter have been carried out by initializing a system in a state with a sharp number of particles in the subsystem of interest [39], letting the system evolve, and then extracting the particlenumber distribution from many snapshots. Despite its straightforwardness, this technique has intrinsic limitations. First, it is limited to initial states with a sharp number of particles in the subsystem of interest, a restriction that eliminates most thermal equilibrium states. Second, in the presence of noise, its performance scales exponentially poorly with subsystem size, since the particle number count is vulnerable to decoherence deep inside the subsystem. Third, the FCS involves arbitrary-order correlation functions, which are difficult to retain faithfully in approximate numerical algorithms based on tensor networks.

In this Letter, we propose and explore a method that circumvents these obstacles. This method is essentially a concrete realization of the thought experiment in Ref. [16], an implementation of which has been proposed earlier in the context of mesoscopic transport [46]. Here, we exploit a further simplification due to the brickwork circuit structure. Our protocol involves an ancilla that acts as a turnstile: the net flux of particles across the turnstile is imprinted on the phase of the ancilla, and the FCS can be read off from this phase information. By tracking the flux instead of the full particle number, this approach evidently overcomes the drawbacks of the direct-counting approach: it is robust against decoherence deep inside the subsystem, and also does not rely on using number-definite initial states. While interferometry using ancillas is a familiar technique for measuring nontrivial many-body observables [47], including FCS [48], our proposal is distinctive in using only local two-site gates, and it can therefore be realized with minimal overhead in near-term experiments. This feature also makes it practical to implement in approximate numerical algorithms: mapping the FCS onto a single-site expectation value opens the door to its computation using tensornetwork methods, such as the density matrix truncation algorithm [49], which are designed to preserve few-site expectation values but not high-order correlations.

Protocol-We consider the brickwork unitary circuit geometry illustrated in Fig. 1(a) acting on a chain of N qubits with conserved particle number (interpreting the states j0i and j1i as empty and occupied, respectively). We are interested in particle transport and fluctuations between two halves of the system, starting from an arbitrary (possibly non-number-sharp) initial state. One can regard this, alternatively, as a Trotterized nearest-neighbor Hamiltonian but the unitary gates can have arbitrary time dependence. Each layer of gates is bracketed by two "turnstile" gates involving the central qubit (labeled in red) and an ancilla, which is initialized in the state jþi ≡ 2 -1=2 ðj0i þ j1iÞ. The turnstile gates apply rotations by an angle AEλ about the z axis on the ancilla, controlled on the central qubit being in the state j1i. The intuition behind the protocol is as follows. The net effect of the turnstile gates before and after the layers depicted in Fig. 1(a) is to imprint a phase λð-λÞ if the occupation of the central site went down (up) during that cycle. Crucially, because of the brickwork geometry, any particle that hops on (off) the central site during this layer must have come from (or gone to) sites above it. During the next layer, likewise, any change in the occupation of the central site corresponds to transport between the central site and those below it. To measure the current, it suffices for the ancilla to keep track of the change in the occupation of the central site. A mathematical description of this protocol is detailed in Sec. SI of the Supplemental Material (SM) [50].

Closed-system numerics-Having established the theoretical framework behind our ancilla-assisted measurement of FCS, we now showcase this method in the context of a particular model, namely, the one-dimensional XXZ chain described by the Hamiltonian H ¼ -

, where S ν ¼ σ ν =2 are spin-1=2 operators. We will focus on spin transport in the isotropic case of Δ ¼ 1, whereupon H reduces to the Heisenberg (or XXX) model. While the XXZ chain is known to be integrable-and its eigenstates and spectrum can be solved for exactly using the coordinate Bethe ansatzcalculating the dynamical correlation functions of this model at finite temperatures is a nontrivial problem. Recently, it was shown that a particular Trotterized version of this model preserves integrability [52,53], which allows us to faithfully study XXZ dynamics via discrete time evolution. In this formulation, the discrete-time Floquet dynamics are governed by the two-step propagator

2kþ1 ðθ; ϕÞ, where the halfstep U even (U odd ) updates all odd-even (even-odd) numbered pairs of spins with the unitary

We choose θ ¼ 0.4π and ϕ ¼ 0.8π, which accordingly sets Δ ¼ sinðϕ=2Þ= sinðθÞ ¼ 1. Of course, our turnstile-based method is completely general and can be used to evaluate the FCS for any value of Δ, including the XY limit (Δ ¼ 0). Adding in the ancilla qubit, we sandwich controlled phase gates with alternating angles AEλ in between the U even and U odd layers.

As mentioned earlier, one of the central advantages of our method is that allows us to study spin transport and FCS starting with initial states that are not necessarily number sharp. Specifically, we consider a family of states with a domain wall between the left and right halves of the system, which are described by the density matrix ρðt ¼ 0Þ ∝ ðe μS z Þ ⊗N=2 ⊗ ðe -μS z Þ ⊗N=2 , with the height of the domain wall given by tanh μ. As μ → ∞, all the spins in the left (right) half point up (down), and the system is furthest from equilibrium. On the other hand, at μ ¼ 0, the initial state is the equilibrium thermal state at infinite temperature, and the SU(2) symmetry of the Heisenberg model is restored. This limit is particularly interesting from the viewpoint of quantum transport as the corresponding dynamics have been conjectured to be in the celebrated Kardar-Parisi-Zhang universality class [54][55][56][57][58][59][60][61].

In Fig. 2, we plot the generating function fðλÞ for μ ¼ 0 obtained using the ancilla; we fit this function to a low-order polynomial near the origin in order to compute its derivatives. This approach is especially well-suited to extract the mean of the total transferred magnetization, hQðtÞi, defined here as the net number of spin excitations that have crossed between the two halves of the chain by time t, QðtÞ ¼ P i>N=2 ½S z i ðtÞ -S z i ð0Þ -P i≤N=2 ½S z i ðtÞ -S z i ð0Þ as well as its variance hðQðtÞ -hQðtÞiÞ 2 i, which we can calculate up to t ¼ 40 using the time-evolving block decimation (TEBD) algorithm with a truncation error < 10 -7 .

Adding noise-Another key feature of the ancilla method is that it remains capable of measuring FCS in the presence of noise. By contrast, direct microscopy is vulnerable to any noise that violates total particle number conservation, since microscopy cannot distinguish changes in particle number due to transport across the midpoint from those due to number-changing noise deep inside the subsystem. In some cases, e.g., with amplitude-damping noise, the effects of noise can be detected and eliminated through postselection (albeit at the cost of exponential overhead in the subsystem size). In other cases, such as depolarizing noise, the noise cannot be corrected for even using postselection. As a result, for large subsystems in the presence of noise, microscopy can measure the FCS only on timescales when essentially no noise has acted anywhere in the subsystem.

The turnstile method is insensitive, by construction, to noise occurring deep in the subsystems. Hence, it faithfully captures the dynamics of charge transfer in noisy systems. When the noise is strong enough, the dynamics of charge transfer are themselves affected by the noise; the turnstile method captures this physical crossover between unitary and noisy dynamics as we now show.

To demonstrate this point, we study transport in the Heisenberg model, in the presence of noise, for the mean, the variance, and the skewness hðQðtÞ -hQðtÞiÞ 3 i as a function of the circuit depth. The third and higher moments can be obtained from the subdominant terms in the expansions for hXi, hYi, and consequently, are very sensitive to the details of the fitting polynomial. It is therefore advantageous to instead construct the full probability distribution p Q by examining λ ∈ ½0; π, which avoids any fitting errors, with the trade-off being that larger counting fields require a higher bond dimension to simulate.

Here, we focus on depolarizing noise, described by the quantum channel ρ → P 3 p¼0 K p ρK † p , where K p are the Kraus operators

Our results for the noisy simulations are shown in Fig. 3. On timescales that are short compared to 1=γ, we reproduce the noiseless behavior, with the mean and variance scaling as t 2=3 [62], seemingly in consistency with the Kardar-Parisi-Zhang prediction. On longer timescales, we find that the mean saturates and the variance begins to grow linearly. This is a consequence of the noise breaking the conservation law: on long timescales, the occupation of the central site fluctuates randomly, so the phase accumulation is random. Thus, its mean stops changing and its variance grows linearly in time. To the best of our knowledge, the turnstile method is the only way to access this crossover; while our predictions are converged in truncation error, the fact that they reproduce theoretical expectations in the large-noise limit is a helpful sanity check.

Besides the mean and the variance, the turnstile protocol also naturally facilitates access to the higher moments, which have been shown to be essential for determining the dynamical universality class in recent experiments on a superconducting quantum processor [63]. Interestingly, Fig. 3(c) indicates that the skewness is much less affected by the noise than the first two moments. These conclusions continue to hold for other types of noise, such as amplitude damping, as detailed in the SM.

Finally, we remark on the effects of noise on the ancilla itself. In general, we can use the ancilla to measure FCS until we reach the T 2 time of the ancilla. In an optimally calibrated protocol, λt would be chosen on a uniform grid from 0 to 1=Q typ , where Q typ is the typical scale of charge transfer over timescale t ≤ T 2 . Estimating the number of shots needed to get a precise histogram of charge transfer is a standard exercise in metrology [64]: this quantity scales as NðλÞ ∼ min½1=λ 2 ; ðt=T 2 Þ 2 =λ 4 . This rapid increase in the number of needed shots limits the resolution with which QðtÞ can be reconstructed at late times. The sensitivity limits of FCS measurements, as well as ways to improve their sensitivity using ideas from error correction [65], will be explored in future work.

Application to approximate algorithms-In addition to its experimental relevance, the turnstile approach has the benefit of being compatible with approximate quantum simulation algorithms. For purely unitary dynamics, the turnstile approach is equivalent to the introduction of a fictitious counting field [16,17,19,[66][67][68][69][70] that modifies the forward and backward evolution in different ways so that the evolution of a systems density operator is no longer given by a quantum channel. The turnstile method avoids this-the evolution of the full system (including the ancilla) can be treated as a quantum channel-making it compatible with quantum simulation algorithms that execute a quantum channel by construction.

An example of such an algorithm, and the one that we will focus our attention on in this section, is the density matrix truncation (DMT) algorithm [49]. By protecting carefully chosen components during a singular-value decomposition and truncation step, DMT preserves the local structure of the system's density matrix-after a full truncation sweep over every bond, the algorithm preserves all reduced densities matrices on three neighboring sites [71]. However, by treating the ancilla as a bona fide site in the chain, three-site operators on the physical sites in the neighborhood of the ancilla are no longer preserved. We can restore the preservation of all physical three-site operators by combining the ancilla and a physical site into a super site and modifying DMT to accommodate the enlarged site.

DMT performs similarly to other methods near equilibrium, but has shown a clear advantage far from equilibrium [49,72]. Combining the turnstile method and the DMT algorithm should enable more efficient evaluation of nonequilibrium FCS, closing the gap to ultracold atomic and trapped ion experiments where initial states are typically far from equilibrium [26,28,39,[73][74][75][76][77][78][79][80][81]. In the rest of this section, we will study the nonequilibrium charge transfer statistics of a U(1)-charge-conserving random circuit on a spin-1=2 chain, using both DMT and its closest competitor, TEBD. We will focus on a semipolarized Néel state as an initial state, from which the system undergoes a global quench. This initial state is a (mixed) product state defined on even and odd sites by ρ even=odd ∝ expðAEμσ z Þ, respectively. In the following, we choose to work with polarization p ≡ tanhðμÞ ¼ 0.9. Both the global quench and the rapid entanglement generation in random-circuit evolution make this environment very challenging for existing numerical tools. We simulate the system, including the ancilla, up to time t ¼ 40 with a selection of counting fields and bond dimensions (d ¼ 100, 200, 400, 800, 1200) for a single circuit realization. Doing so for both TEBD and DMT, we extract the variance of the charge transfer using a low-order polynomial fitting of the cumulant generating function, χðλÞ ¼ log fðλÞ, near the origin. The (scaled) variance is shown in Fig. 4. With regular TEBD, the data quickly becomes unconverged in bond dimension (by times t ∼ 10). DMT significantly outperforms this, remaining converged for much later times (t ∼ 26). We have verified that this separation in performance is generic for random circuits by repeating our analysis for multiple circuit realizations, as FIG. 4. The variance of the charge transfer, as a function of time, across the central bond of a random unitary circuit prepared in the semipolarized (p ¼ 0.9) Néel state. Both TEBD (dashed) and DMT (solid) algorithms are used with bond dimensions up to d ¼ 1200. The FCS is expected to be asymptotically described by a symmetric exclusion process (SSEP) (black dashed).

well as for a different nonequilibrium state (a semipolarized domain-wall state), which can be found in the SM [82].

In addition to improved convergence characteristics, the DMT results are consistent with the FCS in a symmetric simple exclusion process, which is expected to describe the asymptotic FCS of U(1)-charge-conserving random unitary circuits [36]. The diffusivity D ¼ 1 and the susceptibility χ are the same for a symmetric exclusion process and the random circuits considered here. Therefore, the dynamics of charge, including at the level of fluctuations, is given by the same fluctuating hydrodynamic theory [83]. It follows that a symmetric exclusion process and the random circuit considered here should share the same long-time charge transfer statistics. This is borne out by Fig. 4, which overlays the charge transfer variance for the exclusion process initialized with a semipolarized Néel state over the random circuit results.

We now briefly compare the performance of DMT using the turnstile with that of other numerical approaches (beyond the comparison with mixed-state TEBD above). For pure-state evolution, DMT can be compared with either direct matrix-vector multiplication or pure-state TEBD [84]. For comparable resources, those methods are accurate until times in the range t ≤ 24-32. For mixed states, exact density-matrix evolution and TEBD with quantum trajectories are limited to much shorter times t ≤ 16 (in the latter case, because of sampling overhead) [85]. In the equilibrium limit, as noted above, the performance of DMT is comparable to that of mixed-state TEBD. Thus, in the purestate and equilibrium limits, DMT is competitive with algorithms that exploit the special structure present in those limits. Away from those limits-for evolving generic nonequilibrium mixed states as we have done here-DMT outperforms these alternative methods by approximately a factor of 2.

Discussion and outlook-To date, computing FCS has been a challenging task-both experimentally and theoretically-as it requires knowledge of arbitrary-order correlations. In this Letter, we overcome these obstacles and present a new method for measuring FCS, based on a quantum turnstile, that is robust against noise and only involves local operators, thus allowing for its easy implementation in near-term experimental platforms. We showed the versatility of this approach in the context of studying the FCS of magnetization transfer in a Heisenberg spin chain. Our results in this regard directly apply to the recent experiments on superconducting qubit arrays [63], which rely on the higher moments, such as the skewness and kurtosis, to determine the dynamical universality class. We also illustrated how the quantum turnstile protocol can be paired with tensor-network algorithms to study numerically demanding transport problems such as charge transfer in random circuits. Our results pave the way for using FCS to resolve the underlying dynamics of a range of many-body systems that are realizable in noisy quantum simulators (e.g., disordered systems). The turnstile method can be further adapted to explore higher-order temporal fluctuations of more general (i.e., nonconserved) operators in many-body systems; the nature of these fluctuations is largely an open question, although exact predictions have been derived for certain solvable models [86].