{'fraction_non_alphanumeric': 0.049303277832327075, 'fraction_numerical': 0.010800760440386905, 'mean_word_length': 4.024984253621667, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical Friedel oscillations of a Fermi sea

We study the scenario of quenching an interaction-free Fermi sea on a one-dimensional lattice ring by suddenly changing the potential of a site. From the point-of-view of the conventional Friedel oscillation, which is a static or equilibrium problem, it is of interest what temporal and spatial oscillations the local sudden quench will induce. Numerically, the primary observation is that for a generic site, the local particle density switches between two plateaus periodically in time. Making use of the proximity of the realistic model to an exactly solvable model and employing the {Abel regularization} to assign a definite value to a divergent series, we obtain an analytical formula for the heights of the plateaus, which turns out to be very accurate for sites not too close to the quench site. The unexpected relevance and the incredible accuracy of the Abel regularization are yet to be understood. Eventually, when the contribution of the defect mode is also taken into account, the plateaus for those sites close to or on the quench site can also be accurately predicted. We have also studied the infinite lattice case. In this case, ensuing the quench, the out-going wave fronts leave behind a stable density oscillation pattern. Because of some interesting single-particle property, this dynamically generated Friedel oscillation differs from its conventional static counterpart only by the defect mode.


I. INTRODUCTION
Quantum nonequilibrium dynamics is now under extensive study [1][2][3], as theorists are quenching all kinds of models in all kinds of protocols. As a fair observation, in most cases, the model taken up is a many-body system with interaction. Although models of this kind have to be the choice for addressing problems like thermalization, their actual wide employment might to some extent result from the prejudice that only interacting many-body systems could yield novel, interesting dynamics. But actually, when it comes to dynamics, even the simplest model can be nontrivial [4,5]. The point is that, there is a big gap from the static property of a system to its dynamics, as the extra dimension of time makes the latter much more richer and more complicated than the former. As examples of interesting dynamics exhibited by simple models, we have the Landau-Zener-Stueckelberg tunneling of a two-level system [6][7][8], the nonspreading wave packet of a free particle [9,10], the dynamical localization of a particle on a one-dimensional tight binding chain [11], the coherent destruction of tunneling of a particle in a double-well potential [12], etc. Moreover, in the field of atomic physics, many intriguing phenomena, such as above-threshold ionization [13,14] and stability under super-intense lasers [15], can be displayed by the simplest atom, namely, the hydrogen atom. Hence, in search of interesting dynamics, one does not need to resort to complicated many-particle systems-single-particle models or * liu yu@iapcm.ac.cn many-particle models without interaction could suffice.
In this spirit, we had studied the transition dynamics of a Bloch state on a one-dimensional tight binding chain [16,17]. The scenario is extremely simple. We just take a one-dimensional tight binding chain with the periodic boundary condition and put a particle in an arbitrary Bloch state [see Fig. 1(a)]. Then we quench it by suddenly changing the potential of an arbitrary site. Because of the newly introduced local barrier, the particle will be scattered into all other Bloch states, in particular, the one which is the mirror of the initial Bloch state, and its wave function in the real space will be deformed. This is all one can anticipate from general principles. However, down-to-earth numerical simulation reveals that the probability of the particle remaining in its initial state (the survival probability) and the probability of it being reflected into the mirror state (the reflection probability) both show cusps periodically in time [see Fig. 2(a)]. This kind of nonsmooth behavior is totally unexpected. Yet more is in store. In the real space, the evolution trajectory of the probability density on a generic site shows plateau structures [see Fig. 2(b)]. All these unexpected, structured dynamical phenomena can be explained by identifying and studying an ideal model, which we shall sketch below.
It is natural to generalize the single-particle problem to a many-particle one. After all, to single out a particle and prepare it in a certain Bloch state other than the ground one, is no easy task. Experimentally, people more often deal with a collection of particles simultaneously. Hence, in this paper we consider the same scenario as before but with the initial state being a Fermi sea, in which all the Bloch states below the Fermi energy are Initially an arbitrary Bloch state with wave vector qi is occupied. (b) The many-particle case. Initially all the Bloch states with wave vector |q| ≤ q f are occupied.
filled [see Fig. 1(b)]. From the point of view of Friedel oscillation, this is a violent way of introducing the defect into the Fermi sea. It is of interest what temporal and spatial oscillation it will induce. To make things tractable and to adhere to our idea mentioned above, we assume that there is no interaction between the particles. Despite this simplicity, the problem is still very challenging analytically. But at least numerically, it is simple. As there is no interaction, we can just evolve each particle independently, and then for calculating single-particle quantities like the particle number at a site, we just need to sum up the contributions of each particle. It turns out that the particle number at a generic site switches between plateaus of two different heights, with the times of jumping determined by the Fermi wave vector q f . Determining the heights of the plateaus is then our primary aim.
In the following, we shall first describe in detail the quench scenario and the numerically observed plateaus in Sec. II. Then in Sec. III, we shall try to develop an analytic formula for the plateaus. This will be done in three steps. First, we shall review the single-particle case in Sec. III A, and derive many results in a new way. Then, in Sec. III B, by approximating the realistic Fermi sea by a fictitious one, we employ the Abel summation trick to get an analytic formula for the heights of the plateaus, which turns out to be very accurate for sites not too close to the quench site. Motivated by this success, we do some formal extension in Sec. III C, and obtain another analytic formula which reduces to the Abel one for sites far away from the quench site and has the advantage that it is defined also for the quench site. Finally, in Sec. III D, we take into account the contribution of the defect mode, and get a formula accurate for all sites, including the quench site. This is the ultimate result for the plateau heights. On the basis of this result, in Sec. IV, we consider the infinite lattice case. In this case, there is no repeated plateau switching, and for each site, the particle number will settle down to a constant value. The picture is that the sudden quench generates two out-going wave fronts, which leave in their wakes a stable density pattern. This is Friedel oscillation generated dynamically. We shall discuss its relation with the conventional static Friedel oscillation.

II. THE SCENARIO AND THE PLATEAUS
The setting is simply the one dimensional tight binding chain with the periodic boundary condition, or a lattice ring. The unperturbed single particle Hamiltonian is ( = 1 in this paper) where |l denotes the basis function (the Wannier function) at site l. Here for notational simplicity, it is assumed without loss of generality that the lattice size N = 2L + 1 is an odd number. The periodic boundary condition means that |l ≡ |l + N for arbitrary l. As is well known, the unperturbed eigenstates of H 0 are the so-called Bloch states, which are essentially plane waves on the lattice ring. Explicitly, they have the expression Here k is an integer indexing the Bloch state and q = 2πk/N is the corresponding wave vector. Apparently, |k ≡ |k + N . It is easily verified that the eigenenergy associated with the Bloch state |k is ε(q) = −2 cos q. This dispersion relation is illustrated in Fig. 1. In Refs. [16,17], the following quench scenario [see Fig. 1(a)] was studied. Initially, a particle is in some Bloch state |k i with wave vector q i = 2πk i /N . This is an eigenstate of H 0 , and nothing interesting happens. Then at t = 0, the system is quenched suddenly by changing the potential of some site, which, because all the sites in a ring are on equal footing, can be assumed to be the l = 0 site. The subsequent evolution of the single-particle system is then controlled by the final Hamiltonian H f = H 0 + H 1 , with the perturbation where U is the strength of the defect potential, or the quench amplitude. It turns out that this scenario yields very interesting and unexpected dynamics. Let us denote the wave function of the particle as |k i (t) ≡ exp(−iH f t)|k i , where k i denotes the initial state and the argument t indicates how long it has been evolved by H f . This way of notation is used throughout the paper. The state |k i (t) can be easily determined once the Hamiltonian H f is diagonalized numerically. It is then observed that, in the momentum space, both the probability of the particle remaining in its initial state |k i , . That the survival probability shows cusps has also been observed in some other models [18][19][20][21][22][23][24]. and the probability of it being reflected into the momentum-reversed Bloch state | − k i , show cusps periodically; while in the real space, the local probability density, at a generic site l, jumps constantly from plateau to plateau. These two phenomena are illustrated in Fig. 2(a) and 2(b), respectively. An important question is at what times the cusps or the sudden jumps occur. The answer is that the cusps are located at t = rT , where r is an integer and T = N/v i , with v i = 2 sin q i being the group velocity of a wave packet with wave vector q i , and the sudden jumps are located at t = rT ± s c , with s c = |l|/v i . These results are predicted by the ideal model in Refs. [16,17] and agree with numerics very well.
In the present paper, we study the scenario of locally quenching a Fermi sea, as illustrated in Fig. 1(b). Now instead of a single Bloch state, all the Bloch states with wave vector |q| ≤ q f = 2πk f /N are occupied by spinless, noninteracting fermions. The total number of particles is then M = 2k f + 1, and the particle number at each site isn = M/N . The total wave function of the system is always a Slater determinant. Formally, we can denote it as Before proceeding, let us introduce another orthonormal basis. Both H 0 and H f are invariant under the reflection with respect to the defect site l = 0. Hence, it is convenient to have a set of basis states with definite parities. We thus form the even-parity standing waves and the odd-parity standing waves The odd-parity standing waves vanish at the quench site l = 0, and thus are not affected by the defect potential. Therefore, they are still eigenstates of H f . But the evenparity standing waves are coupled to each other by the perturbation Hamiltonian H 1 . In terms of the new basis, the total wave function (7) can be rewritten as The quantity of most interest is the local particle number on an arbitrary site l. It is the sum of the contributions of all the particles, Here in the second line we have used that fact that the odd-parity standing waves are eigenstates of H f and thus they only accumulate a global phase in time. The time variation of the local particle number then comes solely from the even-parity states. Equation (10) suggests studying the time evolution of each even-parity standing wave first. Figure 3 shows the time evolution of the probability density D + l (t) ≡ | l|k + (t) | 2 for three different values of k. In Fig. 3(a), we see that D + l alternates between plateaus on two different heights. This should be compared with Fig. 2(b), where the initial state is not an even-parity state but a Bloch state. We see that while D + l jumps back and forth between just two plateaus, D l has a sinusoidal profile and shows more plateaus. The difference comes from the fact that |k(t) contains also an odd-parity part |k − (t) , which interferes with the even-parity part |k + (t) . The remaining two panels of Fig. 3 demonstrate a point less emphasized in Refs. [16,17]. There we see that as the wave vector q = 2πk/N recedes from the inflection point q = π/2, the two-plateau structure deteriorates. This phenomenon has its root in the fact that the dispersion curve ε(q) is best approximated by a straight line at q = π/2 (or the middle of the energy band) and least at q = 0 or π (or the edges of the band). A detailed theory will be presented in Sec. III A below. Now apart from some constant contribution from the odd-parity standing waves, n(l, t) is the superposition of all the curves | l|k + (t) | 2 . In view of the three panels in Fig. 3, it is uncertain what will result from the superposition; in particular, it is uncertain whether n(l, t) will display any plateau structure. Superficially, there are two unfavorable factors. First, for different k, the switch times rT ±s c are different, as both T and s c depend on k; second, for those k ≃ 0, the plateau structure is not well developed at all. However, as Fig. 4(b) shows, n(l, t) does show a well-shaped two-plateau structure. Moreover, its switching between the two plateaus is synchronized with that of | l|k + f (t) | 2 , as is evident by comparing Fig. 4(a) with 4(b). Here it is interesting that in n(l, t), the features of the low-lying states are smeared out-only that of the states on the top like |k + f (t) is still visible.

III. ANALYTIC FORMULAS
We face the problem to account for the unexpected two-plateau structure displayed by n(l, t). In view of the irregular behavior of the low-lying states, to develop a comprehensive theory seems a task of challenge and is yet to be fulfilled. Here, we confine ourselves to a relatively humble goal. That is, we admit that n(l, t) alternates between two plateaus and that it does so in the same pace as | l|k + f (t) | 2 , and focus only on the problem of predicting the heights of the plateaus. Once this is done, the skeleton of the curve of n(l, t) is determined. As the first plateau corresponds to the pre-quench valuen of n(l), the real problem is to predict the height of the second plateau. This is what we try to do in the following subsections. As we shall see, this humble target is still too high, and in many cases we have to resort to bold approximations or just formal procedures to do what we can only do. However, interestingly, the formulas so obtained, although not really rigorously justified, turn out to be very accurate predictions of the height of the second plateau. In particular, in the end, we will get an analytic formula accurate for an arbitrary site l.
A. An exactly solvable ideal model As Eq. (10) shows, to calculate n(l, t), we can calculate D + l (t) = | l|k + (t) | 2 first. In this subsection, we try to find an approximate value of the latter by using the proximity of the realistic model to an ideal model, which is exactly solvable. This has essentially been done in Refs. [16,17] already, so below we just sketch the construction of the ideal model and present the most relevant result, with which the phenomena in Figs. 2, 3, and 4(a) can be understood. For the sake of completeness, the detailed derivation is retained but deferred to the Appendix. The derivation there is more straightforward than that in Refs. [16,17]. We shall also discuss more about the limitation or validity of the approximation than in Refs. [16,17].
The state |k + i (t) evolves by the equation with the initial value |k + i (t = 0) = |k + i . As both H 0 and H 1 are diagonal with respect to the even-odd partition of the basis functions, |k + i (t) always belongs to the even-parity subspace, and is of the form where a k (t) are time-dependent coefficients. With respect to the even-parity basis functions {|k + }, H 0 is diagonal while H 1 is close to a matrix whose entries are all equal. Specifically, we have and (g ≡ U/N ) Numerically, it is observed that for given U , if the lattice size N is large enough, and if the wave vector q i is not so close to 0 or π, in the time evolution, only those few states |k + with k ≃ k i are significantly populated. Other states are far off resonance and hence are suppressed. This fact means that in the full spectrum of H 0 , only a small segment is effective in the dynamics in question. Within this segment, the spectrum of H 0 is nearly equally spaced, with the gap between two adjacent levels being approximately where η = 2π/N is the difference between two adjacent wave vectors in the Brillouin zone. Also, by (14), the coupling between two arbitrary levels in the segment is of the constant value 2g.

Construction and solution of the ideal model
We now extend this segment, together with the two features, to the whole axis (−∞, +∞). We consider an ideal model consisting of infinitely many levels {|j | − ∞ ≤ j ≤ +∞}, whose energies are equally spaced, i.e., and between any pair of which, the coupling is constant, i.e., for ∀j 1 , j 2 ∈ Z. HereH 0,1 are counter-parts of H 0,1 . By construction, |j is meant to be the counter-part of |(k i + j) + , and its wave function in the real space is postulated to be The correspondence between |j and |(k i + j) + of course fails for j large enough, either because |(k i + j) + no longer belongs to the linear segment, or simply because there is no one-to-one correspondence between a finite set and an infinite one. However, as long as those fictitious levels |j are not significantly populated, the error introduced could be negligible. Also note that in (16), the energy of |j is defined as j∆ but not ε(q i ) + j∆. This is a mere shift of the origin of energy, and does not lead to any observable effect.
Corresponding to the original dynamical problem of the realistic model, in the ideal model, initially the system is in the state |ψ(t = 0) = |0 . It then evolves as At any time t, it is of the form |ψ(t) = j∈Zã j (t)|j . As |j is the counterpart of |(k i + j) + , so isã j (t) that of a ki+j (t)e iε(qi)t . Onceã j (t) is calculated, we get an approximation of l|k + i (t) as Now we have two steps to go. First, we have to solvẽ a j (t) and then we can substitute it into (20) to calculate l|k + i (t) and in turn D + l (t) approximately. The detailed calculation is in the Appendix. Here we just describe the ultimate result. For the current problem, a characteristic time scale is the so-called Heisenberg time T = 2π/∆. For a given site l satisfying −L ≤ l ≤ L, there are two critical times, s + c = T |l|/N and s − c = T (N − |l|)/N . The physical meaning of these times will be clear in the next paragraph. Now for an arbitrary time t = rT + s, where 0 < s < T and r is a nonnegative integer, we have where the angle variable θ i = θ(q i ) is defined as Equation (21) explains the two-plateau structure in Figs. 3 and 4(a). Its content can be best illustrated by presenting the snapshots of the probability distribution | l|k + i (t) | 2 at consecutive times, as is done in Fig. 5. There, 17 snapshots of | l|k + i (t) | 2 , which are sampled equal-distantly in the time interval [0, 2T ], are stacked in a pile. The red dashed lines indicate the fronts of the waves generated by the sudden quench. They move at the constant velocity v i = 2 sin q i . The time for them to finish a loop on the lattice ring is then N/v i , which is exactly T . Similarly, rT + s ± c are just the times when the wave fronts (backwards or forwards) pass the site l. One can recognize that the curves show a pattern inside the two rhombuses and another pattern outside the two rhombuses. On the interfaces between the two regions, the dislocation is quite clear. This observation embodies the meaning of (21) that | l|k + i (t) | 2 alternates between two values in time.
Another thing noticeable in Fig. 5 is that at t = rT , the probability distribution | l|k + i (t) | 2 reconstructs itself to its initial pattern approximately. It would be an exact reconstruction if it were the ideal model-the ideal model has an equal-distant spectrum of (A.7), hence the revival would be periodic and complete.
Some new light will be shed on (21) when we consider the eigenstates of H 0 + H 1 in Sec. IV.

Validity of the approximation
Now we discuss the validity of the approximation. The ideal model agrees with the realistic model only locally. Hence, for the dynamics occurring in the ideal model to be a faithful reflection of that occurring in the realistic model, the levels actively participating (significantly populated) in the ideal model should be as few as possible, and the linearization of the spectrum ε(q) at q = q i should be as good as possible. By (A.13), we see that the population on |j (j = 0) is up bounded by the value where F (x) ≡ 16x 2 /(1 + 4π 2 x 2 ). As the series j≥1 1/j 2 converges absolutely, for any ǫ > 0, there exists a n max , such that |j|>nmax |ã j (t)| 2 < ǫ. Moreover, as F (x) is a monotonically increasing function of |x|, the smaller the ratio |g/∆| is, the smaller n max can be chosen. As g/∆ = U/(4π sin q i ), for fixed U and N , the minimum of |g/∆| is achieved at q i = π/2. On the other hand, the energy of the state |(k i + j) + has the Taylor expansion At q i = π/2, the quadratic term vanishes identically, and the linearization is the best. Hence, the ideal model as an approximation works the best at q i = π/2. As q i recedes from this optimal point (either towards 0 or π), the ratio |g/∆| increases and at the same time the quadratic term in (24) increases too. Both facts are detrimental to the approximation, and the analytic predictions above become less accurate. This is demonstrated in Fig. 3.
Although π/2 is the most favorable value of q i , at other levels of q i , the approximation can work well too as long as the lattice size N is large enough. Actually, as g/∆ = U/(4π sin q i ) is independent of N , n max is independent of N . However, for j in the range −n max ≤ j ≤ n max , the quadratic term will become negligible in comparison with the linear term as N → ∞. Therefore, for given U , the range of q i where the ideal model is a good approximation will become larger and larger as N increases, and will cover (0, π) eventually in the limit of N = ∞.

B. Magic of the Abel summation method
Having treated the single-particle problem, we now proceed to the many-particle problem. By (10), we face the problem of summing over the |k + (t) states. But the problem is that, as pointed out in the proceeding subsection, for a finite lattice and for states |k + lying close to the bottom of the energy band, we do not have an accurate estimation of | l|k + (t) | 2 .
This seems an insurmountable difficulty. But there is an important fact worth noting. In Fig. 6, we study the population on each even-parity state |k at some arbitrary times. We see that the population distribution is close to its initial value Hence, although on the single-particle level, each |k + state mixes with its neighbors, on the many-particle level, population transfer occurs only between the states in the vicinity of the Fermi energy. Those states at the bottom of the Fermi sea, which pose a difficulty for the summation, are dormant. Actually, that the states deep beneath the surface of the Fermi sea are in a sense dormant is also reflected in Fig. 4. Originally, n(l, t) is the superposition of all | l|k + (t) | 2 with 0 ≤ k ≤ k f . By (21) and as demonstrated in Fig. 3, these quantities switch between plateaus at different times. However, the net result of the superposition is that n(l, t) switches between the two plateaus with the same rhythm as | l|k + f (t) | 2 . The observations above motivate us to approximate the realistic Fermi sea by a fictitious one. Suppose that in the fictitious model in the proceeding subsection, all the levels |j with j ≤ 0 are filled initially. After the quench, the levels mix with each other, but overall, population transfer occurs only at the surface layer of the Fermi sea. For this fictitious model, the analytic formula (21) works exactly for all the single-particle levels, and the plateau switchings are all synchronized. The fictitious Fermi sea is infinitely deep, and hence the total number of particles on a local site l is ill-defined, but its variation could be well defined. By (21), the height of the second plateau relative to the first plateau is Here θ f ≡ θ(q f ). Apparently, this series is not convergent but oscillatory. But one can use the trick of the Abel summation method to extract a finite and definite value from it [25]. The idea is that while the series diverges, the series might converge for 0 ≤ x < 1, and the limit might exist. If so, then the original summation (28) is assigned the value S. Carrying out these procedures for the series (27), which involve just some geometric series, we get readily Here the subscript indicates that the formula is from the Abel summation method. It turns out that this formula, which is obtained by a purely formal procedure, is a very accurate prediction of the second plateau in the time evolution of n(l, t).
In Fig. 7, we show the evolution trajectories (solid lines) of the variation of local particle number δn(l, t) = n(l, t) −n for two different sets of parameters. The analytic prediction of (31) is indicated by the horizontal dashed lines. In each panel, we see that in the time intervals where the second plateau is expected, the solid line either almost coincides with the horizontal dashed line, or oscillates around it.
In Fig. 8, we check the accuracy of the formula (31) in a more systematic way. For the two sets of parameters (N, k f , U ) in Fig. 7, we compare the plateaus exhibited  Fig. 7(a) [Fig. 7(b)]. As δnA(l) decays like 1/l, for clarity, we present δn × l and δnA × l.
by δn(l, t) and the analytic prediction δn A (l) in (31) for 1 ≤ l ≤ 20. The height of the plateau is defined as where t 1,2 are chosen in the region where δn has settled down. In our calculation, we take t 1 = 0.5s + c + 0.5s − c and t 2 = 0.1s + c + 0.9s − c . This choice of t 1,2 is of course arbitrary, but other choices yield almost the same result. As δn A (l) decays like 1/|l|, for clarity, instead of δn or δn A , we present δn × l and δn A × l in Fig. 8. We see that the two agree with each other very well in the whole range of l. The difference is only discernible at the lower limit l = 1. But even there, the relative error is smaller than 20%.
The quantity δn A depends on the lattice size N . But in the thermodynamic limit of (N, k f ) → ∞ with q f = 2πk f /N fixed, it converges to which is a cosine function modulated by the 1/|l| function. This limiting behavior is illustrated in Fig. 9. In Figs. 9(a) and 9(b), we study two cases with the Fermi wave vectors q f being the same while the lattice sizes N differing by a factor of three. In each panel, the red dash dotted line indicates the value δn ∞ A , which is independent of N . From Fig. 9(a) to 9(b), we see clearly that as the lattice size increases, the plateau, which is accurately predicted by δn A , tends towards the dash-dotted line predicted by δn ∞ A . This means that the system shows some finite-size effect and that is accurately captured by δn A .  (31), and the dash-dotted line its thermodynamic limit (33). A more systematic study is presented in (c). Fixing (U, l, N/k f ) and starting with (N, k f ) = (401, 100), (N, k f ) are enlarged consecutively by 1, 3, 5, . . . , 15 times. The circles denote the predictions of δnA(l). Apparently, they fall on a straight line with respect to 1/N and, as 1/N → 0, converge towards δn ∞ A (l), which is indicated by the horizontal line.
The scaling behavior of the finite-size effect is investigated more systematically in Fig. 9(c), where the lattice size N in Fig. 9(a) is further increased by 5, 7, . . . , 15 times, and δn A is plotted against 1/N . It is apparent that δn A (the circles) approaches δn ∞ A (the horizontal line) by the 1/N power law.
We have thus seen that the formal result (31) is a very accurate prediction of the height of the second plateau. Its success is the most dramatic thing in this paper. To make sense of the magic Abel summation method, we note that it is essentially a way of regularization. The exponentially decaying factor x m (when 0 < x < 1) helps to suppress the contribution of the low-lying states (with m ≫ 1), which are fictitious in the fictitious Fermi sea and dormant in the realistic Fermi sea.
Of course, formula (31) also has its drawbacks. In the first place, it is not defined at the quench site l = 0. This corresponds to the fact that the series m≥0 1 is not Abel-summable. This shortcoming will be overcome below by other formulas.

C. A formal summation
The success of the formal procedure in the proceeding subsection is very encouraging. We have thus considered applying (21) to the realistic Fermi sea, regardless of the validity of this formula for small k. Substituting (21) into (10), we get another prediction of the relative height of the second plateau to the first one as where the last term corresponds to the state |0 + . We have used the fact that θ(0) = ±π, and that the k = 0 term has to be scaled down by 1/2 in view of (8). Here the subscript stands for summation. Compared with δn A in (31), δn s has the advantage that it is also defined at l = 0.
The summation above can be rewritten more compactly as δn s = (S 1 − S 2 )/(2N ), with The second term was easily calculated as the arguments of the cosine function forms an arithmetic sequence. The first term is more complicated-Because of the θ function, the arguments of the cosine function are not an arithmetic sequence. But they are not far away from an arithmetic sequence, especially when |l| is large. Actually, the θ function defined in (22) is a very regular function in the sense that its derivative is bounded. It can be easily verified that |dθ/dq| ≤ 4/|U |. Therefore, if |U | or |l| is large enough, the gap between two neighboring arguments is approximately 2η|l|.
We thus can employ the same idea as for S 2 , and approximate S 1 as It is easy to verify that the error introduced in the first and the last line is on the order of 1/N , and the approximation in the third line is legitimate if |U l| ≫ 1.
Collecting (35) and (36), we see that the formal summation (34) can reduce to (31) in the large l limit. This is verified in Fig. 10. This is an interesting result, as (31) and (34) come from the fictitious Fermi sea and the realistic Fermi sea, respectively. That they agree with each other asymptotically adds to the plausibility of both. Now, for l large enough, δn s is close to δn A , while the  Fig. 7(a) [Fig. 7(b)]. Except for l = 0, we present δn × l and δns × l instead of δn and δns.
latter in turn is close to δn. Hence, for l large enough, δn s should also be a good approximation or prediction of δn. In Fig. 11, we compare δn s with δn for the two sets of parameters in Fig. 7 or Fig. 8. Indeed, as l increases, the two curves converge together. The two differ significantly only for l small.

D. Influence of the defect mode
So far, we have obtained two analytic formulas for the height of the second plateau. They both work very well for l large, but become inadequate for l small. The Abel formula (31) is not defined at all for l = 0, while the summation formula (34) deviates from the exact value significantly for l around zero.
This local failure forces us to examine the small-l cases more carefully. It turns out that the small l's need extra attention already on the single-particle level. Actually, an attentive reader should have noticed that in Fig. 5, a spike pulses at l = 0, while the pattern of | l|k + i (t) | 2 in its vicinity is mostly stationary. In Figs. 12(a) and 12(b), the local probability density | l|k + (t) | 2 as a function of time is shown for l = 0 and l = 2, respectively, with all other parameters the same as in Fig. 3(a). We see that unlike in Fig. 3(a), where l = 50, for these two small values of l, | l|k + (t) | 2 does not show any sign of plateau, but just oscillates very quickly. More systematic investigation reveals that as l increases, the oscillation amplitude shrinks and the two-plateau structure emerges gradually. The former fact is already visible by comparing Fig. 12(a) with 12(b).
All these observations point to the defect mode induced by the defect potential. It is well-known that on an infinite one-dimensional tight binding lattice, a defect po-tential like H 1 will induce a defect mode |φ d localized around the defect [26]. Its wave function is of the form The corresponding eigenenergy is ǫ d = (sgn U ) √ U 2 + 4, and the normalization factor A = U/ǫ d . In reality, we are dealing with a finite lattice. But as its wave function decays exponentially, as long as the lattice size is much larger than its characteristic size, the defect mode will always be there with its energy and wave function perturbed only slightly.
Besides the defect mode, the post-quench Hamiltonian H 0 +H 1 has also many extended states, whose wave functions spread out on the whole lattice. Let us denote them as |φ m and their energies as ǫ m (1 ≤ m ≤ N − 1). Formally, we have where a d,m (k) ≡ φ d,m |k + , φ d,m (l) ≡ l|φ d,m , and in the second line, as the subscripts indicate, C d and C e denote the contributions of the bound state and the extended states, respectively. Now for l large enough, the term C d is exponentially small and dominated by C e , and | l|k + (t) | 2 ≃ |C e (l, t)| 2 . However, for l small, the term C d will be competitive with C e . Although the overlap a d between the localized mode |φ d and the extended state |k + is on the order of 1/ √ N , the amplitude of the mode function at the site l, i.e., φ d (l), is on the order of unity. Hence, the term C d , like C e , is also on the order of 1/ √ N . In Fig. 3(a), it is verified that for appropriate parameters, and when l is large, | l|k + (t) | 2 shows the predicted plateau structure. As in this region, l|k + (t) ≃ C e (l, t), the observed behavior is actually that of |C e (l, t)| 2 . That is, in the time interval (s + c , s − c ), By analytic continuation, this relation should hold also for l small. This is confirmed in Figs. 12(c) and 12(d).
There, C e is calculated by subtracting C d from l|k + (t) . We see that the well-developed plateau agrees with the prediction of (39) exactly. Making use of (39), and the fact that C d and C e contain no common Fourier component, we have on average in the time interval (s + c , s − c ), Here the first term on the right hand side is the contribution of the defect mode, which is missed by the ideal model. It can be and was ignored for large l. Retrieving this term, we now have a modified prediction of the plateau of δn(l, t) as where δn s is defined in (34) and P d is the total population on the defect mode, i.e., Note that P d is l independent.
With the contribution of the defect mode taken into account, finally (the subscript of δn f stands for final) we have an accurate formula for the height of the second plateau for all l. In Fig. 11, we see that δn S agrees with δn very well for large l, but significantly underestimates it for small l. It turns out that the contribution of the defect mode, which is positive-definite, and is peaked at l = 0, fills the gap exactly. In Fig. 13, the time evolution of δn(l, t) for two small values of l is shown. We see that the plateaus ensuing the quench are located exactly at the height predicted by δn f (l) in (41). For large l, δn f reduces to δn s , and is again accurate.

IV. THE INFINITE LATTICE CASE
So far, we have always assumed a finite lattice. On a finite lattice, as we have seen in the figures above, the local particle number n(l) never really settles down to a constant value-again and again, it bursts up suddenly after remaining quiet for a long time. Actually, the out-going waves generated by the sudden quench will sweep across the lattice ring repeatedly, and each time a wave front passes by a site l, the local particle number δn(l) experiences an abrupt and violent change. An overall picture can be well provided by a series of snapshots of δn(l, t), as we do in Fig. 14. There, like in Fig. 5, the trajectories of the wave fronts divide the space-time into two different regions, where δn show qualitatively different patterns.
But on an infinite lattice, the out-going waves will not come back again. The out-going waves will just pass by each site once and after that the local particle number n(l) will converge to a constant value. That is, the outgoing wave fronts leave behind a stationary pattern of n(l), the region of which expands linearly with time. This phenomenon is actually visible in the lower few curves in Fig. 14. If superimposed together, one can see that they agree very well in the inner region.
This stationary pattern of n(l) can be fairly called dynamically generated Friedel oscillation. As discussed in Sec. III A 2, as N → ∞, the ideal model is a good approximation for all q ∈ (0, π). Hence, in this limit, formula (34) and in turn (41) can be justified almost rigorously. By (41), the stationary value of δn(l, t → ∞) is where the two terms on the right hand side corresponds to the contribution of the defect mode and the extended states, respectively. By straightforward calculation and replacing summation by integral, we have It is natural to compare the dynamically generated Friedel oscillation with the conventional static (or equilibrium) Friedel oscillation [27][28][29][30][31]. It is well known that a defect in a Fermi sea will induce density oscillation around it. To calculate the density oscillation, a convenient approach is using Green functions [32]. But here we shall take a more straightforward way, which will also help us gain some insight of the problem.
Imagine that for the initial Hamiltonian H 0 , all the Bloch states with wave vectors |q| ≤ q f are filled. Now ramp up the defect potential adiabatically. The oddparity eigenstates are not affected and their contribution to the particle density is unchanged. But the even-parity eigenstates do feel the potential and deform accordingly. For H 0 , a generic even-parity eigenstate is of the form where the wave vector q satisfies qN = 2πm for some integer m. For the final Hamiltonian H f = H 0 + H 1 , its even-parity eigenstates fall into two categories. The first category contains only one state, i.e., the defect mode, which in the N → ∞ limit, reduces to |φ d in (37). The second category are the extended states. By Bethe's ansatz, their wave functions are in the form is some phase shift, and B is a normalization factor, which in the thermodynamic limit should be close to 2/N . The eigenvalue equation at l = 0 determines the relation between θ and q, which is exactly Eq. (22). The only thing is that the q's are to take different values than in (46). Specifically, the periodic boundary condition requires ψ (f ) (−L) = ψ (f ) (L + 1), which leads to the condition with m being an integer. By continuity, ψ  (47), it suggests that when we quench an even-parity standing wave like (46) on a sufficiently large lattice ring, the stationary pattern predicted by (21) inside the light cone actually corresponds to the eigenstate (47) of the final Hamiltonian. That is, after the abrupt quench, apart from a component proportional to the defect mode, the initial state is phase shifted to its adiabatic correspondence inside the light cone. The resulting probability distribution is largely the same as if the defect potential is turned on adiabatically. This is confirmed in Fig. 15. Now we can sum over the even-parity states and determine the static Friedel oscillation. Again, it contains two parts, i.e., where the defect mode's contribution is dyn results naturally from the singleparticle behavior discussed above. Hence, we see that the dynamically generated Friedel oscillation differs from its static counterpart only by an exponentially localized part.
It is instructive to study the asymptotic behavior of δn st (l), or that of δn dyn (l) for large |l|. As both δ   dyn (l). The second part of (45) can be easily integrated to sin(2q f |l|)/4π|l|. As for the first part, it is 1 2π with f (q) = cos θ = 4 sin 2 q − U 2 4 sin 2 q + U 2 , g(q) = sin θ = 4U sin q 4 sin 2 q + U 2 . Both f and g are slowly varying functions of q. Integrating by parts, we have 1 2π Combing the two parts, we recover the formula (33).

V. CONCLUSIONS AND DISCUSSIONS
We have investigated the nonequilibrium dynamics of a Fermi sea when a defect potential is introduced suddenly. On the one hand, this is a natural generalization of our previous study of the quench dynamics of a single Bloch state, which is featured by the cusps or plateaus in the time evolution of some physical quantities. The concern is whether these singular behaviors will survive in the many-particle case. On the other hand, from the pointof-view of the well-known (equilibrium) phenomenon of Friedel oscillation, this is a problem of interest: what if the defect potential is introduced not adiabatically but abruptly?
As for the first problem, the primary observation is that the local particle number n(l) alternates between two plateaus. This two-plateau behavior is actually exhibited on the single-particle level by an even-parity eigenstate. But its existence on the many-particle level is far from obvious. Let alone other issues, the plateauswitching times for different single-particle states disperse continuously in a wide range, and hence it could be totally reasonable if the abrupt plateau-switching were smeared out.
Leaving aside the problem of the very existence of the plateaus, we have focused on predicting the heights of the plateaus. This seems like a mere summation job. But the problem is that we do not really know all the summands, as the theory we had developed for the corresponding single-particle problem holds not for all the single-particle states involved. We thus had to resort to arguments and approximations, and when confronted with a series oscillatory instead of convergent, we simply took on the Abel summation method to assign a finite value to it. Magically, this finite value is a very accurate prediction of the numerically exact one for l not very small. We are thus in a (possibly fortunate) dilemma that we get the numbers accurate without knowing why. The unreasonable relevance and accuracy of the formalism is yet to be understood. Anyway, the success of the formal approach encouraged us to try other formal expressions, and ultimately an expression accurate for all sites is obtained. We have thus successfully determined the backbone of the evolution trajectory of n(l). This is possibly all that we can wish for, as the erratic details exhibited by n(l) are definitely beyond any simple analytic formula.
As for the second problem, the interesting thing is on the single-particle level. It is found that on an infinite lattice ring, after the quench, the state adjusts itself towards the eigenstate of the final Hamiltonian with equal energy. Specifically, apart from a localized component corresponding to the defect mode, the time evolving state agrees with the phase-shifted eigenstate of the final Hamiltonian inside the light cone. Because of this nice single-particle behavior, the dynamically generated Friedel oscillation differs from the conventional, static Friedel oscillation only by an exponentially decaying term. It is conjectured that this single-particle property holds regardless of the quench potential or the quench protocol, as long as the quench potential is a local one and it stays unchanged in the end.
As stressed as our point, simple models can yield rich dynamics. Here, in the study of a non-interacting manyparticle model, we have actually noticed many interesting phenomena which are not accounted for yet. For example, in Fig. 3(c), around t = 50, the solid curve oscillates in a very elegant way, which is typical of a state close to the band edge. It reminds us of the Airy function actually. Can we find an analytic formula for it? It might be related to the shape of the out-going wave fronts. Similarly, in both panels of Fig. 4, in the first period and on the second plateau, the oscillation amplitude shrinks with time. But with what rate does it shrink? Is it power law or exponential? To answer these questions, likely we have to go beyond the ideal model. Anyway, more problems are posed than answered, and more effort is to be expended on this simple model. The two terms containing ∆ and g come fromH 0 and H 1 , respectively. The important thing is that the latter is independent of j. Hence, we can introduce the auxiliary quantity Where T ≡ 2π/∆ is the Heisenberg time associated with the linear spectrum in (16). This name apparently comes from the time-energy uncertainty relation. As we shall see below, it is the most important time scale in the model, as it determines the period of the curves in Fig. 4. We want to decompose the meromorphic function L(p) according to its poles and the corresponding residues in the form L(p) = m∈Z b m /(p + iE m ), so that we can recover S(t) from L(p) as S(t) = m∈Z b m e −iEmt . The poles −iE m of L(p) are roots of the equation 1 + gT cot(−ipT /2) = 0. We solve easily with γ = θ(q i )/T , where the angle θ is defined as a function of the wave vector q as θ ≡ 2 arctan gT = 2 arctan U 2 sin q . (A.8) The corresponding residues are b m = 1/p gT d dp cot(−ipT /2) p=−iEm We thus know that S(t) is of the form S(t) = 2g 1 + g 2 T 2 m∈Z e −i(m∆+γ)t m∆ + γ . (A.10) The summation can be carried out by noting that The result is that, for t = rT + s, where 0 < s < T , and r is a nonnegative integer, which is a piece-wise constant function. Substituting this result into (A.3), we get a j (t) = δ j,0 + 2g(e −ij∆s − 1) ∆(1 + igT )j e −irθ . (A.13) Note that the factor in the brackets depends only on j and s, while the exponential factor outside depends only on r. In the special case of j = 0,ã j is a piece-wise linear function of t. This nonsmoothness results in the cusps in Fig. 2(a), as P i,r = |1 ±ã 0 | 2 /4.
Substituting (A.13) into (20) and ignoring the global phase e −irθ−iε(qi)t which is independent of l, we get In the derivation above, we see that it is the two defining features of the ideal model, i.e., (16) and (17), that are responsible for its exact solvability. We argue that these two nice properties are shared by many other onedimensional models [22]. For a generic one-dimensional system, in the midst of its spectrum, both the eigenenergies and the eigenfunctions should be smooth functions of the state index. Hence, at least locally, the two properties should hold and the ideal model could be relevant for the dynamics of the realistic model.
The problem is that its relevance is confined to one dimension. In higher dimensions, the eigenfunctions are indexed with multiple variables (say k x and k y ) and thus generally the eigenenergies are not equally spaced, i.e., (16) does not hold. Therefore, for higher-dimensional generalizations of the local quench scenario considered here, we need new analytic approaches. While this sounds like a challenge, it also means that the dynamic Friedel oscillation in higher dimensions could be fertile of new features.