From Non-Hermitian Linear Response to Dynamical Correlations and Fluctuation-Dissipation Relations in Quantum Many-Body Systems

Quantum many-body systems are characterized by their correlations. While equal-time correlators and unequal-time commutators between operators are standard observables, the direct access to unequal-time anti-commutators poses a formidable experimental challenge. Here, we propose a general technique for measuring unequal-time anti-commutators using the linear response of a system to a non-Hermitian perturbation. We illustrate the protocol at the example of a Bose-Hubbard model, where the approach to thermal equilibrium in a closed quantum system can be tracked by measuring both sides of the fluctuation-dissipation relation. We relate the scheme to the quantum Zeno effect and weak measurements, and illustrate possible implementations at the example of a cold-atom system. Our proposal provides a way of characterizing dynamical correlations in quantum many-body systems with potential applications in understanding strongly correlated matter as well as for novel quantum technologies.

Notwithstanding its fundamental importance, both sides of the FDR have thus far only been measured for classical systems [36][37][38]. For quantum systems, only one side, the unequal-time commutator, is easily accessible thanks to Kubo's celebrated linear response theory [4,21], which has been extensively used to characterize quantum systems out of equilibrium [23][24][25][26][27][28][29][30]. A * kevinthomas.geier@unitn.it main difficulty regarding the measurement of dynamical correlations stems from the fact that a projective von Neumann measurement at a particular time collapses the quantum state [39], which prevents an unperturbed measurement at a later time and thus hinders a measurement of the time-time correlation with respect to the initial state. Various pioneering proposals for measuring unequal-time correlations on various platforms exist [30,[40][41][42][43][44][45][46][47][48], but attempts to overcome the inherent difficulties of such a measurement are often specific to certain setups or apply only to a limited set of observables. As of today, an experimental observation of the unequal-time anti-commutator in a quantum many-body system remains elusive.
Here, we discuss how a linear response to a non-Hermitian perturbation [49,50] permits direct experimental observation of the unequal-time anti-commutator. Combined with a traditional method for measuring the corresponding unequal-time commutator, e.g., standard linear response, this scheme gives access to both sides of the FDR independently, allowing one to track a system's evolution towards thermal equilibrium. We illustrate this possibility by means of numerical simulations at an example motivated by a ground-breaking cold-atom experiment [51] -a Bose-Hubbard system that is quenched from a Mott-insulating initial state to the superfluid phase (see Fig. 1). This analysis provides a blueprint for revealing the FDR using experimental abilities that are common in state-of-the-art engineered quantum systems.
The key to measuring unequal-time anti-commutators is the ability to engineer (effective) non-Hermitian perturbations. In recent years, a tremendous interest in non-Hermitian physics has emerged [52,53], stimulated by the rapid progress in the experimental generation and control of non-Hermitian systems [54][55][56][57][58][59]. Indeed, non-Hermiticity gives rise to a wealth of new physics with arXiv:2104.03983v2 [cond-mat.quant-gas] 28  (b) Thermalization dynamics of the dissipative part of the "Hermitian" dynamic susceptibility χ nn (tw, ω) ("commutator") and the reactive part of the "non-Hermitian" dynamic susceptibility χ (NH) nn (tw, ω) ("anti-commutator"). (c) Dynamic susceptibilities, rescaled according to the FDR (13) at early and late waiting times. The effective temperatures k B T / J = {4.5, 4.2} for Jtw = {0.1, 10}, respectively, are determined by Eq. (16) using the least-squares method. The FDR is clearly violated at early times, but it is restored at late times when the system has thermalized. novel (topological) phases and unconventional critical behavior [60][61][62][63][64][65][66], bearing a vast potential for applications, e.g., in strongly enhanced quantum sensing [67,68] or adiabatic quantum optimization [69,70]. Leveraging on this development, we design a specific protocol to generate effective non-Hermitian dynamics in a system of cold atoms, enabling access to the fluctuation side of the FDR (i.e., to the unequal-time anti-commutator). Our scheme is most conveniently phrased as an application of the quantum Zeno effect [71,72], combining outcoupling to an ancillary system with a projection on the Zeno subspace given by the empty ancilla. In a coldatom implementation, this can be realized through a co-herent or dissipative perturbation in the linear regime, together with the ability of distinguishing zero from nonvanishing ancilla population in post-selection. While a single step in the Zeno evolution yields the unequal-time anti-commutator in time domain, an extended Zeno evolution, which we propose to implement harnessing engineered dissipation [73,74], allows one to probe frequencyresolved responses in the same way as in standard linear response experiments. To demonstrate the feasibility of our proposal, we benchmark our protocol by numerically solving the full quantum evolution, including the stochastic dynamics underlying the dissipative scheme, and discuss experimental error sources. We also examine formal relations to dissipative quantum systems, where non-Hermitian dynamics can be generated by postselecting individual quantum trajectories on the absence of quantum jumps [58,66,[75][76][77], and establish general cross-connections between (non-)Hermitian linear response and ancilla-based weak measurements of dynamical correlations [43,44] (see Ref. [78] for a comprehensive review on weak measurements). Our proposed realization of non-Hermitian linear response is feasible even when existing weak measurement protocols are difficult to engineer experimentally, and it excels in regimes where projective protocols fail as a consequence of their restriction to observables with two eigenvalues [30,41,[43][44][45], as we demonstrate through numerical benchmarks. Our approach thus opens the door to probing the FDR in quantum many-body systems in an unbiased way and for a broad range of observables.

II. THE FLUCTUATION-DISSIPATION RELATION
For a quantum many-body system in thermal equilibrium, the FDR [21] links the symmetrized correlation spectrum S BA (ω) of any two operators A and B across the entire frequency spectrum ω to the dissipative part of the dynamic susceptibility χ BA (ω) via where is the reduced Planck constant and k B is the Boltzmann constant. This elegant relation requires only a single parameter as input, the global temperature T . The ease of accessing χ BA (ω) can then be exploited to obtain S BA (ω). However, when the system is far from equilibrium, the two sides of Eq. (1) become non-stationary and the FDR can be broken [8,29], making it necessary to devise independent handles on both sides of the relation, as has been proposed in Ref. [30]. One can generalize the definitions of S and χ to such a non-equilibrium situation by introducing the response function and the symmetrized dynamic correlation function defined, respectively, in terms of the unequal-time commutator and anti-commutator of the Heisenberg operators A(t) and B(t). Here, θ(t) is the Heaviside step function, ensuring causality of the response, and the subscript in the expectation value · · · 0 signifies that the Heisenberg operators evolve under the (unperturbed) Hamiltonian H 0 . In the context of non-equilibrium quantum field theory, Eqs. (2) and (3) are also known as the spectral function ρ and the statistical function F , respectively [79]. Equation (2) is the non-equilibrium version of Kubo's well-known linear response function [21], which determines the evolution of the expectation value B(t) under the perturbed Hamiltonian H(t) = H 0 + H 1 (t), to linear order in the perturbation H 1 (t) = −f (t)A, according to In contrast to the usual equilibrium linear response scenario, the initial state is not necessarily stationary with respect to H 0 . In this situation, it is common to define the non-equilibrium generalization of the dynamic susceptibility χ as the Fourier transform of the response function φ with respect to the relative time ∆t = t − t at fixed central time τ = (t + t )/2 [79], This quantity is commonly decomposed as χ BA = χ BA + iχ BA into a reactive part χ BA (ω) = [χ BA (ω) + χ AB (−ω)]/2 and a dissipative (or absorptive) part χ BA (ω) = [χ BA (ω) − χ AB (−ω)]/2i [3], the latter entering the right-hand side of the FDR (1). The correlation spectrum S(τ, ω) on its left-hand side can be defined analogously to Eq. (5) as the Fourier transform of Eq. (3). For a thermalizing system, we expect χ BA (τ → ∞, ω) and S BA (τ → ∞, ω) to reach steady values that fulfill the FDR. As such, the restoration of the FDR provides an excellent probe for how and when a quantum manybody system approaches thermal equilibrium [30]. On top of that, the FDR yields the effective temperature at which the system thermalizes [23][24][25][26][27][28][29][30]. Remarkably, this independent way of defining temperature does not require any a priori assumptions other than the FDR. In our numerical benchmarks, we find good agreement between the effective temperature extracted from the FDR and the expected temperature of a thermal ensemble at the equivalent energy density (see Section IV).
While the commutator in Eq. (2) can be accessed rather straightforwardly, for example, by studying how energy is absorbed or how an observable deviates from its equilibrium value following a time-dependent perturbation [21], the determination of the unequal-time anticommutator in Eq. (3) is, unfortunately, considerably more challenging. We now employ a recent extension of linear response theory to non-Hermitian Hamiltonians [49,50] as a general way of gaining access to the left-hand side of Eq. (1), which enables direct probes of the FDR.

III. NON-HERMITIAN LINEAR RESPONSE THEORY
Though long established in the context of open quantum systems [80,81], recent years have seen a surge of interest in quantum systems with non-Hermitian Hamiltonians [53]. Here, we tap into this development by exploiting the linear response to a non-Hermitian perturbation [49,50] in order to extract unequal-time anticommutators.
In contrast to usual linear response theory, we assume that the system is effectively described by a non-Hermitian Hamiltonian H(t) = H 0 + H 1 (t), where H 0 is the unperturbed (Hermitian) Hamiltonian and H 1 (t) = −if (t)A is an anti-Hermitian perturbation with a positive semi-definite operator A and a non-negative timedependent function f (t). For example, such a scenario arises in the quantum trajectories approach to dissipative quantum systems [75][76][77] if the evolution is conditioned on the absence of quantum jumps [56,58,66] (see also Section VI B). In addition, we show in Section VI A that existing ancilla-based weak measurement protocols for the unequal-time anti-commutator [43,44] can also be rephrased in the framework of non-Hermitian linear response. In Section V, we present a scheme for engineering effective non-Hermitian Hamiltonians based on the quantum Zeno effect to probe such responses, even frequency-resolved, for a wide range of observables.
A quantum state described by the density operator ρ(t) evolves in time under the non-Hermitian Hamiltonian H(t) according to the von Neumann equation (6) with initial condition ρ(0) = ρ 0 . Using time-dependent perturbation theory, a straightforward calculation (reported in Appendix A) shows that, to linear order in the perturbation, the unnormalized expecation value of a (Hermitian) observable B is given by The non-Hermiticity of the perturbed Hamiltonian has the important consequence that the state ρ(t) is no longer normalized: as can be seen by inserting the identity operator for B in Eq. (7), its norm decreases with time, to linear order, as Physically, this decrease can be interpreted as the leakage of the wave function into a complementary state space (see also Section V). To account for this loss of probability, we consider the normalized expectation value , describing a post-selected measurement [78]. Combining Eqs. (7) and (8), the disconnected correlations drop out to linear order, and we can write the response as with the "non-Hermitian" response function Here, we insert the Heaviside step function θ to ensure causality of the response. Remarkably, the non-Hermitian response function in Eq. (10) is the desired measurable quantity that gives direct access to the unequal-time anti-commutator (3) by virtue of the relation φ (NH) To establish a link between the response function (10) and the correlation spectrum appearing on the left-hand side of the FDR (1), we define the "non-Hermitian" dynamic susceptibility, similarly to Eq. (5), as the Fourier transform We can split this quantity as χ into the components (for conciseness, we remove the τ argument from the following formulas) which we refer to, in analogy to their Hermitian counterparts, as the reactive and dissipative parts of the non-Hermitian susceptibility, respectively. As shown in Appendix A, the reactive part, Eq. (12a), gives access to the correlation spectrum via the identity . This allows us to rewrite the FDR (1) in thermal equilibrium as which is expressed entirely in terms of the susceptibilities χ (NH) BA and χ BA , accessible using non-Hermitian and standard (Hermitian) linear response, respectively. As such, linear response theory provides an elegant and general framework for independently probing both sides of the FDR (13) out of equilibrium, which works for arbitrary observables in any quantum many-body system. Compared to projective protocols for measuring unequal-time anti-commutators [30,41,[43][44][45], which are restricted to observables with two eigenvalues (see also discussion in Section VI B 2), one of the main assets of linear response theory is its broad applicability.
It is worthwhile emphasizing that the outlined derivation of the response to a non-Hermitian perturbation is by no means restricted to the linear regime only, but, as well known in standard response theory [20], can be extended to non-linear responses. In fact, by calculating the expansions in Eqs. (7) and (8) to the desired nonlinear order, one can in principle access an infinite hierarchy of unequal-time correlations, order by order. By perturbing the system at multiple sites simultaneously, non-linear responses could therefore also enable access to (global) many-body operators which are expected not to thermalize and consequently violate the FDR.
Approaching the problem of measuring dynamical correlations from the (non-)Hermitian linear response perspective turns out to be fruitful for a number of reasons. For one, non-Hermitian linear response is completely agnostic to the way the non-Hermitian perturbation is implemented and therefore directly benefits from any advancements in the field of non-Hermitian physics regarding the generation and control of non-Hermitian Hamiltonians. Furthermore, it provides an ancilla-free interpretation of common ancilla-based weak measurement schemes for the unequal-time anti-commutator [43,44]. So far, it has not been clear whether ancilla-free formulations of such protocols allow for a meaningful physical interpretation [44], but, as we show in Section VI A, this is indeed possible in the light of non-Hermitian linear response. Conversely, any non-Hermitian perturbation can in principle be realized with the help of an ancilla using only unitary evolution and standard projective measurements, although the required couplings may not always be straightforward to implement. In Section V, we present specific ancilla-based schemes with experimentally feasible system-ancilla couplings, providing access to dynamical correlations for a wide range of observables. Finally, from a linear response point of view, it is natural to study responses to periodic perturbations that directly give access to frequency-resolved susceptibilities. As explained in Sections V B and V C, this becomes practical within our framework also for non-Hermitian perturbations by exploiting the quantum Zeno effect.

IV. ILLUSTRATION: QUENCH IN A BOSE-HUBBARD SYSTEM
In this section, we demonstrate how (non-)Hermitian linear response allows one to access both sides of the FDR (13) independently. Such measurements can be used to either probe thermalization or the absence thereof [30]. If a system of interest is coupled to a large thermal bath, it will sooner or later always approach thermal equilibrium with the bath temperature [81], and the FDR will eventually hold. In contrast, the question whether and how a closed quantum system thermalizes once it is brought out of equilibrium is much more subtle. Remarkably, an isolated system can act as its own bath [16,17]: a thermalizing subsystem behaves after long times as if it was in thermal equilibrium with the rest of the system at an effective temperature set by the initial state. According to the eigenstate thermalization hypothesis [14][15][16][17][18], this process occurs on the level of individual eigenstates. Although the precise conditions for its validity are not yet entirely understood, it is believed that (eigenstate) thermalization generally holds for generic states of interacting quantum many-body systems in the bulk of the spectrum. Important scenarios where thermalization fails (with concomitant violation of the FDR) include integrable models [17,26], many-body localization [82,83], as well as Hilbert space fragmentation and the related phenomenon of quantum many-body scars [84][85][86]. On top of that, breaking FDRs is a characteristic signature of far-from-equilibrium systems near a non-thermal fixed point [8,29]. All of these settings represent promising targets for our (non-)Hermitian linear response scheme to reveal either the validity or the breakdown of the FDR.
For illustrative purposes, we focus here on the generic case where the system does thermalize and the FDR is expected to hold. In ground-breaking cold-atom experiments, it has been shown that even in very small interacting quantum systems, expectation values can reach steady states that are consistent with thermal equilibrium [51]. We now illustrate how such an experiment could go one step further by demonstrating the validity of the FDR. To this end, we numerically solve the full quantum evolution for the Bose-Hubbard chain describing the experiment in Ref. [51] (we emphasize that our approach does not depend on such a model choice and can be applied to general quantum systems). The Bose-Hubbard Hamiltonian is given by Here, the optical lattice sites are denoted by = 1 . . . L with associated bosonic annihilation, creation, and number (density) operators a , a † , and n , respectively. J is the strength of the nearest-neighbor hopping, for which we assume periodic boundary conditions, and U the onsite interaction rate. In our numerics, we do not truncate the local Hilbert-space dimension and employ an adaptive Krylov subspace method for time evolution [87][88][89].
While previous numerical studies of FDRs in this model have focused on density autocorrelations at large U/J and low fillings [30], here we consider quenches into the superfluid regime (U/J ∼ 1) at unit filling and also explore off-site correlations as a function of distance. The linear response protocol is illustrated at the top of Fig. 1. We initialize the system of N = L bosons in a Mott-insulating state at U/J → ∞ and then quench it at time t = 0 into the superfluid phase at U/J = 1.5625 [90], chosen consistent with the experiment in Ref. [51]. This quench throws the system heavily out of equilibrium. After a variable waiting time t w , we either apply a Hermitian or an anti-Hermitian perturbation in order access the desired response functions in Eqs. (2) or (10), respectively. The perturbation is applied as a rectangular pulse of strength s and duration δt, The exact shape of the pulse is unimportant as long as the pulse duration is sufficiently short compared to the characteristic time scales of the system (cf. Appendix C). In this case, the pulse can be approximated by a δ function as f (t) ≈ sδ(t − t w ). Figure 1a shows the time trace of the response to a (non-)Hermitian perturbation giving access to density autocorrelations (B = A = n), computed in a system of L = 12 sites for a perturbation of strength s = 0.05 and duration Jδt = 0.01 [91]. The thermalization dynamics of the corresponding dynamic susceptibilities χ (NH) and χ is depicted in Fig. 1b. For the purposes of this section, we evaluate the susceptibilities at fixed waiting time t w [29], i.e., (15) using an exponential filter of characteristic frequency γ/J = 0.2 to ensure convergence of the Fourier integrals (see Appendix B 1 for technical details). χ (NH) is symmetric in ω and grows as a broad central peak with small wings of opposite sign that gradually disappear as t w increases, while χ is anti-symmetric and develops characteristic peaks around non-zero frequencies.
To assess whether the two functions satisfy the FDR, we use the least-squares method to find the effective temperature T that best relates the susceptibilities via the FDR (13), i.e., for a fixed waiting time t w , In Fig. 1c, one can see that the FDR is clearly violated at early times, i.e., there exists no global value of T such that Eq. (13) holds, but at later times, the agreement is remarkable and supports the interpretation that the system undergoes thermalization.
A distinct feature of the FDR in equilibrium is that it holds for any pair of observables A and B. To confirm this prediction for our model system, we have computed χ (NH) and χ for off-site density correlations corresponding to A = n and B = n +d as a function of the distance d. The results are shown in Fig. 2a, where χ (NH) is rescaled according to Eq. (13) with the best-fitting effective temperature T obtained from Eq. (16) for each configuration (t w , d). Qualitatively, it can be seen that the two quantities deviate for early times, but agree well for late times. To make this statement more quantitative, in Figs. 2c and 2d, we show the relative and abso-  n n +d , rescaled according to the FDR (13), for different waiting times tw as a function of the spatial distance d. The effective temperature T is determined for each configuration according to Eq. (16) using the leastsquares method. At early times, clear deviations are visible, but for late times, the two quantities agree and the FDR (13) is fulfilled. (b) Least-squares value of the effective temperature T , (c) relative error, and (d) absolute error of the FDR as a function of the waiting time tw for several distances d. At small distances, after times on the order of J −1 , the effective temperature approaches a constant value consistent with the prediction H0 T = E0 for a thermal state (gray dashed line), and the relative error becomes small. As the distance increases, the relative error grows, but the absolute deviation becomes small (see also Appendix B).
lute error of the FDR, i.e., the L 2 norm of the difference between the left-and right-hand side of Eq. (13) (see Appendix B 1), as a figure of merit measuring how well the FDR (13) is fulfilled at a particular instance of time. For small distances, the relative error becomes vanishingly small after waiting times on the order of J −1 , while for larger distances, the error tends to drop later and fluc-tuates around a non-zero offset. A similar behavior is exhibited by the effective temperature (see Fig. 2b): at small distances and late times, the temperatures are approximately constant and agree with each other, while this is no longer true for larger distances. Only in the former case, where the relative error is small, the effective temperature can be attributed the physical meaning of the temperature at which the subsystem degrees of freedom thermalize. This temperature is consistent with the one obtained for a thermal state at the equivalent energy density using the condition H 0 T = E 0 (gray dashed line at k B T / J = 4.27 in Fig. 2b, calculated for L = 8 using exact diagonalization), where · · · T denotes the expectation value with respect to a canonical ensemble at temperature T , and E 0 is the energy of the initial state after the quench [51].
While global many-body observables are expected to violate the FDR due to the purity of the global quantum state, one would expect two-site observables like the off-site density correlations shown in Fig. 2 to thermalize and thus satisfy the FDR for sufficiently large systems and late times. Although the absolute error gradually decreases with increasing distance due to the lower signal strength (see Figs. 2a and 2d), relative discrepancies persist even after very long times. In Appendix B, we investigate the behavior of the error as a function of system size and study a similar quench scenario for a two-dimensional (2D) Bose-Hubbard system of 4 × 4 lattice sites and N = 16 bosons (unit filling). Our analysis reveals that the relative error at long waiting times decreases with increasing system size. Furthermore, the larger 2D system exhibits only a minor trend towards larger relative errors as the distance increases, and the FDR is overall better fulfilled than in the smaller onedimensional (1D) chain. This points to the conclusion that the observed discrepancies of the FDR in Fig. 2 for large distances are likely due to finite-size effects. Thus, our numerical results indicate that for sufficiently large systems, off-site density correlations fulfill the FDR even at long distances, confirming the expectation that subsystems consisting of few degrees of freedom thermalize.
Having illustrated how the FDR becomes accessible via (non-)Hermitian linear response, we now turn to the question of how to realize the corresponding non-Hermitian perturbations experimentally.

V. REALIZATION OF NON-HERMITIAN LINEAR RESPONSE
There exists a growing body of work that describes how non-Hermitian physics can be generated in quantum many-body systems [52,53]. Non-Hermitian Hamiltonians naturally arise in the context of dissipative quantum systems [80,81], where they govern the evolution of individual quantum trajectories conditioned on the absence of quantum jumps [75][76][77]. This way, it is possible to harness natural sources of dissipation in order to explore  Figure 3. Realization of an effective non-Hermitian Hamiltonian using the quantum Zeno effect, illustrated for an optical lattice. (a) Coupling a single lattice site to an ancilla gives rise to a perturbation by the density operator A = n at that site. (b) A perturbation by the hopping operator a † 1 a 2 + a † 2 a 1 can be achieved by coupling two sites 1 and 2 simultaneously to an ancilla. (c) Single step in the quantum Zeno evolution. The probability p(0) of detecting no particles in the ancilla gradually decreases over time (red). A measurement of the ancilla population, post-selected on the condition that the ancilla is empty (inset), projects the system on the empty-ancilla subspace. The coupled evolution plus projection corresponds to an effective non-Hermitian perturbation (NHH, black). (d) When the projective measurement is performed frequently as compared to the strength of the coherent coupling Ω, the system plus ancilla is kept in the quantum Zeno regime for a prolonged period of time. The resulting pulsed Zeno evolution (red) is interpolated by the evolution under an effective non-Hermitian Hamiltonian (NHH, black). Alternatively, the repeated measurements can be substituted by strong engineered dissipation on the ancilla. The light gray lines show 20 trajectories corresponding to different realizations of engineered classical dephasing noise ξ(t) on the ancilla, whose ensemble average (gray dashed line) approximates an effective non-Hermitian evolution.
novel non-Hermitian physics [56,58,66]. Over the years, ever better techniques of screening experiments as much as possible from any sources of dissipation have been developed, with the goal of observing clean unitary dynamics in isolated quantum systems. This bears the potential to re-introduce channels of engineered dissipation using specifically designed control schemes.
In this section, we propose an ancilla-based protocol that relies entirely on synthetic sources of dissipation in order to realize an effective non-Hermitian Hamiltonian. The perturbation can selectively be applied as a short pulse or under continuous modulation of its strength, allowing one to probe frequency-dependent responses in the same way as in standard linear response scenarios. Moreover, our flexible and experimentally feasible choice of system-ancilla coupling gives access to a wide range of observables. Figure 3 gives an overview of the scheme, which is most conveniently phrased as an application of the quantum Zeno effect [71,72]. Depending on the desired perturbation operator A, the relevant subsystem, e.g., a single site or two neighboring sites in an optical lattice, is coherently coupled to an initially empty ancilla, as depicted in Figs. 3a and 3b. A measurement of the ancilla population projects the system on the subspace with a definite number of particles in the ancilla. As will become clear further below, non-Hermitian dynamics is realized by post-selecting those measurement outcomes where the ancilla remains empty (see Fig. 3c). Repeating this measurement frequently gives rise to a quantum Zeno effect: as the measurement frequency tends to infinity, the probability of populating the ancilla vanishes. If, instead, the measurement frequency is finite, there is a finite probability of populating the ancilla. As illustrated in Fig. 3d, this leads to a "pulsed" leakage of probability from the subspace where the ancilla is empty to a complementary subspace with non-vanishing ancilla population. Instead of the pulsed Zeno effect, we can also use the continuous Zeno effect [72], which can be realized by substituting the repeated measurements with strong engineered dissipation on the ancilla [73]. This has the advantage of not requiring any non-destructive measurements during the evolution, but only a single post-selected measurement at the final evolution time. Both the pulsed Zeno evolution and the ensemble average over many noise realizations in the continuous case can be described by an effective non-Hermitian Hamiltonian [72,92,93] (see Fig. 3d), which realizes the desired anti-Hermitian perturbation for measuring the unequal-time anti-commutator.
While our scheme can be implemented on various platforms, for the sake of concreteness, we focus here on bosons in optical lattices, where the ancilla may correspond to an auxiliary lattice site or an additional internal state. A crucial experimental requirement is the ability to distinguish an empty ancilla from one with non-zero population, which enables the projection on the emptyancilla Zeno subspace. This requirement is met, for instance, by modern quantum gas microscopes, which reach both single-site and single-particle resolution [94,95].
It is instructive to first consider a single step in the Zeno evolution consisting of a short coupling pulse followed by a projection, as depicted in Fig. 3c. It turns out that this scenario corresponds to applying a δ-like perturbation suitable for measuring the time trace of the non-Hermitian response function (10) like in Section IV. Subsequently, we explain how the quantum Zeno effect enables a prolonged evolution under a non-Hermitian Hamiltonian, focusing on the scenario with strong engineered dephasing noise that induces a continuous Zeno effect (cf. Fig. 3d). We benchmark variants of our scheme for measurements in both time and frequency domain at the example of the Bose-Hubbard chain introduced in Section IV. Bose-Hubbard systems subject to dissipation have been extensively studied with the goal of exploring the rich dynamics of open quantum systems [66,96,97], whereas here, we use engineered dissipation as a tool [74] to probe dynamical correlations in closed systems via non-Hermitian linear response. In Section VI, we compare our approach with other protocols for measuring unequal-time anti-commutators, including ancilla-based weak measurement schemes [43,44] and projective protocols [30,41,[43][44][45], and discuss potential sources of errors as well as strategies on how to mitigate them.

A. Non-Hermitian linear response as a single step in the quantum Zeno evolution
In this subsection, we discuss a single step in the Zeno evolution, which corresponds to applying an effective non-Hermitian δ-like perturbation as in Section IV and allows one to access the unequal-time anti-commutator in Eq. (10) in time domain.

Outline of the scheme
We consider a system-ancilla coupling Hamiltonian of the form where a (a † ) and b (b † ) represent the bosonic annihilation (creation) operators of the system mode to be probed and the ancilla, respectively, and Ω is the coupling strength. In the coupling scheme depicted in Fig. 3a, the operator a represents a single lattice site , giving rise to an effective anti-Hermitian perturbation by the density (number) operator A = n = a † a , as becomes clear below. The scheme in Fig. 3b couples two lattice sites 1 and 2 , which may or may not be nearest neighbors, simultaneously to the ancilla. This corresponds to the replacement a → a 1 + a 2 in Eq. (17) and produces a non-Hermitian perturbation by the operator A = n 1 + n 2 + a † 1 a 2 + a † 2 a 1 . This type of perturbation can therefore be used to access FDRs for the hopping operator a † 1 a 2 + a † 2 a 1 , as we demonstrate below for nearest neighbors. It is possible to consider even more general setups [98], e.g., by adding a relative phase between the two couplings in Fig. 3b using laser-assisted tunneling [99], or by coupling a multitude of sites to one or more ancillas, enabling global perturbations by sums of local operators. The general form of the accessible perturbations is given in Appendix C 2.
A single Zeno step of duration δt corresponds to a unitary evolution described by the time evolution operator U (δt) = exp(−iHδt/ ), followed by a projection on the Zeno subspace defined by the performed measurement [72]. During the coupling, the total Hamiltonian is given by H = H 0 + H cpl , but for sufficiently short δt, it is permissible to neglect the evolution under H 0 (for simplicity, we assume the ancilla to have no internal dynamics). A measurement of the ancilla population projects the system on one of the Zeno subspaces with a fixed number of particles in the ancilla, which can be realized experimentally via post-selection. Prior to the coupling, we require the ancilla to be in the vacuum state. Let P denote the projection operator on the empty-ancilla subspace. Then, during one Zeno step, the state ρ(t w ) at the waiting time t w changes, up to a normalization, as As shown in Appendix C 1, to leading order in the effective coupling strength s = (Ωδt) 2 /2, this process corresponds to the evolution under an effective non-Hermitian Hamiltonian, For the purpose of measuring the non-Hermitian linear response in time domain, a single Zeno step is sufficient.
The system subsequently evolves unitarily under the unperturbed Hamiltonian H 0 up to the final observation time t > t w + δt. The unnormalized expectation value of an observable B is then given by In the linear regime, the probability of detecting no particles in the ancilla after the coupling reads p(0) = 1−2s A(t w ) 0 , which can be found by inserting the identity operator for B in Eq. (20). Normalizing Eq. (20) by this probability yields, to leading order in s, the conditional expectation value representing a post-selected measurement conditioned on the empty ancilla. As anticipated, a comparison with Eqs. (4) and (10) shows that this result effectively corresponds to a linear response after applying the anti- As required in Section III, the perturbation operator A = a † a is indeed positive semi-definite, in line with the physical intuition that the norm of the state can only decrease through outcoupling followed by a projection. It is instructive to compare the result in Eq. (21) with the one obtained if no projection on the empty-ancilla subspace is performed, e.g., if the measurement apparatus is unable to distinguish an empty ancilla from one with non-vanishing population or the result of the ancilla measurement is ignored. In this case, a simple average over all ancilla populations is obtained, where, to leading order in the effective coupling strength s, only single occupancies of the ancilla contribute. The unconditional response then reads (see Appendix C) The last term stems from a process where a single particle ends up in the ancilla after the coupling. Post-selecting on the empty ancilla eliminates this undesired contribution, yielding a pure non-Hermitian evolution that gives access to the unequal-time anti-commutator.  Figure 4. Simulation of the linear response to a non-Hermitian perturbation generated by a single step in the Zeno evolution of coupling to an ancilla followed by a projection on the empty-ancilla subspace (see Fig. 3c). (a) Time trace of the density after applying the perturbation to a single site as in Fig. 3a. The unnormalized and normalized responses correspond to Eqs. (20) and (21), respectively. The result agrees well with the response to a non-Hermitian perturbation by the density operator A = n (NHH, green dashed line). (b) The correlation spectrum extracted from the response in (a) agrees well with the exact result SBA. The FDR (1) between χ BA and SBA, calculated using the known temperature k B T / J = 4.27 of the thermal state, is shown for comparison. (c) Time trace of the nearest-neighbor correlator after coupling two neighboring sites simultaneously (black, see Fig. 3b) or individually (gray, Fig. 3a) to the ancilla. Subtracting the latter quantity from the former yields the response to a perturbation by the hopping operator A = a † a +1 + a † +1 a . The respective responses agree well with their effective descriptions in terms of non-Hermitian Hamiltonians (NHH, green dashed lines). (d) The extracted correlation spectrum reproduces the exact one to good accuracy.

Numerical benchmark: non-Hermitian linear response in time domain
To benchmark our scheme, we numerically solve the full quantum evolution describing a measurement of S BA (t, t w ) for a thermal state ρ T = exp(−H 0 /k B T )/Z(T ) in a Bose-Hubbard chain of L = 8 sites at unit filling and with periodic boundary conditions. Here, Z(T ) = Tr[exp(−H 0 /k B T )] is the canonical partition sum, and the temperature T is chosen such that the mean energy H 0 T = Tr(H 0 ρ T ) corresponds to that of a Mott-insulating state. A thermal state is an ideal benchmark for our purposes since the temperature T is known and the FDR is satisfied exactly, so any deviations from the FDR indicate deficiencies of the method.
In Figs. 4a and 4c, we show the time traces of the responses to perturbations corresponding to the coupling configurations in Figs. 3a and 3b, respectively, i.e., for on-site densities (A = B = n ) and nearestneighbor . From the latter measurement, the response for the combination A = B = a † a +1 + a † +1 a can be obtained by subtracting the response of the same observable B to perturbations A involving only the densities at the relevant sites. Experimentally, the nearestneighbor correlator B = a † a +1 + h.c. can be measured, e.g., by projecting the system on non-interacting double wells and monitoring the double-well occupancy as a function of time [100,101]. The coupling to the ancilla is applied as a rectangular pulse of duration Jδt = 0.01 and its strength is chosen such that the effective coupling becomes s = 0.05 for the density and s = 0.02 for the correlator, corresponding to a decay of the norm by about 10 % in both cases. As can be seen in Figs. 4a and 4c, the simulated ancilla measurement agrees well with the description in terms of the effective non-Hermitian Hamiltonian in Eq. (19). In Figs. 4b and 4d, we compare the correlation spectra extracted from the responses in Figs. 4a and 4c, respectively, with the exact result. The Fourier integrals have been calculated using exponential filters of characteristic frequencies γ/J = 0.1 for the density and γ/J = 0.05 for the correlator. Due to the sizable static contribution to the response in the latter case, the height of the central peak in Fig. 4d strongly depends on the choice of γ, but this is irrelevant for probing FDRs because the value of the correlation spectrum at ω = 0 is not constrained by the FDR (1). Up to small deviations resulting from non-linear effects, which can be reduced at the cost of a lower signal-to-noise ratio (see discussion in Section VI C), our scheme provides an accurate measurement of the correlation spectrum for both densities and correlators.

B. Non-Hermitian linear response via the pulsed quantum Zeno effect
We now explain how to realize a prolonged evolution under a (possibly time-dependent) effective non-Hermitian Hamiltonian, suitable for probing frequencyresolved responses as is common in standard linear response experiments. To this end, we generalize the coupling Hamiltonian in Eq. (17) by allowing for an arbitrary modulation g(t) of the coupling strength, i.e., We first note that the naive approach of extending the coupling duration in the previous scheme, consisting of a single Zeno step of coupling plus projection, does not yield the desired result. If the coupling duration in Eq. (18) is prolonged up to the final measurement time t > t w , instead of Eq. (20), we obtain to leading order the response The three-time correlations in the integrand appear because the leading perturbative contribution to the response is of quadratic order in the coupling Hamiltonian (see Appendix C 3 for details). If g(t) is properly normalized and has compact support on the interval [t w , t w + δt] with δt sufficiently short as compared to the characteristic time scales of H 0 , Eq. (24) reduces to Eq. (20), but in general does not yield the desired two-time anticommutator.
The key to obtaining a response as in Eq. (7) is to iterate the Zeno step presented in the previous subsection as depicted in Fig. 3d. Such a repeated series of measurements is the common scenario for the pulsed quantum Zeno effect [71,72,93]. To this end, we split the inter- (18), corresponding to an individual Zeno step of unitary evolution under the Hamiltonian H(t) = H 0 +H cpl (t), followed by a measurement of the ancilla population that projects the system on the subspace with empty ancilla (realizations where one or more particles are detected in the ancilla are discarded). Thus, the state evolves, up to a normalization, as denotes the time evolution operator from time t i−1 to t i . This equation describes the evolution under a continuously applied system-ancilla coupling with intermittent measurements of the ancilla population. The role of the measurements is to destroy the coherences between the relevant Zeno subspaces, giving rise to a different evolution than in Eq. (24), where a measurement is performed only once at the final time. In Appendix C 3, we show that, to leading order in the coupling and for δt sufficiently short as compared to the characteristic time scales of H 0 and g(t), the (unnormalized) expectation value of an observable B after the Zeno evolution is given by with A = a † a and s = (Ωδt) 2 /2. Approximating the sum by an integral, this result coincides with a linear response to the anti-Hermitian perturbation H 1 (t) = −if (t)A according to Eq. (7), where f (t) = g 2 (t)Ω 2 δt/2. Since the operator f (t)A is positive semi-definite (cf. Section III), this effective non-Hermitian Hamiltonian describes a gradual leakage of probability out of the emptyancilla Zeno subspace (see Fig. 3d).

C. Non-Hermitian linear response via the continuous quantum Zeno effect
Unfortunately, implementing the pulsed Zeno effect without destroying the sample during the intermittent measurements poses a prohibitive layer of complexity for many experiments. For this reason, we instead exploit the continuous Zeno effect in what follows (see Fig. 3d). This formulation of the Zeno effect arises in the presence of a strong coupling to an external system, which plays the role of a measurement apparatus and leads to wildly fluctuating phases between the relevant Zeno subspaces [72]. One way of generating such a continuous Zeno effect is by adding engineered classical noise to the system, which has been proposed, e.g., in Ref. [73] to constrain the dynamics of quantum simulators for lattice gauge theories. Here, we apply this idea to realize a time-dependent effective non-Hermitian perturbation.

From engineered dissipation to non-Hermitian dynamics
We consider the ancilla to be subject to classical dephasing noise, as indicated in Fig. 3. Such a source of noise can be engineered via a rapidly fluctuating effective detuning, e.g., in form of a Zeeman or ac Stark shift, acting on the ancilla only. We assume that the fluctuations are sufficiently fast compared to all relevant physical time scales, such that their effect can be approximated by a Gaussian white-noise process ξ(t) satisfying ξ(t) = 0 and ξ(t)ξ(t ) = δ(t − t ), where · · · denotes the ensemble average over all noise realizations. For example, using lasers to generate an ac Stark shift, this technical requirement can be fulfilled using acousto-optical devices [102]. The evolution of the density operator ρ(t) can then be described by the stochastic von Neumann equation [73] with dephasing rate κ > 0 and Wiener increments dW (t) = ξ(t) dt, subject to the Stratonovich interpretation of stochastic calculus [103,104] (see also Appendix D). The deterministic part of Eq. (27) is governed by the Hamiltonian H(t) = H 0 + H cpl (t), where the coupling Hamiltonian H cpl (t) is given by Eq. (23). By virtue of stochastic calculus, it can be shown (see Appendix D 1) that the noise-averaged density opera-tor σ(t) ≡ ρ(t) satisfies the Lindblad master equation with the Hermitian Lindblad operator L = b † b. The stochastic differential equation (27) represents a diffusive unraveling [105] of the master equation (28). Such diffusive unravelings typically arise in the theory of continuous measurements, where a quantum system is continuously monitored and the resulting measurement back action gives rise to diffusive quantum trajectories [106,107]. By contrast, in our case, there are no actual measurements involved and Eq. (27) describes a random unitary evolution with pure dephasing [108,109]. In fact, there exists an infinite number of stochastic unravelings, both diffusive and jump-like, whose ensemble average is described by Eq. (28) [75,81,110]. As an alternative to the approach in Eq. (27) using engineered dephasing, we could also start from Eq. (28) with the Lindblad operator L = b, describing a spontaneous decay of particles in the ancilla at a decay rate κ. As shown in Appendix D, such a setting gives rise to the same effective non-Hermitian Hamiltonian as considered below.
The quantum Zeno effect is realized in the strong-noise limit κ → ∞ [73]. The strong dissipation leads to an exponential decay of coherences between Zeno subspaces, in analogy to the effect of repeated measurements, and thereby suppresses the build-up of population in the ancilla. As shown in Appendix D 2, to leading order in perturbation theory, the density operator σ P (t) = Pσ(t)P, projected on the subspace with no particles in the ancilla, obeys the evolution equation generated by the effective non-Hermitian Hamiltonian As required in Section III, the perturbation operator A is positive semi-definite and f (t) is nonnegative, describing a leakage of probability out of the empty-ancilla subspace. In Fig. 3d, we illustrate that this effective non-Hermitian dynamics arises as the ensemble average over stochastic trajectories governed by Eq. (27). The crucial advantage of the continuous Zeno effect over the pulsed formulation, where repeated non-destructive measurements are required, is that a single projection at the final measurement time is sufficient, which can conveniently be realized as a post-selection on measurement outcomes where no particles are detected in the ancilla.

Numerical benchmark: non-Hermitian linear response in frequency domain
We now demonstrate how our scheme enables access to the FDR directly in frequency domain. From the structure of the general linear response formula (4) it becomes clear that by applying a non-Hermitian perturbation under a suitable periodic modulation f (t) continuously until the final observation time t f , it is possible to directly measure non-Hermitian dynamic susceptibilities of the form (30) It is our goal to extract the reactive part χ (NH) of this quantity from the linear response to the effective non-Hermitian Hamiltonian in Eq. (29). For simplicity, we focus on the common case where χ (NH) corresponds to the real part of Eq. (30), and we consider Hermitian operators A and B such that the response function (10) is real. Due to the non-negativity constraint on f (t), it is not possible to modulate the effective coupling around zero. Instead, we choose the modulation in Eq. (27) as ] Ω 2 /κ. According to Eq. (9), the response is then given by where ] is the conditional expectation value obtained from post-selection. The first two terms on the right-hand side of Eq. (31) represent the response to a static non-Hermitian perturbation with g(t) ≡ 1 and the last term is proportional to the desired real part of Eq. (30), which can be seen after changing the integration variable to ∆t = t f − t.
Thus, it is possible to extract the quantity χ (NH) BA (t f , ω) for a given probe frequency ω from two linear response measurements, one with a periodic modulation and one with a constant perturbation, the latter being subtracted from the former.
To benchmark our protocol, we resort to the previous example of the density autocorrelation spectrum in a periodic 1D Bose-Hubbard chain. For this purpose, we numerically solve the stochastic von Neumann equation (27) -the most fundamental equation in our approach -using stochastic Magnus integration [111][112][113]. By comparing the results of the stochastic simulation to those obtained based on Eqs. (28) and (29), we demonstrate the validity of the approximations underlying the effective description in terms of a non-Hermitian Hamiltonian.
When choosing the final evolution time t f , which determines the cutoff of the integrals in Eq. (31), it is important to keep in mind the trade-off between signalto-noise ratio and accuracy: while a longer propagation time t f can yield a more accurate approximation of the Fourier integral in Eq. (30), the strength of the perturbation must typically be reduced accordingly in order  Figure 5. Simulation of the non-Hermitian linear response scheme for measuring the density autocorrelation spectrum in frequency domain. The effective non-Hermitian perturbation is generated through coupling to an ancilla subject to strong engineered dephasing, exploiting the continuous quantum Zeno effect. The simulations are based on the stochastic von Neumann equation (27) (SvNE), the master equation (28) (ME), and the effective non-Hermitian Hamiltonian (29) (NHH). The stochastic simulation has been averaged over 200 realizations, and the error bars show the ensemble standard deviation of the mean. (a) Decrease of the norm resulting from the projection on the empty-ancilla subspace for the probe frequency ω = 0 and norm decay q = 0.15 up to the final evolution time Jt f = 2. A stronger dephasing rate κ improves the agreement between the ME (SvNE) and NHH descriptions at early times, while the deviations at later times are due to non-linear effects. (b) Unnormalized and normalized (conditional) responses, corresponding to Eqs. (7) and (31), respectively, as a function of frequency for a fixed final time Jt f = 2 and norm decay q = 0.15. (c) Correlation spectra Snn extracted from the responses in (b) according to Eq. (31), in comparison with the exact result and the FDR (1) between S and χ (evaluated for the same truncation Jt f = 2). The different combinations of the parameters κ and q for the ME simulation show that the agreement with the exact result can be improved by going deeper into the limit of strong dissipation and weak perturbations.
to stay within the linear response regime, which lowers the signal-to-noise ratio. The optimal balance between these effects depends on several conditions, such as the targeted frequency range, the properties of the response function, and the resolution of the measurement apparatus. For concreteness, in the following benchmark example we choose Jt f = 2 for all probed frequencies ω. While this truncation affects the form of the extracted spectrum at low frequencies [114], it yields an adequate approximation of the Fourier integral at higher frequencies most relevant for probing the FDR.
To account for the different sensitivities of the responses at different probe frequencies, we parametrize the perturbation strength in terms of the norm decay q due to the effective non-Hermitian Hamiltonian. That is, for each frequency ω, given the fixed final observation time t f and the dephasing rate κ, we adjust the coupling strength Ω such that according to Eq. (8) the norm of the state has decreased by the amount q at the end of the evolution. For a translationally invariant system at unit filling, we have A(t) 0 = n(t) 0 = 1, and therefore Ω = [κq/2t f (1 + sinc(ωt f /π))] 1/2 , where sinc(x) = sin(πx)/πx.
In Fig. 5, we compare various simulations of the scheme based on Eqs. (27) to (29) for a system of L = 4 lattice sites. Figure 5a illustrates a typical decay of the norm over time due to the effective non-Hermitian perturbation. The perturbation strength is adjusted such that by the final evolution time Jt f = 2 the norm has dropped approximately by an amount q = 0.15, which results in a good signal-to-noise ratio, but lies slightly beyond the onset of the non-linear regime (longer propagation times may require balancing with a reduced perturbation strength). The simulated unnormalized and normalized (conditional) responses as a function of frequency are shown in Fig. 5b, from which we extract the correlation spectra presented in Fig. 5c. The results are compared to the exact correlation spectrum, which is evaluated for the same truncation Jt f = 2 of the Fourier integral in Eq. (30) to allow for a consistent benchmark. For the stochastic simulation, we choose the accessible dephasing rate κ/J = 10. As can be seen in Fig. 5, the stochastic simulation based on Eq. (27) agrees with the simulation based on the master equation (28) within the statistical error bars that show the ensemble standard deviation of the mean for an accessible number of 200 realizations. Moreover, Fig. 5c shows that these parameters already yield the correlation spectrum at a reasonable accuracy suitable for certifying the validity of the FDR (1). The description in terms of the effective non-Hermitian Hamiltonian in Eq. (29) is closer to the exact result than the description in terms of the master equation (28) for the same parameters, revealing that the linear regime is wider for the former than for the latter, which could be remedied through extrapolation. In the effective non-Hermitian description (29), the coupling strength Ω and the dephasing rate κ enter only via the ratio Ω 2 /κ, which is proportional to the norm decay q, while these two parameters enter Eqs. (27) and (28) individually. Going deeper into the Zeno limit of large κ improves the validity of Eq. (29) at early times and at higher frequencies, shown in Fig. 5 for κ/J = 100. Decreasing at the same time the effective coupling strength, as illustrated in Fig. 5c for q = 0.05, the agreement between the extracted correlation spectrum and the exact result improves further, which shows that, at the cost of decreasing the signal-to-noise ratio, the exact correlation spectrum can in principle be approximated to arbitrary accuracy.

VI. DISCUSSION
In this section, we put our non-Hermitian linear response approach for measuring dynamical correlations and FDRs in perspective with other schemes. We first demonstrate that common ancilla-based weak measurement protocols [43,44] fit into this general framework since their ancilla-free formulations can be interpreted as a non-Hermitian linear response. In addition, we compare the ancilla-based technique for realizing non-Hermitian linear response presented in Section V with other schemes for accessing non-Hermitian dynamics or measuring dynamical correlations, including non-invasive and projective protocols [30,41,[43][44][45]. We conclude with a discussion of experimental aspects and potential error sources.

A. General relation between non-Hermitian linear response and ancilla-based weak measurements
To reveal the close connection between non-Hermitian linear response and ancilla-based weak (or non-invasive) measurements of the unequal-time anti-commutator, we first briefly review common protocols of the latter kind. While these weak measurement protocols have originally been developed for spin systems [43,44], here we formulate them for general quantum systems and allow for arbitrary durations of the system-ancilla coupling (typically, only short coupling pulses are considered). For further details on the following points, see Appendix E.
System and ancilla are assumed to be initially in a product state, ρ 0 = ρ S ⊗ ρ A . The non-invasive protocol of Refs. [43,44] starts by evolving the system under the unperturbed Hamiltonian H 0 up to a certain waiting time t w , while the ancilla does not participate in the dynamics. System and ancilla are then coupled by the Hamiltonian where A and X are Hermitian operators acting on system and ancilla, respectively, and f (t) represents an arbitrary time-dependent modulation. The form of the coupling Hamiltonian is one of the main differences to our protocol in Section V A (see discussion in the next subsection). After a coupled evolution up to time t > t w , one measures projectively the observables B on the system and Y on the ancilla, respectively. Instead of directly correlating the measurement outcomes as proposed in Refs. [43,44], we consider here conditional expectation values in order to reveal the connection to non-Hermitian linear response. As derived in Appendix E, the expectation value of B under the condition that the ancilla measurement of Y yields the outcome y is given, to linear order in the coupling, by (33) with λ y = P y X 0 / P y 0 ∈ C. Here, P y is the projector on the eigenspace of eigenvalue y and c.c. denotes the complex conjugate. The key to access the unequaltime anti-commutator is to choose the ancilla state as well as the operators X and Y such that the expectation value P y X 0 = Tr[ρ A P y X] becomes purely imaginary (choosing P y X 0 real instead yields the commutator). The conditional expectation value of B is then formally equivalent to the non-Hermitian linear response in Eq. (9). In fact, by tracing out the ancilla, it can be shown (see Appendix E) that the coupled evolution of system and ancilla corresponds to the evolution under the effective non-Hermitian Hamiltonian H eff = H 0 −isf (t)A with s = iλ y ∈ R. Interestingly, this leads to the insight that any such weak measurement protocol for the unequal-time anti-commutator can be interpreted as the linear response to a non-Hermitian perturbation.
The general connection between non-Hermitian linear response and ancilla-based weak measurements is beneficial for both disciplines: particular observables previously accessible only via ancilla-based schemes may be obtainable more efficiently in an ancilla-free way using the tools of non-Hermitian physics, while certain non-Hermitian Hamiltonians difficult to engineer directly may be realized with the help of an ancilla.

Ancilla-based weak measurements
One of the main challenges of the ancilla-based weak measurement scheme for the unequal-time anticommutator discussed above is to engineer the ancilla state as well as the observables X and Y in such a way that P y X 0 becomes purely imaginary. While Refs. [43,44] discuss suitable configurations for spin systems, it is far less obvious how to choose the setup in an experimentally feasible way on other platforms such as bosons in optical lattices (for example, as discussed in Appendix E, to access density correlations, number non-conserving coupling Hamiltonians may be required, which cannot be realized with massive particles).
By contrast, our ancilla-based scheme in Section V relies on the system-ancilla coupling in Eq. (17), which is quadratic in the creation and annihilation operators. Despite its simple form, the coupling can flexibly be adapted to measure the unequal-time anti-commutator of a wide range of previously inaccessible observables such as nearest-neighbor correlators, as discussed in Section V A and Appendix C 2. In addition, our choice of the initial ancilla state in the form of the vacuum is particularly easy to prepare experimentally. On the formal level, an important difference to the common weak measurement approach is that for our choice of the coupling Hamiltonian (17), the linear order in perturbation theory vanishes and the leading contribution to the response stems from the quadratic order (see Appendix C), where the anti-commutator naturally arises and can be isolated by post-selection on realizations without any particles in the ancilla.

Projective protocols
Projective protocols allow one to probe dynamical correlations of dichotomic observables (observables with two eigenvalues) by performing consecutive projective measurements directly on the system and correlating the outcomes in a suitable way [30,41,[43][44][45]. As such, compared to schemes based on weak perturbations such as linear response, projective protocols are backaction-free and feature a higher signal-to-noise ratio. Despite these advantages, the fact that projective protocols work only for dichotomic observables restricts their general applicability.
In Ref. [30], it has been analyzed how projective protocols can be applied to approximately dichotomic observables, in particular densities in Bose-Hubbard systems close to the hard-core limit. In Appendix F, we present numerical benchmarks comparing the performance of projective protocols and non-Hermitian linear response for measuring the unequal-time anticommutator in Bose-Hubbard systems at various fillings and on-site interactions. Our analysis shows that projective protocols perform well at low fillings and large onsite interactions, but yield unsatisfactory results when applied beyond this regime, i.e., as soon as multiple occupancies can no longer be neglected. In particular, the relevant scenario of Bose-Hubbard systems at unit filling and moderate values of U/J, which we study in Section IV inspired by the experiment of Ref. [51], remains beyond the scope of projective protocols. By contrast, our non-Hermitian linear response approach does not have restrictions on observables regarding the number of eigenvalues and performs well across the entire parameter space explored in Appendix F. Thus, our scheme allows one to reliably access unequal-time anticommutators and the associated FDRs also in regimes outside the range of projective protocols, e.g., in the pair superfluid phase of dipolar bosons [115] or other phases where multiple occupancies play an essential role.

Dissipative dynamics and post-selected quantum trajectories
The structure of the second term on the right-hand side of Eq. (22), characterizing the unconditional response af-ter coupling to the ancilla without measuring the ancilla population, resembles to the "recycling term" in Lindblad master equations [77,81]. In fact, the short coupling pulse to the ancilla can be viewed as an effective dissipative perturbation, ρ} is the Lindblad dissipator with dissipation rate γ = s/δt. This yields Eq. (22) for the expectation value of an observable B after a unitary evolution up to time t. In the quantum trajectories picture [75][76][77], the Lindblad dissipator D[ρ] generates an evolution under the non-Hermitian Hamiltonian H 1 = −i γa † a, subject to quantum jumps described by the "recycling term" 2γaρa † . By post-selecting on the absence of quantum jumps, it is possible to isolate the pure non-Hermitian evolution [56,58,66]. From this point of view, the projection on the empty-ancilla subspace in Eq. (18) can be interpreted as a post-selection on the absence of quantum jumps, i.e., particles hopping to the ancilla. This allows us to eliminate the undesired contribution in Eq. (22) due to the "recycling term" and obtain instead the result in Eq. (21), reflecting a purely non-Hermitian perturbation that gives access to the unequal-time anti-commutator.

C. Experimental considerations and error sources
Many experimental setups such as quantum gas microscopes permit the simultaneous readout of all site populations in a single shot [94,95]. This is convenient for simultaneously measuring the responses of different observables B, e.g., B = n for = 1 . . . L, to a fixed perturbation A determined by the coupling scheme. In addition, for the single Zeno step and the continuous Zeno evolution in Sections V A and V C, respectively, the measurement of the ancilla population can be deferred up to the final observation time t and measured along with the other site populations (cf. Ref. [43]). The projection on the empty-ancilla subspace is then achieved by post-selecting those realizations where no particles are detected in the ancilla. Since the effective coupling s needs to be chosen sufficiently weak to stay within the regime of linear response, the fidelity of the post-selection is typically high (see Fig. 3c). However, there is the usual linear response trade-off between maximizing the measurement signal (large s) and staying within the perturbative regime where the linear approximation is valid (small s).
One can distinguish two types of detection errors: false positives, i.e., at least one particle is detected in the ancilla, but there is actually none, and false negatives, i.e., no particles are detected, but there is at least one. Let α be the false positive rate and let β be the false negative rate. If the measurement is post-selected on the condition that no particles are detected in the ancilla, which may in some cases be erroneous, the conditional state in Eq. (18) is replaced by ρ = (1 − α)PρP + βQρQ, where ρ is the state right after the coupling and before the pro-jection, and Q = 1 − P is the projector on the subspace with a non-vanishing ancilla population (for simplicity, we do not distinguish different false-negative probabilities within the Q subspace since the error due to single occupancies dominates in the linear regime). This shows that false positives lower the measurement fidelity, while false negatives contribute a systematic error to the result, arising from the inadvertent projection on a complementary subspace (see the discussion in the context of Eq. (22) and Appendix C 1).

VII. CONCLUSION
In this work, we have demonstrated that non-Hermitian linear response enables access to the unequaltime anti-commutator as the missing piece for the direct observation of the FDR in quantum systems. As an illustration, we have discussed how a Bose-Hubbard system after a global quench reaches thermal equilibrium, and we have derived techniques to generate the required non-Hermitian dynamics in cold-atom systems coupled to an ancillary mode by exploiting the quantum Zeno effect. This proposal provides a concrete scenario for the direct observation of the FDR and an unbiased way of probing thermalization dynamics in state-of-the-art experiments on synthetic quantum matter.
Our non-Hermitian linear response approach is completely agnostic to specific platforms and implementations, and as such can be applied in any non-Hermitian system. It is independent of microscopic details such as interactions, geometry, or particle statistics, and can thus be used with bosons, fermions, or spins alike. Higher orders in the response may be used to access nested unequal-time anti-commutators of increasing order. Moreover, we have shown that common ancillabased weak measurement protocols for dynamical correlations fit in the same framework, as these can be interpreted in the light of (non-)Hermitian linear response. Our proposed ancilla-based realization of non-Hermitian linear response permits the extraction of dynamical correlations for a wide range of previously inaccessible observables beyond density correlations, even frequencyresolved. While we have focused on lattice systems, our protocol can immediately be applied to continuous systems, e.g., via spatially focused laser beams, giving access to dynamical correlations of the field operator coarsegrained over a small region in space [116]. The discussed framework thus provides an array of possibilities to experimentally -and also numerically [8,29] -characterize quantum systems in and out of thermal equilibrium.

Derivation of the non-Hermitian linear response formula
To derive Eq. (7), we first transform to the interaction picture with respect to the unperturbed (Hermitian) Hamiltonian H 0 . The von Neumann equation (6) then reads is the anti-Hermitian perturbation in the interaction picture. This equation can equivalently be expressed in integral form as withρ(0) = ρ 0 . To linear order in the perturbation, we can replaceρ(t ) in the integrand by ρ 0 , yielding The expectation value of an observable B can be computed in the interaction picture as B(t) = Tr[B(t)ρ(t)], whereB(t) = e iH0t/ Be −iH0t/ . Inserting Eq. (A3) into this expression and using the cyclic property of the trace leads to the result in Eq. (7).

Connection between correlation spectrum and non-Hermitian dynamic susceptibility
To show that S BA (τ, ω) = − χ (NH) BA (τ, ω), we first note that the symmetric correlation function (3) obeys the symmetry relation S BA (t, t ) = S AB (t , t). In what follows, we use the short-hand notation S BA (τ, ∆t) = S BA (t = τ + ∆t/2, t = τ − ∆t/2). Then, the aforementioned identity reads S BA (τ, ∆t) = S AB (τ, −∆t). This allows us to express the correlation spectrum as (A4) Using further Eq. (10), noting that the Heaviside step function allows us to extend the integration domain to negative ∆t, as well as the definition of the generalized susceptibility, Eq. (5), we arrive at By contrast, if we consider the Fourier transform at fixed waiting time t w as in Eq. (15), the integrand S BA (t w + ∆t, t w ) does not possess a symmetry with respect to the relative time ∆t in general. Thus, out of equilibrium, we generally have S BA (t w , ω) = − χ (NH) BA (t w , ω), but this relation is restored once the system reaches a stationary state.

Appendix B: Fluctuation-dissipation relations after a quench in a Bose-Hubbard system
In this Appendix, we provide details on our analysis in Section IV of FDRs following a quench in a Bose-Hubbard system. In particular, we investigate in more detail the behavior of the deviation from the FDR for offsite density correlations as a function of distance. To this end, we study the error as a function of system size and compare our results in Section IV to a similar analysis for a Bose-Hubbard system in 2D.

Technical details on the computation of dynamical susceptibilities
The susceptibilities appearing in the FDR (13) can be found by varying both the waiting time t w and the observation time t = t w +∆t, and computing the Fourier transform with respect to the relative time ∆t at fixed central time τ according to Eqs. (5) and (11). However, since at early times τ the integration domains of the Fourier integrals in Eqs. (5) and (11) are limited to only a short range of ∆t, we consider instead the non-equilibrium generalization of the susceptibility as the Fourier transform of the response function at fixed waiting time t w [29] according to Eq. (15) (and analogously for its Hermitian counterpart). This form has the advantage that the integration domain is unbounded, and it is often more efficient to compute as the integrand is directly obtained from the linear response at fixed waiting times. In practice, the integral in Eq. (15) needs to be regulated appropriately, for example, by truncating it once the correlations have decayed sufficiently or by means of a frequency filter accounting for a finite spectral resolution in experiments. Unless stated otherwise, we follow the latter approach, using an exponential filter of characteristic frequency γ/J = 0.2, which amounts to the replacement e iω∆t → e (iω−γ)∆t in the Fourier integral (15).
Note that out of equilibrium, the susceptibilities in Eqs. (11) and (15) generally disagree, and, in particular, S BA (t w , ω) = − χ (NH) BA (t w , ω) in general (cf. Appendix A 2). Once the system has reached a stationary state, the response functions depend only on the relative time ∆t and the different conventions become equivalent, provided the integration domains are chosen appropriately. For our purposes, we probe the FDR out of equilibrium in the form of Eq. (13), expressed in terms of the susceptibilities obtainable from the (non-)Hermitian linear response at fixed waiting time as in Eq. (15).
Moreover, the susceptibility components in Eq. (12) are in general complex [3] and require measurements of both the response of B to a perturbation by A and vice versa. Sufficient conditions for them to be real include the case B = A † , or, if A and B are Hermitian, the property B(t)A(t ) = A(t)B(t ) . The latter is fulfilled, for instance, if A and B are on-site observables in a system that is invariant under both translations and reflections, as it is the case for density-density correlations in a Bose-Hubbard chain with periodic boundary conditions. Thus, for our model system and our choice of observables, the reactive and dissipative parts of the dynamic susceptibility are real and correspond, respectively, to the real and imaginary parts of the susceptibility (11) (and similarly for the Hermitian susceptibility in Eq. (5)).
We quantify deviations from the FDR (13) by the absolute error where · 2 denotes the L 2 norm, which we define by Frequency ω/J Figure 6. FDRs for off-site density correlations at early and late waiting times tw for several distances d. The dynamic susceptibility χ (blue) is compared to χ (NH) (red), rescaled according to the FDR (13) using the least-squares result for the effective temperature indicated in the plots. At small distances and late times, the curves overlap well, while at larger distances discrepancies persist even after long times.

One-dimensional Bose-Hubbard system
To better understand the significance of the increase of the relative error as a function of distance in Fig. 2c, we show in Fig. 6 the extracted FDRs at early and late waiting times for the individual distances. The data are the same as in Fig. 2 and the susceptibility χ (NH) has been rescaled according to the FDR (13) using the indicated effective temperature T obtained from the least-squares fit in Eq. (16). At early waiting times, there is no global value of T to make χ (NH) and χ overlap, and the FDR is clearly violated (in some cases, an attempted fit can even yield unphysical negative temperatures). By contrast, at late waiting times, χ (NH) and χ fulfill the FDR and at small distances, the extracted effective temperatures are consistent with the temperature k B T / J = 4.27 of a thermal state at the same energy density as the initial state (calculated for L = 8 using exact diagonalization). At larger distances, some peaks in Fig. 6 exhibit clear deviations which persist even after very long times and contribute to the increased relative error in Fig. 2c.
For an ergodic system in the thermodynamic limit, it is generally expected that a two-site subsystem, regardless of the distance between the two sites in real space, eventually thermalizes and thus satisfies the FDR. The observed deviations in Figs. 2 and 6 may therefore be an artifact of the finite system size. Apart from that, numerical errors induced in the course of the data analysis, such as integration and truncation errors in the evaluation of Fourier integrals or distortions caused by the frequency filter, may contribute to the deviation. We have checked that improving on the latter points does not alter the picture qualitatively.
To study the influence of finite-size effects, we have calculated the error as a function of the particle number N (corresponding to the number of lattice sites L at unit filling) up to N = L = 16. Figure 7 shows the relative and absolute errors at the moderate waiting time Jt w = 10 for all possible distances d in the respective systems. Note that for a periodic chain of length L, the maximum distance is d = L/2 . The Fourier integrals have been truncated at J∆t = 30 using an exponential filter of characteristic frequency γ/J = 0.1. Both the relative and the absolute errors for all distances decrease as the system size increases until the relative error saturates at a value close to zero. Although the exponential growth of the Hilbert-space dimension makes an exact numerical treatment of even larger systems inaccessible, the clear trend in Fig. 7 suggests that the deviations from the FDR in Figs. 2 and 6 at large distances for L = 12 are likely due to finite-size effects. Our analysis thus confirms the expectation that the two-site subsystems relevant for offsite density correlations thermalize and thus fulfill the FDR, provided the system is not too small.

Two-dimensional Bose-Hubbard system
To show that our results for the 1D Bose-Hubbard chain are generic, we study the analogous quench scenario in a 2D Bose-Hubbard system. We consider a system of 4 × 4 lattice sites with N = 16 particles (unit filling) and periodic boundary conditions in each direction. The larger system size compared to the 1D setting above al- lows us to support the conjecture that the FDR is better fulfilled as the system size increases. As before, we initialize the system in a Mott-insulating state and quench at time t = 0 into the superfluid phase at U/J = 1.5625. In Fig. 8, we present the same analysis for the 2D system as carried out in Fig. 2 for the 1D chain (cf. Section IV). Due to the periodic boundary conditions and the isotropy of the hopping, the curves fall into five classes corresponding to distances d = 0 . . . 4 between the perturbed and probed lattice site. Note that these distances do not correspond to the physical distances in the 4 × 4 lattice, but to the minimum number of hopping events connecting the two sites. Similarly to the 1D setting, the FDR is violated at short waiting times, indicated by the large relative error in Fig. 8c. After times on the order of J −1 , the errors decrease dramatically and the FDR is fulfilled for all accessible distances with only a minor trend towards larger relative errors for larger distances. The effective temperatures in Fig. 8b for the individual distances reach approximately constant values Frequency ω/J Figure 9. Same as Fig. 6, but for the 2D Bose-Hubbard system. The data are the same as in Fig. 8. At long waiting times, the dynamic susceptibility χ (blue) agrees well with χ (NH) (red) after the latter is rescaled according to the FDR (13).
that mutually agree up to deviations of about ten percent or less. Furthermore, the effective temperatures are close to the temperature k B T / J = 8.91 of a thermal state at the same energy density as the initial state (calculated for a 3 × 3 system using exact diagonalization). The fulfillment of the FDR at late times is further illustrated in Fig. 9 (analogously to Fig. 6), where the agreement between χ and the rescaled χ (NH) is remarkable. This supports our conclusion in Appendix B 2 that the deviations observed for smaller systems are likely due to finite-size effects, while sufficiently large systems thermalize as expected.

Appendix C: Derivation of non-Hermitian linear response via the pulsed quantum Zeno effect
In this Appendix, we derive how an effective non-Hermitian perturbation can be generated by coupling the system to an ancilla and exploiting the quantum Zeno effect (see Fig. 3). We first discuss in detail a single step in the Zeno evolution, which corresponds to a δlike non-Hermitian perturbation that allows one to access the unequal-time anti-commutator in time domain. Furthermore, we derive the general form of perturbation operators that can be realized this way. Finally, we explain how the pulsed Zeno effect, generated by repeatedly projecting the coupled system on the subspace with no particles in the ancilla, allows one to realize an effective evolution under a non-Hermitian Hamiltonian for an extended period of time.

Single Zeno step
In what follows, we use time-dependent perturbation theory to derive Eqs. (21) and (22), describing, respectively, the conditional and unconditional response after a single step in the Zeno evolution of coupling to the ancilla followed by a projection on the empty-ancilla subspace, as depicted in Fig. 3c.
The protocol starts by evolving the initial state ρ 0 under the Hamiltonian H 0 up to the waiting time t w , at which the perturbation is applied. Before the coupling, the state is given by ρ(t w ) = e −iH0tw/ ρ 0 e iH0tw/ . In Section V A, we have approximated a δ-like perturbation as a rectangular pulse of duration δt. Here, we consider a slightly more general scenario where we allow for an arbitrarily shaped pulse g(t) as in Eq. (23). The corresponding total Hamiltonian reads H(t) = H 0 + g(t)H cpl with H cpl given by Eq. (17). It is convenient to work in the interaction picture with respect to the unperturbed Hamiltonian H 0 . Time evolution is then governed by the von Neumann equation whereρ(t) = e iH0t/ ρ(t)e −iH0t/ andH cpl (t) = Ω ã † (t)b + b †ã (t) withã(t) = e iH0t/ ae −iH0t/ denote, respectively, the density operator and the coupling Hamiltonian in the interaction picture. Rewriting Eq. (C1) as an integral equation and substituting the left-hand side into the right-hand side, we arrive at For the discussion of a single step in the quantum Zeno evolution, we consider a pulse g(t) with compact support on the interval [t w , t w + δt], normalized such that tw+δt tw dt g(t) = δt. If the pulse duration δt is much shorter than the characteristic time scales of H 0 , we can approximateρ(t ) ≈ρ(t w ) andH cpl (t ) ≈H cpl (t ) ≈ H cpl (t w ) in the integrands, yielding, up to second order in δt, the result We require the ancilla to be empty before the coupling. More specifically, we assume that the combined state of system and ancilla at the waiting time t w is given by the product stateρ(t w ) =ρ S (t w ) ⊗ρ A , where the ancilla is in the pure vacuum stateρ A = |0 0|. (For concreteness, we focus here on bosonic systems, but the derivation for fermions proceeds analogously and yields the same result for the unequal-time anti-commutator.) Inserting this state into Eq. (C3), we obtaiñ whereñ(t w ) =ã † (t w )ã(t w ) is the number operator and h.c. denotes the Hermitian conjugate. After coupling the system to the ancilla, a single step in the Zeno evolution is completed by measuring the population of the ancilla, projecting the state on a subspace with a definite number of particles in the ancilla. Let P n = 1 ⊗ |n n| be the projection operator on the subspace with n particles in the ancilla. Since [P n , H 0 ] = 0, the measurement can optionally be deferred up to the final observation time (cf. Ref. [43]). The projected states read where s = (Ωδt) 2 /2 is the effective coupling strength, and P nρ (t w + δt)P n = 0 for n ≥ 2, up to second order in δt. The probability of detecting n particles in the ancilla is then given by p(n) = Tr[P nρ (t w + δt)P n ], which yields and p(n ≥ 2) = 0, up to second order in δt. Here, we have used Tr[ρ(t w )ñ(t w )] = Tr[ρ 0 n(t w )] = n(t w ) 0 . Remarkably, the result in Eq. (C5a) can, to leading order in the coupling, be expressed as the evolution under an effective non-Hermitian Hamiltonian, withH eff (t) = −i sÃ(t)/δt and perturbation operator A = n = a † a. According to Lüders' rule [117], the conditional state, given that n particles have been detected in the ancilla, is obtained by normalizing the projected states (C5) by the respective probabilities (C6), i.e.,ρ(t w + δt|n) = P nρ (t w + δt)P n /p(n). Up to leading order in s, we find where we have discarded (or traced out) the ancilla and omitted the subscripts indicating system density operators. By contrast, if the ancilla population is not measured or if the measurement outcomes are ignored, the state after the coupling is instead described by the unconditional density operator which can be obtained directly from Eq. (C4) after tracing out the ancilla. For times t > t w + δt, the coupling is switched off and the system evolves solely under the Hamiltonian H 0 . According to Eq. (C1), this evolution is trivial in the interaction picture, such thatρ(t) =ρ(t w + δt). The (unnormalized) expectation value of an observable B with respect to the state (C5a) projected on the empty-ancilla subspace reads where, in the second step, we have transformed from the interaction picture to the Heisenberg picture, and in the last step, we have used the cyclic property of the trace. This is the result reported in Eq. (20) in the main text. For the conditional state (C8a), we recover Eq. (21), which describes a post-selected measurement conditioned on the empty ancilla.
If we were to post-select on the condition that a single particle is detected in the ancilla, corresponding to the conditional state Eq. (C8b), we would instead obtain This quantity contributes a systematic error in the case of faulty detection with false negatives (see the discussion in Section VI B). From the derivation in this Appendix it becomes clear that, up to second order in δt, i.e., up to linear order in the effective coupling strength s, the error is dominated by single occupancies of the ancilla site. Finally, the unconditional expectation value (22) follows from Eq. (C9). This shows that post-selection is essential in order to remove the undesired contribution in form of the "recycling term" due to Eq. (C8b), enabling access to the unequal-time anti-commutator.

General system-ancilla coupling
The coupling schemes in Figs. 3a and 3b are designed to realize non-Hermitian perturbations by the density operator and the hopping operator, respectively. We now consider the general situation where an arbitrary number of system modes is coupled to up to M ancillary modes. This scenario is described by the general coupling Hamiltonian where the operator is a linear combination of system modes a with coefficients λ m ∈ C, coupled to the m-th ancilla with coupling strength Ω m ≥ 0. The configuration in Fig. 3a, a single lattice site * coupled to a single ancilla, is recovered for M = 1 and λ 1 = δ * , while Fig. 3b, two sites 1 and 2 simultaneously coupled to a single ancilla, corresponds to M = 1 and λ 1 = δ 1 + δ 2 . As before, we consider a short coupling pulse of duration δt such that the state after the coupling is given by Eq. (C3). Subsequently, a measurement of the individual ancilla occupancies is performed and the state is conditioned on the outcome of that measurement (as mentioned above, the measurement may also be deferred up to the final observation time). Given the outcome (n 1 , . . . , n M ), the post-measurement state, up to a normalization, reads P n1...n Mρ (t w + δt)P n1...n M , where P n1...n M = 1⊗|n 1 · · · n M n 1 · · · n M | is the projection operator on the subspace with a definite ancilla population corresponding to the measurement outcome.
Up to leading order in the coupling, only processes where at most a single particle ends up in one of the ancillas contribute. Let P 0 = P 0...0 denote the projector on the subspace with all ancillas empty. The projector on the subspace with a single particle in the m-th ancilla and all others empty can then be expressed as P we find the (unnormalized) post-measurement states where we have traced out the ancillas and introduced the effective coupling strengths s m = (Ω m δt) 2 /2. We note that this result holds for fermions as well, where instead of Eq. (C14) the corresponding fermionic anticommutation relations apply. The respective probabilities of finding no particles in any ancilla or a single particle in the m-th ancilla read A comparison of Eqs. (C15) and (C16) with Eqs. (C5) and (C6) shows that the coupling to the m-th ancilla in the general coupling Hamiltonian (C12) generates an effective non-Hermitian perturbation by the operator Coupling to multiple ancillas simultaneously can be used to realize perturbation by (arbitrarily weighted) sums of the operators A m . This demonstrates that our scheme enables flexible access to unequal-time correlations and FDRs for a wide range of observables, two specific examples of which, namely densities and nearest-neighbor correlators, we have illustrated in Section V.

Prolonged non-Hermitian evolution via the pulsed quantum Zeno effect
To gain a deeper understanding of how the pulsed quantum Zeno effect enables a prolonged evolution under an effective non-Hermitian Hamiltonian, we consider in detail two consecutive steps in the Zeno evolution, using a similar formalism as in Ref. [72], and investigate the role of the projective measurement after the first step.
To this end, let P denote the projection operator on the empty-ancilla subspace H P and Q = 1 − P the projector on the complementary subspace H Q = H ⊥ P with at least one particle in the ancilla. It is convenient to write the density operator ρ on the total Hilbert space H = H P ⊕ H Q in the form where ρ PP = PρP and ρ QQ = QρQ are the populations of H P and H Q , respectively, and ρ PQ = PρQ = ρ † QP are the coherences between these two subspaces. Similarly, the time evolution operator from time t 0 to t corresponding to Eq. (C1) can be expressed as Let us denote the initial stateρ(t w ) of the Zeno evolution by ρ 0 and set t 0 = t w . To keep the notation simple, for the purposes of this subsection, we omit the tilde indicating interaction picture operators and use the abbreviation H(t) ≡ g(t)H cpl (t). Since the ancilla is initially empty, we have ρ 0 PQ = ρ 0 QQ = 0. The unitary evolution from time t 0 to t 1 in the presence of the system-ancilla coupling changes the state as where U i PQ = PU (t i , t i−1 )Q. From Eq. (C2), using P 2 = P, Q 2 = Q, PQ = QP = 0, and PH(t)P = 0, we obtain the populations and coherences of ρ 1 , up to quadratic order in the coupling, as ΩPã † (t)bQ. Measuring the ancilla population projects the state on the subspace with a definite number of particles in the ancilla. Without registering the measurement outcome, this yields the unconditional state Pρ 1 P + Qρ 1 Q. Crucially, the measurement process destroys any coherences ρ PQ and ρ QP between the Zeno subspaces H P and H Q . We are interested in measurement outcomes where no particles are detected in the ancilla. Conditioning the state on this outcome corresponds to a projection on the emptyancilla subspace H P , The second Zeno step proceeds analogously to Eqs. (C20) and (C24): the state first evolves unitarily from time t 1 to t 2 in the presence of the system-ancilla coupling and is then projected on the empty-ancilla subspace H P , with up to leading order in the coupling. It is instructive to compare this result to the one obtained if no measurement is performed after the first step. The state then receives additional contributions from the coherences, yielding, to leading order in the coupling, The result in the last line could have been directly obtained from Eq. (C2) for t = t 2 by applying the projector P on both sides. This is evident because without the projection after the first step, the system plus ancilla evolves unitarily from time t 0 to t 2 . However, Eq. (C27) explicitly exposes the crucial effect of the measurement after the first step: the last term in the second-to-last line is precisely the contribution from the coherences ρ 1 PQ and ρ 1 QP , which is missing in Eq. (C26) since the coherences have been destroyed by the measurement.
Iterating the Zeno evolution for n steps up to time t n (including projections after each step), the resulting state is given, to leading order in the coupling, by If the duration t i+1 − t i of each Zeno step is sufficiently short as compared to the time scales of the unperturbed Hamiltonian as well as the modulation g(t), the integrand in each integral is approximately constant, yielding Eq. (26). As discussed in Section V B, this result can in turn be interpolated by a continuous evolution under an effective non-Hermitian Hamiltonian (see Fig. 3c). By contrast, if the state is only projected at the final observation time, but no projections are performed during the evolution as in Eq. (C27), we obtain, to leading order in the coupling, which corresponds to Eq. (24) in the main text. Since the evolution time t n − t 0 may be on the same order or longer than the characteristic time scales of the unperturbed Hamiltonian, it is not possible to approximate the integrand as constant here. Consequently, this procedure does not yield the desired two-time anti-commutator in general.
As these discussions show, exploiting the Zeno effect allows us to apply effective non-Hermitian perturbations for an extended period of time. The essential mechanism is the destruction of the coherences between the Zeno subspaces due to the intermittent measurements. As explained in Section V C and Appendix D, this effect can be mimicked if the ancilla is exposed to strong (engineered) dissipation, which represents an alternative way of realizing non-Hermitian linear response via the quantum Zeno effect.
is only guaranteed if the Stratonovich interpretation is used [108].

Derivation of the effective non-Hermitian Hamiltonian
To derive the effective non-Hermitian Hamiltonian governing the evolution in Eq. (29), following Ref. [73], we consider the strong noise limit of Eq. (28) projected on the empty-ancilla subspace. It is convenient to work in the interaction picture, i.e, in a rotating frame with respect to the unperturbed Hamiltonian H 0 . Equation (28) then reads The operators b and b † as well as the Lindblad operators remain unchanged as they act on the ancilla only and therefore commute with H 0 . We now use the projection operator on the emptyancilla subspace P = P 0 as well as its complement Q = 1 − P to derive coupled equations for the populationsσ PP = PσP andσ QQ = QσQ of the two subspaces, as well as for their coherencesσ PQ = PσQ andσ QP = QσP. The projection operators are Hermitian and satisfy the properties P 2 = P, Q 2 = Q, PQ = QP = 0, as well as [P, H 0 ] = [Q, H 0 ] = 0, the latter following from the fact that H 0 does not change the number of particles in the ancilla. Furthermore, since P projects on the empty-ancilla subspace, we have bP = Pb † = 0. Applying the projectors P and Q to Eq. (D2) from the left and from the right yields the coupled system of equations where the operatorsH PQ (t) = g(t) ΩPã † (t)bQ and H QP (t) = g(t) ΩQb †ã (t)P mix the two subspaces. In deriving Eq. (D2), we have considered the engineered dephasing scenario described by the stochastic von Neumann equation (27), in which case the Lindblad operator L = b † b is Hermitian and the projectors commute with L. In the alternative setting, where the ancilla is subject to spontaneous decay, the Lindblad operator is given by L = b and does not commute with the projectors. In this case, Eqs. (D3a) to (D3c) receive an additional contribution from the "recycling terms" 2κbσ QQ b † , whose effect is to incoherently remove particles from the ancilla. Since these terms are proportional toσ QQ , which is initially zero and whose growth is suppressed by the Zeno effect, their presence does not change the following line of arguments. Nonetheless, it is possible to get rid of these terms completely by keeping track of all the modes the ancilla decays to and post-selecting on the condition that the ancilla plus these additional modes are empty. To see this, we can assume that the ancilla decays only to a single mode with associated annihilation and creation operators c and c † . The corresponding Lindblad operator L = c † b now conserves the number of particles in the ancilla plus the extra mode. Consequently, the contribution from the "recycling terms" to Eqs. (D3a) to (D3c) vanishes due to the action of the projector P.
We now consider the strong noise limit of Eq. (D3). The terms on the right-hand side of the equations for the coherences (D3b) and (D3c) rotate at characteristic frequencies of the unperturbed Hamiltonian H 0 viã a(t) = e iH0t/ ae −iH0t/ as well as via the modulation function g(t), whose role is to probe dynamic correlations in the system at a given frequency. In contrast, the terms proportional to the dissipation rate κ cause a damping of the coherences. If κ is sufficiently large, in particular, if it is much larger than the characteristic frequencies of H 0 , we can make the approximation that the coherences are instantaneously damped to a momentary equilibrium state given by dσ PQ / dt ≈ 0 (and analogously forσ QP ). This allows us to adiabatically eliminate the fast incoherent dynamics and to solve Eqs. (D3b) and (D3c) for the coherences. To leading order in Ω/κ, we find where (· · · ) −1 denotes the Moore-Penrose pseudoinverse. To leading order in Ω/κ, we can furthermore neglect the terms proportional toσ QQ , which is initially zero and grows, according to Eqs. (D3d) and (D4), only slowly at a rate Ω 2 /κ. This suppression of the growth of population in the ancilla is precisely a manifestation of the Zeno effect. Thus, plugging Eq. (D4) into Eq. (D3a), we obtain with the effective non-Hermitian Hamiltoniañ (D6) Due to the action of the projector P in this expression, the pseudoinverse acts only on states with exactly one particle in the ancilla, where it reduces to a multiplication by unity. Thus, the effective non-Hermitian Hamiltonian takes the simple formH eff (t) = −ig 2 (t) Ω 2 Pã † (t)ã(t)P/κ. Finally, Eq. (29) follows after transforming back to the non-rotating frame.

Appendix E: Connection between ancilla-based weak measurements of dynamical correlations and (non-)Hermitian linear response
Ancilla-based weak measurement schemes for dynamical correlations can be adapted to probe either the unequal-time commutator or anti-commutator through a suitable choice of the ancilla state, the system-ancilla coupling, and the projective measurement performed on the ancilla [43,44]. It has been shown that those variants that probe the unequal-time commutator can be cast into an ancilla-free formulation [44], giving rise, e.g., to rotation-based protocols [30,41,[43][44][45]. For weak perturbations, e.g., small rotation angles, these ancillafree schemes correspond in fact to (standard) linear response. By contrast, the interpretation of ancilla-based weak measurement protocols that target the unequaltime anti-commutator is far less obvious. For instance, Ref. [44] poses the question of whether an ancilla-free measurement of this quantity is possible in general. Here, we show that, indeed, any ancilla-based weak measurement protocol for the unequal-time anti-commutator can be described in an ancilla-free way as a non-Hermitian linear response, exposing the close connection between these frameworks.
To this end, we consider a general ancilla-based weak measurement that uses only projective measurements of standard (Hermitian) operators on the ancilla. The following derivation proceeds in analogy to the one for spin systems presented in Refs. [43,44], but here we consider a more general scenario: we do not specify the type of system, work with general mixed states, and consider arbitrary durations of the system-ancilla coupling. Let us denote the initial state of system and ancilla by ρ S and ρ A , respectively, and assume the combined system to be in a product state initially, ρ 0 = ρ S ⊗ ρ A . The target system evolves under the Hamiltonian H 0 , while we assume the ancilla to have no internal dynamics. System and ancilla are coupled via the general coupling Hamiltonian with a time-dependent function f (t) and Hermitian operators A and X acting on system and ancilla, respectively. The total Hamiltonian of the combined system then reads It is convenient to work in the interaction picture,ρ(t) = e iH0t/ ρ(t)e −iH0t/ . The von Neumann equation can equivalently be expressed in integral form as whereH cpl (t) = f (t)Ã(t) ⊗ X is the interaction-picture coupling Hamiltonian withÃ(t) = e iH0t/ Ae −iH0t/ . In the last line, we have assumed the coupling to be sufficiently weak such that we can replaceρ(t ) in the integral, to linear order in H cpl , byρ(0). Note that the validity of this linear approximation is not necessarily restricted to short times t, but can also be ensured for longer times by a sufficiently weak coupling strength f (t). After a coupled evolution up to time t, during which system and ancilla become entangled, we measure projectively the observable B ⊗ Y , where B and Y are Hermitian operators acting on system and ancilla, respectively, and post-select on the outcome of the ancilla measurement. Although in practice system and ancilla are often measured simultaneously, it is instructive to treat this process as a consecutive measurement of the ancilla first and the system second. Without loss of generality, we assume the observable Y to have a discrete spectrum of (real) eigenvalues { y }. Let P y denote the projector on the eigenspace of the eigenvalue y. After obtaining this outcome, according to Lüders' rule [117], the state collapses toρ where p(y) = Tr [P yρ (t)P y ] is the probability of measuring the outcome y. For the coupling Hamiltonian (E1), the unnormalized post-measurement state reads where h.c. denotes the Hermitian conjugate, while the probability of measuring y becomes (E6) Here, O(t) 0 denotes the expectation value of the Heisenberg operator O(t), evolving under the unperturbed Hamiltonian H 0 , with respect to the initial state ρ 0 = ρ S ⊗ ρ A . Note that expectation values involving only ancilla operators are time independent since we assumed the ancilla to have no internal dynamics. Using , we obtain the normalized, conditional post-measurement state, to linear order in the coupling, as Next, we are interested in the conditional expectation value of the system observable B, given that the measurement of Y on the ancilla yields the outcome y. In a first step, we trace out the ancilla, This yields the conditional expectation value where c.c. denotes the complex conjugate. With this result at hand, we can choose the ancilla state ρ A as well as the ancilla operators X and Y such that the integrand contains either the unequal-time commutator or the anti-commutator of the system observables A and B.
If P y X 0 = Tr[P y Xρ A ] is real, i.e., P y X 0 / P y 0 = −s with s ∈ R, Eq. (E9) gives access to the unequal-time commutator, (E10) This expression coincides with Kubo's linear response formula (cf. Eqs. (2) and (4)) up to a constant factor in the response function. There are two special cases worth discussing. First, if X = 1, P y X 0 = P y is always real and the scheme always yields the unequal-time commutator. This is not surprising: for X = 1, system and ancilla always remain in a product state and the coupling Hamiltonian (E1) corresponds to a Hermitian perturbation on the target system only, which is exactly the linear response scenario. Second, it is instructive to consider the unconditional expectation value B(t) = y B(t) y p(y), which corresponds to not measuring the ancilla at all or disregarding the outcome of the ancilla measurement. By combining Eqs. (E6) and (E9), and using the completeness relation y P y = 1, we find, to linear order, which again always yields the unequal-time commutator. These two examples illustrate two essential ingredients for extracting the unequal-time anti-commutator from ancilla-based weak measurements: firstly, the coupling must entangle system and ancilla, and secondly, it is necessary to correlate the measurement on the target system in some way with the outcome of the ancilla measurement, e.g., through post-selection.
In order to extract the unequal-time anti-commutator from Eq. (E9), P y X 0 must be purely imaginary, i.e., This expression corresponds, up to a constant factor in the response function, directly to the non-Hermitian linear response scenario in Eqs. (9) and (10). Equations (E10) and (E12) demonstrate the fact that any ancilla-based weak measurement designed to probe the (anti-)commutator can effectively be described as a (non-)Hermitian linear response. To make this connection even more explicit, we trace out the ancilla in the unnormalized post-measurement state (E5), In analogy to Ref. [44], to linear order, this result can be re-written in terms of generalized measurement (or Kraus) operators [107,118] as (E15) It is easy to verify that, to linear order, these operators fulfill the completeness relation y M † y M y = 1. The measurement operator M y describes the effect of the system-ancilla coupling, conditioned on the outcome y of the ancilla measurement, without explicitly referencing the ancilla. This ancilla-free description corresponds to the evolution under the effective Hamilto- For real P y X , this evolution is unitary [44] and corresponds to standard (Hermitian) linear response, giving access to the unequal-time commutator. Remarkably, the case of purely imaginary P y X , which according to Eq. (E12) probes the unequal-time anti-commutator, corresponds to an anti-Hermitian perturbation and maps directly to the non-Hermitian linear response scenario described in Section III. This shows that for any ancillabased weak measurement of dynamical correlations there is a corresponding ancilla-free linear response description. Conversely, any linear response protocol can, at least in principle, be realized via an ancilla-based weak measurement. While obvious for standard (Hermitian) linear response, in the non-Hermitian case, the challenge is to choose the ancilla state as well as the ancilla operators X and Y appropriately such that P y X ∈ iR. To see that this is always possible in general, let Y be any Hermitian operator on a (complex) Hilbert space of dimension two or higher with at least two distinct eigenvalues y 1 and y 2 . Then, take the ancilla state to be the equal superposition of the corresponding eigenstates, |φ = (|y 1 + |y 2 )/ √ 2, with ρ A = |φ φ|. Now, let P y = P y1 be the projector on the eigenspace with eigenvalue y 1 and set X = −i(|y 1 y 2 | − |y 2 y 1 |). Then, we have P y X 0 / P y 0 = −i, as desired. Clearly, for a given anti-Hermitian perturbation, the choice of ρ A , X, and Y is not unique, and the challenge consists in finding the configuration that is most convenient for the desired application. Appropriate choices for spin systems have been discussed, for instance, in Refs. [43,44], but their experimental realization on other platforms, e.g., bosons in optical lattices, is unfortunately not straightforward. To illustrate this, assume we are interested in perturbations by the density operator A = n and consider the above scenario of achieving purely imaginary P y X 0 for a bosonic ancilla, which translates to and |φ = (|0 + |1 )/ √ 2. However, neither the superposition of Fock states |φ nor the coupling Hamiltonian H cpl ∝ n ⊗ X, which would be cubic in boson operators, can be realized with massive, non-relativistic particles. More generally, in order to probe unequal-time anti-commutators involving the density A = n through a particle number measurement Y = b † b on a bosonic ancilla, the operator X cannot be diagonal in the Fock basis, as this would imply P y X ∈ R, regardless of the ancilla state. In other words, a particle number nonconserving coupling Hamiltonian would be required in such a setting. In our proposal of Section V, this difficulty does not arise because the leading contribution to the response is quadratic in the coupling Hamiltonian, which enables non-Hermitian perturbations for a wide range of observables including densities and correlators with experimentally feasible system-ancilla couplings. All in all, we have established a general connection between ancilla-based weak measurement protocols for dynamical correlations and (non-)Hermitian linear response theory. In particular, our results pave the road to measuring the left-hand side of the FDR (1), i.e., the unequal-time anti-commutator, via non-Hermitian linear response in an ancilla-free fashion, harnessing the rapidly developing toolbox of non-Hermitian physics [52,53]. mitian operator A with two distinct eigenvalues a 1 , a 2 ∈ R. Let P 1 and P 2 be the projection operators on the corresponding eigenspaces such that A = a 1 P 1 + a 2 P 2 with P 1 + P 2 = 1. This allows us to express both projectors entirely through the operator A and the known eigenvalues, which would not be possible if A had more than two eigenvalues. The protocol starts by evolving the initial state ρ 0 (under the target Hamiltonian H 0 ) to the waiting time t w . Then, the observable A is measured projectively and the state is conditioned on the outcome a 1 or a 2 of this measurement, yielding the conditional postmeasurement states with (i, j) ∈ { (1, 2), (2, 1) }, where we have switched to the Heisenberg picture. Solving for the unequal-time anti-commutator, we obtain The probabilities p(a i ) can be expressed through the expectation value A(t w ) with the help of Eq. (F1), yielding This result states that the desired unequal-time anticommutator of A and B can be extracted from a measurement of the unconditional expectation value B(t) (without previous projective measurement) as well as the conditional expectation values B(t) a1 and B(t) a2 , given that the outcomes a 1 and a 2 have been obtained from the projective measurement of A at the waiting time t w , respectively. A few remarks are in order. The projective measurement of A at time t w can be deferred up to the final observation time t with the help of an ancilla using shelving techniques [30]. This way, the need for non-destructive measurements can be avoided. Furthermore, it is worth emphasizing that there is no restriction on the number of eigenvalues of the operator B, i.e., the dichotomic constraint applies only to A.

Numerical benchmark: projective protocols versus non-Hermitian linear response
We now specialize the projective protocol in Eq. (F5) to density correlations in a Bose-Hubbard system. In the hard-core limit U/J → ∞, multiple occupancies of the same lattice site are prohibited. The density n at site then becomes a dichotomic observable with only two eigenvalues 0 and 1. We thus recover the protocol reported in Ref. [30], − n 2 (t) 0 (1 − n 1 (t w ) ) .
For soft-core bosons, Eq. (F6) does not hold in general since the density operator A = n 1 can take more than two eigenvalues. However, the projective protocols in Eqs. (F4) and (F5) can still be used to measure the exact unequal-time anti-commutator for an arbitrary observable B and any dichotomic observable A. For instance, a possible choice of A is the parity Π of the particle number at site , which in conventional quantum gas microscopes is even more easily accessible than the density itself due to pairwise atom loss caused by the nearresonant imaging light [74]. If we associate the eigenvalues a even = 0 and a odd = 1 with even and odd parity, respectively, the operator Π = a even P ,even + a odd P ,odd coincides with the density in the hard-core limit. Thus, in the regime where multiple occupancies can be neglected, we can approximate the density-density anticommutator in Eq. (F6) by the (exactly obtainable) quantity {n 2 (t), Π 1 (t w )} .
An alternative strategy to approximate {n 2 (t), n 1 (t w )} for soft-core bosons is to take Eq. (F6) literally and compute the conditional expectation values n 2 (t) 0 and n 2 (t) 1 from only those realizations where the projective measurement of n 1 (t w ) yields the outcomes 0 and 1, respectively, discarding realizations with higher occupancies. By contrast, n 2 (t) still represents the (full) unperturbed expectation value. This way, the asymptotic behavior of the unequal-time anti-commutator for t t w is correctly reproduced: two local observables A and B typically become uncorrelated in an ergodic system after sufficiently long times and the anti-commutator reduces to the disconnected product 2 B(t) A(t w ) . For the Bose-Hubbard model, the conditional expectation values of the local densities in Eq. (F6) are expected to eventually re-equilibrate to their unperturbed value n 2 (t) , such that the righthand side indeed becomes 2 n 2 (t) n 1 (t w ) . As long as the system is sufficiently close to the hardcore limit, we can expect Eq. (F6) to reproduce the unequal-time anti-commutator for any t ≥ t w to good accuracy. In what follows, we analyze how well this approximation works for on-site densities (B = A = n ) as a function of the filling and the on-site interaction. In Fig. 10, we compare the performance of the projective protocol in Eq. (F6) to that of the non-Hermitian linear response scheme discussed in Section IV for a 2D Bose-Hubbard system as a function of the filling. To this end, we vary the number of particles N on a square lattice with open boundary conditions consisting of 4 × 4 sites, labeled by a pair of indices ( x , y ) with x , y ∈ { 1, . . . , 4 }. We initialize the system in a single Fock state where the particles are distributed to maximize their mutual distances without initially occupying the interior site (2,2), at which we probe the density correlations. Figure 10a shows the probability p(n) of finding zero, one, or more than one particle at the probe site for U/J = 5 as a function of time. The initial oscillations quickly damp and the probabilities become approximately stationary. In Fig. 10b, we show the probability p(n) = t −1 t 0 dt p(n, t ), time-averaged up to time Jt = 10, as a function of the particle number N . For small N , higher occupancies n > 1 can be neglected and the density operator at the probe site is approximately dichotomic. Figures 10c and 10d show, respectively, the time trace of the unequal-time anticommutator (B = A = n 2,2 ) and the reactive part of the non-Hermitian dynamic susceptibility χ (NH) (correlation spectrum) at the waiting time Jt w = 10 for several values of N . The exact results are compared to those extracted using the projective protocol and the non-Hermitian linear response scheme. For the latter, we have used a rectangular pulse of duration Jδt = 0.01 and a perturbation strength s = 0.05 as in Section IV. The Fourier integral in Eq. (11) has been computed using an exponential filter of characteristic frequency γ/J = 0.2 (see Appendix B 1). While for small N the projective protocol correctly reproduces both the exact time trace and the exact spectrum, there are sizable deviations as the number of particles N (and thus the contribution of higher occupancies) grows. By contrast, the non-Hermitian linear response scheme reproduces the exact results to good accuracy regardless of the filling. In Figs. 10e and 10f, we show, respectively, the L 2 norm of the relative error χ  Fig. 10d. For the projective protocol, both errors increase with increasing particle number, while the errors remain small for the non-Hermitian linear response scheme. For larger onsite interactions U , higher occupancies are suppressed more strongly, which delays the rise of the error curve for the projective protocol as the filling increases: given a certain acceptable tolerance for the relative error of, say, less than 20 %, the projective protocol for U/J = 5 (U/J = 10) yields acceptable results for up to N = 4 (N = 9) particles.
We now investigate the performance of projective protocols for the scenario in Section IV, i.e., a quench in a 1D Bose-Hubbard chain of length L = 12 with periodic boundary conditions at unit filling, initially prepared in a Mott-insulating state. Since n (t) ≡ 1 in this case, the projective protocol in Eq. (F6) reduces to {n 2 (t), n 1 (t w )} ≈ 1 + n 2 (t) 1 . If we evaluate this expression at t = t w , the right-hand side takes the value 2 and therefore the connected anti-commutator extracted from the projective protocol vanishes. This behavior is qualitatively different from that of the true anti-commutator, which is maximal at t = t w . Consequently, Eq. (F6) represents a rather poor approximation of the unequal-time anti-commutator in this scenario, especially at small U/J. To obtain a slightly better approximation, we resort to the projective protocol in Eq. (F4), which is no longer equivalent to Eq. (F5)  Figure 11. Same as Fig. 10, but for a 1D Bose-Hubbard system at unit filling as a function of the on-site interaction U . Despite the suppression of higher occupancies at large U (b), the relative error (e) of the projective protocol (P) remains sizable. By contrast, the non-Hermitian linear response scheme (LR) yields good results for any value of U .
if A has more than two eigenvalues. However, unlike Eq. (F5), the protocol in Eq. (F4) does not reproduce the correct asymptotic behavior of the anti-commutator for t t w if A is not dichotomic. This can be fixed by replacing B(t) (a 1 + a 2 ) on the right-hand side by α B(t) with α = 2 A(t w ) − (a 1 − a 2 )[p(a 1 ) − p(a 2 )]. Since B(t) is usually stationary in the regime of interest, this replacement merely contributes a constant offset to the time trace of the anti-commutator, which ensures {B(t), A(t w )} → 2 B(t) A(t w ) for t t w and avoids spurious static peaks in the correlation spectrum.
In Fig. 11, we present a similar analysis as in Fig. 10 for the 1D system at unit filling as a function of the on-site interaction U . In Fig. 11b, it can be seen be seen that there is a significant contribution from states with higher occupancies at small on-site interactions U/J. As expected, the projective protocol does not perform well in this regime, while the non-Hermitian linear response scheme yields good results. As we move to larger U/J, we enter the Mott-insulating regime where single occupancies dominate and the dynamics is governed by particlehole excitations [119]. Although the probability of higher occupancies p(n > 1) diminishes with increasing U/J, its contribution remains on the same order as that of the probability for vacancies p(0). Thus, the density is nowhere well approximated by a dichotomic observable since the neglected higher occupancies (particle excitations on top of the Mott insulator) are of equal importance as vacancies (hole excitations). This explains why the absolute error of the projective protocol in Fig. 11f decreases substantially with increasing U/J, while the relative error in Fig. 11e decreases only slowly and remains comparatively large even at large U/J. We have checked that the error behaves similarly if we approximate the density-density unequal-time anti-commutator by replacing A = n 1 with the parity Π 1 , as discussed above. Thus, as opposed to non-Hermitian linear response, projective protocols are not well suited for probing unequal-time anti-commutators and the associated FDRs for densities at unit filling.
Our numerical benchmarks suggest that projective protocols generally work well at low fillings and large onsite interactions, where multiple occupancies can be neglected. However, they do not represent a good alternative to measure unequal-time anti-commutators and FDRs in regimes where the relevant observables are not approximately dichotomic. In our example of the Bose-Hubbard model, this limitation unfortunately applies to a major part of the physical parameter space, including the relevant setting of a system at unit filling and moderate on-site interaction strengths. In order to explore these regimes of interest, we must therefore resort to alternative methods like non-Hermitian linear response, which performs well across the entire parameter space.