Radiation dynamics in homogeneous plasma

Physica D 126 (1999) 236–260

Radiation dynamics in homogeneous plasma M. Escobedoa,1, ∗ , M.A. Herrerob,2 , J.J.L. Velazquezb,2 ´ Vasco, Apartado 644, E-48080 Bilbao, Spain Departamento de Matem´aticas, Universidad del Pais Departamento de Matem´atica Aplicada, Facultad de Matem´aticas, Universidad Complutense, E-28040 Madrid, Spain a

b

Received 5 January 1998; received in revised form 10 July 1998; accepted 14 September 1998 Communicated by U. Frisch

Abstract We study in this paper the asymptotic behaviour of solutions of a nonlinear Fokker–Plank equation. Such an equation describes the evolution of radiation for a gas of photons, which interacts with electrons by means of Compton scattering and Bremsstrahlung radiation. Assuming that a suitable adimensional parameter ε (which measures the strength of the Bremsstrahlung effect) is small enough, we show that the problem considered has two natural timescales. For times t = O(1), the dynamics is conducted by that of a reduced problem, corresponding to setting ε = 0 in the original equations. Solutions of that problem may blow up in finite time, and the total number of photons is no longer preserved after the singularity formation. Nevertheless, solutions of this problem can be continued for all times, if defined in a suitable sense. When t → ∞, solutions of such a modified problem converge towards a Bose–Einstein distribution with a suitable (in general nonzero) chemical potential. However, at times of order t = O((ε|log ε|)−2/3 ), the Bremsstrahlung term becomes dominant at low frequencies, and c 1999 Elsevier Science B.V. All drives the photon distribution to approach to a Planck distribution as time goes to infinity. rights reserved. Keywords: Radiation; Fokker–Plank equation; Bremsstrahlung effect

1. Introduction The purpose of this paper is to describe the dynamics of a low energy, homogeneous, isotropic photon gas, that interacts with a low energy, low temperature electron gas with a Maxwellian distribution of velocities. This is a simplified but relevant problem in the study of nonequilibrium radiation transport in highly ionized gases, and as such it has deserved considerable attention in the study of kinetic theory of plasmas (see Refs. [2,4,7,9] etc.), including in particular the analysis of the interaction between galactic dust and X-ray and cosmic background radiation at the early stages of the Universe [8,10]. The derivation of the relevant kinetic equation can be seen in [4], whereas a ∗

Corresponding author. Partially supported by DGICYT Grant PB96-0663 and EEC Contract ERB FMRX CT960033 HCL. 2 Partially supported by DGICYT Grant PB96-0614. 1

c 0167-2789/99/$ – see front matter 1999 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 2 7 8 9 ( 9 8 ) 0 0 2 5 5 - 3

M. Escobedo et al. / Physica D 126 (1999) 236–260

237

Fokker–Planck approximation in a suitable asymptotic limit, which is more convenient for an analytic treatment of the problem, has been obtained in [2,7,9]. More precisely, assume that (i) photons have low energy compared with the rest energy of the electron (hν/me c2  1 where me is the rest mass of the electron), (ii) temperature T is in the range 1  KB T  me c2 , KB being the Boltzmann constant. Then the corresponding Fokker–Planck equation reads as follows:    ∂f −2 ∂ 4 ∂f 2 =x x +f +f (1.1) + σ (x)(f0 − f ) for 0 < x < +∞, t > 0. ∂t ∂x ∂x In (1.1) the first term on the right is responsible for the approach to equilibrium through Compton scattering, whereas the second one corresponds to emission–absorption by Bremsstrahlung radiation. More precisely, f (x, t) is the distribution radiation, where x is the frequency of photons in units of kB T / h, h is the Planck constant, and the time t is measured in units of (KB T ne /mc)σc . Here ne is the electron density, and α(x) = σc x 4 is the scattering cross-section due to Compton scattering, where σc = (8π/3)(e2 /me c2 )2 , e being the charge of the electron. The second term on the right in (1.1) represents Bremsstrahlung at a rate given by the Born approximation x  , (1.2) σ (x) = εx −3 e−x/2 K0 2 where K0 is the zero-order modified Bessel function, and ε is the inverse of the corresponding relaxation time. In many practical cases it can be considered as a small parameter. Finally, f0 (x) is the Planck distribution of energy, which in our current units reads f0 (x) =

ex

1 −1

for x > 0.

(1.3)

The reader is referred to Refs. [7,10] for additional details concerning the derivation of (1.1). This last is to be considered together with initial data f (x, 0) = g(x) ≥ 0

for 0 ≤ x < +∞,

as well as the natural boundary conditions   ∂f + f + f 2 → 0 when x → 0 and x → +∞. x4 ∂x

(1.4)

(1.5)

Notice that (1.5) states that for the equation obtained by setting ε = 0 in (1.1) (which shall be referred to in the sequel as the reduced equation), there is no flux either at x = 0 or at x = ∞, so that one formally has that the quantity Z +∞ x 2 f (x, t) dx (1.6) M(t) = 0

R∞ should be constant in time, M(t) = M(0) − 0 x 2 g(x) dx in such a case. In physical terms, M(t) represents the total number of photons, and will often be referred to as the total mass for convenience. It is worth noticing that, while f0 (x) in (1.3) is an explicit stationary solution of (1.1), the Bose–Einstein distributions fµ (x) =

1 ex+µ

−1

for x > 0,

(1.7)

where µ > 0, are stationary solutions of the reduced equation. Our concern here is the analysis of problem (1.1)–(1.5) (henceforth termed as problem (Pε )), as well as the relation, as ε → 0, to the reduced problem obtained by setting ε = 0 in (Pε ) (problem (P0 )). These questions have been addressed in [2,3,7,9,10], and more recently in [5]. Roughly speaking, one expects that solutions of (Pε ) will converge for large times to the Planck distribution f0 in (1.3). When the reduced problem (P0 ) is considered,

238

M. Escobedo et al. / Physica D 126 (1999) 236–260

however, the fact that the mass M(t) seems to be preserved has raised the question of whether a transient Bose– Einstein condensation phenomenon might appear [3,9]. More precisely, assume that we start from an initial value g(x) such that Z +∞ Z +∞ x2 2 x g(x) dx > dx = M. M(g) = ex − 1 0 0 Since the only stationary solutions to    ∂f ∂ ∂f = x −2 x4 + f + f2 ∂t ∂x ∂x

for 0 < x < +∞, t > 0,

(1.8)

compatible with (1.5), are the Bose–Einstein and Planck distributions, it is natural to guess that global solutions f (x, t) will concentrate the excess of mass at x = 0 when t → +∞, so that f (x, t) → (ex − 1)−1 + (M(g) − M)x −2 δ(x),

(1.9)

where δ(x) denotes the Dirac mass at the origin. The analysis done in [5] shows, however, that this is not the case. Indeed, it was proved in that paper that solutions of (P0 ) may develop singularities in a finite time T > 0. At the onset of such singularities (which may appear for arbitrarily small values of the initial mass M(g)), the flux condition (1.5) is lost at x = 0, where one merely has the estimate 0 ≤ f (x, T ) ≤

C x2

as x → 0,

(1.10)

for some C > 0. Notice that (1.10) precludes the condensation phenomenon corresponding to (1.9). Nevertheless, if we replace (1.5) at x = 0 by (1.10), it turns out that the corresponding modified problem, that we shall denote by (MP0 ), has a unique solution for all times t > T . In another words, solutions can be continued after a singularity sets in, provided that we replace the flux condition at the origin by estimate (1.10). In general, mass is not conserved anymore, and one has instead that  Z +∞ d 2 x f (x, t) dx ≤ 0 for t > T . (1.11) dt 0 Bearing in mind these facts, our goal here is to discuss the dynamical behaviour of problem (Pε ) in the case where 0 < ε  1. To begin with, and in contrast with the situation for (P0 ), for any fixed ε > 0, solutions of (Pε ) are defined for all times. Since the proof of such result is a bit technical, we shall postpone it to Appendix A at the end of the paper. On the other hand, we shall concern ourselves with the study of solutions of (Pε ) when ε → 0. As it turns out, the limit problem obtained in this case is not (P0 ), but (MP0 ). Actually, the effect of the Bremsstrahlung term prevents the formation of singularities at x = 0, although solutions can become very large near the origin for times t = O(1). This question will be examined in detail in Section 4, where a boundary layer analysis of the case ε → 0 will be presented. Before this, we shall provide in Section 2 a review of the mechanism of formation of singularities (blowup) for problem (P0 ), that has been studied in [5]. Then, in Section 3, we analyse the long-time behaviour of solutions of the modified problem (MP0 ). In particular, it is shown that the radiation distribution f (x, t) will converge, as t → +∞, to either f0 (x) in (1.3), or to some of the Bose–Einstein distributions fµ (x) in (1.7). However, the precise value of the limit profile seems to depend strongly on the actual form of the initial value g (notice that M(t) is not necessarily preserved in time, so that we can only assert that (1.11) holds), and we cannot determine at this stage what asymptotic limit will correspond to any initially given distribution of radiation f (x, 0) = g(x). Finally, the long-time asymptotics of problem (Pε ) when 0 < ε  1 will be considered in Section 5. The presence of the small parameter ε leads then to a two-timescales problem. For times t = O(1), the problem is dominated by (MP0 ), whose asymptotics is discussed in Section 3. Then, in a timescale t = O((ε|log ε|)−2/3 ), solutions equilibrate

M. Escobedo et al. / Physica D 126 (1999) 236–260

239

towards the Planck distribution f0 , a fact that we examine in detail for a particular class of initial values g. Notice that such a timescale is much faster than that corresponding to t = O(1/ε), that could be obtained from a naive look at (1.1) (cf. to this respect the observations before (5.2) in Section 5). A few concluding remarks are then gathered in a short Section 6. The paper ends with Appendix A where a global existence result for problem (Pε ) is provided. 2. Singularity formation for ε = 0: a review Consider the problem    ∂ ∂f ∂f = x −2 x4 + f + f2 , ∂t ∂x ∂x

when 0 < x < +∞, t > 0,

(2.1)

together with the following conditions: f (x, 0) = g(x) ≥ 0 for x > 0,   4 ∂f 2 +f +f → 0 as x → 0, x ∂x   ∂f + f + f 2 → 0 as x → +∞, x4 ∂x Z +∞ x 2 g(x) dx < +∞,

(2.2) (2.3a) (2.3b) (2.3c)

0

where g is a continuous and non-negative function. We have proved in [5] that, for a large class of initial values g(x), solutions of (2.1),(2.2),(2.3a)–(2.3c) will develop singularities near x = 0 in a finite time. In particular, it was shown in that paper that condition (2.3a) fails to hold when the singularity develops. In this section we shall briefly recall the blow-up mechanism analysed in [5], since it will be relevant for later purposes. We first remark that the change of variables v(y, t) = xf (x, t),

x = ey ,

(2.4)

transforms (2.1) into vt = vyy + 2vvy + 2v 2 + vy − 2v + ey (vy + 3v),

−∞ < y < +∞, t > 0.

(2.5)

Let T > 0 be the time of formation of singularities for a given solution to (2.1),(2.2),(2.3a)–(2.3c). We introduce new variables as follows: v(y, t) = (T − t)−1 Φ(y − log(T − t), −log(T − t)),

(2.6a)

ξ = y − log(T − t),

(2.6b)

τ = −log(T − t).

Then Φ satisfies the equation Φτ = e−τ Φξ ξ + (2Φ − 1)(Φξ + Φ) + e−τ (Φξ − 2Φ) + e−2τ eξ (Φξ + 3Φ).

(2.7)

The solutions constructed in [5] are such that, as τ → +∞, they converge to solutions of (Φ ξ + Φ)(2Φ − 1) = 0, which consists on arrays of exponential functions Φ = Ce−ξ having jumps in a finite number of shocks, where Rankine–Hugoniot and entropy conditions hold (see Fig. 1) .

240

M. Escobedo et al. / Physica D 126 (1999) 236–260

Fig. 1. A limit blow-up profile of (2.1),(2.2),(2.3a)–(2.3c) in rescaled variables.

More precisely, the following condition must be satisfied at any shock: Φ(s + ) + Φ(s − ) = Ce−s , where ξ = s is the position of the shock, and, as usual, s + = limr→s,r>s r, s − = limr→s,r
C , x2

as x → 0+ ,

(2.8)

for some constant C > 0. In particular, (2.8) implies that (2.3a) cannot hold at t = T . As a matter of fact, once f (x, t) becomes singular near x = 0 at some time, it remains so for all later times. This follows from the fact that, for any real constant ξ0 , the function v(y, t) = (e−(y+αt−ξ0 ) − C)+ , where (s)+ = max{s, 0}, is a subsolution of (2.5) provided that α > 2 and 0 < C < (α − 2)/2. In view of (2.8), we also have that v(y, 0) ≤ v(y, T )

for − ∞ < y < +∞,

provided that ξ0 < 0 and |ξ0 | be large enough. It then turns out by comparison that lim x 2 f (x, t) ≥ eξ0 −αt

x→0

for any t > T .

It is worth noticing that solutions can be continued after a singularity sets in. To this end, we have to replace (2.3a) by a suitable condition at x = 0. More precisely, the following result has been shown in [5]. Assume that 0 ≤ g(x) ≤

C x2

for 0 < x < 1 and some C > 0,

lim x 4 g(x) = 0.

(2.9) (2.10)

x→+∞

Then there exists a unique solution f (x, t) of the problem given by (2.1),(2.2),(2.3b),(2.3c),(2.9) and (2.10). Moreover, such solution satisfies for any t > 0 the following estimates: 0 ≤ f (x, t) ≤

C x2

for 0 < x < 1 and some C > 0 depending on g,

(2.11a)

M. Escobedo et al. / Physica D 126 (1999) 236–260

and

 lim x 4

x→+∞

 ∂f (x, t) + f (x, t) + f 2 (x, t) = 0. ∂x

241

(2.11b)

3. Long-time asymptotics for the problem without Bremsstrahlung From now on, we shall define the modified problem (MP0 ) as that consisting of (2.1),(2.2),(2.3a)–(2.3c) for times t before the time T of singularity formation (if any), and by (2.1),(2.2),(2.3b),(2.9) and (2.10) for times t > T . The goal of this section consists in proving that any solution of (MP0 ) converges, as t → ∞, towards one of the Bose–Einstein distributions given in (1.7). More precisely, the following result holds: Theorem 1. Let f be the solution of (MP0 ), and assume that g does not vanish identically. Then there exists µ ≥ 0 such that lim f (x, t) = fµ (x) =

t→+∞

and

Z

+∞

Z x f (x, t) dx → 2

0

0

1 , ex+µ − 1

+∞

unif ormly on compact sets of (0, +∞),

x 2 fµ (x) dx,

as t → +∞.

(3.1)

(3.2)

Moreover, if (2.3a) holds for all t > 0, then the solution f (·, t) converges in the same sense described above towards the unique stationary solution fµ such that Z +∞ Z +∞ 2 x g(x) dx = x 2 fµ (x) dx. (3.3) 0

0

Before starting with the proof of Theorem 1, a few remarks are in order. To begin with, in general we do not know how to characterize the parameter µ in the limit distribution. For instance, even when µ ≥ 0 exists such that (3.3) holds, we cannot conclude that such a condition determines the limit profile fµ . The reason is that singularities may R +∞ appear in finite time for which (2.3a) no longer holds, and in this case the quantity M(t) = 0 x 2 f (x, t) dx need not be constant in time anymore. If, however, (2.3a) holds for any time, then the asymptotic profile is determined by (3.3). There are a number of cases when the flux condition (2.3a) remains true for all times. For instance, this happens when g is such arguments for that 0 < g(x) < f0 (x). This can be seen by means of classical regularity R +∞ theory and Rcomparison +∞ parabolic equations. Incidentally, when the initial value g is such that 0 x 2 g(x) dx > 0 x 2 f0 (x) dx, it follows that singularities necessarily develop in finite time. However, as was observed in [5], even solutions for which M(0) is arbitrarily small may blow up in finite time. Once the flux condition is lost, M(t) may decrease with time, but we always have that M(t) ≥ C > 0, for some positive constant C depending on g. This can be seen by comparing f (x, t) with f (x, t), solution of the same problem with a compactly supported initial value g such that 0 ≤ g ≤ g and 0 ≤ g ≤ f0 . Proof of Theorem 1. For convenience of the reader, the proof will be split into a number of lemmas. Let us begin by introducing a definition. We say that two continuous functions f1 (x, t) and f2 (x, t) intersect at time t, if there are two open intervals in the real line, I1 and I2 , such that f1 (·, t) > f2 (·, t) in I1 and f1 (·, t) < f2 (·, t) in I2 . The following result then holds: Lemma 2. Let f1 (x, t) and f2 (x, t) be two classical solutions of (2.1) and (2.3b) corresponding to two initial values g1 , g2 that satisfy (2.2),(2.9) and (2.10). If f1 and f2 intersect in a time interval t1 < t < t2 , then the quantity

242

M. Escobedo et al. / Physica D 126 (1999) 236–260

Z φ(t) =

+∞

(f1 − f2 )+ (x, t)x 2 dx

(3.4)

0

is strictly decreasing for t1 < t < t2 . Proof. Let us denote by 0(t) = {x : f1 (x, t) = f2 (x, t)}. Set sgn+ (s) = 1 for s > 0, sgn+ (s) = 0 otherwise. If we write the equation satisfied by f1 − f2 , multiply both sides by sgn+ (f1 − f2 ) and then integrate over the whole line, we eventually obtain   Z d +∞ + 2 4 ∂ + ((f1 − f2 ) ) (f1 − f2 ) (x, t)x dx ≤ x ≡ I ≤ 0, (3.5) dt 0 ∂x |0(t) where in any of the intervals (ai , bi ) which make up 0(t), [F (x)]|0(t) = limx→bi ,xai F (x). Let now τ1 , τ2 , be such that t1 < τ1 < τ2 < t2 , and suppose that φ(τ1 ) = φ(τ2 ). We should then have that Z τ2 Z τ2 0 φ (s) ds ≤ [x 4 ((f1 − f2 )+ )x ]|0(t) dt. 0 = φ(τ2 ) − φ(τ1 ) = τ1

τ1

Since f1 and f2 intersect for τ ∈ (τ1 , τ2 ), by the strong maximum principle (cf. [6], p. 52) we deduce at once that [((f1 − f2 )+ )x ]|0(s) < 0 for τ < s < τ + δ and for some δ > 0. Therefore I < 0 and the result follows.  We next derive the following ODE result. Lemma 3. For any α > 0 there exists a unique positive function f α which satisfies 2

x 4 (f x + f + f ) = α,

f or 0 < x < +∞.

(3.6)

0

Moreover, f α < 0 for all x > 0, and √ α α f α (x) ∼ 2 as x → 0, f α (x) ∼ 4 as x → +∞. x x

(3.7)

On the other hand |f α (α 1/4 y) − h(y)| = o(h(y))

as α → +∞, unif ormly on 0 < y < +∞,

(3.8)

where h is the unique positive solution of 2

h+h =

1 . y4

Finally, f α < f β if α < β, and f α (x) → f0 (x) =

1 ex − 1

as α → 0, unif ormly on compact subsets of (0, +∞).

(3.9)

Proof. The existence and uniqueness of a positive solution of (3.6) follows from standard ODE techniques. As a matter of fact, local non-negative solutions of Eq. (3.6) are approximately described in Fig. 2. The asymptotic behaviours given in (3.7) are easily obtained by means of classical dominated balance techniques (cf. for instance [1]). To check (3.8), we set y = xα −1/4 , hα (y) = f α (x), and observe that hα satisfies y 4 (α −1/4 (hα )y + hα + h2α ) = 1. Finally, to derive (3.9) we merely observe that f0 satisfies (3.6) with α = 0 there. The corresponding convergence statement then follows from standard continuous dependence results for ODEs. 

M. Escobedo et al. / Physica D 126 (1999) 236–260

243

Fig. 2. Non-negative trajectories of (3.6) defined for x > 0.

Let f (x, t) be the solution of (MP0 ) for a given initial data g(x). In view of (2.11a) and (2.11b) it turns out that 0 ≤ f (x, t) ≤ f α (x)

for all x > 0, t > 0, provided that α is selected large enough.

(3.10)

Inequality (3.10) enables us to derive local bounds on f and its derivatives by means of classical parabolic estimates (cf. [6]). It then turns out that for any sequence (tj ) with tj → +∞, there exists a subsequence (still denoted by (tj )), and a function f ∗ (x) (possibly depending on that subsequence), such that lim f (x, tj ) = f ∗ (x) uniformly on compact subsets of (0, +∞).

j →+∞

(3.11)

By (3.10), we certainly have that f ∗ (x) ≤ f α (x),

for all x > 0.

(3.12)

Moreover, the following result holds: Lemma 4. The function f ∗ defined above is such that f ∗ (x) ≤ f0 (x) =

1 ex − 1

f or all x > 0.

(3.13)

Proof. We claim that the following alternative takes place. Either: (i) f ∗ (x) ≤ f0 (x) for all x > 0, or else: (ii) f ∗ (x) intersects f α (x) for any α > 0 sufficiently small. If (ii) does not happen, then there exists a sequence (αj ), αj → 0, such that f ∗ does not intersect f αj . Thus, by (3.7), f ∗ ≤ f αj and therefore f ∗ ≤ limj →+∞ f αj ≡ f0 . Let us show next that case (ii) never occurs. To this end we argue by contradiction, and assume that f ∗ intersects f α for all α ≤ α0 . By (3.11) we then have that for any α such that 0 < α < α0 /2 there exists t0 = t0 (α) with the propertyRthat f (·, t) intersects f α for t > t0 . Arguing as in the proof of Lemma 2, we then derive that the function +∞ ψ(t) = 0 x 2 (f (·, t) − f α )+ dx is strictly decreasing for t > t0 . Moreover, limt→+∞ ψ(t) = 0. If this were not the case, i.e., if limt→+∞ ψ(t) = a > 0, we would have that Z +∞ x 2 (f ∗ (x) − f α (x))+ dx = a. (3.14) 0

244

M. Escobedo et al. / Physica D 126 (1999) 236–260

Let now F (x, t) be the solution of (P0 ) with initial value F (x, 0) = f ∗ (x). By Lemma 2 we should have that Z +∞ x 2 (F (x, 1) − f α (x))+ dx < a. 0

On the other hand Z Z +∞ x 2 (f (x, tj + 1) − F (x, 1))+ dx ≤ 0

0

+∞

x 2 (f (x, tj ) − f ∗ (x))+ dx

and this last quantity goes to zero as j → +∞. Hence, ψ(tj + 1) < a for large j , which contradicts (3.14) since ψ(t) is strictly decreasing. Thus a = 0, and by (3.14) this in turn implies that f ∗ ≤ f α . This concludes the proof. As a next step, we show: Lemma 5. The following alternative holds. Either: (i) f ∗ = fµ for some µ ≥ 0, or else: (ii) f ∗ intersects fµ for some µ > 0. Proof. Suppose that f ∗ does not intersect any of the fµ with µ > 0. Then either there exist µ0 > 0 for which f ∗ ≥ fµ0 , or otherwise f ∗ ≤ fµ for any µ ≥ 0. The second possiblity would imply that f ∗ = 0, which contradicts R +∞ the fact (already remarked at the beginning of this section) that 0 x 2 f (x, t) dx ≥ C(g) > 0 for all t > 0. We now define µ∗ = inf{µ > 0, f ∗ (x) ≥ fµ (x)}. Clearly, µ0 ≥ µ∗ ≥ 0. Take any sequence (δj ) with δj > 0 such that δj → 0 as j → +∞. By our choice of µ∗ , for any j there exists an interval of the positive half-line where fµ∗ ≤ f ∗ < fµ∗ −δj . On the other hand, since f ∗ does not intersect fµ∗ −δj for any j , it follows that f ∗ < fµ∗ −δj for all x > 0. Letting now j → +∞, we obtain  that f ∗ = fµ∗ , and the proof is concluded. As a last step in the argument, we now show that case (ii) in Lemma 5 never occurs. Lemma 6. The function f ∗ does not intersect any distribution fµ with µ > 0. Proof. We first point out that, if f ∗ intersects fµ for some µ > 0, then it also intersects fµ+δ for some δ > 0 small enough. In that case, there exists t0 ≡ t0 (µ, δ) such that f (·, t) intersects fµ+δ for all t ≥ t0 . We now may argue as R +∞ in Lemma 4 to show that the functional ζ (t) = 0 x 2 (f (x, t) − fµ+δ (x))+ dx is such that ζ (t) → 0 as t → +∞,  whence f ∗ = fµ+δ , which contradicts the assumption that f ∗ intersects fµ . End of the Proof of Theorem 1. Putting together Lemmas 2–6, we have already shown that for any sequence (tj ) → +∞ as j → +∞, there exists a subsequence (still denoted by (tj )) and a number µ ≥ 0 such that lim f (x, tj ) = fµ (x)

j →+∞

uniformly on compact sets of (0, +∞).

(3.15)

It merely remains to show that the value of µ in (3.15) is independent of the sequence (tj ). To this end, we observe R +∞ that by Lemma 2 (with f2 = 0), function M(t) = 0 x 2 f (x, t)dx is decreasing in time. Assume that there exist two different sequences of time (tn1 ) and (tn2 ) with limn→+∞ tn1 = limn→+∞ tn2 = +∞, and two non-negative numbers µ1 , µ2 with µ1 6= µ2 and such that f (·, tn1 ) → fµ1 ,

f (·, tn2 ) → fµ2 ,

as n → +∞.

Then, letting n → +∞, we would have that Z i fi∗ (x)x 2 dx ≡ M(fµi ) for i = 1, 2. M(tn ) → R+

Since µ1 6= µ2 , this would contradict the monotonicity of M(t). The proof is now concluded.

M. Escobedo et al. / Physica D 126 (1999) 236–260

245

4. The complete problem. Influence of Bremsstrahlung In this and the following sections we shall consider the complete problem (P) consisting of (1.1)–(1.5). Our basic startpoint is the fact that there exists a unique solution of (P) which exists for all times t > 0. To keep the flow of the main arguments here, the proof of that fact will be postponed to Appendix A at the end of the paper. Since we already know that blow-up may occur at the origin in absence of Bremsstrahlung, it is apparent that such term should play a crucial role near x = 0. In this section we shall describe the smoothing effect of the term σ (x)(f0 − f ) near the origin, for values 0 < ε  1, by means of matched asymptotic expansions methods. To proceed, we assume that, for some given initial value g(x), the solution of the reduced problem (P0 ) develops a singularity at x = 0, t = T > 0, according to the mechanism recalled in Section 2. Keeping to the notation introduced in that section, v(y, t) will become large for y → −∞ and t ∼ T . In such regions, we can approximate f0 to the lower order as follows: f0 (x) ∼

1 = e−y . x

(4.1)

On neglecting smaller terms, in the region where v  1 we may then replace (2.5) by the equation vt = vyy + 2vvy + 2v 2 − εe−3y y(1 − v),

as y → −∞.

(4.2)

In terms of the variables Φ, ξ, τ introduced in (2.6a),(2.6b) and (4.2) reads Φτ ∼ e−τ Φξ ξ + (2Φ − 1)(Φξ + Φ) − εeτ e−3ξ (ξ − τ )(1 − eτ Φ) as ξ + log(T − t) → −∞.

(4.3)

Let us assume that, when ε = 0, the corresponding solution of (P0 ) approaches asymptotically to a function Φ(ξ ) of the type described in Section 2. For simplicity, we shall restrict ourselves to the case where Φ has a single shock at ξ = 0. One then has that Φ(ξ ) = e−ξ

for ξ > 0,

Φ(ξ ) = 0

for ξ < 0.

In Eq. (4.3), one readily sees that the term eτ Φ  1 for large times when Φ = O(1). On the other hand, the contribution from e−τ Φξ ξ is expected to be negligible as τ → +∞, everywhere except near the discontinuity of the limit profile Φ. Far from the limit shock, we may then approximate (4.3) as follows: Φτ = (2Φ − 1)(Φξ + Φ) − εeτ e−3ξ (ξ − τ )(1 − eτ Φ) as ξ + log(T − t) → −∞.

(4.4)

Near the shock ξ = λ(τ ), we impose entropy and Rankine–Hugoniot conditions, namely Φ(λ(τ )− , τ ) < Φ(λ(τ )+ , τ ),

(4.5)

and −λ0 (τ ) =

(Φ 2 − Φ)(λ(τ )+ ) − (Φ 2 − Φ)(λ(τ )− ) . Φ(λ(τ )+ ) − Φ(λ(τ )− )

(4.6)

When ε > 0 is small, the effect of the last term on the right in (4.4) becomes important as soon as ετ e2τ ∼ 1. This suggests introducing the new timescale s = τ + 21 log ε + 21 log(|log ε|),

(4.7)

which transforms (4.4) into Φs = (2Φ − 1)(Φξ + Φ) −

  1 e2s e−3ξ − 21 log ε − 21 log(|log ε|) + s − ξ Φ. |log ε|

(4.8)

246

M. Escobedo et al. / Physica D 126 (1999) 236–260

For 0 < ε  1, (4.8) is in turn dominated by the reduced equation Φs = (2Φ − 1)(Φξ + Φ) − 21 e2s e−3ξ Φ.

(4.9)

Eq. (4.9) describes a transition layer between an aysmptotic behaviour Φ for ε = 0 and a dynamics dominated by the absorption rate σ (x)(f0 − f ) in Eq. (1.1). Such an equation has to be complemented with the matching conditions Φ(ξ, s) → Φ(ξ ) away from ξ = 0, as s → −∞.

(4.10)

Problem (4.9)–(4.10) can be explicitly solved. To do that, we define new variables G, u, θ given by Φ = e−s G(ξ − s, s),

u = ξ − s,

θ = −e−s .

Eq. (4.9) is then recast in the form   ∂G 1 ∂G = 2G + G − e−3u G. ∂λ ∂u 2

(4.11)

The equations for the characteristics of (4.11) are du = −2G, dθ

dG 1 = 2G2 − e−3u G, dθ 2

which give 1 dG = −G + e−3u . du 4 Integrating this last equation yields G = − 18 e−3u + Ce−u , whence Φ(ξ, s) = − 18 e−3ξ e2s + Ce−ξ ,

(4.12)

where C is an arbitrary integration constant. Actually, (4.12) has to be truncated by zero at the position of the shock ξ = λ(τ ), where condition (4.6) should hold. The actual function Φ is therefore given by (4.12) for ξ > λ(s), and Φ = 0 for ξ < λ(s). Taking into account our matching condition (4.10) and our choice of Φ(ξ ), it turns out that C = 1 in (4.2), and lim λ(s) = 0.

(4.13)

s→−∞

The Rankine–Hugoniot condition (4.6) gives now the following ODE that has to be satisfied by λ(s) λ0 = 1 + 18 e−3λ+2s − e−λ .

(4.14)

Eq. (4.14) can be analysed by means of classical ODE techniques. As a matter of fact, trajectories are approximately described in Fig. 3. It can be shown that (4.13) and (4.14) admits a monoparametric family of solutions satisfying   1 K 2 2s (4.15a) − e , as s → −∞, λ(s) ∼ Kes + 8 2 λ(s) ∼ s + A(K),

as s → +∞,

(4.15b)

M. Escobedo et al. / Physica D 126 (1999) 236–260

247

Fig. 3. Trajectories of Eq. (4.14).

where A(K) is a uniquely determined function of parameter K in (4.15a). We notice that (4.12),(4.15a) and (4.15b) provide the following asymptotics for Φ: Φ(ξ, s) ∼ e−s Q(ξ − s)

as s → +∞,

(4.16a)

Q(z) = − 18 e−3x + e−z for z > A(K), Q(z) = 0

for z < A(K).

(4.16b)

Observe that the shock amplitude vanishes exponentially as s → +∞. It remains yet to select the value of K in (4.15a) and (4.15b). To this end, we derive an expansion in terms of the parameter ε of the solution of Eq. (4.3) in the region where ετ e2τ  1. Namely, we try Φ = Φ0 + εΦ1 + · · · , and assume that the initial data is being perturbed at most by terms of order ε as ε → 0. Function Φ0 will then solve Eq. (2.7), whereas Φ1 should satisfy Φ1,τ = e−τ Φ1,ξ ξ + (2Φ0 − 1)(Φ1,ξ + Φ1 ) + 2Φ1 (Φ0,ξ + Φ0 ) − eτ e−3ξ (ξ − τ )(1 − eτ Φ0 ), where the initial value for Φ1 would include any possible correction on the initial value of Φ0 up to order O(ε). When τ  1 and ξ is far away from the shock, the relevant terms in that equation are Φ1,τ = (2Φ0 − 1)(Φ1,ξ + Φ1 ) − τ e2τ e−3ξ Φ0 . In the region under consideration, one has that Φ1 ∼ τ e2τ W (ξ ), where W satisfies 2W = (2Φ0 − 1)(W 0 + W ) − e−3ξ Φ0 ∼ (2e−ξ − 1)(W 0 + W ) − e−4ξ . Integrating this last equation, we obtain that C0 e−3ξ 1 , W (ξ ) = − e−3ξ + 4 (2e−ξ − 1)2

248

M. Escobedo et al. / Physica D 126 (1999) 236–260

where in order to avoid singularities we have to set C0 = 0 in the expression above. This yields the following expression for Φ1 : Φ1 ∼ − 41 τ e2τ e−3ξ , when τ  1. Concerning Eq. (4.6), we now have that λ0 (τ ) = λ(τ ) − εΦ1 (0, τ ) ∼ λ + 41 ετ e2τ = λ + 18 e2s , where we have made use of (4.7) to derive the last formula on the right above. Bearing on mind (4.15a) and (4.15b), we readily see that matching is now achieved by setting K = 0 there. To conclude this section, we just rewrite the previous results in terms of the variables v and y. For simplicity, we shall set δε = − 21 (log ε + log(|log ε|)). Taking into account (4.16a),(4.16b),(4.7) and (2.7), we obtain 1 Q(y + δε ), v(y, T ) ∼ p ε|logε|

for 1  −y ≤ δε + C, with y < 0 and C arbitrary.

(4.17)

Estimate (4.17) corresponds to a boundary layer near x = 0, where Bremsstrahlung is dominant. Notice that (4.17) gives v(y, T ) ∼ e−y for |y|  1, which coincides with the profile obtained in [5] for the case ε = 0. Notice also that the final profile (4.17) exhibits a shock at y = −δε + A(0), where A(K) is given in (4.15a) and (4.15b). The intensity of the shock is of order (ε|log ε|)−1/2 . As a matter of fact, this discontinuity will be smoothed out by a Burgers-like diffusive boundary layer of width (ε|log ε|)1/2 , that will not be described in detail here (cf. [5] for a similar analysis). For completeness, we now write (4.17) in the original variables f and x   p 1 ε|logε| 1 + 1 , for ε|logε|eA(0) ≤ x  1, f (x, T ) ∼ 2 − 2 8 x x   p 1 (4.18) , for 0 < x ≤ ε|logε|eA(0) . f (x, T ) ∼ o 2 x p It is worth noticing that the discontinuity in (4.17) and (4.18) disappears for times t − T ∼ ε|logε|. This can be seen by means of the following change of variables: 1 w, v=p ε|logε|

y = −δε + ξ,

t =T +

p ε|logε|ρ.

(4.19)

To the lowest order, Eq. (4.2) now becomes wρ = 2wwξ + 2w 2 − e−3ξ w,

(4.20a)

which is to be supplemented with the initial value w(ξ, 0) = Q(ξ ),

(4.20b)

together with entropy and Rankine–Hugoniot conditions near the shock, that are needed since diffusive terms have been dropped. Since f (ξ ) = −e−3ξ /2 + e−ξ solves (4.20a), it turns out that the solution of (4.20a) and (4.20b) is explicit, namely w(ξ, ρ) = − 18 e−3ξ + e−ξ ,

for ξ ≥ λ(ρ),

w(ξ, ρ) = 0

for

ξ ≤ λ(ρ).

(4.21)

The modified shock λ(ρ) moves now according to the equation λ0 = 18 e−3λ − e−λ ,

λ(0) = A(0).

(4.22)

M. Escobedo et al. / Physica D 126 (1999) 236–260

249

p Fig. 4. Early stage of the process ( ε|log ε|  T − t  1) corresponding to dominance of the second term on the right in (4.16b) when z > A(0) there. This situation closely resembles the case ε = 0.

Fig. 5. Regularization of the profile in the hyperbolic equation (4.9) after Bremsstrahlung becomes dominant, i.e. at times 0 < T − t  (cf. (4.17)). The jump observed in the picture will actually be smoothed out by the effect of the neglected diffusive effects.

p

ε|log ε|

Observe that λ(ρ) approaches exponentially fast towards the constant −(log8)/2 as ρ  1. A sketch of the solution profiles at different stages (before and after Bremsstrahlung takes over) is provided in Figs. 4–6. We finally remark that it is possible to describe the subsequent evolution of f (x, t) at times t ∗ = T + O(1), provided that some information on the asymptotics of f (x, t) near x = 0 is available. As a matter of fact, estimate (2.8) suggests that solutions of (MP0 ) will behave in the form f (x, t) ∼

a(t) , x2

as x → 0, for t ≥ T ,

(4.23)

where a(t) is a continuous function, except perhaps at a discrete set of points (tj ) where jump discontinuities of the form a(tj+ ) > a(tj− ) may appear. Actually, Eq. (4.2) with ε = 0 admits travelling wave solutions of the form v(y, t) = e−(y−t) , which in the original (x, t) variables yields the asymptotics f (x, t) ∼

et x2

as x → 0.

250

M. Escobedo et al. / Physica D 126 (1999) 236–260

Fig. 6. Radiation profile for times 1  t − T  smoothed by diffusion.

p

ε|log ε| (cf. (4.21) and (4.22)). Again, the corner appearing at the left of the picture will be

Fig. 7. A limit profile corresponding to a jump discontinuity in (4.23).

The fact that jump discontinuities as described above may occur can be illustrated by considering solutions as those remarked upon in Section 2, where the limit profile is given by a function as that in Fig. 7. (Compare with Fig. 4). At the shock position of Φ(ξ ) in Fig. 7, Rankine–Hugoniot condition is to be satisfied. It is readily seen that, at the time t = T where a singularity is formed, function a(t) in (4.23) jumps from C1 to C2 . It is important to remark that a(t) is related to the loss of mass of function f (x, t) as time evolves. Indeed, from (2.1),(2.3b), and (4.25) we deduce that Z +∞  d 2 x f (x, t)dx = −(a(t))2 , dt 0 which is strictly negative if a(t) 6= 0. We finally observe that the singular behaviour of f (x, t) as x → 0 described in (4.23) is smoothed out by the effect of Bremsstrahlung radiation as ε → 0. Indeed, suppose that at any time t ∗ > T where a(t) is continuous we make use of the change of variables (4.19) with T replaced by t ∗ there. The function w thus defined then satisfies (4.20a), and is such that w(ξ, ρ) → a(t ∗ )e−ξ

as ξ → +∞.

M. Escobedo et al. / Physica D 126 (1999) 236–260

251

A simple analysis by characteristics of (4.20a),(4.20b),(4.22) shows that, as ρ → +∞, w behaves as the stationary solution  + (4.24) ws (ξ ) = − 18 e−3ξ + a(t ∗ )e−ξ . The asymptotics of f near the origin is then described by function ws in (4.24). The shape of the corresponding v is similar to that in Fig. 6, where the only difference is the change in time of the term a(t ∗ ) in (4.24), which produces modifications in the tail of the figure for y  1. 5. Long time asymptotics for the problem with 0 < ε  1 The purpose of this section is to describe the asymptotics of solutions of problem (P) when t → +∞ and ε > 0 is sufficiently small. On physical grounds, it is reasonable to expect that solutions will approach to the Planck function f0 (x) = (ex − 1)−1 , which corresponds to an equilibrium distribution for photons. On the other hand, we have shown in Section 3 that, when the Bremsstrahlung term is neglected, solutions will converge (away from the origin) to a distribution fµ (x) = (ex+µ − 1)−1 , for some chemical potential µ ≥ 0 whose precise value depends on the details of the initial distribution g(x). Therefore, a crucial question consists in determining the asymptotic behaviour of a solution f (x, t) corresponding to an initial value f (x, 0) = g(x) ≡ (ex+µ − 1)−1 with µ > 0. We shall address this issue here under the assumption that 0 < ε  1, and will show that f (x, t) will slowly move towards f0 as t → +∞. Namely, we have that f (x, t) ∼ fµ(t) (x),

where µ(t) → 0 and µ(t) = 0 at a time t = t0 + O((ε|log ε|)−2/3 ).

(5.1)

Concerning (5.1) a few remarks are in order. To begin with, it is stated there that µ(t) = 0 at some finite time t. As a matter of fact, (5.1) will be valid only until µ(t) becomes small enough. Then a slower motion will set in, which is not described in detail here. This fact will become apparent in the discussion to follow, and will be stressed when appropriate (cf. the remark following (5.8a) and (5.8b)). A second observation is that (5.1) provides a stabilization rate towards Planck distribution which is faster than it might appear at first glance. Indeed, since f0 is a solution of (1.1), subtracting the equations satisfied by f (x, t) and f0 (x), multiplying both sides by x 2 sgn(f − f0 ) and integrating on space yields at once the estimate: Z Z +∞ d +∞ 2 x |f (x, t) − f0 (x)| dx ≤ −ε σ (x)|f (x, 0) − f0 (x)| dx (5.2) dt 0 0 (cf. (A.5) in Appendix A). Note that (5.2) provides stabilization to f0 (x) in times at most of order O(1/ε). To proceed further, it will be convenient to analyse in detail the evolution in time of the quantity M(t) = R +∞ x 2 f (x, t)dx instead of µ(t). The corresponding estimates for this last function can be readily recovered 0 R +∞ through the relation M(t) = 0 x 2 (ex+µ(t) − 1)−1 dx. To begin with, we note that a direct integration of (1.1)– (1.4) gives Z Z ∞   d ∞ 2 x 1 −x/2 dM = e x fε (x, t) dx = ε K0 (f0 (x) − f (x, t))dx. (5.3) dt dt 0 2 x 0 If we now try (5.1) into (5.2), we see that the right-hand side will diverge. This strongly suggests that (5.1) cannot be uniformly valid near x = 0, where a boundary layer should appear. In order to examine that region, we will use the v and y variables (cf. (2.4)) and approximate (1.1) by (4.2). The Bremsstrahlung term becomes relevant when εe−3y |y| = O(1), which suggests introducing a new space variable as follows: ξ = y − 13 log ε − 13 log(|logε|).

(5.4)

252

M. Escobedo et al. / Physica D 126 (1999) 236–260

Eq. (4.2) transforms then into vt ∼ vξ ξ + 2vvξ + 2v 2 + 13 e−3ξ (1 − v),

(5.5)

which has to be complemented with matching conditions: v→0

as ξ → +∞,

(5.6a)

v→1

as ξ → −∞,

(5.6b)

Indeed, (5.6a) follows from (5.1), whereas (5.6b) comes from assuming fast approximation of f (x, t) to f0 (x) for small x (i.e., at low frequencies of the spectrum). Actually, a pointwise estimate of the difference |f − f0 | which is in agreement with this assumption will be given in (A.33) in Appendix A. On taking now δ > 0 arbitrary, we may approximate dM/dt in (5.3) as follows:   Z δ Z log δ |log x| 1 dM ∼ε − f (x, t) dx + O(ε) = ε |y|e−y (1 − v)dy + O(ε) dt x x 0 −∞ Z +∞ 1 e−z (1 − v)dz + O(ε). (5.7) = ε2/3 |logε|2/3 3 −∞ In view of (5.7) and (5.1) will follow provided that we show that Z +∞ e−z (1 − v(z)) dz = O(1) for all t > 0, I (t) ≡ −∞

I (t) > 0

remains strictly positive as t → +∞.

(5.8a) (5.8b)

If conditions (5.8a) and (5.8b) hold, (5.7) can be recast in the form dM/dt ∼ C(t)ε2/3 |logε|2/3 for some C(t) ≥ θ > 0 and some positive constant θ. This would imply that M(t) → +∞ as t → +∞, which is obviously impossible by (5.2). This is why (5.1) has to be understood as t → +∞, which is obviously impossible by (5.2). This is why (5.1) has to be understood as providing convergence to f0 (x) up to times t ∗ > 0 for which Z +∞ x 2 (ex − 1) dx, M(t ∗ ) ∼ M0 = 0

where a new boundary layer should be introduced. To conclude our result, we have yet to show that (5.8a) and (5.8b) hold. This we shall do by means of a comparison argument as follows. Since f (x, 0) = (ex+µ − 1)−1 < (ex − 1)−1 , and this last is an explicit solution of (1.1), it turns out that f (x, t) ≤ (ex − 1)−1 for all t > 0. In terms of v and ξ , this reads v(ξ, t) ≤ 1

for − ∞ < ξ < +∞, t > 0.

(5.9)

Set now w = v − 1. We shall look for stationary solutions w = w(ξ ) of the form w = σ Φ, where 0 < σ  1 and Φ = O(1). Then Φ will essentially satisfy the stationary equation obtained from (5.5) by setting w = v − 1 there, when only linear terms are retained. Namely, the equation for Φ(ξ ) will be Φ 00 + 3Φ 0 + 2Φ − 13 e−3ξ Φ = 0.

(5.10)

We now claim that (5.10) has solutions Φ(ξ ) such that Φ(ξ ) < 0 for all ξ , and Φ(ξ ) → 0 as ξ → ±∞. Actually, the fact that Φ(ξ ) → 0 algebraically as ξ → +∞ follows since the dominant terms in (5.10) when ξ  1 are Φ 00 + 3Φ 0 + 2Φ = 0. A similar argument shows that Φ(ξ ) → 0 at an exponential rate when ξ → −∞. To check that Φ < 0 everywhere it is convenient to change variables as follows: S=

Φ0 , Φ

(5.11)

M. Escobedo et al. / Physica D 126 (1999) 236–260

253

whereby (5.10) transforms into S 0 + S 2 + 3S + 2 − 13 e−3ξ = 0.

(5.12)

We now select S(ξ ) as the solution of (5.12) such that S(0) = 0. A standard ODE analysis reveals that S(ξ ) > 0 for ξ < 0, and S(ξ ) < 0 for ξ > 0. In view of the matching conditions (5.6a) and (5.6b) to be imposed on v, it follows that we may select σ small enough so that v(ξ, t) ≤ 1 + σ Φ(ξ ) ≤ 1,

for − ∞ < ξ < +∞, t > 0,

(5.13)

where Φ(ξ ) is obtained from Φ(ξ ) through relation (5.11). As a next step, we proceed to obtain a suitable subsolution for v. To this end, we consider the equation 9 00 (η) + 9 0 (η) − 16 e−3η 9(η) = 0

for − ∞ < η < +∞.

We claim that, for any constant K > 0, there exists a solution of (5.14) such that:   2 −3η2 0 9(η) ∼ Kexp − √ e 9(η) > 0 and 9 (η) > 0 for all η > 0, 3 6

(5.14)

when η → −∞.

(5.15)

To check this, it is convenient to make use again of the transformation (5.11). Actually, on setting S = 9 0 /9, (5.14) gives S 0 + S 2 + S = 16 e−3η ,

(5.16)

a Riccati-type equation alike to (5.12). An analysis of the trajectories of (5.16) reveals that there exists a solution S(η) such that S(η) > 0 for all η, and S(η) ∼

e−3η/2 √ 6

for η → −∞,

S(η) ∼ Ae−η

for η → +∞,

(5.17)

for some A > 0. From (5.17),(5.15) follows at once by integration. Consider now the auxiliary function g(t) = (logt)/3, which satisfies g 0 = e−3g

for t > 0, g(0) = −∞.

(5.18)

Let us finally take R0 large enough (but fixed). We now claim that, if we select K sufficiently large in (5.15), the function w(ξ, t) = −min{9(ξ − g(t)), 1}

(5.19)

is a subsolution of the equation obtained by setting w = v − 1 in (5.5), namely wt = wξ ξ + 3wξ + 2wwξ + 2w(1 + w) − 13 e−3ξ w

(5.20)

in the region where ξ < −R0 , with initial and side conditions w(ξ, 0) = −1

for ξ ≤ −R0 ,

w(−R0 , t) = −1

for t > 0.

(5.21)

To show this, it actually suffices to observe that z(ξ, t) = −9(ξ − g(t)) is a subsolution for (5.20) and (5.21), whenever w = z. Notice that L(z) ≡ zt − zξ ξ − 3zξ − 2zzξ − 2z(1 + z) + 13 e−3ξ z = g 0 9 0 (η) + 9 00 + 39 0 − 299 0 + 29(1 − 9) − 13 e−3ξ 9   = e−3g 9 0 + e3g 9 00 − 16 e−3η 9 − 16 e−3ξ 9 + 39 0 − 299 0 + 29(1 − 9).

254

M. Escobedo et al. / Physica D 126 (1999) 236–260

As we are actually interested in what will happen up to times t ∗ = O((ε|logε|)−2/3 ), we may select K in (5.15) such that 9 00 > 0 for η ≤ −R0 , and such that g(t) < 0 for t ≤ t ∗ , so that e3g 9 00 ≤ 9 00 in (5.21). On the other hand, since 9(η) will be of order O(1) where w = z, we may assume that 39 0 − 299 0 + 29(1 − 9) ≤ C9 0 for some positive constant C. Bearing these remarks on mind, we deduce from (5.21) and (5.14) that, setting again S = 9 0 /9  −3ξ   −3g −3η  1 e 1 e e e−3ξ 9 + C9 0 = − + CS = − + CS . (5.22) L(z) ≤ − 6 9 6 9 6 Recalling (5.17), it now follows that for t ≤ t ∗ , L(z) ≤ 0, and we may apply the maximum principle to derive that w(ξ, t) + 1 ≤ v(ξ, t)

for − ∞ < ξ < ∞, 0 < t < t ∗ .

(5.23)

Putting together Eqs. (5.13) and (5.23), estimate (5.8a) and (5.8b) follows. We conclude this section by observing that assumption (5.9) is crucial in our previous argument. One may certainly wonder what would occur if such hypothesis is dropped, or, in another words, if we select an initial value f (x, 0) = g(x) which does not coincide with a Bose–Einstein distribution fµ (x) with µ > 0. While this case is not discussed in detail here, a rather simple argument may be advanced that suggests that the corresponding dynamics will be similar to the one just considered. To begin with, we know that when ε = 0, f (x, t) → fµ (x) for some µ ≥ 0 as t → +∞, so that it makes sense to expect that for 0 < ε  1, f (x, t) will be close to such fµ in times t = O(1) < t ∗ = O((ε|log ε|)−2/3 ). However, at such times (that could be reset as t = 0 by means of a translation) there could be regions where f (x, 0) > fµ , and we might have v > 1 there. In such a case, the effect of the last term in the right-hand side of (5.5) will consist in propagating an increasing perturbation to the region where ξ < 0 and |ξ |  1, which should be taken care of by introducing a boundary layer near ξ = −∞ (i.e., near x = 0). Since the analysis of such a region is rather technical, we shall omit further details. 6. Concluding remarks We have analysed in this paper the evolution of the distribution of radiation for a gas of photons interacting with a background of electrons by means of Compton scattering, under the assumption that we remain in a regime in which the photon–electron interaction (in particular Bremsstrahlung radiation) are small. To this end, we have studied a Fokker–Planck equation which has been derived as a suitable approximation to the corresponding Boltzmann equation in such regime. We have shown that the evolution of the distribution of photons is governed by a two-scale problem. In the shorter scale (t = O(1) in suitable units), the dynamics is initially driven by an equation obtained from (1.1) by setting ε = 0 there, which is known to develop singularities in finite time at x = 0, i.e., in the limit of low frequencies. However, before such singularities appear, the complete equation produces a meaningful decrease in the number of photons at low frequencies. In particular, this result is in sharp contrast with some customary assumptions in the literature (cf. for instance [3,9] etc.), that admit that at this scale the number of photons is preserved, a fact that could yield transient behaviours alike to the Bose–Einstein condensation for sufficiently large initial distributions of photons. The analysis in this paper shows that in general this is not the case. On the other hand, for very large times the distribution of photons approaches towards the Planck equilibrium distribution in a timescale of order t = O((ε|log ε|)−2/3 ). The precise size of this scale is due to the influence of Bremsstrahlung at low frequencies, that plays a key role in the long-time dynamics of solutions. Appendix A We shall prove herein a global well posedness result for problem (1.1)–(1.6) with ε > 0. This signals a difference with the case ε = 0, where solutions which develop singularities in finite time are known to exist (cf. [5]). More precisely, we show the following:

M. Escobedo et al. / Physica D 126 (1999) 236–260

255

Theorem A.1. Let g be a continuous and nonnegative function such that 0 ≤ g(x) ≤

C f or some C > 0, whenever 0 < x < +∞, x2

(A.1)

lim x 4 g(x) = 0.

(A.2)

x→+∞

Then, for every given ε > 0 there exists a unique solution f (x, t) of (1.1)–(1.6) which is defined for all times t > 0, and is such that lim x 2 f (x, t) = 0

f or all t ∈ (0, +∞),

x→0

lim x 4 f (x, t) = 0

(A.3)

f or all t ∈ (0, +∞).

(A.4)

Moreover, the following inequality holds: Z ∞ Z +∞ x 2 |f (x, t) − f0 (x)| dx ≤ x 2 |g(x) − f0 (x)| dx f or all t ≥ 0,

(A.5)

x→+∞

0

0

where f0 = (ex − 1)−1 . Proof. We proceed in several steps. As a startpoint, we shall first approximate problem Eqs. (1.1)–(1.6) by a sequence of suitable regularized problems. To this end we observe that, on making use of the change of variables (2.4), Eq. (1.1) reads vt = vyy + 2vvy + 2v 2 + vy − 2v + ey (vy + 3v) + ελ(y)(f0 ey − v), where λ(y) = e−3y K0



ey 2



−∞ < y < +∞, t > 0,

(A.6)

y

e−(1/2)e .

For any given integer n, we now consider the following auxiliary problem: vnt = vnyy − 2vn + 2vn vny + 2ψn (vn ) + vny + ey (vny + 3vn ) + ελ(y)(f0 ey − vn ), when − n < y < n t > 0,

(A.7)

where ψn (v) = min(v 2 , n2 ), together with side conditions vn (y, t) = 0

for y = ±n,

(A.8)

and initial data vn (y, 0) = ey g(ey )

for − n ≤ y ≤ n.

(A.9)

It then follows from classical theory of parabolic equations (cf. for instance [6]) that there exists a unique solution vn to Eqs. (A.7)–(A.9) which is defined for all times t > 0. For any α > 0, consider now the function v α (y) = xf α (x),

(A.10)

256

M. Escobedo et al. / Physica D 126 (1999) 236–260

where f α is given in Lemma 2. We readily see that, on selecting α large enough (depending on constant C in (A.1)), there holds vn (y, t) ≤ v α (y)

for − n < y < n and t > 0.

(A.11)

Indeed, one readily sees that vn (y, 0) ≤ v α (y) in the interval −n ≤ y ≤ n. On the other hand, v α is a stationary solution of (A.6) with ε = 0. Since the term ελ(y)(f0 ey − v α ) has a negative sign, inequality (A.11) follows from classical comparision results for solutions of parabolic equations (cf. [6]). On letting now n → +∞, we obtain a subsequence, still denoted by (vn ), for which (A.11) holds. From such a bound, classical regularity theory for parabolic equations, as described in [6], provides enough estimates on the derivatives of vn appearing in (A.7), so as to allow to pass to the limit in that equation when n → +∞. In this way, a classical solution v(y, t) of (A.6) is obtained which is defined for all −∞ < y < +∞ and t > 0, and satisfies 0 ≤ v(y, t) ≤ v α (y)

when − ∞ < y < +∞, t > 0.

(A.12)

We point out that, as a consequence of (A.12), we have that, for some constant C > 0, and all t > 0 0 ≤ f (x, t) ≤

C x2

for 0 ≤ x ≤ 1,

0 ≤ f (x, t) ≤

C x4

for x > 1.

(A.13)

To proceed further we check now that (1.3) holds. To this end, we examine first the region x  1, and for given η > 0 and C > 0 as in (A.13), we consider the auxiliary function h(x, t), solution of the following problem:    1 ∂ ∂h ∂h = 2 x4 + h + h2 + σ (x)(f0 − h), for x ≥ R, t > 0, ∂t ∂x x ∂x h(R, t) =

C R4

for t > 0, h(x, 0) =

η x4

for x ≥ R.

(A.14)

The existence of a global solution of (A.14) is obtained in a similar way as before, by means of approximation by suitable problems with bounded coefficients. By (A.2), for any η > 0, there exists R = R(η) > 0 such that f (x, 0) ≤

η x4

for x ≥ R.

(A.15)

We now set w(x, t) = x 4 h(x, t) and observe that w solves    w ∂ ∂w 4 w2 ∂w = x2 − w + w + 4 + σ (x) f0 − 5 , ∂t ∂x ∂x x x x w(R, t) = C,

for all t > 0,

w(x, 0) = η,

for x ≥ R, t > 0,

for x ≥ R.

We now claim that: For R > 0 large enough, (A.16) admits a stationary supersolution Z +∞  x 4 w R (x) = Ce−(x−R) + 2ηe−x x 4 s −4 es ds. R R

(A.16)

(A.17)

To check (A.17), we first observe that the result is true when ε = 0. To see this, we drop the subscript R in w R and argue as follows. Function w is actually the solution of the linear ODE problem   4 w = 2σ for x > R, w(R) = C. w0 + 1 − x

M. Escobedo et al. / Physica D 126 (1999) 236–260

257 0

Clearly, limx→+∞ w(x) = 2σ . We now set w∗ (x) = 2σ + 2w(x)/x, and observe that w∗ + (1 − 4/x)w∗ < 2σ for x > R. Since w ∗ (R) = 2σ + 2C/R < C = w(R) if R  1, it turns out that w(x) > 2σ + 2w(x)/x. From this inequality, we readily see that ! ! !   2 2 0 2 ∂w 4w 2ww ∂ w ∂ w w w w 2 2 2 0 − +w+ 4 =x 2σ + 4 = x −4 5 =2 2 w −2 x ∂x ∂x x ∂x x x x x4 x x   w w (A.18) = 2 2 2σ + 2 − w < 0, x x whence (A.17) follows under our current assumptions. When ε > 0, we merely observe that the term σ (x)(f0 − w/x 5 ) decays exponentially when x → +∞. Since the last term on the right-hand side in (A.18) is algebraic, (A.17) also holds in that case. We have thus obtained that lim sup x 4 f (x, t) ≤ 2η.

x→+∞

Since η > 0 may be chosen arbitrarily small, (A.4) holds. In order to obtain estimates for ∂f/∂x as x → +∞, we rescale as follows: gR (y, t) = R 4 f (x + R, R 2 t),

y = x + R.

Function gR then satisfies the equation     1 ∂g ∂g ∂ = (y + R)−2 (y + R)4 + g + 4 g 2 + σ (y + R)(f0 − g) , R2 ∂s ∂y ∂y R

for y ≥ −R, t > 0.

The bound (A.4) implies that gR (y, t) → 0 uniformly on |y| ≤ 1, t > 0 as R → +∞. Since the source term σ (y + R)(f0 − g) decays exponentially as R → +∞ on the same region, by classical parabolic theory we obtain that |∂gR /∂y| → 0 as R → +∞, uniformly on sets |y| < 1/2 and R 2 t ≥ 1. This implies in particular that   ∂f =o 1 as x → +∞, (A.19) ∂x x4 uniformly for any t > 0. Putting together (A.4) and (A.19) we obtain (1.3) when x → +∞. As a next step we proceed to prove (1.3) as x → 0. The argument on this region is quite different to that done for x  1. This is due to the fact that the key term to be reckoned with is now σ (x)(f0 − f ). As a matter of fact, when that term is dropped, condition (1.3) may fail to hold at x = 0 after a finite time, as it has been recalled before. Let us consider now the function G(t) given by the solution of G0 (t) = 2G2 − G − AεG4 |log g|,

t > 0,

lim G(t) = +∞,

t→0

(A.20)

where A is a positive constant to be selected presently. A dominated balance argument involving the first term on the left and the third on the right in (A.20) gives that G(t) ∼

1 (Aεt)1/3 |logt|1/3

as t → 0.

(A.21)

Moreover, if A is small enough, function H (G) = 2G2 − G − AεG4 |logG| has three non-negative roots, x3 > x2 > x1 = 0. Analysis of the phase diagram associated to (A.20) shows that limt→+∞ G(t) = x2 provided that G(t) 6= 0 and G(t) 6= x3 . From this and (A.21), we obtain that   1 for all t > 0 (A.22) G(t) ≤ C 1 + (Aεt)1/3 |logt|1/3

258

M. Escobedo et al. / Physica D 126 (1999) 236–260

and some positive constant C = C(A, ε). Furthermore, estimating the value of G(x2 ) by means of a dominated balance between the first and third terms on the right of (A.20), provides the lower bound G(t) ≥

1 1 10 (Aε)1/2 |logε|1/2

for t > 0.

(A.23)

Consider now the set 6 = {(y, t) : G(t) ≤ Av α (y)}. Clearly, we have that ey f0 (ey ) ≤ 2 on 6. Selecting A = A(ε) small enough, it follows from (A.23) that G(t) ≥ 2 ≥ ey ,

for all t > 0.

(A.24)

We now observe that G(t) is a supersolution to (A.6) in the set 6 provided that  y e y y −3y e−(1/2)e (G − 2) ≥ Aε|log(G)|G4 . (1 − 3e )G + εe K0 2

(A.25)

Recalling the asymptotic behaviour of v α = xf α (cf. (3.7)), it turns out that (A.25) holds if we select A < α −3/2 /2. To obtain a global supersolution to (A.6), we now define the auxiliary function V (y, t) = min{v α (y), G(t)}.

(A.26)

Actually, the function V defined above is also a supersolution for the approximating problems Eqs. (A.7)–(A.9). Back to the original variables, we have thus derived the bound √ α C , (A.27) f (x, t) ≤ (1 + (Aεt|log t|)−1/3 ), for 0 < x < x 2G(t) which implies (A.3). To fulfill the flux condition at x = 0, we need to improve estimate (A.27). This we shall do as follows. Consider the function h(x, t) = f (x, t) − f0 (x) ≡ f (x, t) − (ex − 1)−1 , which satisfies ∂ ∂h = x −2 ∂t ∂x



 x

4



∂h + h + h2 + 2f0 h ∂x

− σ (x)h.

From (A.27) and the definition of f0 , we see at once that 2C (1 + (Aεt|logt|)−1/3 ), |h(x, t)| ≤ x

√ α . for 0 < x < 2G(t)

(A.28)

For any given t0 > 0, we now select R > 0 such that √ 1 α . R≤ 4 4G(t0 ) Let us define GR (x, t) = Rh(Rx, t). Function GR satisfies    ∂GR ∂ ∂GR = x −2 x4 + RGR + G2R + 2Rf0 GR − σ (Rx)GR , ∂t ∂x ∂x

for − ∞ < x < +∞,

(A.29)

M. Escobedo et al. / Physica D 126 (1999) 236–260

259

whereas by (A.28) and our choice of R, we also have that GR (x, t) ≤ 4C(1 + (Aεt|log t|)−1/3 ) ≡ M(ε, t0 ) for

1 2

≤ y ≤ 2, t ≥ t0 .

(A.30)

We next multiply both sides of (A.29) by sgn G, and make use of Kato’s inequality i.e., the fact that for any smooth function F (y), ∂ 2 F /∂y 2 · sgn(F ) ≤ ∂ 2 |F |/∂y 2 in a suitable weak sense (for instance, in the sense of distributions). Recalling the definition of σ (x), we thus obtain that GR satisfies:    ∂|GR | −2 ∂ 4 ∂|GR | 2 ≤x x + R|gR | + gR + 2Rf0 |gR | ∂t ∂x ∂x |logR| |gR |, for t > t0 and 21 ≤ x ≤ 2, (A.31) −B R3 for some B > 0 independent of R. We now remark that the function   ch(β(|logR|1/2 /R 3/2 )(x − 1)) |logR| + M(ε, t0 )exp −B (t − t0 ) wR (x, t) = M(ε, t0 ) ch(β(logR|1/2 /2R 3/2 )) R3 ≡ wR1 (x, t) + wR2 (x, t),

(A.32)

where β is a fixed constant independent of R, is a supersolution of (A.31) in the region x ∈ [1/2, 3/2], t ≥ t0 . Actually, wR1 is obtained after examining the ODE x2

∂ 2 |G| B|logR| − |G| = 0, ∂x 2 R3

which contains the most relevant contribution of the first term on the right of (A.31) as well as the last one there, whereas wR2 is derived after considering the first order ODE B|logR| ∂|G| =− |G|. ∂t R3 In view of our bound (A.30), we readily see that wR (x, t) ≥ |GR (x, t)| when x = 1/2, x = 3/2 and t ≥ t0 , as well as that wR (x, t0 ) ≥ |GR (x, t0 )| for 1/2 ≤ x ≤ 3/2. By comparison, we then have that |GR (x, t)| ≤ wR (x, t) whenever t ≥ t0 and 1/2 < y < 3/2. In particular, this last inequality implies that, for any fixed δ ∈ (0, 1) with 0 < δ < 1, we have that, for some C1 > 0:    1/2  ∂GR ≤C1 M(ε, t0 )exp −β |logR| δ +Mexp −β |logR| δ , for 3 ≤ y ≤ 5 and t − t0 ≥δ. (A.33) ∂y 4 4 4R 3/2 4R 3 Notice that the second term on the right above is negligible with respect to the first one there whenever 0 < R  1. By (A.33) and classical regularity theory for (A.31), we now arrive at  1/2  ∂GR ≤ CMexp −β |logR| δ ∂x R 3/2 for some C > 0, provided that 4 5

≤ y ≤ 56 ,

t ≥ t0 + 2δ

and R > 0 is small enough, which in the original variables gives   ∂h CM(ε, t0 ) |logx|1/2 ≤ exp −β δ , for t ≥ t0 + δ, 0 < x ≤ R0 ≡ R0 (ε, t0 , δ). ∂x |x|2 x 3/2 From (A.28) and (A.34), condition (1.3) is shown to hold at x = 0.

(A.34)

260

M. Escobedo et al. / Physica D 126 (1999) 236–260

Finally, in order to conclude the proof of Theorem A.1, it remains to show (A.5) (which follows at once upon integrating (5.2) in time), as well as the uniqueness statement. This last is achieved by subtracting the equations satisfied by two possible solutions f1 , f2 corresponding to the same initial value f1 (x, 0) = f2 (x, 0) = g(x). Multiplying the equation thus obtained by sgn+ (f1 − f2 ) and integrating in space gives then Z +∞ x 2 (f1 − f2 )+ dx = 0, 0

whence f1 ≤ f2 . By interchanging then the roles of f1 and f2 , we eventually deduce that f1 = f2 .



References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]

C.M. Bender, S.A. Orszag, Advanced Mathematical Methods for Scientists and Engineers, Mc Graw Hill, New York, 1978. G. Cooper, Compton Fokker–Planck equation for hot plasmas, Phys. Rev. D 3 (1971) 231–22316. R.E. Caflisch, C.D. Levermore, Equilibrium for radiation in a homogeneous plasma, Phys. Fluids 29 (1986) 748–752. H. Dreicer, Kinetic theory of an electron–photon gas, Phys. Fluids 7 (1964) 735–753. M. Escobedo, M.A. Herrero, J.J.L. Velazquez, A nonlinear Fokker–Planck equation modelling the approach to thermal equilibrium in a homogeneous plasma, Trans. Amer. Math. Soc. 350 (1998) 3837–3901. A. Friedman, Partial differential equations of parabolic type, Krieger, New York, 1983. A.S. Kompaneets, The establishment of thermal equilibrium between quanta and electrons, Soviet Phys. JETP 4 (1957) 730–737. P.J.E. Peebles, Principles of Physical Cosmology, Princeton University Press, Princeton, NT, 1993. R. Weymann, Diffusion approximation for a photon gas interacting with a plasma via the Compton effect, Phys. Fluids 8 (1965) 2112–2114. Ya.B. Zel’dovich, D. Novikov, Relativistic Astrophysics, vol. 2: The Structure and Evolution of the Universe, The University of Chicago Press, Chicago, 1983.