Total Energy Dynamics and Asymptotics of the Momentum Distribution Following an Interaction Quench in a Two-component Fermi Gas

TL;DR

This paper investigates the energy dynamics and momentum distribution asymptotics in a two-component Fermi gas after an interaction quench, addressing divergence issues and analyzing the evolution of Tan's contact across different quench protocols.

Contribution

It introduces a method to remove divergences in energy calculations and derives the time evolution of Tan's contact during interaction ramps.

Findings

Divergences in energy dynamics can be regularized for arbitrary interaction protocols.

Tan's contact exhibits rapid oscillations and reaches a stationary value after a sudden quench.

In the adiabatic limit, the contact grows quadratically and saturates to equilibrium.

Abstract

The absence of a characteristic momentum scale in the pseudo-potential description of atomic interaction in ultracold (two-component Fermi) gases is known to lead to divergence in perturbation theory. Here we show that they also plague the calculation of the dynamics of the total energy following a quantum quench. A procedure to remove the divergence is devised, which provides finite answers for the time-evolution of the total energy after a quench in which the interaction strength is ramped up according to an arbitrary protocol. An important result of this analysis is the time evolution of the asymptotic tail of the momentum distribution (related to Tan's contact) and the contact for a linear interaction ramp are obtained, as a function of the interaction ramp time in the crossover from the sudden quench to the adiabatic limit are reported. In sudden quench limit, the contact, following a rapid oscillation, reaches a stationary value which is different from the equilibrium one. In the adiabatic limit, the contact grows quadratically in time and later saturates to its equilibrium value for the final value of the scattering length.

Equations

$$\hat{V}=\frac{g}{2\Omega }\sum _{\mathbf{pkqr}}\sum _{\sigma \ne \alpha }{c}_{\mathbf{p}\sigma }^{†}{c}_{\mathbf{k}\alpha }^{†}{c}_{\mathbf{q}\alpha }{c}_{\mathbf{r}\sigma }{\delta }_{\mathbf{p}+\mathbf{k},\mathbf{q}+\mathbf{r}},$$

$$K\cot {\delta }_s(K)=-\frac{1}{a_s}+\frac{K^{2}R}{2}+\cdots$$

$$T_{\sigma \alpha }(\mathbf{k},\mathbf{p})=g+\frac{2g^2}{\Omega }\sum _{\mathbf{qr}}\frac{{\delta }_{\mathbf{p}+\mathbf{k},\mathbf{q}+\mathbf{r}}}{E_{pkqr}}+O(g^3)=\frac{4\pi a_s}{m}.$$

$$H(t)=H_0+S(t)V,$$

$$S(t)=\theta (t)\left[\theta (T_r-t)\frac{t}{T_r}+\theta (t-T_r)\right],$$

$$n_{p\sigma }(t)=\mathrm{Tr}\left[\rho (t)c_{\mathbf{p}\sigma }^{†}{c}_{\mathbf{p}\sigma }\right],$$

$$E_{\text{tot}}(t)$$

$$=\mathrm{Tr}\left[\rho (t)H_0\right]+S(t)\mathrm{Tr}\left[\rho (t)V\right],$$

$$E_{\mathrm{kin}}(t)=\mathrm{Tr}\left[\rho (t)H_0\right]=\sum _{\mathbf{p}\sigma }{ϵ}_p{n}_{p\sigma }(t).$$

$$⟨O(t)⟩$$

$$O(t)$$

$$+\frac{(-i{)}^2}{2!}{\int }_{C}d{t}_{1}d{t}_2⟨\mathcal{T}\left[\tilde{V}(t_1)\tilde{V}(t_2)O(t)\right]{⟩}_c$$

$$+\cdots$$

$$i{G}_{p\sigma }(t_1,t_2)=⟨\mathcal{T}\left[c_{\mathbf{p}\sigma }(t_1)c_{\mathbf{p}\sigma }^{†}(t_2)\right]⟩,$$

$${\mathcal{G}}_{p\sigma }(a,b)=\left(\begin{array}{cc}i{G}_{p\sigma }^T(a,b)& i{G}_{p\sigma }^{<}(a,\bar{b})\\ i{G}_{p\sigma }^{>}(\bar{a},b)& i{G}_{p\sigma }^{\tilde{T}}(\bar{a},\bar{b})\end{array}\right).$$

$$i{G}_{p\sigma }^{<}(t_1,{\bar{t}}_2)$$

$$i{G}_{p\sigma }^{>}({\bar{t}}_1,t_2)$$

$$i{G}_{p\sigma }^T(t_1,t_2)$$

$$+\theta (t_2-t_1)i{G}_{p\sigma }^{<}(t_1,{\bar{t}}_2),$$

$$i{G}_{p\sigma }^{\tilde{T}}({\bar{t}}_1,{\bar{t}}_2)$$

$$+\theta (t_1-t_2)i{G}_{p\sigma }^{<}(t_1,{\bar{t}}_2).$$

$$\begin{gathered} n_{p\sigma }(t)=⟨{\hat{n}}_{p\sigma }⟩-i{\int }_{C}d{t}_1⟨\mathcal{T}\left[\tilde{V}(t_1){\hat{n}}_{p\sigma }(t)\right]{⟩}_c \\ +\frac{(-i{)}^2}{2!}{\int }_{C}d{t}_{1}d{t}_2⟨\mathcal{T}\left[\tilde{V}(t_1)\tilde{V}(t_2){\hat{n}}_{p\sigma }(t)\right]{⟩}_c+\cdots \end{gathered}$$

$$\begin{gathered} n_{p\sigma }(t)=⟨{\hat{n}}_{p\sigma }⟩+{\int }_{C}d{t}_1{\mathcal{G}}_{p\sigma }(t,t_1){\Sigma }_{p\sigma }^{(1)}(t_1){\mathcal{G}}_{p\sigma }(t_1,t) \\ +{\int }_{C}d{t}_1{\int }_{C}d{t}_2{\mathcal{G}}_{p\sigma }(t,t_2){\Sigma }_{p\sigma }^{(2)}(t_2,t_1){\mathcal{G}}_{p\sigma }(t_1,t)+\dots \end{gathered}$$

$${\mathcal{G}}_{p\sigma }(t,b)=e^{-i{ϵ}_p(t-b)}\left(\begin{array}{cc}1-n_{p\sigma }^0& -n_{p\sigma }^0\\ 0& 0\end{array}\right),$$

$${\mathcal{G}}_{p\sigma }(a,t)=e^{-i{ϵ}_p(a-t)}\left(\begin{array}{cc}-n_{p\sigma }^0& 0\\ 1-n_{p\sigma }^0& 0\end{array}\right).$$

$${\Sigma }_{\sigma }^{(1)}=\frac{g}{2\Omega }\sum _{\mathbf{k}}{n}_{k,-\sigma }^0.$$

$${\Sigma }_{p\sigma }^{(2)}(b,a)=-\frac{2g^2}{{\Omega }^2}S(t_1)S(t_2)\sum _{\mathbf{kqr}}{\delta }_{\mathbf{p}+\mathbf{k},\mathbf{q}+\mathbf{r}}$$

$$\times \left(\begin{array}{cc}{\bar{\Sigma }}_{\sigma }^{(2,T)}(b,a)& {\bar{\Sigma }}_{\sigma }^{(2,>)}(\bar{b},a)\\ {\bar{\Sigma }}_{\sigma }^{(2,<)}(b,\bar{a})& {\bar{\Sigma }}_{\sigma }^{(2,\tilde{T})}(\bar{b},\bar{a})\end{array}\right).$$

$${\bar{\Sigma }}_{\sigma }^{(2,<)}(t_2,{\bar{t}}_1)$$

$$=-(1-n_{q\sigma }^0)(1-n_{r\alpha }^0)n_{k\alpha }^0{e}^{i(t_1-t_2)({ϵ}_q+{ϵ}_r-{ϵ}_k)},$$

$${\bar{\Sigma }}^{(2,>)}({\bar{t}}_2,t_1)$$

$$=n_{q\sigma }^0{n}_{r\alpha }^0(1-n_{k\alpha }^0)e^{i(t_1-t_2)({ϵ}_q+{ϵ}_r-{ϵ}_k)},$$

$${\bar{\Sigma }}^{(2,T)}(t_2,t_1)$$

$$=\theta (t_2$$

$${\bar{\Sigma }}^{(2,\tilde{T})}({\bar{t}}_2,{\bar{t}}_1)$$

$$=\theta (t_1$$

$$n_{p\sigma }^{(2)}(g,t)=-\frac{2g^2}{{\Omega }^2}\sum _{\mathbf{kqr}}{A}_{pkqr}^{\sigma }{F}^{(2)}(E_{pkqr},t){\delta }_{\mathbf{p}+\mathbf{k},\mathbf{q}+\mathbf{r}},$$

$$\begin{gathered} A_{pkqr}^{\sigma }=\sum _{\alpha \ne \sigma }[n_{p\sigma }^0{n}_{k\alpha }^0(1-n_{q\sigma }^0)(1-n_{r\alpha }^0) \\ -(1-n_{p\sigma }^0)(1-n_{k\alpha }^0)n_{q\sigma }^0{n}_{r\alpha }^0]. \end{gathered}$$

Full text

I Introduction

The single-channel model Lee et al. (1957) provides a compact, single-parameter description of interactions in ultracold gases. Within this model, interactions in a two-component Fermi gas are described by the following term in the Hamiltonian:

[TABLE]

where $c_{\boldsymbol{p}\sigma}$ ($c^{{\dagger}}_{\boldsymbol{p}\sigma}$) are fermion destruction (creation) operators obeying $\{c_{\boldsymbol{p}\sigma},c^{{\dagger}}_{\boldsymbol{k}\alpha}\}=\delta_{\boldsymbol{pk}}\delta_{\sigma\alpha}$ ($\sigma,\alpha=\uparrow,\downarrow$) and $\{c_{\boldsymbol{p}\sigma},c_{\boldsymbol{k}\alpha}\}=0$. However, in terms of the parameter $g$, it appears as if Eq. (1) holds for arbitrarily large particle momentum exchange, $\boldsymbol{K}=\boldsymbol{p}-\boldsymbol{r}$. In other words, the interaction in (1) lacks of a characteristic momentum scale, which nevertheless exists for real interactions but depends on the microscopic details of the two-atom potential.

The lack of a characteristic momentum scale has important consequences for the calculation of physical properties using the single-channel model. For instance, the perturbation series for the ground state energy is plagued with divergences arising from the behavior of integrals at high momenta (see e.g.  Abrikosov et al. (1975a); Pathria and Beale (1996a)). In addition, for momenta smaller than the inverse effective range (i.e. $p\ll R^{-1}$ 111This condition obeyed by most ultracold gases even in the vicinity of broad Feshbach resonances where $|a_{s}|$ diverges), the momentum distribution, $n_{p\sigma}=\langle c^{{\dagger}}_{\boldsymbol{p}\sigma}c_{\boldsymbol{p}\sigma}\rangle$ exhibits a $p^{-4}$ tail at large momenta $p\gg k_{F}$ ($k_{F}$ is the Fermi momentum). This renders divergent the kinetic energy, $E_{\text{kin}}=\sum_{p\sigma}n_{p\sigma}\epsilon_{p}$ Viverit et al. (2004); Tan (2008) ($\epsilon_{p}=p^{2}/2m$ is the single-particle dispersion and $m$ is the atom mass). In connection to this problem, Tan Tan (2008) has recently shown that the total energy can be entirely written in terms of the momentum distribution $n_{k\sigma}$ and a parameter, $C$, which he termed ‘contact’. The latter controls the asymptotic behavior $\sim p^{-4}$ of the momentum distribution $n_{p\sigma}$.

In many-body perturbation theory, the removal (renormalization) of the divergences described above begins by recognizing that the coupling $g$ in Eq. (1) is not physical and must be replaced by the measurable $s$-wave scattering length, $a_{s}$. Thus, Eq. (1) can be used to compute the two-particle scattering amplitude (see Eq. 3 below). We can parametrize the latter in terms of the $s$-wave phase shift $\delta_{s}(K)$. In the limit of small momentum exchange of the colliding particles $\boldsymbol{K}=\boldsymbol{p}-\boldsymbol{r}$, the phase shift obeys:

[TABLE]

This expression allows to relate the coupling $g$ in Eq. (1) to the scattering length $a_{s}$ ($R$ is the effective range). An equivalent treatment relies on perturbation theory to relate $g$ to $a_{s}$ through the following expression for the two-particle $T$-matrix Abrikosov et al. (1975a); Pathria and Beale (1996a):

[TABLE]

In the above equation, $\sigma\neq\alpha$ and $E_{pkqr}=\epsilon_{p}+\epsilon_{k}-\epsilon_{q}-\epsilon_{r}$.

In this paper, we study how to remove the divergences that appear in the expressions describing the nonequilibrium dynamics of the total energy following an interaction quantum quench. We consider a fairly general quench in which the system evolution is dictated by the Hamiltonian:

[TABLE]

where $H_{0}=\sum_{\boldsymbol{p}\sigma}\epsilon_{p}c^{{\dagger}}_{\boldsymbol{p}\sigma}c_{\boldsymbol{p}\sigma}$ is the kinetic energy and $V$ is given in Eq. (1). The function $S(t)$ describes the experimental protocol followed to turn on the interaction starting from the non-interacting system (accounting for the effects of a weak initial interaction can be also done perturbatively and will be reported elsewhere Huang and Cazalilla (2018)). We will see that the divergences that plague the calculation of the ground state energy in equilibrium are also present in the perturbative calculation of the total energy dynamics. Specializing to the case of a linear ramp, where

[TABLE]

we obtain the time-evolution of the total energy and the momentum distribution as a function of ramp time $T_{r}$. Tuning $T_{r}$ allows us to study the crossover from the sudden (i.e. $T_{r}\to 0$) to the quasi-adiabatic (i.e. $T_{r}\to+\infty$) limit.

In order to compute the total energy dynamics, we also compute the instantaneous momentum distribution, i.e.

[TABLE]

where $\rho(t)$ describes the state of the system at time $t$ (see Sec. II for details). The perturbative expansion for this quantity does not contain any divergent integrals. However, the leading perturbative corrections appearing at $O(g^{2})$ behave as $n_{p\sigma}\sim p^{-4}$ at large momentum $p$. Likewise, the dynamics of the total energy is obtained from:

[TABLE]

where the first term in the second line is the instantaneous kinetic energy, i.e.

[TABLE]

The high momentum tail $\sim p^{-4}$ of the instantaneous momentum distribution $n_{p\sigma}(t)$ at $p\gg k_{F}$ renders $E_{\mathrm{kin}}(t)$ divergent, as in the equilibrium case. Interestingly, when the coupling $g$ is replaced by the scattering length $a_{s}$ by the renormalization method described in Sec. III, the divergence $E_{\mathrm{kin}}(t)$ is cancelled by an additional contribution from the interaction energy $E_{\mathrm{int}}(t)=S(t)\mathrm{Tr}\>\left[\rho(t)V\right]$. This cancellation parallels the similar cancellation happening in equilibrium and allows us to isolate the leading term in perturbative expansion of the tail of the instantaneous momentum distribution. For a linear interaction ramp, we find that the definition of Tan’s contact requires some care in order to fully capture the crossover behavior of the tail from the sudden to the adiabatic limit.

Our results also allow to explore the dynamics of the two-component Fermi gas following linear ramp in the interaction strength. A subject of particular interest in this situation is existence of a pre-thermalized regime  Moeckel and Kehrein (2008, 2009); Nessi et al. (2014); Hamerla and Uhrig (2014); Eckstein et al. (2009); Nessi and Iucci (2015); Marcuzzi et al. (2013); Essler et al. (2014); Langen et al. (2016); Moeckel and Kehrein (2010); Rademaker (2017); Mitra (2013); Cazalilla (2006); Buchhold et al. (2016); Gong and Duan (2013); Kollar et al. (2011); Cazalilla and Chung (2016); Babadi et al. (2015); Alba and Fagotti (2017); ping zou and Zhang (2018). Thus, we have explored the existence of a pre-thermalized regime as a function of the ramp time. Previous studies on pre-thermalization in ultracold Fermi gases have focused on the behavior of the momentum distribution near the Fermi momentum $k_{F}$ at zero temperature Moeckel and Kehrein (2008, 2009); Nessi et al. (2014); Hamerla and Uhrig (2014); Eckstein et al. (2009); Nessi and Iucci (2015); Kollar et al. (2011); Cazalilla and Chung (2016); Babadi et al. (2015). It was concluded that the persistence of a discontinuity at $k_{F}$ at the same time that the total energy has reached its final (thermalized) value Berges et al. (2004) characterizes the pre-thermalized regime. Here we report results for the dynamics of the full momentum distribution at finite temperatures as well as large momenta, which can provide a more experimentally accessible way to characterize the pre-thermalized regime. This is because in realistic systems, the discontinuity of the momentum distribution at $k_{F}$ is absent due to finite temperature effects and trap confinement. On the other hand, as we argue below, the dynamics of the full momentum distribution at finite temperature and its asymptotic behavior at high momenta contains a great deal of useful information about the pre-thermalized regime.

The rest of this article is organized as follows. In section II, we give the derivation of observables. In section III, we discuss the renormalization procedure in and out of equilibrium. In section IV, we specialize to the a linear ramp quench and study the dynamics of total energy and the momentum distribution to show the emergence of pre-thermalization in short time. In section V, we derive the contact for a ramp of interaction in the interaction strength and discuss its dynamics and relation to pre-thermalization. In Section VI, we summarize our results and present the conclusions of this work. Throughout, we use units where $\hbar=1$ and $k_{B}=1$.

II Evolution of Observables

In what follows, we use perturbation theory to obtain the short to intermediate time dynamics of the system as the interaction is quenched. To leading order in the scattering length, $a_{s}$, this treatment is valid when the strength of interaction being quenched is small. Indeed, earlier work Moeckel and Kehrein (2008); Nessi et al. (2014); Cazalilla and Chung (2016); Moeckel and Kehrein (2009, 2010); Chiocchetta et al. (2016); Gramsch and Potthoff (2015); Kollar et al. (2011); Tsuji and Werner (2013); Stark and Kollar (2013), it has been established that pre-thermalization is accessible through perturbation theory in the quenched interaction. We calculated the observables to second order in the quench interaction. The valid time scale therefore fulfills the condition $E_{F}t\sim(k_{F}a_{s})^{-2}\ll(k_{F}a_{s})^{-3}$, wheres $E_{F}=k^{2}_{F}/2m$ is the Fermi energy. This requires $k_{F}a_{s}\ll 1$ in order to apply to hold at long times.

In this section, we provide the details of the derivation of the time evolution of the observables considered in this work. In the interaction picture, the evolution of an observable can be calculated from the following expression, which is amenable to perturbative expansion:

[TABLE]

In Eq. (9), $\tilde{V}(t)=e^{iH_{0}t}Ve^{-iH_{0}t}S(t)$ is the quench interaction in interaction picture, and $O=c^{{\dagger}}_{\boldsymbol{p}\sigma}c_{\boldsymbol{p}\sigma},V,\ldots$ stands for the observable of interest. Time $t$ lies on the closed contour $C$ shown in Fig. 1 and $\mathcal{T}$ is the time-ordering symbol on $C$. Since we are working with a closed time contour, trictly speaking, the denominator of Eq. (9) equals unity. However, it is needed when expanding in powers of $\tilde{V}(t)$ in order to cancel disconnected Feynman graphs. The expectation values are computed according to $\langle...\rangle=\text{Tr}[e^{-\beta H_{0}}...]$, which is the thermal average with respect to an non-interacting initial Hamiltonian at absolute temperature $T=\beta^{-1}$.

Expanding Eq. (9) in powers of the quenched interaction $\tilde{V}(t)$ yields:

[TABLE]

As mentioned above, we shall only take into account the fully connected contributions (denoted by $\langle\ldots\rangle_{c}$ in the expression above) resulting from the application of Wick’s theorem. When using the latter, there are four possible choices of time arguments for the fermion propagator:

[TABLE]

where $t_{1}$ and $t_{2}$ can be either in the $\tau$ or $\bar{\tau}$ branches of $C$. The free fermion propagator can be written in matrix form as follows:

[TABLE]

Using $c_{\boldsymbol{p}\sigma}(t)=c_{\boldsymbol{p}\sigma}e^{-i\epsilon_{p}t}$, the elements of the above matrix can be evaluated to yield:

[TABLE]

II.1 Instantaneous momentum distribution

Let us first consider the evolution of the momentum distribution. The latter is obtained from the Eq. (9) by setting $O=\hat{n}_{p\sigma}=c^{{\dagger}}_{p\sigma}c_{p\sigma}$. Hence, expanding the evolution operator up to second order in the quenched interaction, the following expression is obtained:

[TABLE]

The expectation values in the above expressions can be represented diagrammatically. In the closed time contour, we can express the expectation value using the self-energy and propagator matrices, which to second order in the quenched interaction yield:

[TABLE]

The propagator $\mathcal{G}_{p\sigma}(a,b)$ is defined in Eq. (12), and $\Sigma^{(1)}_{p\sigma}(t_{1})$, $\Sigma^{(2)}_{p\sigma}(t_{2},t_{1})$ can be found using diagrams (see below). For the calculation of equal-time expectation values, we choose the time argument of the observable (i.e. $t$) to lie slightly before the turning point of the contour $C$, which is on the time ordered (i.e. $\tau$) branch. In this case, the fermion propagators must be obtained from Eq. (11), which yields:

[TABLE]

where the non-vanishing entries correspond to either $b$ lying before or after $t$ on the contour $C$. Similarly,

[TABLE]

The two non-zero entries in the above matrix correspond to $a$ lying before or after $t$ on the contour $C$.

The self-energy can be calculated from the diagrams shown in Fig. 2 and the propagators, Eqs. (13) to (16). Thus, to first order in $V(t)$, we obtain:

[TABLE]

Combining the matrix from Eq. (LABEL:eq:n), the propagators, Eq. (19) and Eq. (20), and the self-energy for the first order correction, Eq. (21), we obtain that first order correction to the instantaneous momentum distribution vanishes, i.e. $n_{p\sigma}^{(1)}(t)=0$. Ultimately, this is a consequence of the initial state being an eigenstate of the occupation operator $\hat{n}_{p\sigma}=c^{{\dagger}}_{\boldsymbol{p}\sigma}c_{\boldsymbol{p}\sigma}$.

At second order in the interaction, we need to use the following self-energy matrix, which contains four different combinations of the time arguments $(t_{1},t_{2})$ on the two branches of the contour $C$:

[TABLE]

Evaluating the diagrams in Fig. 2 using the propagators in Eqs. (13) to (16), we obtain the following expressions for the elements of the self-energy matrix:

[TABLE]

Combining Eq. (LABEL:eq:n), the propagators (cf. Eq.  13 to Eq. 16) and the second order corrections to the self-energy, Eqs. (23) to (26), we arrive at:

[TABLE]

where $A^{\sigma}_{pkqr}$ denotes the function:

[TABLE]

In Eq. (27), $E_{pkqr}=\epsilon_{p}+\epsilon_{k}-\epsilon_{q}-\epsilon_{r}$, and the function

[TABLE]

has been introduced. Note that $F^{(2)}(E,t)$ depends on the explicit form of $S(t)$ (see Sec. IV and Appendix A for the form of this function in a number of important limiting cases).

For $p>k_{F}$, the above expression, Eq. (27), reduces to

[TABLE]

Furthermore, momentum conservation requires that $\boldsymbol{k}=\boldsymbol{q}+\boldsymbol{r}-\boldsymbol{p}$. Since the occupation factors in this expression force $q,r\leq k_{F}$, it follows that $|\boldsymbol{q}+\boldsymbol{r}|\leq 2k_{F}$. Thus, for $p\gg k_{F}$, $k\gg k_{F}$, and the above expression simplifies to

[TABLE]

It will be shown in Sec. V that this expression leads to a $\sim p^{-4}$ tail of the momentum distribution, which renders the kinetic energy divergent.

II.2 Unrenormalized kinetic energy

From the above result for the instantaneous momentum distribution, we can derive the second order correction to the kinetic energy:

[TABLE]

where, in the last line, we have used that $\sum_{\sigma}A^{\sigma}_{pkqr}=\sum_{\sigma}A^{\sigma}_{kpqr}=-\sum_{\sigma}A^{\sigma}_{qrpk}=-\sum_{\sigma}A^{\sigma}_{rqpk}$ and $F^{(2)}(-E,t)=F^{(2)}(E,t)$. Notice that in two extreme limits, i.e. the sudden and adiabatic limits $F^{(2)}(E,t)\sim 1/E$ for $E\to+\infty$ (see Appendix A). Thus, the above expression is divergent, as in equilibrium (the expression of the adiabatic limit coincides with the equilibrium one). Generically, for any other quench protocol function $S(t)$ between these two extreme limits, we expect the same behavior and the above expression to remain divergent, as it is also confirmed for a linear ramp quench, see Sec. IV.

II.3 Unrenormalized interaction energy

The interaction energy can be obtained from Eq. (9) by setting $O=\tilde{V}(t)$. Expanding this expression to second order in $g$ yields

[TABLE]

The calculation in this case proceeds in a similar fashion to the one in previous section. Evaluating the Feynman diagrams in Fig. 3, we obtain:

[TABLE]

Note that the first order term, $V^{(1)}(t)$, Eq. (35), is proportional to the function $S(t)$. The second order term, $V^{(2)}(t)$ can be obtained from the matrix:

[TABLE]

In the last line, we have used that $t$ lies on the time-ordered branch of the contour $C$. Using the free propagator (cf. Eqs. 13 to 16), the second order term becomes:

[TABLE]

From the above expression, using Eq. (36), we obtain the second order correction to the interaction energy:

[TABLE]

In the last line, we have introduced the function

[TABLE]

which depends on the function $S(t)$ that defines the quench protocol. The sine function in Eq. (41) appears after swapping the dummy momenta around in order to symmetrize the integrand, that is, after using:

[TABLE]

Again, notice that that in the sudden and adiabatic limits, $F^{(1)}(E,t)\sim 1/E$ for $E\to+\infty$ (see Appendix A), which means that Eq. (40) is divergent, as in equilibrium (the expression of the adiabatic limit coincides with the equilibrium result, see Appendix. A). Generically, for any other quench protocol function $S(t)$ between these two extreme limits, we expect the above expression to remain divergent. This is explicitly confirmed for a linear ramp protocol in Sec. IV.

II.4 Unrenormalized total energy

From the results obtained in previous sections for the interaction and kinetic energy, we can obtain the dynamics of the total energy for a quench in which the interaction is switched on according to an arbitrary protocol described by $S(t)$. The first order correction is:

[TABLE]

and the second order correction reads:

[TABLE]

where, in the last line, we have introduced the function:

[TABLE]

and used that $F_{\mathrm{tot}}(-E,t)=-F_{\mathrm{tot}}(E,t)$ to obtain the expression in the last line of (44). Furthermore, notice that, as shown in Appendix A, $F_{\mathrm{tot}}(E,t)$ vanishes in the limit of a sudden quench. This ensures that the second order correction to the total energy vanishes, as required by energy conservation. See Eq. (94) in Appendix A and Sec. IV for a more in depth discussion of this point.

III Elimination of divergences

As we have mentioned in the Introduction, perturbation theory in the powers of the bare interaction $V\propto g$ (cf. Eq. 1) yields an expression for the equilibrium total energy containing divergent integrals Abrikosov et al. (1975b); Pathria and Beale (1996b). The divergences appear at second in the coupling $g$ and their elimination, i.e. their ‘renormalization’, is possible by realizing that the unphysical $g$ must be replaced by the physical $s$-wave scattering length, $a_{s}$. The latter is related to $g$ via the two-particle scattering amplitude (i.e. the $T$-matrix), see Eq. (3).

The same divergences reappear in the perturbative expansion for the dynamics of the total energy obtained in previous sections. Indeed, as we show by explicit calculation below, the kinetic energy is divergent because the leading order [i.e. $O(g^{2})$] correction to instantaneous momentum distribution behaves as $\sim p^{-4}$ for $p\gg k_{F}$. Hence, $E_{\mathrm{kin}}(t)=\sum_{p\sigma}\epsilon_{p}n_{p\sigma}(t)$ is divergent at all times $t$. Similar divergences appear in the expressions for the interaction energy at $O(g^{2})$ as well.

However, in the nonequilibrium case, it is difficult to rely on the $T$-matrix because defining the latter requires the introduction of asymptotic scattering states. Those states are well defined for interactions that are switched on and off adiabatically, as it is assumed in equilibrium, but this becomes difficult when the interaction changes (sometimes, very rapidly) in time as it is the case of our study. Thus, the generalization of the renormalization procedure employed in equilibrium to the problem of interest here is not straightforward. Instead, we show below that the renormalization procedure can be carried out by computing the perturbative corrections to the evolution of the total two-particle energy.

III.1 Renormalization in equilibrium

Let us begin by reviewing how the divergences are eliminated in equilibrium case. As mentioned above, we shall compute the shift of the total energy for two particles to relate $g$ to the $s$-wave scattering length, $a_{s}$ Pathria and Beale (1996a). This approach reproduces the well-known equilibrium results based on the scattering matrix approach Pathria and Beale (1996a); Abrikosov et al. (1975b). The shift to the total energy for two-particles is defined by:

[TABLE]

where $|\Psi^{0}_{\boldsymbol{kp},\sigma\alpha}\rangle=c^{{\dagger}}_{\boldsymbol{p}\sigma}c^{{\dagger}}_{\boldsymbol{k}\alpha}|0\rangle$ describes a state of two free particles with $\sigma\neq\alpha$, and $|\Psi_{\boldsymbol{pk},\sigma\alpha}\rangle$ is the two-particle state perturbed by the interactions. Using time-independent perturbation theory,

[TABLE]

where $E_{kpqr}=\epsilon_{k}+\epsilon_{p}-\epsilon_{q}-\epsilon_{r}$. Hence, the shift of the kinetic energy is

[TABLE]

To the same order in $g$, the shift to the interaction energy reads:

[TABLE]

Thus, the total energy shift is

[TABLE]

Let us define the physical scattering amplitude by requiring that:

[TABLE]

after equating it to Eq. (50), we arrive at Eq. (3). Inverting the series gives the coupling $g$ in terms of the scattering length

[TABLE]

which allows to remove the divergences in the expressions for the many-particle ground state energy Abrikosov et al. (1975a); Pathria and Beale (1996a). Note that the same result can be obtained directly from the perturbative expression for the total energy shift Pathria and Beale (1996a). However, we have chosen this more cumbersome method in order to emphasize that using the interaction energy instead would introduce in Eq. (52) a spurious factor of two (compare Eqs. (49) and (50)), which would not remove the divergences.

III.2 Nonequilibrium renormalization

Let us now turn to the nonequilibrium situation. In the calculation described in Sec. II, the evolution of the total energy shift is obtained in a perturbative series in the coupling $g$:

[TABLE]

where $E^{(n)}_{\mathrm{tot}}(g,t)=O(g^{n})$. However, this expansion contains divergent integrals. Our goal in this section is to remove the divergences, which can be carried out after replacing $g$ by $a_{s}$, in a procedure similar to the one outlined above for the equilibrium case.

Adding up the first and second order corrections obtained in Sec. II yields:

[TABLE]

For two-particles, the evolution of the total energy shift can be obtained in a similar fashion, which leads to:

[TABLE]

Next, we require that:

[TABLE]

Inverting the series will gives the coupling $g$ in terms of the scattering length 222Strictly speaking the coupling $g$ should be (weakly) time dependent. This is reminiscent of the equilibrium situation where Eq. 52 would require $g$ to be (weakly) energy dependent.:

[TABLE]

In the limit where the interaction is switched adiabatically, notice that $F_{\text{tot}}(E)=2E^{-1}$ (see Appendix A) and thus we recovers the equilibrium result, Eq. (52).

Substituting Eq. (57) into the first order correction for total energy shift, Eq. (54), yields the following correction:

[TABLE]

which is second order in the scattering length $a_{s}$ and divergent. However, this divergence exactly cancels the one already present in the second order correction to the energy, i.e. $E_{\mathrm{tot}}^{(2)}(g^{2},t)$ (after replacing $g$ by $(4\pi a_{s}/m)$). Thus, the finite expression for the dynamics of the total energy to second order in $a_{s}$ is :

[TABLE]

Furthermore, according to the dependence on $S(t)$ and time, we can split $\Delta_{\mathrm{tot}}(t)$ into the sum of $\Delta_{\mathrm{kin}}(t)\sim F^{(2)}(E,t)$ and $\Delta_{\mathrm{int}}(t)\sim F^{(1)}(E,t)$. Explicitly,

[TABLE]

Symmetrizing the dependence on momenta $\boldsymbol{p},\boldsymbol{k},\boldsymbol{q}$ and $\boldsymbol{r}$ of the integrand with the help of $F^{(2)}(-E,t)=F^{(2)}(E,t)$, $\Delta_{\text{kin}}(t)$ can be written as

[TABLE]

where

[TABLE]

After setting $g=4\pi a_{s}/m$ in Eq. (32) and subtracting $\Delta_{\text{kin}}(t)$, the divergence in the kinetic energy is cancelled and it is possible to obtain a finite result for the dynamics of the kinetic energy. Furthermore, for $p\gg k_{F}$ (compare Eq. (63) and Eq. (31) after setting $g=4\pi a_{s}/m$), we obtain:

[TABLE]

This explains why adding $-\Delta_{\mathrm{kin}}(t)$ to the kinetic energy cancels the divergence which arises from the behavior of the instantaneous momentum distribution at large momenta. Indeed, this result provides the basis for our analysis of the dynamics of Tan’s contact following the interaction quench, which is presented in Sec. V. In the following section, we shall evaluate the above expressions for a linear ramp of the interaction strength, for which $S(t)$ is given by Eq. (5).

IV Results for a linear ramp quench

In this section, we obtain explicit results for quenches in which the interaction strength is ramped up linearly in time and $S(t)$ is given by Eq. (5). As we have seen above, the first order correction to the total energy is just the expectation value of $V$ in the initial state (cf. Eq. 1) multiplied by $S(t)$. At second order in the scattering length $a_{s}$, using Eq. (59) and evaluating the time integrals for the functions $F^{(1)}(t)$ (cf. Eq. 41) and $F^{(2)}(t)$ (cf. Eq. 29) yields the following expression for the time dependence function, $F_{\mathrm{tot}}(E,t)=F_{\mathrm{tot}}(E,T_{r},t)$:

[TABLE]

For the purpose of the discussion below, we have introduced a scaling function, $G(x,y)$, along with $F_{\mathrm{tot}}^{\text{eq}}(E)=E^{-1}$. The latter is the form that $F_{\mathrm{tot}}(E,t)$ takes in the adiabatic limit (see Appendix A).

Let us next consider the form of $F_{\mathrm{tot}}(E,T_{r},t)$ in specific limiting cases. For a sudden quench, $T_{r}\to 0$ and $G(Et,ET_{r})=0$. Therefore,

[TABLE]

which implies that, to second order, the shift to total energy vanishes. Indeed, this is expected from energy conservation in the limit of a sudden quench. To this this, let us compute the total energy directly for a sudden quench:

[TABLE]

where $E^{(0)}_{\text{tot}}=\langle H_{0}\rangle$ is the total energy before the quench (and $\langle\ldots\rangle$ is the average over the initial state). In other words, the dynamics of the total energy in the sudden quench limit simply amounts to a constant energy shift happening at $t=0$. The shift is the expectation value of the interaction in the initial state, which is first order in $a_{s}$. Hence, not only the $O(a^{2}_{s})$ term vanishes (in agreement with our findings) but so do all higher order corrections.

However, in the adiabatic limit where $T_{r}\to+\infty$, we have

[TABLE]

This means that at times $t<T_{s}$, there is a transient during which energy grows quadratically in time. The growth saturates to the equilibrium value for $t>T_{r}$, i.e. the interaction strength reaches its full value.

We have numerically evaluated the momentum integrals and obtained the evolution the second-order correction to the total energy and the instantaneous momentum distribution following a linear ramp in the interaction strength for a uniform Fermi gas in three dimensions. The left panel of Fig. 4 shows the time evolution of second order correction to the total energy for different ramp times $T_{r}$. As the sudden quench limit where $T_{r}\to 0$ is approached, the second order energy shift vanishes, as explained above. In the opposite limit, as the ramp time $T_{r}$ increases, second-order energy shift saturates at a value that approaches the equilibrium energy shift after initially growing quadratically at short times.

Concerning the dynamics of the instantaneous momentum distribution, as discussed in Sec. II.1, it is described by the function $F^{(2)}(E,T_{r},t)$. For a linear ramp in the interaction strength, this function can be written as follows (see Appendix A for the details) :

[TABLE]

where

[TABLE]

is a scaling function of the dimensionless variables $x=ET_{r}$ and $y=Et$ and $F^{(2)}_{\text{eq}}(E)=1/E^{2}$ is the equilibrium value of $F^{(2)}(E,t)$. Thus, as the sudden quench limit where $T_{r}\to 0$ is approached

[TABLE]

Therefore, we see that the second order correction to instantaneous momentum distribution reaches a stationary value that obeys the following relation:

[TABLE]

where $n^{(2,\mathrm{eq})}_{k\sigma}$ is the second order correction to the equilibrium momentum distribution of the Fermi gas with interaction strength equal to $g=4\pi a_{s}/m$. This result can be regarded as a generalization of the relationship obtained in Ref. Moeckel and Kehrein (2008, 2009) between the discontinuity of the momentum distribution at Fermi surface i.e. $Z(t)=\lim_{\delta\to 0}\left[n_{k_{F}+\delta,\sigma}(t)-n_{k_{F}-\delta,\sigma}(t)\right]$ in the pre-thermalized regime and in equilibrium:

[TABLE]

where $Z_{\mathrm{st}}$ and $Z_{\mathrm{eq}}$ stand for the stationary (i.e. pre-thermalized) and equilibrium values of $Z$, respectively.

On the other hand, in the adiabatic limit where $E_{F}T_{r}\gg 1$, we have

[TABLE]

Hence, the long time behavior of the instantaneous momentum distribution approaches the equilibrium momentum distribution with an interaction strength $g=4\pi a_{s}/m$, i.e.

[TABLE]

as expected.

In the right-hand panel of Fig. 4, we have plotted the second-order correction to the instantaneous momentum distribution $n^{(2)}_{p,\sigma}(t)$ for $\sigma=\uparrow$. For reasons of experimental interest, we consider an initial (non-interacting) state at finite temperature $T=0.1T_{F}$, where $T_{F}=E_{F}=k^{2}_{F}/2m$ is the Fermi temperature (in $\hbar=k_{B}=1$ units), rather than the zero-temperature ground state, as considered in previous studies Moeckel and Kehrein (2008); Kollar et al. (2011); Nessi et al. (2014). Notice that the momentum distribution reaches a stationary value for $t/T_{r}\gg 1$. This follows from the asymptotic behavior of $F^{(2)}(E,T_{r},t)$ for $t/T_{r}\gg 1$:

[TABLE]

Thus, for any finite $T_{r}\neq 0$, the long time limit of instantaneous momentum distribution for $k\sim k_{F}$ (see dashed line in Fig. 4) differs from the equilibrium value $n^{(2,\textrm{eq})}_{k}$.

V Dynamics of the contact

In the previous section, we have studied the behavior of the perturbative corrections to the total energy and the instantaneous momentum distribution $n_{k\sigma}(t)$ for $k\sim k_{F}$ ($k_{F}$ being the Fermi momentum). In this section, we turn our attention to the asymptotic behavior of the latter for $k\gg k_{F}$. In this regime, $n_{k\sigma}(t)$ exhibits a $k^{-4}$ decay, similar to its behavior in equilibrium Tan (2008). However, the pre-factor (the non-equilibrium version of Tan’s contact) depends both on time, $t$ and the energy, $\epsilon_{k}=k^{2}/2m$.

In connection with the asymptotic behavior of the instantaneous momentum distribution, in Sec. III.2 we have shown that the correction that removes of the divergences from the perturbative expression for the total energy factorizes into a kinetic and an interaction-energy contribution. The correction to the kinetic energy, $\Delta_{\mathrm{kin}}(t)$ can be written in terms of the function $\delta n_{k}(t)$ (see Eq. (62) and following equations). As shown in Sec. III.2, this function also yields the asymptotic behavior of the instantaneous momentum distribution, $\sum_{\sigma}n_{k\sigma}(t)$ for $k\gg p_{F}$ (see Eq. 64). Hence, it is possible to extract the dynamics of Tan’s contact by studying the asymptotic behavior of Eq. (63). To this end, let us rewrite $\delta n_{k}(t)$ as follows:

[TABLE]

where have implemented momentum conservation by requiring that $\boldsymbol{p}=\boldsymbol{q}+\boldsymbol{r}-\boldsymbol{k}$. In order to evaluate the above integrals numerically in the large $k$ limit, it is convenient to parametrize Doggen and Kinnunen (2015) $\boldsymbol{q}=(\boldsymbol{p}+\boldsymbol{k})/2+\boldsymbol{s}$, $\boldsymbol{r}=(\boldsymbol{p}+\boldsymbol{k})/2-\boldsymbol{s}$. The result of the numerical evaluation of $\delta n_{k}(T_{r},t)$ is shown in Figure 5 as a function of $k/k_{F}$ for $t=20E^{-1}_{\mathrm{F}}$. For this time, $\delta n_{k}(T_{r},t)$ has reached a stationary value. At large $k/k_{F}$, we observe that it exhibits $1/k^{4}$ dependence. However, the detailed asymptotic behavior of $\delta n_{k}(T_{r},t)$ is quite rich and also depends on the ramp time $T_{r}$, as shown in Fig. 5. Indeed, a careful comparison the results for $T_{r}=6E^{-1}_{\mathrm{F}}$ and $T_{r}\to 0$ is worth here. For the former value of $T_{r}$, the asymptotic $\sim k^{-4}$ behavior is reached for smaller values of $k/k_{F}$ than for the $T_{r}\to 0$ results, which approach the sudden quench. In addition, close inspection of the curves shows that $\delta n_{k}(k,T_{r}\to 0)$ first approaches the $k^{-4}$ with a different coefficient, and only at large $k/k_{F}$ finally converges to the same asymptote as $\delta n_{k}(T_{r}=6E^{-1}_{F},t)$. We explain this behavior below.

Since the asymptotic behavior of the $\delta n_{k}(T_{r},t)$ appears to be rather complex, a naïve generalization of Tan’s contact from the equilibrium case, i.e.

[TABLE]

does not suffice to fully capture it. Thus, in order to illustrate this point, let us analyze the asymptotic behavior of the instantaneous momentum distribution in the sudden quench and adiabatic limits.

V.1 Sudden quench limit

In the sudden quench limit, using the following property (cf. Eqs. 71 and Eq. 72):

[TABLE]

we obtain the following behavior for the $k^{-4}$ asymptote:

[TABLE]

where $C_{\mathrm{eq}}$ is the equilibrium contact calculated to leading order in perturbation theory (see Ref. Doggen and Kinnunen (2015) and below). The above relation is similar to the relation obeyed in the pre-thermalized regime by the discontinuity of the momentum distribution at the Fermi momentum, Eq. (74).

V.2 Adiabatic limit

In the adiabatic limit where $\epsilon_{k}T_{r}\to\infty$, we obtain the following expression for the asymptote:

[TABLE]

Thus, like the total energy energy there is a transient for $t<T_{r}$ where it grows quadratically in time. However, after the interaction reaches its full value at $t\geq T_{r}$, it saturates to the equilibrium value, $C_{\mathrm{eq}}$, in the steady state.

V.3 In between

Finally, let us consider the general case of a finite ramp time $T_{r}>0$. The left panel in Fig. 6 shows the time evolution of the asymptote of $k^{4}\delta n(k,T_{r},t)$ normalized to the equilibrium contact for different ramp times, $T_{r}$ (time is measured in $\epsilon^{-1}_{k}$ units). At short times, the asymptote of $k^{4}\delta n(k,T_{r},t)$ oscillates with a period $t_{\mathrm{os}}\sim 2\left(k/k_{F}\right)\epsilon^{-1}_{k}$. As $k\to\infty$ the oscillation dies out and, at large $\epsilon_{k}t$, the asymptote of $\delta n_{k}(T_{r},t)/C_{\mathrm{eq}}$ approaches a constant. The value of the latter varies between $1$ and $2$, depending on the value of $T_{r}$. The dependence on $T_{r}$ is shown on the right panel of Fig. 6, which illustrates how the asymptotic behavior of $k^{4}\delta n_{k}(T_{r},t)$ crosses over from the adiabatic limit where

[TABLE]

to the sudden-quench limit where

[TABLE]

The crossover happens for $T_{r}\sim\epsilon^{-1}_{k}$. The explanation for this behavior is as follows: For large but fixed $k\gg k_{F}$, any ramp of the interaction strength in a time much longer than $\epsilon^{-1}_{k}$ is regarded by the fermions at momentum $\boldsymbol{k}$ as adiabatic and therefore the behavior of asymptote approaches the adiabatic limit. However, for $T_{r}\ll\epsilon^{-1}_{k}$, the fermions at momentum $\boldsymbol{k}$ experience the quench as sudden, and therefore the asymptote approaches the sudden limit.

VI Conclusions

In conclusion, we have discussed the renormalization of total energy obtained from the single-channel model in the context of interaction quenches. The resulting expressions are finite and, when evaluated for a linear ramp in the interaction strength starting from an non-interacting state in a two-component Fermi gas, allow us to study the crossover from the sudden-quench to the adiabatic limit. In addition, we have also studied the behavior of instantaneous momentum distribution. Thus, we found signatures of the pre-thermalization emerging at short to intermediate times after a linear ramp of the interaction. We have shown that the pre-thermalization signatures persist even at finite temperatures and they are also visible in the high-momentum tail of the distribution function. These results are important for the experimental characterization of this nonequilibirum state.

We have analyzed the dynamics of the high momentum tail of the instantaneous momentum distribution, which is related to the non-equilibrium dynamics of the Tan’s contact. Thus, we have uncovered an interesting crossover from adiabatic to sudden-quench dynamics in the asymptotic behavior of the momentum distribution as a function of the ramp time and the momentum scale of the fermions at the tail. Although explicit results were obtained only for a specific quench protocol that assumes a linear ramp of the interaction strength, they should also hold for more general quenches provided the ramp time $T_{r}$ is replaced by the relevant switching-on time scale of the quench protocol.

The results reported in this work can be experimentally verified by preparing a three-dimensional two-component Fermi gas with an accessible broad Feshbach resonance in a non-interacting state, and quenching it to an interacting state. The total energy dynamics can be accurately measured from the time-of-flight images obtained by turning the scattering length to zero before releasing the gas from the trap. In addition, the renormalization method described here can be readily applied to the computation of other interesting nonequilibrium properties such as the full work distribution. In addition, in future work it would be interesting to extending it beyond the second order in the scattering length.

Acknowledgements.

We thank Y. Takahashi and Y. Takasu, and S.-Z. Zhang for useful discussions and commennts. We acknowledge support from the Ministry of Science and Technology of Taiwan through grants 102-2112-M-007-024-MY5 and 107-2112-M-007-021-MY5, as well as from the National Center for Theoretical Sciences (NCTS, Taiwan).

Appendix A Time dependence of $F^{(1)},F^{(2)}$ and $F_{\mathrm{tot}}$

In the previous sections, we have used a number of results for the functions $F^{(1,2)}(E,t)$, $F_{\mathrm{tot}}(E,t)$. In this section, we obtain their form for a linear ramp as well as several other important limiting cases, namely the sudden-quench and adiabatic limits. Let us recall that

[TABLE]

describes the time dependence for the total energy, and

[TABLE]

describe the time dependence for interaction and kinetic energy (momentum), respectively. For a ramp quench, setting $S(t,T_{r})=\theta(t)[\theta(t-T_{r})+\theta(T_{r}-t)t/T_{r}]$, we can derive:

[TABLE]

They describe the time dependence of total energy shift:

[TABLE]

For the the sudden quench $S(t)=\theta(t)$. We found the two time dependent functions $F^{(1)}(E,t)$ and $F^{(2)}(E,t)$:

[TABLE]

from which, the function that controls the time dependence of total energy in the sudden quench limit:

[TABLE]

as required by the conservation of total energy (see discussion in Sec. IV, Eq. 68). In the equilibrium limit, using $S(t)=e^{-\eta|t|}$ with $\eta\to 0^{+}$ to describe the adiabatic switching of the interaction, we obtain:

[TABLE]

which lead to the equilibrium results. Hence,

[TABLE]