Quantum wakes in lattice fermions

The wake following a vessel in water is a signature interference effect of moving bodies, and, as described by Lord Kelvin, is contained within a constant universal angle. However, wakes may accompany different kinds of moving disturbances in other situations and even in lattice systems. Here, we investigate the effect of moving disturbances on a Fermi lattice gas of ultracold atoms and analyze the novel types of wake patterns that may occur. We show how at half-filling, the wake angles are dominated by the ratio of the hopping energy to the velocity of the disturbance and on the angle of motion relative to the lattice direction. Moreover, we study the difference between wakes left behind a moving particle detector versus that of a moving potential or a moving particle extractor. We show that these scenarios exhibit dramatically different behavior at half-filling, with the"measurement wake"following an idealized detector vanishing, though the motion of the detector does still leaves a trace through a"fluctuation wake."Finally, we discuss the experimental requirements to observe our predictions in ultracold fermionic atoms in optical lattices.


INTRODUCTION
Many signature effects of classical hydrodynamics have counterparts in quantum systems and serve to provide intuition as well as a spectacular source for interesting new physical situations. Due to the absence of internal scale in hydrodynamics, it can be applied for physical scenarios of vastly different scales. For example, relativistic hydrodynamics has been successfully used to explain collective effects in heavy-ion collisions at RHIC and LHC [1]. On the other hand, studies of hydrodynamic-like effects in strongly interacting electron systems show unexpected effects due to their similarity to viscous fluids. For example, ref. [2] shows that in certain situations, conductance may exceed Landauer's ballistic limit due to viscous effects, while ref. [3] demonstrates that slow "swimming" in a Fermi gas is of a topological nature, and can be fine-tuned to be done without dissipation.
Another interesting example of a hydrodynamics inspired study is the investigation of wake waves produced as a response to a moving potential interacting with a two-dimensional electron gas, recently described in [4]. There, it was pointed out that the pattern formed is determined by a Mach number and has similarities to Kelvin wakes in water and to Mach shock waves following a supersonic projectile. This behavior can be traced to the coherent interference between plasma excitations in the medium, with a dispersion which is waterlike (ω(q) 2 ∝ q) at long wavelengths. A related effect, Cerenkov radiation due to a moving charge in a dielectric has also been studied extensively, most recently in photonic crystals where a host of new variations on the effect have been uncovered where, for example, the direction of radiated energy can be flipped see e.g. [5,6]. * mbw5kk@virginia.edu † ik3j@virginia.edu Here, we consider an altogether different system and, with it, a new set of non-equilibrium problems. We examine the discrete-time steady-state generated by the interaction of different types of disturbances, as described below, with fermions on a lattice, as the disturbances move from site to site. Thus, the discrete time, the lattice and the many-body nature of the system play essential roles in the definition of our model. We find that non-classical disturbances may yield a drastically different response. The case in point is that of a moving quantum particle detector, which in this context has no classical counterpart. In addition, we study a moving particle extraction site, in which particles can be ejected out of the system. These two types of disturbances are compared with results we obtain for a moving potential. We note that in recent years, there have been many investigations of measurement-induced dynamics in many-body systems [7][8][9][10][11][12]. Measurement induced dynamics has been also observed in experiments [13][14][15]. Moving detectors have been considered before as well, most notably in describing the Unruh effect [16], where a uniformly accelerated detector observes a thermal radiation in vacuum. However, the question we consider here is, to the best of our knowledge, completely new: what type of a steady state density pattern will a moving detector leave behind when measuring particle densities in a Fermi sea.
The recent progress of quantum simulation with ultracold atoms [17] makes them an ideal platform for studying these effects. Here, we focus on cold Fermi gases which became a valuable tool in recent years to study non-equilibrium dynamics in analogy to electronic systems. Indeed, recently, increasingly sophisticated techniques became accessible leading to the measurement of spin dynamics [18][19][20] and charge transport [21][22][23]. In particular, it also became possible to observe spin charge separation in one-dimensional lattice systems [24] and to study spin-and charge transport in the two-dimensional Fermi-Hubbard model in the regime of low temperatures and strong correlations that challenge current theory cal-culations [25,26]. These experiments demonstrate that ultracold fermionic atoms are an effective platform for quantum simulation of non-equilibrium phenomena even beyond the capabilities of exact calculations.
Classically, the universality of wakes following moving ships has been characterized by Kelvin's seminal result, that a (gravity) wake behind a moving ship in water is delimited within a constant angle 39 • , irrespective of the ship's velocity [27]. Recent results emphasize finite-size effects [28] through the dependence on the Froude number v/ √ gL of a moving pressure source traveling at velocity v of length L, where g is the gravitational acceleration constant [29]. Here, we study wakes created by point disturbances moving through a Fermi lattice gas including the quantum effects (Fig. 1). Our results depend on both the velocity and the angle with respect to the Bravais lattice directions, as well as on the type of disturbance. Concretely, we consider a tip traveling through a lattice of cold fermionic atoms, interacting with an atom on a lattice site, and then during a time τ moving on to the next site.
We find several unexpected results. For example, we observe a dramatic difference between the wakes of a moving particle detector and a traveling potential disturbance. In particular, on the square lattice at half filling the detector wake vanishes identically due to particle hole symmetry at any temperature. Another surprising result is that, at half-filling, the wake formed by a "particle" extractor is independent of temperature. To find an analytic form for the wake left behind a moving potential we use a co-moving steady-state equation and employ a strategy of identifying nodal lines where the (co-moving) disturbance is exactly zero, in contrast to most treatments of water wakes, which seek for extrema, i.e. troughs and crests. Due to the scale inherent in the lattice structure, our wakes depend explicitly on the time τ characterizing the effective speed of the moving tip, compared to the hopping energy t hop of the fermions in a tight-binding lattice.
To describe these effects of dynamics in many-particle quantum systems we use the non-equilibrium framework derived in [30]. This framework allows for the study of a variety of non-equilibrium problems including particle detection and injection/extraction. It was shown in [30] that in certain statistical mechanics problems, which we detail in the Formalism section, it is possible to make a systematic connection between the evolution of n body density functions with n + 1 density functions, similar in spirit to the Bogoliubov-Born-Green-Kirkwood-Yvon (BBKGY) hierarchy, which is the essential structure leading to the Boltzmann equation for single particle densities from higher order correlation functions (see, e.g. [31]). This approach allows the buildup of tractable nonequilibrium problems utilizing combinations of four elementary operations: detection, particle injection, particle extraction and free evolution. While some of these FIG. 1: A lattice of cold atoms interacting with a moving disturbance. The disturbance can be an applied potential, a detector or an extractor. The blue dots represent the fermionic atoms and the red focussed laser beam illustrates the disturbance moving into the direction of the green arrow.
FIG. 2: Snapshots of the wake developing following a moving particle detector at quarter filling. Plotted are the local particle densities of each lattice site as given by the color bar on the right. The location of the disturbance is marked by a red dot. The simulations are done by iteration of free evolution (Eq. (4)) interspersed with interactions with a particle detector (Eq. (6)) beginning from the ground state of free fermions on a 30 × 30 lattice. The detector was initialized at position (8,15) and moves horizontally to the next site during time 1/(3.4 t hop ), where t hop is the hopping parameter of the free evolution.
ideas have been applied to problems in 1D (e.g. driven and dissipative XX spin chains [32] and steady states of a driven hopping model [30]), here we study an essential 2D problem: the emergence of wakes behind moving objects interacting with a Fermi sea. In particular, we discuss the difference between the motion of a detector, particle extractor, and a potential in detail. The approach of [30] allows for an efficient numerical calculation of the dynamics in such problems. An example of the development of a wake a moving detector is described in Fig. 2, while a comparison between a moving detector and potential is provided in Fig. 3 at different filling fractions of a Fermi sea in a 2d hopping model.  (8,15) and moves horizontally to the next site during time 1/(3.4 t hop ), where t hop is the hopping parameter of the free evolution. Snapshots are taken after 18 time steps, when the full wake pattern has been sufficiently developed. At half filling the difference is most dramatic (top) but differences remain at quarter filling (bottom).
The structure of the paper is as follows. We start by briefly introducing the formalism of [30]. Next we study the effect of a potential hopping from site to site, solving for the characteristic angles of the traveling pattern. We then continue to study the motion of a detector and the motion of a particle extractor and compare these with our results for the moving potential. Finally, we suggest an experimental setup to directly observe the wake patterns.

FORMALISM
First, we provide a formal description of the system depicted in figure 1. We will denote by a r the annihilation operator for a fermion at lattice site r. To describe the density distribution we will focus on the two point function, defined as: where ρ(t) is the density matrix at time t. The evolution of G rr depends on the problem at hand. Due to the discrete nature of the lattice, we find it natural to consider the time evolution in discrete time steps τ , pertaining to the time disturbance moves from site to site. After a step in the evolution process, G(t) → K(G) ≡ G(t + τ ), where K(G) is specified for various processes below. In general, if the system undergoes Hamiltonian evo-lution during a time τ , we have where U = T e − i t+τ t H(s)ds is the many-body evolution. For a general interacting Hamiltonian, G rr (t + τ ) is not determined by G(t) alone and would depend on all higher order correlation functions.
For a non-interacting Hamiltonian, however, the evolution of G does not depend on high order correlations. Let us take the Hamiltonian to be H(t) = rr H rr (t)a † r a r , where H is an N × N matrix if there are N fermion sites. The evolution of G from time t to time t + τ is simplified by the fact that for such a Hamiltonian, where H T (s)ds ] qq is an N × N single particle evolution operator [35]. We therefore find for G: In other words, the matrix G undergoes the evolution G(t) → U G(t)U † , independent of higher correlation functions. A few other operations that yield a closed equation for G are possible and described in detail in [30].
We will use two of the aforementioned operations. We start with the elementary particle detection measurement at a site r. It is described by the following Krauss map of the many body density matrix: wheren r = a † r a r is the number operator associated with site r. Note that for fermions,n r is a projection operator, and the Krauss map (5) describes complete decoherence between the number on the site r and other sites. Substituting (5) in Eq. (1), fermion detection induces the following map on G: where P r = |r r| is the (single-particle) projector on site r and P ⊥ r = I − P r . An additional operation is an extraction event of a particle at site r. Such an operation can be described by the Krauss map where 0 ≤ ≤ 1 describes the efficiency of the extraction procedure. Again, substituting this map in the definition (1) we obtain In this paper, we combine the three types of maps above to represent the dynamics of a disturbance interacting with a lattice as it moves from site to site. We prepare the system at an initial state, with a two point function G(t = 0) = G 0 . The disturbance will interact with the system at position r = 0, and moves to act on an adjacent point aw, where it acts again after which it moves to 2aw and so on, with the disturbance at the nth step acting at position r = naw. We only consider motion at angles where the tip hits the actual lattice sites, i.e. the direction of motion w needs to be a lattice vector.
The evolution of the correlation matrix G between successive time steps is described by the maps given in Eqs (4), (6), or (8)-as elaborated in detail in each section. Let the evolution from time τ n to τ (n + 1) be given by K n (for example, detection at point r = naw using (6), followed by non-interacting evolution for a time τ using (4)). The evolved system at time nτ will therefore have the correlation matrix: and in particular, the local density change compared to the initial density is we follow these dynamics numerically, explicitly affecting the iteration procedure for each of the cases, as explained below, and compare the results to the co-moving steady state which we now define. For a moving disturbance, a steady state can only be formed in the co-moving frame. Consider an elementary operation G → K(G) (for example detection, followed by non-interacting evolution for a time τ ), which is then repeated, but shifted in space by the vector aw. Let S be the translation operator along the direction of motion w, via S † r = r + aw. We then define a steady state for the correlation matrix in the co-moving frame via the requirement that namely G steady is invariant under the combination of the operation K and moving to the next site. To identify relevant steady states, we will seek solutions to (11) in the vicinity of states associated with an unperturbed system. Indeed, as we shall see, the nature of the steady states depends both on the form of the dynamics K as well as on the initial background state (for example, a Fermi gas at different filling fractions). We note that, for simplicity, we focus here on initial states suitable for non-interacting systems. We emphasize that the formalism is valid also for systems prepared in an interacting state as long as the subsequent dynamics is well approximated by non-interacting evolution.

A MOVING POTENTIAL
Consider a tip traveling along the lattice, in a direction w taking a time τ to move between two sites. We approximate this process as a discrete process, where a potential V hops from site to site, remaining a time unit τ at each site. For the purposes of this paper, we will focus on the simplest case of a square lattice with nearest neighbor hopping as the free evolution, i.e. H 0 = −t hop r,r a † r a r with single particle energies where a is the lattice spacing. We will take the tip potential at a fixed reference point r 0 to be V = V a † r0 a r0 . We will mostly concentrate on half filling in this section.
We describe below the wake formed behind a point potential moving at a general speed and angle with respect to the lattice. We begin by summarizing the main results of this section before providing full derivations. Fig. 4 shows the simulation of the wake pattern formed by evolving the system in real time following a successive application of the tip along a horizontal line moving at various speeds. Fig. 5 represents the simulation of the wake formed by similarly evolving the system in real time except with the tip moving at several different angles with respect to the lattice. Denoting we use Eq. (11) to find that the angles of lines of zero disturbance are described by These "zero disturbance" lines are represented as red lines in the figures, and delineate the shape of the wake openings. As expected, since we are not in a Kelvin regime, the angle depends on the speed of the disturbance and is discussed below. Note, in contrast to the classic Kelvin wakes and potential wakes in a two-dimensional electron gas [4], here the lattice breaks rotational symmetry and the wake pattern changes as the potential path rotates with respect to the lattice. Before moving on to the derivation, let us comment briefly on the limiting behavior of Eq. (13). Note, that as α → 0, we find ry rx → ±1, i.e. the two main diagonal directions. This result is consistent with the expectation that as the velocity vanishes, the moving potential is almost static and will radiate via the underlying D-wave symmetry of the lattice.
As α → ∞, i.e. the limit of an extremely fast moving tip, we find that ry rx → wy wx , in other words, the wake converges onto a line following the disturbance, as any disturbance would not have time to disperse. Hence, Eq. (13) implies that the wake will essentially vanish for a potential moving at α → ∞. Red lines represent the angles given by Eq. (35). Each line corresponds to the solution for a given quadrant in Fig. 6. Note the forward pointing cone is a result of a forward dwave radiation when the source is moving slowly.
FIG. 5: Density plots for varying the angle of a moving potential compared to the lattice vectors. A potential moving at 0, 23, and 45 degrees with respect to the lattice. Here, α = 1.7. A smeared potential (Gaussian with half a lattice spacing width) is used instead of a point potential to include effects when the tip is not precisely on a lattice site. It is argued in the appendix that the wake geometry is unaffected by this change away from a point potential. Red lines represent the angles given by Eq. (35). Note that for motion at 45 degrees the angles for quadrants 1 and 4 in Eq. (35) coincide reducing the number of lines to 3.
The co-moving steady state to be found for our system is described by, Eq. (11), where K(G) = e iτ (H0+V ) Ge −iτ (H0+V ) where H 0 = −t hop r,r |r r | is the unperturbed single particle Hamiltonian and V = V |r 0 r 0 | is the tip potential at some initial reference point r 0 . Namely: In general, equation (14) admits infinitely many solutions for G. In particular, any correlation matrix G that satisfies: will automatically be a co-moving non-equilibrium steady state. In the physical scenario we are interested in, however, we have an initial reference state, the correlation matrix G 0 of the unperturbed Fermi system, and in the following we consider the wake as a weak perturbation on this state, allowing us to analytically establish the dominant behavior of the wake pattern.
Since we are only perturbing the free evolution by a small potential, we can assume the steady state G will be close to the steady state of free evolution, G 0 , where G 0rr = a † r a r equilibrium . Thus, we write G = G 0 + δG where δG is assumed a small perturbation. Since H 0 is translation-invariant, we write the co-moving nonequilibrium steady state equation in momentum space as: Substituting lowest order perturbation theory, keeping terms up to linear order in V and δG, we find that the real space density disturbance, at zero temperature is given by: where: and ε f is the Fermi energy. Derivation details can be found in the appendix.
Our main objective now is to compute the large-scale features of the resulting pattern, namely the typical angle that appears in the wake pattern. As in the case of the original Kelvin wake, which is typically derived by a stationary phase method, the present treatment requires careful consideration of the dominant contribution to the density variation [Eq. (16)]. The terms A and R in Eq. (16) will provide us with regions that are particularly important for the integral over k and k . Due to the Fermi functions, we can write Eq. (16) as: Note that |A(k, k )| < 1 (this follows from the inequality |e iθ − 1| ≤ |θ| ) and is dominated by k,k near ε(k) − ε(k ) = 0. We thus see that in contrast to the measurement and extraction wakes considered next, the integral is dominated by momenta near the Fermi surface since we can take such momenta to satisfy both conditions τ (ε(k) − ε(k )) 1 and ε(k) > ε f and ε(k ) < ε f . We will henceforth consider the situation at half filling. Looking at the Fermi surface for our system, we break up the expansion around the Fermi surface into four quadrants given in Fig. 6. Close to the Fermi lines, we will use the variables δ y and δ y instead of k y , k y as the small shifts away from the Fermi surface. Explicitly, Let us now concentrate on R in Eq. (16). This term diverges when for n integer. Here we concentrate on the n = 0 contribution which already recovers some basic features of the wake pattern, and leave the analysis of n = 0 contributions for a future work. The equation can also be interpreted as a Mach-Cherenkov-Landau condition [33] for the momenta emitted by the wake due to creating a particle-hole excitation of momentum K = k − k . Perhaps a more familiar way to write the condition is: For the square lattice, we have where α is defined by Eq. (12). Now, combining the restriction τ (ε(k) − ε(k )) 1 with Eq. (23), we find that where δ x is given by δ x ≡ τ (ε(k)−ε(k )) awx and where w ≡ wy wx . Comparing to Eq. (20) we arrive at: where δ x ≡ w(δy−δ y )+δx 1+(−1) b w . Here, b = 1 for quadrant 1 and 4 (Fig. 6). Otherwise, b = 0. Note that our treatment of δ x as small breaks down when w is close to 1. Indeed, when w x = w y , the constraints on energy together with Eq. (23) are insufficient to force k and k to be close, since (k − k ) can be large, with both k, k close to the Fermi surface, and (k − k ) perpendicular to the w vector. While more refined analysis is needed to describe the special point w x = w y exactly, here we simply observe numerically that our treatment works well for w x < w y and w y > w x , and that the wake pattern change is continuous at w x = w y and is well described by our Eq. (13).
We now combine Eqs. (25), (20), and (23), and expand in small δ x , δ y , δ y to second order. Solving those, we can relate δ x and δ y to δ y and solve for k x − k x and k y − k y as If we assume our potential is at site r 0 = (0, 0), the term e i(r0−r)·(k−k ) in Eq. (19) becomes where B =

2[wxα+sin (akx)]
(wx+(−1) b wy)α[(wx+(−1) b+1 wy)α+2 sin (akx)] (29) × r x (w y α + (−1) b+1 sin ak x ) + r y (w x α + sin ak x ) greatly reducing our momentum space integral to two coordinates, k x and δ y . Note that when the denominators in A, R in Eq. (19) vanish, the leading behavior of the combination AR is real, we arrive at: where Q is the set of four quadrants in Fig. 6. We have checked numerically that the integral (30) indeed captures the main wake pattern of the moving potential well. Our next task is to use Eq. (30) to find the main wake angles. We now estimate analytically the main angles involved in the wake pattern left behind the moving potential. In the case of water wakes, we are interested in the wavefronts, which are lines of maximal disturbance. Here, we find that a more direct approach is to look instead for lines of zero disturbance, i.e. r|δG|r = 0. We will begin by looking at the effects of individual quadrants in Eq. (30). Integrating over δ y and looking first at quad-rant 1, we find To find the characteristic wake lines, we now look for directions r such that Eq. (31) vanishes. Assuming that we could treat the equation by a stationary phase method, a condition for Eq. (31) vanishing would be that there exists a k 0 such that k x B ≈ (k x − k 0 ) 2 . In this case, using the stationary phase approximation around k 0 makes the dominant contribution to the integral in Eq. (31) real, and r|δG|r vanishes. Specifically, for this to happen, we need B = 0 and d dkx B = 0 when evaluated at k 0 . Looking at r x , r y 1, i.e. far away from the potential, the dominating behavior of B (Eq. (29)) comes from r x (w y α + (−1) b+1 sin ak x ) + r y (w x α + sin ak x ). Hence, we find the equations: r x (w y α + (−1) b+1 sin ak x ) + r y (w x α + sin ak x ) = 0, (32) and (−1) b+1 r x cos ak x + r y cos ak x = 0.
Therefore, cos ak x = 0 implying k 0 = π 2a . Plugging this into Eq. (32) yields Repeating this calculation for the other three quadrants, we find and hence our main result Eq. (13). Figures 4 and 5 show agreement between the simulations of the potential wakes and our Eq. (35). While the above treatment was perturbative in V , our analytical treatment for the angles should hold asymptotically at large distances from the source. This is because, in that regime, the response to the disturbance is weak regardless of the strength of the perturbation. Close to the source, the density profile will, in general, not be linear in the strength of the perturbing potential. In particular, at half filling, the requirement δG G 0 implies that r|δG|r 1 2 . To see the range of validity of the description close to the source, let us consider r|δG|r as expressed in the integral (31). We note that B grows linearly in |r|, and therefore, asymptotically r|δG|r oscillates and decays at least as fast as |r| −1 (consistent with the numerical observations). The expression (31) shows that the density profile will be perturbative when V τ a |r| 1, or, in other words at distances |r| a V τ .

MOVING PARTICLE EXTRACTOR AND MOVING DETECTORS
We proceed to consider a moving detector or particle extraction from the system. Note, that these processes are non-unitary. In this section we establish the co-moving steady-state for this problem. In particular, we show that in marked contrast with a moving potential, a moving detector at half filling does not generate a wake because of particle-hole symmetry.
We assume that the detection or extraction process is dominant when the tip is at a given site, but quickly weakens as the tip moves away from that site. It is therefore natural to discretize the process in such a way that we have a disturbance at a given site, followed by a free evolution of the system during a time τ that the tip is traveling to the next site on it's trajectory. The appropriate transformation rules K(G) for detection and extraction are given in Eq. (6) and (8) respectively. If we allow for pure detection to happen with probability p (associated with the efficiency of the detector) and similarly extraction protocol with probability q, we can combine them, together with the free evolution U = e −iτ H0 into the general form which can be written as: where ξ = 2p + 2 q, γ = p + q, {G, P } ≡ GP + P G indicates the anti-commutator, and P is the projection onto a site r 0 where the tip acts. In particular, pure detection will be described by q = 0, hence γ = p and ξ = 2p.
In the next sections we work under the assumption that p, q 1 and hence γ 1. The co-moving steady state equation (11) now reads: Written explicitly in momentum space we have: Assuming that γ 1, G ≈ G 0 + δG, with δG being a small correction, and zero temperature, we find that the local density variation is given by: where ρ f is the density of fermions in G 0 (i.e. the diagonal of the G 0 matrix). Note, that like in the potential case, the term R in Eq. (39) implies that Eq. (23) still characterizes a dominant region for the integral. However, Eq. (39) has two added difficulties when compared with the moving potential case Eq. (16). First, we no longer have the helpful constraint that τ [ε(k) − ε(k )] ≈ 0. Secondly, ε(k) and ε(k ) can now be on the same side of ε f as well as on the opposite side. Nonetheless, in numerical experiments we have observed that the geometry for a moving potential, Eq. (35), does appear to also match with the wake patterns of a moving detector and extractor, as can be observed, e.g., by comparing the wake pattern Fig. 4 to the extractor pattern Fig. 8 . By iterating the evolution equation for the two point function G, the wake pattern can be generated numerically. For the case of a particle removal site moving through a half-filled Fermi sea we obtain the images shown in Fig. 7 and 8. The geometry of the wake patterns is similar to the ones described for the moving potential Fig. 4, however the density variation is always negative due to the depleted particles.

Moving detector at half filling
Particle detection at half filling shows marked contrast with the density wake due to a moving potential. Indeed, due to particle hole symmetry it leaves the average density profile, namely the diagonal of G unchanged. On the other hand, a potential perturbation breaks particlehole symmetry and generates the density wake described above.
In fact, we can establish a stronger property, namely: where r|δG|r µ is obtained by successive applications of K from Eq. (36) on an initial state G 0 = 1 1+e β(H−µ) , with q = 0, using an arbitrary choice of measuring site at each step, and in the end subtracting G 0 . In other words, changing the sign of the chemical potential changes the sign of the wake.
An immediate consequence of Eq. (40) is that at the point, where our many-body Hamiltonian has particlehole symmetry, namely µ = 0 i.e. at half filling: i.e. there should be no wake pattern created by a moving detector. This is shown in Fig. 9 by comparing a detector moving through a half-filled versus a quarter-filled Fermi sea. The image also shows how the quarter filled wake is opposite in sign to the wake generated in the Fermi system at three quarter filling. A full non-perturbative proof of the remarkable relation Eq. (40) is presented in the appendix. Here for simplicity we establish Eq. (40) starting from the (zero temperature) perturbative result Eq. (39) with F = µ and ξ = 2γ. Note that the sum r|δG|r µ + r|δG|r −µ is given by Eq. (39) with the term in brackets replaced with 2−Θ(µ−ε(k))−Θ(µ−ε(k ))−Θ(−µ−ε(k))−Θ(−µ−ε(k )) (42) where we used that ρ f ( F ) + ρ f (− F ) = 1. Now consider the following map reflecting points about the Fermi surface: A "fluctuation" wake The above results suggest at first glance that there is no effect of the detector at half filling. In fact, this is not the case! While the moving detector does not affect the average density at half filling, it does perturb correlations, and thus may be observed through fluctuations. For example, such correlations may be observed by looking at the number of particles in a mode FIG. 9: Density plots of a detector moving through a Fermi sea at several filling fractions. From left to right, the detector is moving at α = 1.7 for quarter-filling, half-filling, and three-quarter-filling, respectively. The color bar has the same scale but different offset (centered around the initial filling fraction).
representing an equal weight superposition in a region lattice neighborhood A of a point r. We have: We will focus on A being the set of nearest neighbors: an example of the wake in the density of the n A (r) is then shown in Fig. 10. That this wake may be non-zero can be observed by generalizing Eq. (39) for off-diagonal elements. In this case, the only change to (39) is that e −ir·(k−k ) → e −i(r·k−r ·k ) . Now, combining Eqs. (39) and Eq. (44), we find where κ = {0, k x , k y , −k x , −k y }. Note, Eq. (45) is not symmetric under a reflection of ε(k),ε(k ) about the Fermi energy. This implies that, unlike the diagonal of G rr , the wake generated for the modes like A for a detector are non-zero.
A couple of remarks are in order. (1) While we focused here on the density of the a † A modes as an indicator of correlations, a more natural quantity for an experimental consideration is density-density correlations, and number fluctuations in the region A. A preliminary check shows that such number fluctuations will also exhibit a wake, which can be studied by considering the four-fermi correlation generalization of Eq. (4) and (6), which give a closed hierarchy of 4 point functions. This calculation will be presented in a future work.
(2) A moving detector at half filling is an interesting example of a "hierarchy" of steady states. In this hierarchy, the local average density or "diagonal" of G at half filling is steady for any path a detector makes, and is thus in a steady state. However, the correlations depend on the trajectory of the detector and would, in general, not be in a steady state, moreover the many-body density matrix would not be in a steady state. It is not hard to construct examples where G is in a steady state, while the many body density matrix is time-dependent.

FINITE TEMPERATURE STATES
In this section we analyze the effect of a non-zero temperature of the system on our moving disturbances. We assume that the system is prepared initially at finite temperature, and we neglect thermal dissipation on the time scale of the motion of our disturbances. We find that at a generic filling the amplitude of the wakes are decreased, as may be expected on general grounds, i.e. with increased density fluctuations of the background. These results are shown in Figs. 11 and 12. Furthermore, we find that, at ρ f = 1 2 , a moving detector continues to produce no wake at finite temperature. Perhaps the most surprising effect we find is that the extractor wake at half filling is temperature independent. This behavior is striking when compared to the moving potential source, see Fig. 11.
Thus, the steady state of δG for a moving potential source, Eq. (16), becomes while for detection/extraction at finite temperature, Eq. (39) becomes We can understand the temperature independence of the moving extractor at half filling as follows. Consider the difference between the moving extractor and moving detector steady state equations (Eq. (47) with ξ = 2γ and π/a −π/a dkdk R(k, k , w)e i(r0−r)·(k−k ) ρ f Note, Eq. (48) depends only on the density ρ f . Therefore the difference Eq. (48) is independent of temperature if temperature is varied at a fixed density. Since detection creates no wake when ρ f = 1 2 at any temperature, this implies that a moving particle extractor is temperature independent at ρ f = 1 2 . The result above, Eq. (48), that the difference between the detector and extractor wakes is temperature independent has been done perturbatively in δG for illustration purposes. In fact it is possible to establish that where G det ,G extr are the non-perturbative steady states for a moving detector and extractor respectively, as derived in the appendix. This result matches simulations of the moving particle extractor as shown in Fig. 11. We note, in passing, that numerical checks show that the fluctuation wake is temperature dependent at half filling, even though there is no detector wake.

DISCUSSION OF EXPERIMENTAL REALIZATIONS
Here we consider a setup where the wakes may be explored in experiments with ultracold 6 Li fermions in optical lattices. Quantum gas microscopes with singleparticle and single-site resolution can directly observe the wake structure, as follows. The disturbance can be created by a focused laser beam with a waist on the order of the lattice spacing by employing a high-resolution objective [34]. Experimental system sizes of more than 30x30 lattice sites have been realized for fermions [26]. We have checked by numerical simulations that using a Gaussian smeared potential instead of a point potential does not significantly alter wake geometry. Moreover, the wake pattern is not significantly changed if the initial Gaussian is not centered on a lattice site.
Because of the light mass of 6 Li the timescales for Hubbard physics are still convenient for lattice spacing of approximately 1 µm [24] leading to moderate requirements of the objective (NA > 0.5). A dynamically movable disturbance can be implemented, for example, by piezo-actuated mirror mount or an acousto-optical deflector. An experimental run will then start with the preparation of a two-dimensional Fermi-Hubbard system and the movement of the focused laser beam through this system up to a certain position. Finally, the system is frozen by increasing the lattice depth and then imaged via fluorescence imaging. A single realization does not contain enough information to extract the details of the wake pattern. Reaching a density resolution of about 2% requires averaging about 2500 experimental realizations (1/ √ n) while keeping track of the final position of the disturbance. The required precision and amount of data is comparable to recent experiments at existing quantum gas microscopes [24,34]. The parameter τ is proportional to the duration of beam motion, which should be compared to the hopping energy t hop . It is possible to swipe the beam at different rates to obtain values of α spanning the full range from the slow-moving D-wave like wakes to the disappearance of the disturbance at high speeds. We remark that in a finite optical lattice setting, the large scale wake pattern may not have enough time to develop if the speed is such that the co-moving steady state cannot be effectively reached.
We note that all three types of disturbances can be implemented in experiments. A moving potential can be created by a far-detuned laser beam. A moving detector can be realized by a near-resonant laser beam. Scattering of photons at low intensity leads in good approximation to a measurement of the on-site particle density. Last, a particle extractor can be implemented via a defocussed red-detuned optical dipole trap. Caused by the out-ofplane minimum in the potential, atoms will be sucked out-of-plane and will be lost in the experiment.
Beyond cold atoms, we expect the effects we predict to also hold in other systems which can be well described by non-interacting fermions. We note that, for spinpolarized fermions contact-interactions are suppressed at low temperatures and dynamics is dominated by the effects of Pauli exclusion, making our approach particularly effective in this case. Moreover, we emphasize that our treatment is essentially exact, and thus can provide a benchmark for studies of the perturbative effect of interactions.
It is also important to note that while the present discussion is focused on non-interacting systems, the formalism presented in [30] is valid also for systems prepared in an interacting state, as long as the subsequent step of unitary evolution while the tip is traveling between sites is well approximated by non-interacting evolution. Thus, a system prepared in a strongly correlated state, such as a Mott insulator, for example, that undergoes a quantum quench where interactions are turned off will still be described by the current method.
The discrete-time co-moving steady state equation on the lattice G = S † K(G)S, Eq. (11), states that as the potential (or another type of disturbance) moves to the next lattice site, G remains invariant, up to a shift of the coordinates co-moving with the disturbance. Written explicitly in coordinate representation, Eq. (11) reads G rr = [S † K(G)S] rr , or, equivalently, using that for momentum states, S |k = e −iaw·k |k , we have Let us write G = G 0 + δG, where G 0 is the initial steady state before turning on the traveling perturbation. Using the form (57), taking the Hamiltonian H between steps to be of the form H 0 + V , the steady state equation (58) is explicitly given by: Note that at this no approximation has been made up to this point. Since we are only perturbing the free evolution by a local potential, we can assume the steady state G will be close to the steady state of free evolution, G 0 , and thus δG is assumed a small perturbation. To proceed, we now use perturbation theory in Eq. (59), by expanding to lowest order in V . Namely, we use the expansion : and keep terms up to linear order in V and in δG. The resulting equation is: Next we note that without the perturbation, G 0 is assumed to be a steady state of free evolution. Since H is transnational invariant, G 0 can be taken to be diagonal in momentum space, therefore we can set e iaw·(k−k ) − 1 k|G 0 |k = 0. Doing the s integral, we find At zero temperature where vol is the system volume. Thus, Finally, density variation in the co-moving non-equilibrium steady state is given by Fourier transforming (66) back to real space and taking the diagonal element, In experiments the point potential may be realized by a broader beam. As we are looking for far field wake patterns, we do not expect local structure of the disturbance to have a significant effect. For example, if the point potential was instead Gaussian, the term e i(r0−r)·(k−k ) in Eq. (67) would be replaced by e − σ 2 2 |k−k | 2 +i(r0−r)·(k−k ) where σ 2 is the variance of the Gaussian potential. Note, however, that when calculating the wake geometry we look far away from the potential source (|r 0 − r| large) and also need only to consider terms where |k − k | is small since these k,k dominate the integral. Hence, we can safely neglect the e − σ 2 2 |k−k | 2 term so long as |r 0 − r| >> σ and thus find that the wake geometry of the Gaussian potential is equivalent to that of the point potential. In Figs. 13 and 14 we simulate potential wakes for a selection of Gaussian potentials and find that the wake geometry is indeed equivalent to that of the point potential.
(2) Moving detection/extraction. The co-moving steady state equation reads G = S † K(G). Written explicitly in momentum space, using (36), we have: Assuming that γ 1, G ≈ G 0 + δG, and zero temperature, Hence, Now turning to the P GP term, where ρ f is the density of fermions for G 0 . Plugging Eqs. (70) and (71) into Eq. (68), we get Finally, the local density variation is given by: which is Eq. (39) in the main text.

Non-perturbative Results
In this section, we show that no detection wake is created at ρ f = 1 2 and that the difference between detection and extraction is temperature independent even non-perturbatively.
We start by looking at a series of non-perturbative detections on G 0 . A single detection at site r and evolution for time τ is G = U P r G 0 P r U † + U P ⊥ r G 0 P ⊥ r U † (74) ≡ a={0,1} P a r (τ )G 0 P a r (τ ) since [U, G 0 ] = 0 and where P 0 ≡ P , P 1 ≡ P ⊥ , and U P U † ≡ P (τ ). Hence, after doing a series of m measurements, we have G = a1,a2,...,am n=m,m−1,...,1 P an rn ((m − n + 1)τ ) G 0 n=1,2,...,m P an rn ((m − n + 1)τ ) (75) Looking at the diagonal of G in real space and inserting a resolution of identity, dq n |q n q n |, to the right of every P an rn sitting in the first term in brackets in Eq. (75) and inserting dq n |q n q n | to the left of every P an rn sitting in the second bracketed term in Eq. (75) we find r|G|r ≡ ζ m (µ) = dkdk e −ir·(k−k ) dQdQ a1,a2,...,am k|P am rm (τ )|q m q m |P am−1 rm−1 (2τ )|q m−1 ...
FIG. 14: A Gaussian potential at half-filling, α = 1.7, and 0.5a standard of deviation where a is the lattice spacing. Left is a Gaussian starting directly on lattice site (8,15). Right is a Gaussian starting in-between lattice sites at (8.3,15.6).
i.e. the detection wake for a chemical potential of µ is one minus the detection wake for a chemical potential of −µ. Hence, when µ = 0 there is no detection wake. We emphasize that this result assumed no particular path for the moving detector.