Contents

TIFR/TH/18-16

Quantum quench and thermalization of one-dimensional Fermi gas via phase space hydrodynamics

Manas Kulkarni111manas.kulkarni@icts.res.in, Gautam Mandal222mandal@theory.tifr.res.in, and Takeshi Morita333morita.takeshi@shizuoka.ac.jp

a. International Centre for Theoretical Sciences

Tata Institute of Fundamental Research, Shivakote, Bengaluru 560089, India.

b. Department of Theoretical Physics

Tata Institute of Fundamental Research, Mumbai 400005, India.

c. Department of Physics

Shizuoka University, 836 Ohya, Suruga-ku, Shizuoka 422-8529, JAPAN.

d. Graduate School of Science and Technology

Shizuoka University, 836 Ohya, Suruga-ku, Shizuoka 422-8529, JAPAN.

July 13, 2019

By exploring phase space hydrodynamics description of one-dimensional free Fermi gas, we discuss how systems settle down to steady states described by the generalized Gibbs ensembles through quantum quenches. We investigate time evolutions of the Fermions which are trapped in external potentials or a circle for a variety of initial conditions and quench protocols. We analytically compute local observables such as particle density and show that they always exhibit power law relaxation at late times. We find a simple rule which determines the exponent of the relaxation. Our findings are, in principle, observable in experiments in one dimensional free Fermi gas or Tonk’s gas (Bose gas with infinite repulsion).

## 1 Introduction

It is an important issue to find out if and when a closed system, subjected to some disturbance— such as in a quantum quench, reaches equilibrium after some time. From a theoretical viewpoint, this is a fundamental question since it lies at the heart of quantum statistical mechanics; e.g. subtleties arise in the applicability of thermodynamics in systems showing many-body-localization. The question of thermalization is also important from an experimental viewpoint. In experiments dealing with ultra-cold atoms, the time scale of equilibration as well as the nature of the equilibrium are measurable quantities and help characterize the system in question. In the latter context, it has also been realized, for a little over a decade now, that there is a non-trivial notion of thermalization of local observables in a so-called free systems or integrable systems where the thermal ensemble is often characterized by an infinite number of chemical potentials, corresponding to an infinite number of conserved quantities. Such an ensemble is called a GGE (generalized Gibbs ensemble). These issues have been summarized in a number of articles in the recent years; for a partial list of references, see [1, 2, 3].

The issue of thermalization in integrable models [4, 5, 6, 7, 8, 9, 10] has been discussed in detail some time ago in the paper, Ref [11], where thermalization has been shown to occur under some general assumptions. In spite of such general arguments, it is important to see if one can obtain some explicit results about time evolution of observables, especially at long times. For quantum quenches leading to gapless Hamiltonians, such explicit results were obtained at long times in Ref. [12] which rigorously showed thermalization of the reduced density matrix of subsystems of a large system. Fully exact time-dependence and subsequent results on thermalization were obtained for free scalars and Fermions in 1+1 dimensions for mass quenches ending at zero mass in [13].

In this paper, we will discuss quantum quench in a free Fermi gas in one space dimension, by using a large- technique. This subject was dealt with earlier in a paper [14] by two of the present authors, where it was shown how moments of the Fermion density approached thermalization from particular time-dependent initial conditions. Apart from studies on thermalization, there has also been several works on nonlinear dynamics and certain hallmarks of it (such as onset of shock fronts) in Fermi gases in large- limit [15, 16, 17, 18] and finite- corrections [19].

The importance of the large- limit is that the Fermions are described by a semiclassical fluid in the phase space. As it turns out, for simple phase space configurations, the equations describing such phase space fluids boil down to equations of conventional hydrodynamics. Thus, it is tempting to think that thermalization can be understood somehow in terms of the equations of conventional hydrodynamics, which would be of significant interest. However, conventional hydrodynamics of the Fermions typically breaks down way before we reach asymptotic times because of formation of shock fronts [14, 15, 16, 17, 18, 19]. (See Figure 9 also.) We will show that one can easily proceed beyond such times when one sticks to methods of phase space hydrodynamics, where shock fronts and singularities in real space merely become folds in phase space which are smooth, continuous configurations and do not present any singularity. Thus, using the phase space formulation, we will derive results on thermalization of large- Fermion systems at large times.

The time evolutions of the one-dimensional free Fermi gases which are equivalent to the gases of the hardcore bosons (Tonks-Girardeau gas [20, 21, 22, 23]) are actively being investigated in cold atomic systems and condensed matter physics. If these gases are confined in an external potential or a periodic circle, they evolve to a steady state which is described by GGE. 444An exception is the motion in an external harmonic potential [24]. In this case, every particle moves by the common periodicity and the dephasing [11] does not occur. There, various observables have been investigated in various quench procedures [25, 26, 27, 28, 29, 29, 30, 31, 32]. We will see that our phase space hydrodynamics approach is quite powerful and makes it possible to reveal the general nature of the late time behaviour of the time evolution. Particularly we show that the local observables exhibit power law relaxation where the exponent is fixed through a simple universal rule. This rule works for arbitrary quench procedures and external potentials (except the harmonic one in which the relaxation does not occur). We demonstrate it in several examples which are summarized in Table 1. We will also find the entropy formula for the late time GGE state.

It is worth mentioning that the two main ingredients of our setup, namely external potentials and quenches are both experimentally realistic. Potentials and confinements of various kinds such as harmonic trap, quartic traps [33, 34, 35], ring-shaped traps [36, 37], box-like traps [38, 39, 40, 41, 42, 43], sinusoidal potentials have been realized. Various quench protocols have been successfully demonstrated [44, 45, 46].

The summary and organization of our paper are as follows. In Section 2, the phase space hydrodynamics method is introduced in the limit of large number () of Fermions. A simple parametric solution of the phase space density is presented in case where the confining potential vanishes. In Section 3, we continue the case of for motion of Fermions on a circle. In Subsection 3.1 we find explicit formulae for the Fermion density and exhibit power law relaxation where the exponent depends on the initial profile. In Subsection 3.1 we consider free motion after releasing the Fermi gas from a confining potential of the form . In this case the Fermion density approaches equilibrium according a power law determined by . In Subsection 3.3, the Fermions are released from a box, and the power law is now universal, viz. . In Section 4, we consider motion of Fermions in a potential at the post-quench stage. We find in Subsection 4.1 that there is power law relaxation even in this case. An explicit example of quenching from to a cosine potential is shown in Subsection 4.2. As claimed in the beginning, the post-quench reduced density matrix is expected to relax to that of a thermal or a GGE state. We explicitly verify this in Section 5 by comparing the time-evolving phase space density after a long time with that in a GGE; we also compute, in Subsection 5.1, the relevant entropy production. A discussion on conventional hydrodynamics and its breakdown is given in Appendix A.

## 2 Time evolution of free Fermions in a thermodynamical limit

In this section, we review the phase space formulation of the dynamics of a one-dimensional Fermi gas (which is equally applicable to a hard-core Bose gas in the Tonks limit).

Let us consider non-relativistic free Fermions in one dimension whose single particle Hamiltonian is given by

 ^h=−12ℏ2∂2x+V(x). (1)

Here is an external potential. The physics of this Fermi gas can be described by the second quantized Fermion field

 ψ(x,t)=∑n^cnχn(x)exp[−iEnt/ℏ],^hφn(x)=Enφn(x), (2) iℏ∂tψ(x,t)=−12ℏ2∂2xψ(x,t)+V(x)ψ(x,t), (3) ∫dx ψ†(x,t)ψ(x,t)=N. (4)

The dynamics of the Fermions can be alternatively be expressed in terms of the Wigner phase space density which is defined by

 u(x,p,t)=∫dη ψ†(x+η/2,t)ψ(x−η/2,t) exp[iηp/ℏ]. (5)

The constraint (4) translates to the following constraints in terms of the -variable:

 u(x,p,t)∗u(x,p,t)=u(x,p,t),∫dxdp2πℏu(x,p,t)=N, (6)

whereas the equation of motion (3) translates to

 ∂∂tu(x,p,t)+{h(x,p),u(x,p,t)}MB=0. (7)

The derivation of the above equations is straightforward, and is given in detail in [47, 48, 49]. Here we have used the notation

 f∗g(x,p) ≡[cosℏ2(∂p∂x′−∂p′∂x)(f(x,p,t)g(p′,x′,t))]p′=p,x′=x, (8) {f,g}MB =f∗g−g∗f. (9)

Indeed, it has been shown in [49] that the Fermion path integral can be entirely rewritten in terms of the variable, subject to the constraints (6).

It is easy to show that observables of the Fermi fluid, such as the density and the fluid velocity can be expressed in terms of the phase space density

 ρ(x,t)=∫dp2πℏu(x,p,t),v(x,t)=1ρ(x,t)∫dp2πℏ pu(x,p,t). (10)

### 2.1 Large N limit of Fermi gas and phase space hydrodynamics

We will define the large limit as

 N→∞,ℏ→0,Nℏ=1. (11)

It is easy to see that in this limit the star product in (8) becomes an ordinary product and the Moyal bracket (9) becomes a Poisson bracket. Thus, the phase space density in the large limit satisfies the equation of motion

 ∂∂tu(x,p,t)+{h(x,p),u(x,p,t)}PB=0, (12)

which is simply Liouville’s equation for the classical phase space density. With , this becomes 555If we substitute into (13), we obtain the classical Hamilton equation and . Thus the points inside the droplet correspond to the single Fermions.

 ∂tu+p∂xu−V′(x)∂pu=0. (13)

The constraints (6) become

 u2=u,∫dpdx2π u(x,p,t)=Nℏ=1. (14)

The first constraint implies that at any given phase space point the phase space density can either be =0 or =1. The regions where are called droplets, representing regions occupied by Fermions, 1 each in every small cell, of area . represents regions without Fermions (see Figure 1). The second constraint implies that the area of a droplet (or in case of multiple disconnected droplets, combined area of all droplets) is . As an example, the Fermi sea for a harmonic trap is represented by a circular droplet centred at the origin and of area 1.

### 2.2 u(x,p,t) in V=0 case

In the case of , the solution to the equation (13) and (14) is very simple. The solution of the constraint is given by using step function

 u(x,p,t)=θ(x+(p,t)−x)θ(x−x−(p,t)), (15)

where describe the boundaries of the droplet as functions of at time . See Figure 1. If the droplet has multiple boundaries for a given , for example see Figure 3, the corresponding step functions should be added. Then the equation (13) is satisfied if

 x±(p,t)=x0±(p)+pt. (16)

Here are the boundaries of the droplet at . This is because the Fermions move obeying the classical equations and . Thus the Fermions with move towards right and those with move towards left as sketched in Figure 1.

We need to choose this initial profile satisfying the second constraint of (14)

 ∫dxdp2π u(x,p,t=0)=∫dp2π(x0+(p)−x0−(p))=1. (17)

Once we impose this constraint at , Liouville’s theorem ensures that satisfies the constraint for any .

We can solve equation (13) and (14) in the case similarly, but the solution is a bit complicated. Hence we first consider the time evolution problem in the case and argue the case later.

## 3 Time evolution of particles on a circle in the V=0 case

We consider the time evolution of the particles in the case. To make the Fermions confined, we put the Fermions on a circle with a period and investigate how the particles settle down to a steady state (if it exists). By regarding the particle motion (16), we can speculate that the system would reach a steady state as shown in Figure 2. We will confirm this picture through detailed computations.

Due to the periodicity, we need to modify the previous result (15) as

 u(x,p,t)= ∑mθ(x+(p,t)−(x+mL))θ(x+mL−x−(p,t)) = ∑mθ(x0+(p)+pt−x−mL)θ(x+mL−x0−(p)−pt), (18)

where the summation represents the contributions of the mirrors. Now we apply Poisson summation formula

 ∑mf(mL)=∑k1L∫dzf(z)e−i2πkzL, (19)

to this equation and obtain

 u(x,p,t)= ∑k1L∫∞−∞dz θ(x0+(p)+pt−x−z)θ(x+z−x0−(p)−pt)e−i2πkzL = ∑k1L∫x0+(p)+pt−xx0−(p)+pt−xdz e−i2πkzL = x0+(p)−x0−(p)L+∑k≠012πik[exp(i2πkLz)]z=x0+(p)+pt−xz=x0−(p)+pt−x. (20)

This result is suggestive. While the second terms depend on time , the first term does not. Particularly, for large , the second terms are highly oscillating and may be irrelevant. Thus the first term may describe the late time steady state and the second terms may represent the damping terms.

To see it more concretely, we apply this formula to the particle density (10) and obtain

 ρ(x,t)= ∫dp2πℏu(x,p,t) = NL+∑k≠0∫dp4π2iℏk(exp(i2πkL(x0+(p)+pt−x))−exp(i2πkL(x0−(p)+pt−x))). (21)

Since the second terms would be suppressed at large through the integral due to the oscillation, we would obtain

 ρ(x,t)=NL,(t→∞). (22)

This result indicates that the particles uniformly spread over the circle independent of the initial profile, and this is exactly what we have expected in Figure 2. Thus it supports our conjecture that the first term in the formula (20) describes the late time steady state.

### 3.1 Power law relaxation

Now we evaluate the time dependent terms of (21) at large to see whether they are really suppressed and, if so, how they damp as the system evolves. We will show that it exhibits a power law relaxation and the exponent is fixed by the extrema of the initial droplet at . In the case of the droplet shown in Figure 1, the behaviors around the minimum and maximum are relevant. We investigate this case as an example.

Suppose that can be expanded around and as

 x0+(p) =x1+γ1(p−p1)α1+⋯,(p∼p1), x0+(p) =x2+γ2(p2−p)α2+⋯,(p∼p2), (23)

where and () are positive constants. If the droplet is smooth, are satisfied. See Figure 1 (Right). We evaluate the integral of the second term of (21) at large , and consider the third term later. We calculate the positive and negative cases separately. First we consider the positive . Since is defined in , the integral becomes

 ∞∑k=1∫p2p1dp4π2iℏkexp(i2πkL(x0+(p)+pt−x)) = ∞∑k=114π2iℏk[∫p1+i∞p1dpexp(i2πkL(x0+(p)+pt−x))−∫p2+i∞p2dpexp(i2πkL(x0+(p)+pt−x))] = ∞∑k=114π2ℏkei2πkL(x1+p1t−x)∫∞0dηe−2πkLtη(1+i2πkLγ1(iη)α1+⋯) −∞∑k=114π2ℏkei2πkL(x2+p2t−x)∫∞0dζe−2πkLtζ(1+i2πkLγ2(−iζ)α2+⋯). (24)

In the second line, we have changed the integral contour as shown in Figure 1 (Right bottom). Since is large, the contributions near the real axis would dominate and we ignore the integral along the horizontal line666Depending on , the integrand of (24) may diverge when . Besides there may be poles or cuts in the Im region. However, for sufficiently large , we can avoid them by taking the horizontal line of the integral contour in a finite region, since only the region contributes to the integral.. In the third line, we have defined the variable and . We have also used (23) and assumed that and are small, since the integral in the large and regions are exponentially suppressed. Note that the first term of the expansion in the integral is canceled by the same term coming from the integral (21) concerning . Hence the second term is the leading contribution. We can perform the integral of the second term by using the gamma function , and obtain

 12πℏL∞∑k=1(γ1Γ(α1+1)(iL2πkt)α1+1ei2πkL(x1+p1t−x)+γ2Γ(α2+1)(−iL2πkt)α2+1ei2πkL(x2+p2t−x)) = γ1Γ(α1+1)2πℏL(iL2πt)α1+1Liα1+1(ei2πL(x1+p1t−x))+γ2Γ(α2+1)2πℏL(−iL2πt)α2+1Liα2+1(ei2πL(x2+p2t−x)), (25)

where we have used the polylogarithm function Li.

For the negative in (21), we change the integral contour similar to Figure 1 but Im region, and obtain

 γ1Γ(α1+1)2πℏL(−iL2πt)α1+1Liα1+1(e−i2πL(x1+p1t−x))+γ2Γ(α2+1)2πℏL(iL2πt)α2+1Liα2+1(e−i2πL(x2+p2t−x)).

By adding them to (25), we finally obtain777These terms have cusp singularities at and . These are the location of the peaks of the droplet and the singularities correspond to the shock fronts at large . See Figure 9..

 γ12πℏL(Lt)α1+1ζ(−α1,x−x1−p1tL−⌊x−x1−p1tL⌋) +γ22πℏL(Lt)α2+1ζ(−α2,x2+p2t−xL−⌊x2+p2t−xL⌋), (26)

where we have used the Hurwitz’s formula

 i−sLis(e2πix)+isLis(e−2πix)=(2π)sΓ(s)ζ(1−s,x−⌊x⌋), (27)

where is the Hurwitz zeta function and is the floor function. Note that is a periodic function with a period 1. This is the leading contribution to from the integral of the terms in (21) at large and is suppressed by the power law as we announced.

Similarly we have the contributions from . If can be expanded near the extrema as

 x0−(p) =x1−γ−1(p−p1)α−1+⋯,(p∼p1), x0−(p) =x2−γ−2(p2−p)α−2+⋯,(p∼p2), (28)

where and () are positive constants888If the shape of the droplet is , and where are positive integers., we obtain the same expression to (26) but and .

Thus the density behaves as

 ρ(x,t)=NL+O(t−(α+1)),α=min(α1,α2,α−1,α−2). (29)

Note that a smaller means a flatter extremum of the initial droplet with respect to . Therefore exponent of late time power law relaxation is fixed by the flattest extremum.

If the initial droplet has more than two extrema as shown in Figure 3, we should divide the boundary at each extrema and define the boundary function between them. (Hence we need four boundary functions in the cases of Figure 3.) Then the calculations are almost same to the case and obtain a similar result to (26) 999If the extremum at is a minimum, we obtain (30) and if it is a maximum, we obtain (31) where and are defined through the expansion around the extremum similar to (23). . Therefore again the flattest extremum is relevant to determine the exponent of the late time relaxation.

#### Other local observables:

We can apply the calculation of to other local observables. We consider the following quantity

 F(x,t)≡∫∞−∞dp2πℏf(p,x)u(x,p,t), (32)

where is a smooth function. Actually various observables in our model are given in this form. For example, if we set and , we obtain the particle density and the velocity field through (10). If we take , we obtain the two point function as

 GF(x,y,t) =∑n=0(x−y)nn!⟨ψ†(y,t)∂nyψ(y,t)⟩ =∑n=0(x−y)nn!∫∞−∞dp2π(−ip)nu(y,p,t)=∫∞−∞dp2πe−i(x−y)pu(y,p,t). (33)

We calculate for the droplet in the case of Figure 1. We substitute (20) to (32) and perform the integral similar to (24). There, since the extrema and dominate in the integrals, we can approximate . Thus we obtain

 F(x,t)= ∫p2p1dp2πℏf(p,x)x0+(p)−x0−(p)L +γ12πℏLf(p1,x)(Lt)α1+1ζ(−α1,x−x1−p1tL−⌊x−x1−p1tL⌋) +γ22πℏLf(p2,x)(Lt)α2+1ζ(−α2,x2+p2t−xL−⌊x2+p2t−xL⌋) +γ−12πℏLf(p1,x)(Lt)α−1+1ζ(−α−1,x−x1−p1tL−⌊x−x1−p1tL⌋) +γ−22πℏLf(p2,x)(Lt)α−2+1ζ(−α−2,x2+p2t−xL−⌊x2+p2t−xL⌋)+⋯, (34)

at large . Again it shows the power law relaxation. Note that the exponent is independent of . Thus, the local observables described by (32) will show the same power law damping. This is an important finding as it demonstrates some sense of universality, i.e, various relevant quantities, i.e., density field, velocity field and two-point correlators have the same long time behaviour.

### 3.2 Particles released from a potential

To see the relaxation process more concretely, we consider the time evolution of the particles released from a potential . Suppose particles are trapped by this potential and stay at the ground state. Then the initial profile are given by,

 ϵF=p22+cx2m⇒x0±(p)=±(p20−p22c)12m,p20≡2ϵF. (35)

Here is the Fermi energy which is determined by the constraint (17) as follows

 1= 2∫p0−p0dp2π(p20−p22c)12m=p0√π(p202c)12mΓ(1+12m)Γ(1+12m+12) ⇒p0=⎛⎜ ⎜⎝√π(2c)12mΓ(1+12m+12)Γ(1+12m)⎞⎟ ⎟⎠mm+1. (36)

At , we suddenly turn off this potential (). Then the particles start evolving according to (18). At large , the particle density becomes

 ρ(x,t)= NL+1πℏ(Lp0c)12m(1t)1+12m[ζ(−12m,p0t−xL−⌊p0t−xL⌋) +ζ(−12m,p0t+xL−⌊p0t+xL⌋)]+⋯, (37)

through (34). Therefore, the particle density shows a power law relaxation . Note that this result is consistent with the previous study [27] which investigated the case and observed the damping of the density as .

### 3.3 Particles released from a box

As the second example, we consider the evolution of the particles released from a box (infinite potential walls). Suppose particles are confined in a region by infinite potential walls and stay at the ground state. Then the profile of is given by

 x0±(p)=±x0θ(p+p0)θ(p0−p),p0=π2x0, (38)

where is determined so that it satisfies the constraint (17). At , we remove the potential and release the particles on the circle.

In this case, we cannot apply the assumption (23). However we can compute the density exactly. By substituting (38) to (21), we obtain

 ρ(x,t)= NL+∞∑k≠0∫p0−p0dp2π2ℏk(e2πikL(x0+pt−x)−e2πikL(−x0+pt−x)) = +L2πℏt{ζ(−1,p0t−x+x0L−⌊p0t−x+x0L⌋)−ζ(−1,p0t−x−x0L−⌊p0t−x−x0L⌋)}. (39)

Here we use the relation between the function and the Bernoulli polynomial , and , and we obtain

 ρ(x,t)= NL+L4πℏt{p0t+x+x0L−⌊p0t+x+x0L⌋−(p0t+x+x0L−⌊p0t+x+x0L⌋)2} (40)

We plot this result in Figure 4. Thus shows a power law relaxation . Here this relation may be interpreted as the case in (29).

## 4 Particle motions in the external potential (the V≠0 case)

We investigate the time evolution in the case. Again we put the Fermions on a periodic circle but the arguments in this section can be applied to the relaxation of the Fermions confined in a potential without introducing a circle.

In the phase space, each particles move along the constant slices, where is energy of the individual particles. By regarding this point, we can solve the phase space hydrodynamic equations (13) and (14) as

 u(x,p,t)= θ(t−∫xx0+(E(x,p))dyp(E(x,p),y))θ(∫xx0−(E(x,p))dyp(E(x,p),y)−t), E(x,p)≡p22+V(x),p(E,y)≡√2(E−V(y)). (41)

Here we have not taken into account the periodicity of the motion and will do it later. In this equation, denote the locations of the boundary of the droplet at at a given energy . Therefore

 ∫xx0±(E(x,p))dyp(E(x,p),y),

represents a ‘traveling time’ which a particle at with energy spends for traveling from to (related to the time-of-flight coordinate). Thus equation (41) simply says that the droplet is 1 only if

 ∫xx0+(E(x,p))dyp(E(x,p),y)≤t≤∫xx0−(E(x,p))dyp(E(x,p),y).

This result is physically reasonable and we can also check that equation (41) satisfies (13) and (14).

Note that there are three types of particle motions in this system: (I) right (left) moving forever, (II) trapped by the potential, and (III) the separatrix of these motions. See Figure 5. (In the case, only the case I appears.) We should treat separately by regarding these three motions.

For simplicity, we assume that the potential has only one minimum and one maximum, and we set for the minimum and for the maximum respectively. Also we set and . Then, if , the particle motion is type I, and, if , the particle motion is type II, and corresponds to type III.

Now we consider the periodicity of the motions. In the case I, because of the periodicity of the circle, the particle with energy moves with a periodicity

 TI(E)=∫L0dyp(E,y). (42)

In the case II, the particle shows a periodic motion with a period

 TII(E)=2∫x2(E)x1(E)dyp(E,y), (43)

where denotes the two turning points at which . See Figure 5. On the other hand, in the case III, the particles just approach to by spending a infinite amount of time, and they do not show periodicity. Related to this, as approaches to , the period diverge logarithmically. Indeed, by assuming , we can estimate the periodicity for as

 TI(E)= ∫L0dyp(E,y)∼∫xc+αxc−αdy√2(E−Ec)+ω20(y−xc)2 = 2 ∫α0dz√2(E−Ec)1√1+ω20z22(E−Ec), (44)

where is a some scale and we have used . Then by defining , and, through the integral, we obtain

 TI(E)∼1ω0log(αω0√E−Ec), (45)

for and it diverges logarithmically. We obtain a similar result for in the case II (). This logarithmic divergence will be important when we discuss the power law relaxation 101010The particle motion near the critical point saturates the bound on chaos [50], once we turn on the quantum effect [51]. .

By taking into account these periodicities, the Wigner distribution (41) for the case I becomes

 u(x,p,t)= ∑mθ(t−mTI(E(x,p))−∫xx0+(E(x,p))dyp(E(x,p),y)) ×θ(∫xx0−(E(x,p))dyp(E(x,p),y)−t+mTI(E(x,p))) = 1TI(E(x,p))∫x0+(E(x,p))