Avalanche Photon Cooling by Induced Compton Scattering: Higher-Order Kompaneets Equation

Induced Compton scattering (ICS) is an interaction between intense electro-magnetic radiations and plasmas, where ICS transfers the energy from photons to plasmas. Although ICS is important for laser plasma interactions in laboratory experiments and for radio emission from pulsars propagating in pulsar wind plasmas, the detail of photon cooling process has not been understood. The problem is that, when ICS dominates, evolution of photon spectra is described as a nonlinear convection equation, which makes photon spectra to be multi-valued. Here, we propose a new approach to treat evolution of photon spectra affected by ICS. Starting from the higher-order Kompaneets equation, we find a new equation that resolves the unphysical behavior of photon spectra. In addition, we find the steady-state analytic solution, which is linearly stable. We also successfully simulate the evolution of photon spectra without artificial viscosity. We find that photons rapidly lose their energy by ICS with continuously forming solitary structures in frequency-space. The solitary structures have the logarithmically same width characterized by an electron temperature. The energy transfer from photons to plasma is more effective for broader spectrum of photons such as expected in astrophysical situations.


INTRODUCTION
Nonlinear interactions between strong electromagnetic waves and plasmas have been studied in the context of laser plasma physics and astrophysical phenomena. Depending on intensities of radiations and conditions of plasmas, plenty of nonlinear interactions exist and are studied in different approaches. For example, quantum electrodynamics predicts vacuumpolarization effects, where photons interact with virtual electron-positron pairs (c.f. Ref [1]). Parametric processes in plasmas, such as induced scattering processes, filamentation instability and so on, have been studied in classical electrodynamics with the two-fluid description of plasmas (c.f. Ref. [2]). There are also studies of parametric processes in the semi-classical formulation, which is used in this paper (c.f. Refs. [3][4][5]). Here, we focus on induced Compton scattering (ICS), which is a parametric instability between photons and each electron, rather than plasmon. The condition for ICS to be dominant process is when both the central frequency and the spectral width of strong electromagnetic radiations are greater than the Langmuir plasma frequency [6]. ICS cools photons toward the Bose-Einstein condensation, while plasmas are heated by this process (c.f. Ref. [7]). ICS has been studied in the context of laser plasma interaction both theoretically [8][9][10] and experimentally [11][12][13][14]. Most of them studied the process with the two-fluid approximation, which accounts only for the linear regime of the instability. On the other hand, ICS has also been studied for the scattering process around high intensity astrophysical sources with the semi-classical formulation, which can partly treat the nonlinear regime (e.g. Refs. [15][16][17][18]). However, they constrained plasma conditions around some astrophysical objects considering only the scattering optical depth to ICS. It is known that there is a difficulty to treat the nonlinear regime of ICS even with the semi-classical formulation as described in this paper.
In the semi-classical formulation, a relaxation process of isotropic photons interacting with a rarefied Maxwellian plasma by Compton scattering is described by the Kompaneets equation [31]. The Kompaneets equation is the Fokker-Plank equation describing evolution of a photon spectrum including the ICS term, which comes from the Boson nature of photons and, which is quadratic in the intensity of the radiation. Therefore, for high intensity radiations, this quadratic term plays a dominant role and the Kompaneets equation is reduced to the nonlinear convection equation (Eq. (22)). It is known that the nonlinear convection equation has the implicit solution, which will be multi-valued after a finite time starting from certain initial conditions. Such a spectrum is unphysical and should be bent to be single-valued, i.e., there should be an appropriate physical process which we have ignored (e.g. Ref. [19]).
Different approaches have been adopted to resolve the unphysical behavior of a photon spectrum when ICS dominates.   [20][21][22] considered heuristically to add a dispersive term (the third derivative term) to the nonlinear convection equation and made the equation Korteweg-de Vries-type, i.e., soliton formation in frequency-space.   [23,24] considered to recover the diffusive term (the second derivative term), which is originally included in the Kompaneets equation, and made the equation Burgers-type, which allows formation of a shock wave structure like hydrodynamics. However, because the original diffusive term is linear in the intensity of radiations, this will not be applicable when the intensity of radiations becomes higher and higher. On the other hand, Zel'dovich and colleagues [25,26] adopted the integral form of the equation, i.e., the Boltzmann equation (Eq. (1)). They predicted formation of solitary structures in a radiation spectrum. It is worth distinguishing which spectral behaviors appear in a given condition between photons and electrons. There are also numerical studies of this problem by Montes (1977Montes ( , 1979 [27,28] and Coppi et al. (1993) [29]. We compare their results with ours in Section 4.3.
In this paper, we propose a new approach to study evolution of a radiation spectrum when ICS dominates. In Section 2, we extend the Kompaneets equation to the higher-order terms. This is a preparation to obtain a new equation and its derivation helps understanding what is the resolution of the problem in the nonlinear convection equation. In Section 3, we find a new equation that describes evolution of a radiation spectrum for intense radiations. We obtain the analytic solution in steady state and also discuss its linear stability. In Section 4, we show some examples of numerical studies of the new equation. We give interpretations of our results and comparison with past studies. Section 5 is devoted to summary. and absorption of photons by plasmas, i.e., the photon number is conserved in this system (c.f. Ref. [7]). Effects of the background magnetic field and polarizations of photons are also neglected. Below, a wavenumber k is in unit of the inverse of the electron Compton wavelength λe ≡ /m e c and a momentum p of an electron is in unit of m e c. We start from the Boltzmann equation for photons including Compton scattering by free electrons of a density n e . We use a normalized momentum distribution of electrons f (p), i.e., d 3 pf (p) = 1. Assuming the spatial homogeneity of the system, we have where the terms 1 + n represent spontaneous and induced scattering terms, respectively. D = 1 − β · Ω and D 1 = 1 − β · Ω 1 represent the factor due to relative velocities between interacting photons and electrons, and (k 1 /k) 2 comes from the difference in phase space volumes between d 3 k 1 and d 3 k, where Ω = k/k, Ω 1 = k 1 /k 1 , and |β| = (1 + p −2 ) −1/2 , respectively. The differential scattering cross section for Compton scattering (the Klein-Nishina cross section) from an initial state k i to a final state k f of a photon is (e.g. Ref. [30]) where quantities with˜represent those in the electron-rest frame, γ = (1 + p 2 ) 1/2 is the Lorentz factor of an electron, µ = Ω i · Ω f is the cosine of the angle between k i and k f , δ is the Dirac's delta function, and σ T is the Thomson cross section. We have used the following Lorentz transformation laws and We execute the integral with respect to k 1 by using the δ function in Eq. (2) and introduce notations y = Θn e σ T ct, Θ ≡ k B T pl /m e c 2 , k ± ≡ γDk/(γD 1 ∓ k(1 − µ)), and n(k ± ) ≡ n(k ± , Ω 1 ), where k B is the Boltzmann's constant and T pl is an electron temperature. Then, Eq. (1) is expressed as ∆S and ∆I correspond to contributions from spontaneous and induced scattering, respectively. We find that net contributions to scattering are differences between an upper state k + and the incident state k for spontaneous scattering, and between k + and a lower state k − for induced scattering. Now, we assume that electrons and photons are non-relativistic, i.e., p ≪ 1 and k ≪ 1. The lowest order in p and k gives Thomson scattering where k ± ≈ k, n(k ± ) ≈ n(k, Ω 1 ), and R ± ≈ 3(1 + µ 2 )/16π. The result is where ∆S (ξ) and ∆I (ξ) represent terms on an order ξ for ∆S and ∆I, respectively. When the photon distribution is isotropic, we obtain ∆S 0 = 0. On the other hand, ∆I 0 = 0 is always satisfied because ICS requires the frequency shift in essence. On the first order with an assumption of isotropy, we obtain the well-known Kompaneets Equation. The right-hand side of Eq. (4) becomes where we have used 4πp 4 f (p)dp = 3Θ, i.e., the Maxwell-Boltzmann distribution. Though the result is simple, the derivation is fairly long. We just note the expansion of k ± on the first order in k and p 2 , where α = Ω p · Ω and α 1 = Ω p · Ω 1 . Odd orders in p will vanish because of isotropy ( dΩ 1 dΩ p α 2m−1 α 2l−1 1 = 0 with natural numbers m and l) and an order of p 2 corresponds to the first order of the plasma energy, i.e., the temperature Θ. Introducing normalized frequency x = k/Θ, Eqs. (4) and (10) give the Kompaneets equation [31] ∂n(x) ∂y In addition, multiplying Equation (12) by x 3 and integrating over dx, we obtain evolution of the photon energy density (ǫ ph ≡ x 3 n(x)dx) as (c.f. Ref. [32]) The first term, which comes from ∆S (k) of Eq. (10), means energy loss of photons by spontaneous Compton scattering. Corresponding induced term (∆I (k) ) is the second term, which 4/15 also implies energy loss of photons. On the other hand, the last term corresponding to ∆S (p 2 ) expresses energy gain by inverse Compton scattering, even though this term comes from the elastic approximation (the Thomson limit, zeroth order in k). It should be noted that there is no induced term associated with inverse Compton scattering, i.e., ∆I (p 2 ) = 0. As is wellknown, this Kompaneets equation is not enough to follow evolution of the photon spectrum for n ≫ 1. Now, we proceed to the next order. We have ∆S (k 2 ) , ∆S (kp 2 ) , and ∆S (p 4 ) for spontaneous scattering but only a cross term ∆I (kp 2 ) exists for induced scattering, because the numerator of Eq. (6) vanishes for even orders in each k and p. The expansion of k ± on the second order in k and p 2 (Θ) is Note that only the terms ∝ k have double sign on the right-hand side of Eq. (14) and this is the reason for ∆I (k 2 ) = ∆I (p 4 ) = 0, i.e., odd orders in k is required for induced scattering.
The integrals are carried out as and for spontaneous scattering and as for induced scattering. We also need the relativistic correction to the Maxwell-Boltzmann distribution f (p) = e −γ/Θ /4πΘK 2 (Θ −1 ), where K 2 is the modified Bessel function of the second kind of order 2. The moments of the distribution are written as 4πp 2 f (p)dp = 1,

5/15
This correction arises another second-order spontaneous term on the order of ∆S (p 4 ) in addition to Eq. (16) because we have where the second term of the last expression is on the order of ∆S (p 4 ) while the first term is already appeared in Eq. (12). Combining all of them, we obtain the higher-order Kompaneets equation that satisfies conservation of photon number, and the Bose-Einstein distribution is the equilibrium solution for this equation (c.f., the study for Eq. (12) is given by Caflisch and Levermore (1986) [33]). Note that the resultant equation is exactly the same as Eqs. (15) and (16) of Challinor and Lasenby (1998) [34], although they took a different approach to derive the equations.

The Case for Large Occupation Number: Steady-State Solution and Its Stability
We consider the case for large occupation number n ≫ 1. When we start from the first order (Eq. (12)), for n ≫ 1, we obtain where N (x) ≡ x 2 n(x). This is the nonlinear convection equation, which leads to an unphysical multi-valued solution. The situation is completely analogous to the shock formation in hydrodynamics so that the presence of a viscosity (the diffusion term in the right-hand side) can be the resolution (c.f. Ref. [35]). The viscous term arises when we extend the ideal fluid (the Euler equations) to the viscous fluid (the Navier-Stokes equations) (c.f. Ref. [36]). On the other hand, in this case of ICS, we will see that the physical resolution is not a diffusion term. Now, we start from the higher-order Kompaneets equation and find a new equation that resolves the difficulty of Eq. (22). Because we expand Eq. (1) assuming k ≪ 1 and Θ ≪ 1 (p ≪ 1), we neglect the second-order spontaneous terms (Eqs. (15) − (17) and (21)) compared with the first-order spontaneous terms and we obtain Efficiencies of each term on the right-hand side of Eq. (23) are characterized by scattering optical depths τ . In order of magnitude calculation, dividing each term by n/y, the right-hand side of Eq. (23) includes four components whose scattering optical depths are described as τ spon = xy (the first term: spontaneous Compton), τ inv = y (the second term: inverse Compton), τ ind = xyn (the third term: induced Compton) and τ sec = xynΘ ≪ τ ind (the other terms proportional to Θ: second-order induced Compton), respectively. ICS becomes dominant process if τ ind ≫ τ spon , τ inv , i.e., n × min(1, x) ≫ 1. If we consider only the τ ind term, the mathematically troublesome Eq. (22)   τ sec / max(τ spon , τ inv ) = nΘ × min(1, x) ≫ 1, the first-order spontaneous terms are neglected compared with the second-order induced terms even though the τ sec terms are higher-order in k or Θ. This is because the spontaneous and induced terms have different dependence on n so that this does not mean an importance of the third and further higher-order terms in k or Θ. The large occupation number satisfying the condition nΘ × min(1, x) = min(nΘ, nk) ≫ 1 is expected in some situations, for example, n ∼ 10 42 (brightness temperature T b (ν) ≡ hνn(ν)/k B ∼ 10 41 K at frequency ν ∼ 10GHz) is implied from an observation of the Crab pulsar in astrophysics [37] and n ∼ 10 22 (petawatt laser with ∼ 3.3 nm spectral width at wavelength ∼ 10.5µm) is obtained from some current laser facilities (e.g. Ref. [38]). In this case, we find that the natural extension of Eq. (22) is written as Note that Eq. (24) exactly corresponds to the isotropic case of Eq. (4) with ∆S = 0 and ∆I = ∆I (k) + ∆I (kp 2 ) . In contrast to Eq. (22), the plasma temperature Θ appears explicitly on the right-hand side that corresponds to ∆I (kp 2 ) (Eq. (18)). Because the right-hand side is the higher-order correction term, the velocity of photon flux (in frequency-space) is basically determined by the second term of the left-hand side, i.e., order N (x) to the negative direction of x. The first term of the right-hand side is important to resolve the problem described in Section 1, while the second term slightly modifies the second term of the left-hand side. Eq.
Eq. (26) shows no characteristic photon energy x, because the first-order spontaneous terms are neglected in Eq. (24), where only τ inv has different dependence on x from τ spon , τ ind and τ sec . This steady-state solution N (x) includes the solution of Eq. (22), as A = 0, i.e., the first term of Eq. (26) is the important contribution from ∆I (kp 2 ) . We require B ≥ A for N (x) > 0. Inside the derivative of the second line of Eq. (25) is written as a constant (7Θ/5)k 2 Θ (B 2 − A 2 ) ≥ 0, i.e., photon flux in frequency-space is negative direction of x. The photon flux is zero for A = B, i.e., contribution from ∆I (kp 2 ) compensates ∆I (k) term in this case and plays as a kind of heating, while ∆I (k) term plays a role of cooling (see Eq. (13)). We require proper boundary conditions at finite frequencies to determine the constants of integration, because the number density N (x)dx and the energy density xN (x)dx of photons diverse for either x → 0 or x → ∞. Since the photon flux is the negative direction, a photon injection at a large x and absorption at a small x are expected to achieve the steady-state (c.f., a similar study for Eq. (12) is given by Dubinov (2009) [39]). We plot Eq. (26) in Fig. 1, for example. One of important features of Eq. (26) is that the Doppler wavenumber k Θ ≈ Θ −1/2 associates with logarithmic scale in frequency x. This arises from the third-derivative term in Eq. (24) and is interpreted as the Doppler width ∆ν ∼ Θ 1/2 ν by induced-inverse Compton process (∆I (kp 2 ) ) already discussed by Zel'dovich and Sunyaev (1972) [25] in the integral form. Note that this Doppler width is a different process from the first-order spontaneous inverse Compton scattering ∆S (p 2 ) . For the later convenience, we define the Doppler wavelength λ Θ ≡ 2π/k Θ . For every integer m, the peaks of N (x) appear at x = Φe mλ Θ with a constant Φ = e −λ Θ φ/2π . Finally, we analyze linear stability for the steady-state solution N (x). We substitute N (x, y) = N (x) + N 1 (x, y) (N 1 /N ≪ 1) into Eq. (24) and linearize the equation, Considering a Fourier component of N 1 (x, y) ∝ e i(ωy−κx) , we obtain the dispersion relation Because Im(ω) > 0, the steady-state solution Eq. (26) is linearly stable. 8/15

Numerical Simulation
Here, we numerically study the initial evolution of photon spectra by solving Eq. (24) with n ≫ 1. We pay particular attention to evaluate the third-derivative in Eq. (24) precisely in order to see developments of solitary structures quoted by Zel'dovich and colleagues [25,26].

Set up
Eq. (24) is solved by fourth-order in x derivative (five-point stencil) and fourth-order in y derivative (fourth-order Runge-Kutta). We use x grid points spaced linearly from 0.1 to 10 and put fixed boundary conditions at both sides of x boundaries. Initial spectra are assumed to be a Gaussian spectrum, including three parameters, a normalization N init , a mean x init and a half width σ init of a Gaussian. The ratio of the initial spectral width σ init to the Doppler wavelength λ Θ ≡ 2π/k Θ , which is a parameter associated with Eq. (24), characterizes the spectral evolution. While sufficiently broad spectra (σ init ≫ λ Θ ) are expected in astrophysical situations, spectral widths of laser beams in laboratories are not always broad. We choose initial conditions that satisfy n init Θ × min(1, x init ) ≫ 1 (x 2 init n init ≡ N init ). Since the optical depth to ICS is written as τ ind = yN/x, we measure time y with y 0 ≡ x init /N init and an adequate time-step is searched to resolve spectral evolution for each calculation.
Eq. (24) is approximated form of Eq. (23) for nΘ × min(1, x) ≫ 1. We checked that both of the equations give almost the same initial evolution until y 0.4y 0 for the same parameter sets, satisfying nΘ × min(1, x) ∼ 10 and 10 2 at the peak of a Gaussian x = x init . For 0.4y 0 y < y 0 , numerically unstable structures start to develop at the portion of N (x, y) ≪ N init in our current scheme. This does not allow us to pursue the spectral evolution for different σ init and λ Θ beyond y ∼ 0.4y 0 . In what follows, we show numerical results of Eq. (24) for y ≤ 0.4y 0 . We fixed (N init , x init ) = (10 7 , 1), since different sets of (N init , x init ) does not change results as long as we measure y with y 0 . For remaining two parameters, we study the case that both σ init and λ Θ are an order of 10 −1 (λ Θ = 0.1 corresponds to an electron temperature of ∼ 10 2 eV).
Eq. (24) has two invariants of motion, the photon number density N (x, y)dx and a quantity ln N (x, y)dx [25]. Conservation of photon number is immediately found when we rewrite x derivatives in Eq. (24) to the second line expression of Eq. (25) and integrate both sides the equation over dx. We obtain conservation of the quantity ln N (x, y)dx from dividing both sides of Eq. (24) by N (x) and integrating over dx. For all the results of our calculation below, we checked that these quantities were conserved for y ≤ 0.4y 0 .  Table 1. Numbering of solitary structures corresponds to the number i on Table 1.  solitary form results from the first term of the right-hand side of Eq. (24). The heights of each solitary structure increase with time without changing their logarithmic width. The smaller value of λ Θ , i.e., the colder plasma, the more and the narrower solitary structures are formed. This behavior is consistent with what is inferred from the steady-state solution (Eq. (26)). On the other hand, numerical simulations without the third-derivative in Eq. (24) gives an unstable shock-like discontinuous structure rather than solitary structures with the use of the same scheme.

Results
To characterize the solitary structures, we consider a zeroth-order log-normal distribution for each solitary structure and fit some of the results at y = 0.4y 0 with a function where we adopt a half width of w = λ Θ /4 considering Eq. (26). Fitted values of n i , m i and the separation between the neighboring peaks ∆m i (for example, ∆m 1a = m 2a − m 1a ) are listed in Table 1. Numbering of solitary structures is found on the insets in Fig. 2; i = 1a − 4a for (σ init , λ Θ ) = (0.20, 0.10), i = 1b − 7b for = (0.20, 0.05) and i = 1c − 4c for = (0.10, 0.05), respectively. We do not fit all the structures because some structures have a smaller amplitude than N 0 (x). For example, we see more structures on the right of the structure numbered '4c' of the inset on the right panel in Fig. 2. Logarithmic widths of the structures are well characterized by the Doppler width λ Θ . However, the spacing of neighboring peaks ∆m i does not follow Eq. (26), i.e., ∆m i = λ Θ (see Table 1), or rather ∆m i is larger for lower frequency peaks. This relates with larger n i at lower frequency, because the velocity of a wave is proportional to N (x) in Eq. (24). Fitted values m 1a < m 1b indicate that the velocity of solitary structures is slower for larger value of λ Θ . Although the velocity relates with energy loss of photons, it is not easy to discuss the dependence on (σ init , λ Θ ) with Fig. 2, because number and height of structures are also different for different sets of (σ init , λ Θ ) (see the discussion about photon energy transfer in Section 4.3).

Discussion
Radio emissions from pulsars (e.g. Ref. [37]) and also fast radio bursts (e.g. Refs. [41]) have extremely large brightness temperature. Although the optical depth to ICS τ ind would attain to unity for these objects (c.f. Refs. [15,16,18]), the past studies did not discuss what kind of signatures are expected to be imprinted in their observed spectra. Our present study predicts a spectral break at the frequency ν ind corresponding to τ ind (ν ind ) 0.2 and solitary structures below ν ind in the photon spectra (red lines in Fig. 2). In this regard, the discrete emission bands in the dynamic spectra of the giant radio pulse occurred at the interpulse phase in the Crab pulsar reported by Hankins & Eilek [37] ('zebra bands') are intriguing phenomena. If we interpret their reported value ∆ ln ν = ∆ν/ν ∼ 0.06 by ICS, the value is not far from ∆m i for λ Θ = 0.05 in Table 1 corresponding to the electron temperature of ∼ a few ×10 eV. However, for applications to realistic astrophysical situations, we should take account for anisotropic photon distributions and for relativistic effects of both bulk and thermal motions of plasmas. These effects will be studied in subsequent papers. While the quantities N (x, y)dx and ln N (x, y)dx are conserved with time, the photon energy density ǫ ph (y) = xN (x, y)dx decreases with time, i.e., the energy of photons is transfered to plasmas by ICS [32,40]. The rate for the energy transfer by ICS is estimated 11/15 from the integration of Eq. (22), dǫ ph (y)/dy = − N 2 (x, y)dx (see Eq. (13)). Considering the initial phase of evolution y ≪ y 0 , i.e., N (x, y) ≈ N 0 (x), we estimate the rate as where we approximate ǫ ph (y) = xN (x, y)dx ≈ 2x init σ init N init and recall y 0 = x init /N init for the last expression. Eq. (31) can be applicable to describe an initial evolutionary phase and gives exponential loss of the photon energy with time ǫ ph (y) ≈ ǫ ph (0) exp(−y/y 0 ). Fig. 3 plots evolution of the normalized photon energy density with time for the results of Fig. 2. We also show the case for an initial width of σ init = 0.05 as dotted lines in Fig. 3. Thin lines (σ init = 0.20) show the exponential energy loss, although the slope is not exactly ln[ǫ ph (y)/ǫ ph (0)] = −y/y 0 . Thick (σ init = 0.10) and dotted (σ init = 0.05) lines deviate from the exponential-law (thin lines) at y ∼ 0.15y 0 and y ∼ 0.05y 0 , respectively. Dotted red line (σ init , λ Θ ) = (0.05, 0.10) and thick black line (σ init , λ Θ ) = (0.10, 0.20) show a bit different slope even at y 0.05y 0 from other lines and these exceptional behaviors are characterized by a relatively large Doppler wavelength λ Θ = 2σ init , which will be discussed later in this section. The energy transfer from photons to plasmas by ICS is effective for the case of broad initial spectra than narrow ones (the blue lines are always below the red lines for the same σ init ). This means that the deviation of N (x, y) from the initial spectrum N 0 (x) is faster for narrower spectra and this behavior is found from the inset on the right panel of Fig.  2 (the blue line is well below the green line compared with those on the left panel of Fig.  2). Although we see a temperature dependence of the energy transfer for the same initial width, Eq. (31) has no information of electron temperature explicitly. The third-derivative in Eq. (24) may play as a heating term because larger λ Θ (higher temperature) gives slower energy transfer in Fig. 3. This is consistent with the discussion at the steady-state solution Eq. (26), i.e., the photon flux is written as ≈ B 2 − A 2 for Θ ≪ 1 and a contribution from ∆I (kp 2 ) (A) resists that from ∆I (k) (B).
We should note the case of λ Θ = 2σ init , i.e., full width of the initial spectrum is almost the same as the Doppler width. For λ Θ > 2σ init , evolution of photon spectra is unstable at least in our numerical scheme. In those calculations, although all of the results shown in this paper have no time variation of N (x > x init , y) (see Fig. 2), there appears unstable features at x x init for λ Θ > 2σ init . This unstable feature is not improved by higher numerical resolutions of both x and y. Currently, it is uncertain whether the exceptional behaviors of dotted red line (σ init , λ Θ ) = (0.05, 0.10) and thick black line (σ init , λ Θ ) = (0.10, 0.20) in Fig. 3 is numerical or physical. However, we should also take care of our formulation (Eq. (1)) for a narrow spectrum σ init ≪ λ Θ . Galeev and Syunyaev (1973) [6] argued that, when σ init ≪ λ Θ , collective behaviors of plasmas are important rather than the interaction with free electrons, i.e., Compton scattering. This is the reason why we do not study the case for λ Θ > 2σ init more in this paper, although study of this regime is important for the application to laboratory experiments.
We should also note that whole of present paper is based on the assumption Θ = const. As seen in Fig. 3, photons clearly lose energy by ICS and transfer their energy to plasmas. The time-scale of energy transfer to electrons is estimated from Eq. (31) as where ǫ pl ≈ n e Θ is the energy density of a plasma divided by m e c 2 . We need to take into account evolution of a plasma temperature for evolution beyond this time-scale (c.f. Refs. [23,24]). Note that Eq. (32) is allowed to be used even for ǫ pl (y) > ǫ ph (y) (c.f. Ref. [7]). Our calculation needs only x init , N init , σ init and λ Θ that determine y 0 , ǫ ph (0), and Θ, i.e., we do not specify n e and t separately. Typical time-scales of photon cooling t ph (Eq. (31) and also Eq. (24)) is found from y ph ∼ y 0 and that of electron heating t pl (Eq. (32)) is found from y pl ∼ y 0 (ǫ pl /ǫ ph ). These are expressed as where the energy densities of photons and electrons are ǫ * ph = cΘ 4 /(π 2 λ -4 e ) xN (x)dx and ǫ * pl ∼ n e Θm e c 2 including dimension, respectively. For given x init , N init , σ init and λ Θ , photons lose energy before heating electrons (t ph ≪ t pl ) for sufficiently large n e . For example, we obtain t ph t pl ∼ 10 ps for n e 10 22 cm −3 (close to conditions of laser experiments) with the parameters which we used in this paper.
We finally discuss differences from past calculations of spectral evolution by ICS [28,29]. Fig. 1 of Coppi et al. (1993) [29] is spectral evolution for isotropic case. Basically, they solved Eq. (22) but including numerical viscosity, i.e., their calculations did not have information of a plasma temperature (λ Θ ). Although their result (the right panel of Fig. 1 of Coppi et al. (1993) [29]) showed two solitary structures that have the same logarithmic width but their behaviors seem different from ours, such as the height of the low frequency structure is smaller 13/15 than that of the high frequency one. We consider that these structures are numerical, i.e., combination of the numerical viscosity and logarithmically-spaced frequency grids in their calculations. Montes (1979) [28] solved the integro-differential equation given by   [26], but with some simplifications. Fig. 2 of Montes (1979) [28] (the calculation for an initially broad Gaussian spectrum) shows solitary structures of the linearly same width. We consider that Eq. (27) of Montes (1979) [28] is over simplified, especially, their assumption ∆ν ≈ Θ 1/2 ν 0 = const. in their Eq. (20) is crucial, while we consider the case ∆ν ≈ Θ 1/2 ν [25].

Conclusions
In this paper, we study evolution of photon spectra when ICS dominates. To get rid of the well-known difficulty of the nonlinear convection equation (Eq. (22)), we consider the higherorder Kompaneets equation. We obtain the new equation (Eq. (24)) that describes evolution of photon spectrum by ICS and that overcomes the difficulty. The second-order induced term ∆I (kp 2 ) (Eq. (18)) obtained from the higher-order Kompaneets equation improves the formulation. In addition, Eq. (24) has the steady-state analytic solution (Eq. (26) and Fig.  1), which is linearly stable. The steady-state analytic solution predicts the formation of solitary structures of logarithmically same width in frequency-space and the ∆I (kp 2 ) term plays as a heating term against the ICS term ∆I (k) , which is the pure cooling term.
We also study evolution of photon spectra by ICS numerically. ICS intermittently forms solitary structures moving toward lower frequency (Fig. 2) and these behaviors are consistent with predictions by some of past studies [25,26]. The solitary structures have the logarithmically same width well characterized by the Doppler width λ Θ and this behavior is also inferred from the steady-state solution. On the other hand, the spacing between peak frequencies of the structures does not follow Eq. (24) (∆m i = λ Θ ). The number and height of solitary structures depend on λ Θ .
The results of our numerical simulation satisfy the two invariants of motion, conservations of photon number N (x)dx and of the quantity ln N (x)dx. On the other hand, the energy density of photons xN (x)dx is transferred to electrons (Fig. 3). Energy density decays exponentially in an initial phase and this behavior is almost consistent with the analytic estimate (Eq. (31)). The energy transfer from photons to plasmas by ICS is effective for the case of broad initial spectra such as expected in astrophysical situations. The numerical results shown in the present paper are only initial phases of evolution y ≤ 0.4y 0 . We need more optimized numerical scheme to study evolution for y > 0.4y 0 .