1. INTRODUCTION
High-gain targets for the inertial confinement fusion (ICF) are meant for the ignition of a small part of the fuel with the following propagation of the thermonuclear burn wave on its main part. As the deuterium–tritium (DT) reaction rate at the temperatures of about 10 keV is much higher than the rates of other thermonuclear reactions, the ignition of the equimolar DT mixture is considered commonly as the primary task.
There are two approaches to ICF. One of the two is based on compression of a spherical layer of the fuel by a single driver providing sufficiently high values of density and temperature (Lindl et al., Reference Lindl, Amendt, Berger, Glendinning, Glenzer, Haan, Kauffman, Landen and Suter2004). At that, the thermonuclear burn wave arises in the target volume. Another approach known as the fast ignition is based on using two drivers (Basov et al., Reference Basov, Gus'kov and Feoktistov1992; Tabak et al., Reference Tabak, Hammer, Glinsky, Kruer, Wilks, Woodworth, Campbell, Perry and Mason1994; Gus'kov, Reference Gus'kov2013). The first driver compresses the target up to necessary value of density, while the second driver provides for fast rise of temperature. As a variant of such approach, one can consider ignition at the target center by a converging shock wave (Shcherbakov, Reference Shcherbakov1983; Betti et al., Reference Betti, Zhou, Anderson, Perkins, Theobald and Solodov2007). A possibility of using multi-shock wave compression for fast ignition is considered by Eliezer and Martinez Val (Reference Eliezer and Martinez Val2011).
The present paper is devoted to the variant of the fast ignition for which the burn wave arises near the target surface and propagates inside the target. Laser-generated beams of electrons (Tabak et al., Reference Tabak, Hammer, Glinsky, Kruer, Wilks, Woodworth, Campbell, Perry and Mason1994), and ions (Roth et al., Reference Roth, Cowan, Key, Hatchett, Brown, Fountain, Johnson, Pennington, Snavely, Wilks, Yasuike, Ruhl, Pegoraro, Bulanov, Campbell, Perry and Powell2001; Gus'kov, Reference Gus'kov2001; Caruso & Strangio, Reference Caruso and Strangio2003; Honrubia et al., Reference Honrubia, Fernandez, Temporal, Hegelich and Meyer-ter-Vehn2010) as well as heavy ion beams (Churazov et al., Reference Churazov, Aksenov and Zabrodina2001; Medin et al., Reference Medin, Churazov, Koshkarev, Sharkov, Orlov and Suslin2002; Gus'kov et al., Reference Gus'kov, Il'in, Limpouch, Klimo and Sherman2010), macroparticles (Caruso & Strangio, Reference Caruso and Strangio2001), and high velocity flows of matter (Gus'kov, Reference Gus'kov2013; Gus'kov & Zmitrenko, Reference Gus'kov and Zmitrenko2012) are considered as a driver for the fast heat of the highly-compressed fuel. Results of many theoretical works show initiation of burn waves propagating inside the target at certain values of beam parameters.
The contemporary concept of fast ignition assumes that the fuel is compressed up to the initial density ρ0 ~ 103ρs, where ρs is the density of the solid-state DT-mixture at atmospheric pressure and temperature of 4 K. Such preliminary compression of the fuel is a serious technical problem due to, particularly, Rayleigh–Taylor instability (Anisimov et al., Reference Anisimov, Drake, Gauthier, Meshkov and Abarzhi2013) arising at the stage of deceleration of a heavy shell while compressing the fuel. In the present paper, we study a possibility to use targets with the fuel density ρ0 ≤ 100ρs. In this connection, the work by Avrorin et al. (Reference Avrorin, Bunatyan, Gadzhiev, Mustafin, Nurbakov, Pisarev, Feoktistov, Frolov and Shibarshov1984) should be mentioned, where propagation of thermonuclear detonation on unevenly compressed and heated non-spherical target is studied numerically, in particular, from a small part of the hot (T = 10 keV) and dense (ρ0 ≈ 900ρs) DT mixture to a cold and much less compressed part (ρ0 ≈ 45ρs).
The first theoretical estimate of the ignition energy threshold for DT fuel of the density ρ0 was obtained by Tabak et al. (Reference Tabak, Hammer, Glinsky, Kruer, Wilks, Woodworth, Campbell, Perry and Mason1994) in the form E ig ~ ρ0−2. The similar estimate with an improved value of the constant of proportionality was obtained by Atzeni (Reference Atzeni1999). In many papers (see, for example, Caruso & Pais, Reference Caruso and Pais1996; Atzeni, Reference Atzeni1999; Churazov et al., Reference Churazov, Aksenov and Zabrodina2001; Caruso & Strangio, Reference Caruso and Strangio2003), the ignition energy threshold is determined from numerical solutions of the two-dimensional (2D) axially symmetric problem with given radius and time dependencies of the particle beam intensity. The improved dependency E ig ~ ρ0−1.85 based on results of many computations for ρ0 ≥ 50 g/cm3 ≈ 230ρs is proposed by Atzeni (Reference Atzeni1999). Note also a number of simple burn models based on ordinary differential equations of evolution, which allow estimating of ignition parameters for uniformly heated volumes of fuel (Nayak & Menon, Reference Nayak and Menon2012).
In the following three works, shell cylindrical targets with pre-compressed DT fuel irradiated from the target end are studied using 2D numerical codes. The initiation of a thermonuclear burn wave in DT fuel of the density about 500ρs by a heavy ion beam was considered by Churazov et al. (Reference Churazov, Aksenov and Zabrodina2001). A possibility to decrease the ignition energy by heating a small part of the fuel near the symmetry axis before the fuel density reaches its maximal value was considered by Caruso and Strangio (Reference Caruso and Strangio2003). Propagation of thermonuclear burn waves within a gold shell for DT fuel density of 50, 100, and 200 g/cm3, as well as the ignition of such a target by heavy ion beams were considered by Ramis and Meyer-ter-Vehn (Reference Ramis and Meyer-ter-Vehn2014).
Symmetrically converging plane thermonuclear burn waves initiated by laser pulses for the initial density ρ0 = ρs and 5ρs have been studied by Khishchenko and Charakhch'yan (2013) using the model of total absorption of laser radiation in the point with the critical density. For such kind of model, epithermal particles generated by laser-plasma interaction and heating the domain of supercritical density are ignored. A slow combustion wave arises only after at least one interaction of the shock wave reflected from the symmetry plane with the ablation front and moves in the fuel preliminary compressed and heated by several shock waves. The slow wave generates before itself a compression velocity profile that increases rapidly the fuel density. As a result, the slow combustion wave can transform into two detonation waves moving in opposite directions. After reflecting the slow or detonation wave from the symmetry plane, the intensive burning continues and increases considerably the burn-up factor of the target.
Further (Charakhch'yan & Khishchenko, 2013), the study of thermonuclear burn waves at ρ0 = ρs and 5ρs was supplemented by the waves generated by proton beams. It was shown that in spite of different ways of ignition, various models of α-particle heat, whether the burn wave remains slow or transforms into the detonation wave, and regardless of way of such a transformation, the final value of the burn-up factor depends essentially on the only parameter Hρ0, as in the known approximate formulas for the fuel expansion in spherical geometry (Basko, Reference Basko2009).
In the present paper, we extend our study to the more wide range of the initial density 5ρs ≤ ρ0 ≤ 100ρs. The one-dimensional (1D) problem to be considered can be treated as a rough model of burning the target shown in Figure 1. The fuel of the initial density ρ0 is in the cylindrical channel with the length of 2H surrounded by a shell of a heavy material. Following, for example, Caruso and Strangio (Reference Caruso and Strangio2003), the monoenergetic beam of protons of kinetic energy 1 MeV is considered as the ignition driver. The beam intensity is 1019 W/cm2 and the duration is 50 ps. The beam of the same energy with the intensity 1018 W/cm2 and the duration 500 ps is also considered. The burn-up factor and the gain obtained from the 1D calculation should be considered as the maximum possible values for the target of the same initial density and length.
Fig. 1. The target scheme.
For the first time such type of target was suggested by Pashinin and Prokhorov (Reference Pashinin and Prokhorov1972), for compression of a gaseous deuterium by laser pulses. Note that the target was experimentally studied by Stöckl and Tsakiris (Reference Stöckl and Tsakiris1991).
The heavy shell can be compressed by a magnetic field (Prut et al., Reference Prut, Khrabrov, Matveev and Shibaev1979) or by heavy ions depositing energy to a certain shell of a multilayer target (Churazov et al., Reference Churazov, Aksenov and Zabrodina2001; Dolgoleva & Zabrodin, Reference Dolgoleva and Zabrodin2006). So we suppose that it is possible to create the configuration close to that shown in Figure 1 with a cylinder of the fuel compressed to the necessary density and with two holes for injecting the proton beams.
Let us assume that the key effect determining the fuel ignition is the heat by α-particles while the lateral expansion of the fuel, the heat transfer from the fuel to the shell and the self-radiation of plasma are insignificant. Such assumption has certain grounds. For the initial density ρ0 = 100ρs ≈ 22 g/cm3, the shell density in its compression can exceed the above value more than ten times. For example, Ramis and Meyer-ter-Vehn (Reference Ramis and Meyer-ter-Vehn2014) considered the gold shell with a density of 400 g/cm3. For the time 150 ps, during which a detonation wave arises in our calculations (see Section 3), the fuel radius in the target can be changed insignificantly. In the case of compression by a strong magnetic field, the heat flux between the fuel and the shell can be essentially diminished (Pashinin & Prokhorov, Reference Pashinin and Prokhorov1972). Finally, we refer to the work by Churazov et al. (Reference Churazov, Aksenov and Zabrodina2001), where the ignition of a similar target by a heavy ion beam is simulated numerically. The authors conclude that the radiative transfer is not determinative when the ignition.
If the above assumption is true, there is an interesting possibility to evaluate the ignition energy of the cylindrical target E ig using 1D calculations. To ignite the target, it is necessary to leave a significant part of the α-particles energy within the fuel. To determine approximately the ignition energy of the cylindrical target in Figure 1, we introduce the parameter R α and solve the 1D problem taking into account the escape of α-particles from a cylinder of the radius R α identifying the latter with the beam and the fuel radiuses. If the solution of the 1D problem contains the thermonuclear burn wave, we set E ig = πR α2I(∞), where I = I(t) is the proton beam energy inserted at the time t per unit of cross-sectional area
$$ I\lpar t\rpar = \int\limits_0^t \, {J_b}\lpar t^{\prime}\rpar {\rm d}t^{\prime}\comma \;$$where J b(t) is the beam intensity given in the 1D problem.
The trajectories of α-particles are relying straight lines. A possibility to use strongly magnetized cylindrical shells confining the trajectories of α-particles and therefore reducing the beam area necessary for the ignition (Kemp et al., Reference Kemp, Basko and Meyer-ter-Vehn2003) is not considered here. For simulation of α-particle transport, we use the track method (Brueckner & Jorna, Reference Brueckner and Jorna1973; Duderstadt & Moses, Reference Duderstadt and Moses1982) and its modification that is a numerical method for the known Cauchy problem for the kinetic steady-state homogeneous equation in Fokker–Plank approximation (Gus'kov & Rozanov, Reference Gus'kov and Rozanov1982).
The argument in favor of the quasi-1D model described above is its qualitative agreement with theoretical estimates and numerical simulations of the fuel ignition by the particle beam of the radius much less than the fuel size. Apart from numerical results presented in Section 5, one can see such agreement from the following simple reasoning. The mean free-path length of α-particles during the isochoric heating is about l α ~ ρ0−1. The energy fraction of α-particles that remains within the channel of the radius R α can be estimated as β ≈ R α/l α ~ R αρ0. Assuming that β is independent of ρ0, obtain R α ~ ρ0−1, and, setting that I(t) is also independent of ρ0, obtain the known estimate E ig ~ ρ0−2 (Tabak et al., Reference Tabak, Hammer, Glinsky, Kruer, Wilks, Woodworth, Campbell, Perry and Mason1994).
A mathematical model is presented in Section 2. The target ignition by proton beams of different intensity and the same energy, as well as generation of a detonation wave are considered in Section 3. The detonation wave reflection from the symmetry plane is studied in Section 4. Integrated characteristics for the initial density 5ρs, 25ρs, and 100ρs are discussed in Section 5. Conclusions are summarized in Section 6. Additionally, the track method for simulation of α-particle transport and its modification using the known Cauchy problem for the kinetic steady-state homogeneous equation in the framework of Fokker–Plank approximation are described in Appendix.
2. PROBLEM STATEMENT AND NUMERICAL METHOD
We use an equation of state (EOS) of hydrogen p = p H(ρ, T), ε = εH(ρ, T) based on the wide-range semiempirical EOS model (Khishchenko, Reference Khishchenko2008). Here p is the pressure, ε is the specific internal energy, ρ is the density, and T is the temperature. The EOS provides for a good agreement with the Thomas–Fermi model with quantum and exchange corrections (Kalitkin, Reference Kalitkin1960; Kalitkin & Kuzmina, Reference Kalitkin and Kuzmina1975) at high densities (ρ≳2 g/cm3 for H) and low temperatures, as well as gets a form of model for ideal-gases mixture of electrons and nuclei at moderate densities and high temperatures.
At the initial point of time t = 0, a motionless plane layer of the equimolar DT mixture occupies the domain 0 ≤ x ≤ H. The one-temperature EOS of the medium is described by formulas
$$ p = p\lpar {\rm \rho} \comma \; T\rpar = {\,p_{\rm H}}\lpar {A^{ - 1}}{\rm \rho} \comma \; T\rpar \comma \; \quad {\rm \varepsilon} = {\rm \varepsilon} \lpar {\rm \rho} \comma \; T\rpar = {A^{ - 1}}{{\rm \varepsilon} _{\rm H}}\lpar {A^{ - 1}}{\rm \rho} \comma \; T\rpar \comma \;$$where A = 2.5 is the atomic weight of the mixture.
The initial density of the fuel is specified as ρ0 = 5ρs, 25ρs, and 100ρs, where ρs ≈ 0.22 g/cm3. The initial temperature whether equals 300 K while the initial pressure p 0 is determined from the EOS, or is chosen (together with p 0) on the isentrope passing through the point (ρs, p a), where p a = 0.1 MPa. Numerical results in both cases are very close to each other.
A free boundary with the pressure p a is initially at the point x = H. The monoenergetic beam of protons of kinetic energy 1 MeV acts upon this boundary during τpb = 50 ps. With the exception of the short time interval δτpb = 0.02τpb, the beam has the constant intensity J 0 = 1019 W/cm2:
$${J_b}\lpar t\rpar = \left\{{\matrix{ {{J_0}t/{\rm \delta} {{\rm \tau} _{\,pb}}\comma \; }&{t \le {\rm \delta} {{\rm \tau} _{\,pb}}\semicolon \; } \cr {{J_0}\comma \; } &{{\rm \delta} {{\rm \tau} _{\,pb}} \lt t \le {{\rm \tau} _{\,pb}}\semicolon \; } \cr {0\comma \; } &{t \gt {{\rm \tau} _{\,pb}}.} \cr } } \right.$$The less intensive beam of the same energy with J 0 = 1018 W/cm2 and τpb = 500 ps is also considered. At the point x = 0, we set the symmetry condition that is equivalent to action of the identical proton beam on a symmetrical layer of the fuel.
Only the primary fusion reaction between deuteron and triton with α-particle and neutron as the reaction products
$$ {\rm D} + {\rm T} \to {\rm \alpha} \lpar 3.5\,{\rm MeV}\rpar + {\rm N}\lpar 14\,{\rm MeV}\rpar \comma \;$$is taken into account. Neutrons are supposed to be escaped from the fuel without interaction.
In computing of all of the coefficients entering into the model, the plasma is supposed to be completely ionized.
A mathematical model is based on the equations of one-fluid two-temperature hydrodynamics. Electron and ion heat conduction, self-radiation of plasma and plasma heating by both the proton beam and α-particles are taken into account (Afanas'ev et al., Reference Afanas'ev, Gamalii and Rozanov1982):
$$\displaystyle{{{\rm d} \rho} \over {{\rm d}t}} = - {\rm \rho} \displaystyle{{\partial u} \over {\partial x}} \comma \;$$
$${\rm \rho} \displaystyle{{{\rm d}u} \over {{\rm d}t}} = - \displaystyle{{\partial p} \over {\partial x}}\comma \;$$
$${\rm \rho} \displaystyle{{{\rm d}{{\rm \varepsilon} _e}} \over {{\rm d}t}} = - {\,p_e}\displaystyle{{\partial u} \over {\partial x}} + \displaystyle{\partial \over {\partial x}}{{{\varkappa} }_e}\displaystyle{{\partial {T_e}} \over {\partial x}} + \displaystyle{3 \over 2}{n_i}{k_{\rm B}}\displaystyle{{{T_i} - {T_e}} \over {{{\rm \tau} _T}}} \!+\! {D_e} + {W_e} + R\comma \,$$
$${\rm \rho} \displaystyle{{{\rm d}{{\rm \varepsilon} _i}} \over {{\rm d}t}} = - {\,p_i}\displaystyle{{\partial u} \over {\partial x}} + \displaystyle{\partial \over {\partial x}}{{{\varkappa} }_i}\displaystyle{{\partial {T_i}} \over {\partial x}} + \displaystyle{3 \over 2}{n_i}{k_{\rm B}}\displaystyle{{{T_e} - {T_i}} \over {{{\rm \tau} _T}}} + {D_i} + {W_i}\comma \;$$where u is the mass velocity, d/dt = ∂/∂t + u∂/∂x is the Lagrangian derivative with respect to time, p e and p i are the electron and ion pressure, p = p e + p i is the total pressure, εe and εi are the electron and ion specific internal energy, T e and T i are the electron and ion temperature, ϰe and ϰi are the electron (Kalitkin & Kostomarov, Reference Kalitkin and Kostomarov2006) and ion (Silin, Reference Silin1971) heat conductivity. The third term in the right parts of Eqs. (2.5) and (2.6) defines the energy exchange between electrons and ions, n i = ρ/(Am u) is the ion number density, A is the atomic weight, m u is the atomic mass unit, k B is the Boltzmann constant, τT is the temperature relaxation time (Kalitkin & Kostomarov, Reference Kalitkin and Kostomarov2006). Apart from Kalitkin & Kostomarov (Reference Kalitkin and Kostomarov2006) and Silin (Reference Silin1971), for ϰe, ϰi and τT, we used formulas from Charakhch'yan et al. (2011). The rest terms in Eqs. (2.5) and (2.6) define the heating of electrons and ions by the proton beam (D e and D i) and by α-particles (W e and W i) as well as the energy exchange between electrons and self-radiation of plasma (R). The radiation pressure and the momentum transfer under deceleration of α-particles are neglected in the equation of motion (2.4). We also leave out of account terms describing the change of plasma composition induced by reaction (2.2) in the equation of discontinuity (2.3).
EOS for electrons in Eqs. (2.3)–(2.6) is taken in a form, which corresponds to thermal contribution of ideal electronic gas in completely ionized plasma of hydrogen isotopes (with taking into account degeneracy):
$${\,p_e}\lpar {\rm \rho} \comma \; T\rpar = \displaystyle{2 \over 3}{\rm \rho} {{\rm \varepsilon} _e}\lpar {\rm \rho} \comma \; T\rpar \comma \;$$
$${{\rm \varepsilon} _e}\lpar {\rm \rho} \comma \; T\rpar = \displaystyle{3 \over 2}{R_A}T\displaystyle{{{{\rm \beta} _e}T} \over {3{{\rm \rho} ^{2/3}} + {{\rm \beta} _e}T}}\comma \;$$where R A = k BN A(Am u)−1, N A is the Avogadro constant,
$$ {{\rm \beta} _e} = {\lpar {\rm \pi} /3\rpar ^{2/3}}\displaystyle{{{m_e}{k_{\rm B}}} \over {{{\rm \hbar }^2}}}{\lpar A{m_u}\rpar ^{2/3}}\comma \;$$m e is the electron mass, ħ is the Planck constant.
EOS for ions in Eqs. (2.3)–(2.6) is taken as follows,
$${\,p_i}\lpar {\rm \rho} \comma \; T\rpar = {\,p_{\rm H}}\lpar {A^{ - 1}}{\rm \rho} \comma \; T\rpar - {\,p_e}\lpar {\rm \rho} \comma \; T\rpar \comma \;$$
$${{\rm \varepsilon} _i}\lpar {\rm \rho} \comma \; T\rpar = {A^{ - 1}}{{\rm \varepsilon} _{\rm H}}\lpar {A^{ - 1}}{\rm \rho} \comma \; T\rpar - {{\rm \varepsilon} _e}\lpar {\rm \rho} \comma \; T\rpar.$$So, in the case of temperature equality T = T e = T i, the EOS (2.1) is satisfied, p e + p i = p(ρ, T), εe + εi = ε(ρ, T).
It should be stressed that, at high temperatures, T≫3βe−1ρ2/3, relations (2.7) and (2.8) have forms of pressure and internal energy of mono-particle ideal gas of Boltzmann,
$$ {\,p_e}\lpar {\rm \rho} \comma \; T\rpar = {\rm \rho} {R_A}T\comma \; \quad {{\rm \varepsilon} _e}\lpar {\rm \rho} \comma \; T\rpar = \displaystyle{3 \over 2}{R_A}T.$$The same relations take place at high temperatures for ionic components (2.9) and (2.10) in the EOS model from Khishchenko (Reference Khishchenko2008):
$$ {\,p_i}\lpar {\rm \rho} \comma \; T\rpar = {\rm \rho} {R_A}T\comma \; \quad {{\rm \varepsilon} _i}\lpar {\rm \rho} \comma \; T\rpar = \displaystyle{3 \over 2}{R_A}T.$$This implies in particular that for the problems in which the temperature difference of electrons and ions occurs at a high degree of plasma heating when the Boltzmann ideal gas approximation is applicable for all kinds of particles, and ionization is complete, EOS of hydrogen isotopes in Eqs. (2.3)–(2.6) can be taken as follows (Charakhch'yan et al., Reference Charakhch'yan, Gryn' and Khishchenko2011),
$$ {\,p_e} = \displaystyle{1 \over 2}p\lpar {\rm \rho} \comma \; {T_e}\rpar \comma \; \quad {{\rm \varepsilon} _e} = \displaystyle{1 \over 2}{\rm \varepsilon} \lpar {\rm \rho} \comma \; {T_e}\rpar \comma \; \quad {\,p_i} = \displaystyle{1 \over 2}p\lpar {\rm \rho} \comma \; {T_i}\rpar \comma \; \quad {{\rm \varepsilon} _i} = \displaystyle{1 \over 2}{\rm \varepsilon} \lpar {\rm \rho} \comma \; {T_i}\rpar \comma \;$$where functions of p and ε correspond to the one-temperature case (2.1). Such an approach to the separation of pressure and internal energy on the electronic and ionic parts was used previously (Khishchenko & Charakhch'yan, Reference Khishchenko and Charakhch'yan2013; Charakhch'yan et al., Reference Charakhch'yan, Gryn', Khishchenko and Fortov2013; Charakhch'yan & Khishchenko, Reference Charakhch'yan and Khishchenko2013) with the EOS model from Khishchenko (Reference Khishchenko2008). Calculations presented further with taking into account the degeneracy of electrons by Eqs. (2.7) and (2.8) did not show significant differences between the obtained results and the case of using Eqs. (2.11).
The number of events of reaction (2.2) per unit time and unit volume is as follows from Afanas'ev et al. (Reference Afanas'ev, Gamalii and Rozanov1982):
$$ F = {n_{\rm D}}{n_{\rm T}}{\left\langle {{\rm \sigma} v} \right\rangle _{{\rm DT}}}\comma \;$$where n D and n T are the deuteron and triton number density respectively, 〈σv〉DT is the ion-temperature dependence of the reaction rate averaged over Maxwellian distribution of ions (Brueckner & Jorna, Reference Brueckner and Jorna1973). The burnout of fuel nuclei is described by the equation
$$ \displaystyle{{{\rm d}{n_j}} \over {{\rm d}t}} = - {n_j}\displaystyle{{\partial u} \over {\partial x}} - F\comma \;$$where subscripts j = D and T correspond to the cases of deuterium and tritium.
Our model of α-particle heat considers the fuel burnout, but ignores the change of the density and EOS due to the change of the plasma composition. The computations for the model of local heat by α-particles performed previously (Khishchenko & Charakhch'yan, Reference Khishchenko and Charakhch'yan2013) show, that the difference between numerical results of simulations with and without taking into account the change of the plasma composition turns out insignificant.
Following Lindl (Reference Lindl1995); Basko (Reference Basko2009), we introduce the local burn-up factor
$${B_{{\rm loc}}} = \displaystyle{{{n_{\rm R}}} \over {{n_{\rm R}} + {n_{\rm D}}}}\comma \;$$where the number density of the fusion reaction events n R is defined by the equation
$$ \displaystyle{{{\rm d}{n_{\rm R}}} \over {{\rm d}t}} = - {n_{\rm R}}\displaystyle{{\partial u} \over {\partial x}} + F\comma \;$$and by the initial condition n R = 0 at t = 0. Excluding the derivative ∂u/∂x from (2.12) and (2.13), one can obtain the following ordinary differential equation along trajectories of Lagrangian particles
$$ \displaystyle{{{\rm d}{B_{{\rm loc}}}} \over {{\rm d}t}} = {{\rm \chi} _B}\lpar 1 - {B_{{\rm loc}}}\rpar \comma \;$$where the function χB = F/n D = n T〈σv〉DT is the burn-up rate. The solution of (2.14) satisfying the initial condition B loc = 0 at t = 0 has the form
$$ {B_{{\rm loc}}}\lpar s\comma \; t\rpar = 1 - \exp\left({ - \int\limits_0^t \, {{\rm \chi} _B}\lpar s\comma \; t^{\prime}\rpar {\rm d}t^{\prime}} \right)\comma \;$$where
$$ s\lpar x\comma \; t\rpar = \int\limits_0^x \, {\rm \rho} \lpar x^{\prime}\comma \; t\rpar {\rm d}x^{\prime}$$is the Lagrangian coordinate.
We use two models of α-particle heat. The first one is the track method (Brueckner & Jorna, Reference Brueckner and Jorna1973; Duderstadt & Moses, Reference Duderstadt and Moses1982) based on simple physical reasons applying to a discrete media. The domain of x is divided into cells by a numerical grid. The number of α-particles per unit cross-section area which are produced within a grid cell is divided into several angular groups basing on the uniform directional distribution of α-particles. The α-particles of one angular group move along the ray intersecting one or more grid cells and are decelerated in compliance with the equation
$$ v\displaystyle{{{\rm d}v} \over {{\rm d}{\rm \xi} }} = a\lpar x\comma \; v\rpar \comma \; \quad v\lpar {x_0}\rpar = {v_0}\comma \; $$where v is the α-particle velocity, v 0 ≈ 13 Mm/s is the initial α-particle velocity, a(x, v) is the deceleration (negative acceleration) of α-particles in plasma, x = x 0 + ξμ, x 0 is the center position of the grid cell where the group of α-particles is produced, ξ is the coordinate along the ray, μ is the cosine of the angle between the x-axis and the ray direction.
The deceleration a(x, v) = a e(T e(x), ρ(x), v) + a i(T i(x), ρ(x), v), where the first term in the right-hand part is related to electrons (Vygovskii et al., Reference Vygovskii, Il'in, Levkovskii, Rozanov and Sherman1990) as well as the second term—to ions (Sivukhin, Reference Sivukhin1964). The functions a e(T e, ρ, v) and a i(T i, ρ, v) from Gus'kov and Rozanov (Reference Gus'kov and Rozanov1982) were also used in our simulations and gave close results.
Solution of Eq. (2.16) together with condition v ≥ v th(x), where v th = (3k BT i/m α)1/2 is the velocity of α-particle thermalization, m α is the α-particle mass, enables one to determine the contributions of the group to the right-hand parts W e and W i of Eqs. (2.5) and (2.6) for the respective grid cells (see Appendix).
If plasma is located within a cylinder, the α-particles escaping from the cylinder do not contribute to the plasma heating. The 1D track method enables us to take into account approximately the above three-dimensional effect. We introduce a length parameter R α and consider the solution of Eq. (2.16) along the bounded interval 0 ≤ ξ ≤ ξmax, where
$$ {{\rm \xi} _{\rm max}} = {R_{\rm \alpha} }{\lpar 1 - {{\rm \mu} ^2}\rpar ^{ - 1/2}} $$corresponds to the intersection point of the ray and the lateral boundary of the cylinder of the radius R α, which symmetry axis coincides with the x-axis.
Let R α → 0. Then ξmax → 0 for all the interval −1 ≤ μ ≤ 1 with the exception of its extreme points μ = ±1. The right parts of Eqs. (2.5) and (2.6), which are integrals along μ, W e → 0, W i → 0, and the ignition is impossible. By increasing R α, one can determine its value, starting from which the burn wave arises, and interpret it as the beam radius necessary for the ignition.
The second model of α-particle heat is based on the kinetic steady-state equation for the distribution function f(x, v, μ) in Fokker–Plank approximation. If all of the produced α-particles have the same initial velocity v 0, and the diffusion of the function f(x, v, μ) in the velocity space can be ignored, the known Cauchy problem for the kinetic homogeneous equation arises (Gus'kov & Rozanov, Reference Gus'kov and Rozanov1982):
$$ {\rm \mu} v\displaystyle{{\partial f} \over {\partial x}} + \displaystyle{{\partial af} \over {\partial v}} = 0\comma \; \quad {v_{{\rm th}}}\lpar x\rpar \le v \le {v_0}\comma \; \quad f\lpar x\comma \; {v_0}\comma \; {\rm \mu} \rpar = - \displaystyle{{\tilde F\lpar x\comma \; {\rm \mu} \rpar } \over {a\lpar x\comma \; {v_0}\rpar }}\comma \; $$
where 
$\tilde F\lpar x\comma \; {\rm \mu}\rpar $ is the distribution function upon μ of the production rate of α-particles in a unit volume near the point x. Since
$$ \int\limits_{ - 1}^1 \, \tilde F\lpar x\comma \; {\rm \mu} \rpar {\rm d}{\rm \mu} = F\lpar x\rpar \comma \;$$
for isotropic distribution of produced α-particles, 
$\tilde F\lpar x\comma \; {\rm \mu}\rpar = F\lpar x\rpar /2$.
The right-hand parts of Eqs. (2.5) and (2.6) have the form
$${W_{e\comma i}}\lpar x\rpar = - {m_{\rm \alpha} } \int\limits_{ - 1}^1 \, \int\limits_{{v_{{\rm th}}}\lpar x\rpar }^{{v_0}} \, f\lpar x\comma \; v\comma \; {\rm \mu} \rpar {a_{e\comma i}}\lpar x\comma \; v\rpar v{\rm d}v{\rm d}{\rm \mu} .$$As shown in Appendix, a minor modification of the track method is a numerical method for the Cauchy problem (2.18). The quasi-1D model is similar to that described above for the track method.
The intensity of a monoenergetic proton beam is determined by the proton velocity v as J = n pvm pv 2/2, where n p is the proton number density, m p is the proton mass. Protons are supposed to be decelerated in the plasma according to Eq. (2.16) for μ = −1 and x 0 = x b. The value of n p is determined by given values of the initial proton velocity and the boundary beam intensity J b(t). Assuming that n p is independent of x, the function J(x) is determined by the solution of Eq. (2.16) v(x). The right-hand parts of Eqs. (2.5) and (2.6) are
$$ {D_e} = \displaystyle{{{a_e}} \over {{a_e} + {a_i}}}\displaystyle{{\partial J} \over {\partial x}}\comma \; \quad {D_i} = \displaystyle{{{a_i}} \over {{a_e} + {a_i}}}\displaystyle{{\partial J} \over {\partial x}}.$$Self-radiation of plasma is described by the steady-state transfer equation in the diffusion approximation by solid angle (Zel'dovich & Raizer, Reference Zel'dovich and Raizer1967). Similarly to Marchuk et al. (Reference Marchuk, Imshennik and Basko2009), we take into account the cooling of electrons by the inverse Compton effect using the known approximate formula (Zel'dovich, Reference Zel'dovich1975; Basko, Reference Basko2009). Resulting equations are as follows,
$$ \displaystyle{{\partial {q_{\rm \nu} }} \over {\partial x}} = {\rm \kappa} \lpar {B_{\rm \nu} }\lpar {T_e}\rpar - {u_{\rm \nu} }\rpar \comma \; \quad \displaystyle{{\partial {u_{\rm \nu} }} \over {\partial x}} = - 3{\rm \kappa} {q_{\rm \nu} }\comma \;$$where ν is the frequency, B ν(T e) is the Planck function, κ = κ(ρ, T e, ν) is the absorption coefficient with accounting for the induced emission. The term in the right-hand part of Eq. (2.5) has the form
$$ R = - \displaystyle{{\partial Q} \over {\partial x}} - \displaystyle{{4{{\rm \sigma} _{\rm T}}{n_e}U} \over {{m_e}{c^2}}}{k_{\rm B}}\lpar {T_e} - {T_r}\rpar \comma \;$$where
$$ Q = \int\limits_0^\infty \, {q_{\rm \nu} }{\rm d}{\rm \nu} \comma \; \quad U = \int\limits_0^\infty \, {u_{\rm \nu} }{\rm d}{\rm \nu} \comma \;$$n e = zn i is the electron number density, m e is the electron mass, c is the speed of light, σT is the Thomson scattering cross-section of photons by free electrons, T r is the photon temperature determined by the equality
$$ \int\limits_0^\infty \, {B_{\rm \nu} }\lpar {T_r}\rpar {\rm d}{\rm \nu} = U.$$Numerical method for the simulations is based on splitting into physical processes. The Godunov first order accurate method in Lagrangian variables (Godunov, Reference Godunov1959) is used for the hydrodynamics equations. For the heat conduction and the energy exchange between electrons and ions, the implicit over time method is used. The uniform Lagrangian grid contains from 350 to 700 nodes. Number of angular groups in the track method varies from 12 to 24. The results of calculations vary insignificantly when changing to double the number of grid points and angular groups.
The grid on the frequency ν for solving of Eqs. (2.19) occupies the range from 4 to 8 decimal exponents. Number of grid nodes per a decimal exponent varies from 5 to 20. The integration with respect to the frequency in Eq. (2.20) is performed by the trapezium method.
3. IGNITION
Let us consider the initial stage of the target ignition, which is limited to the proton beam duration, t ≤ τpb. Suppose the heating by the proton beam is so fast that the motion of the fuel and its density change can be ignored. Such kind of heating will be referred to as isochoric. If the thermonuclear reaction and the plasma self-radiation are also ignored, the sum of Eqs. (2.5) and (2.6) takes the form
$$ {{\rm \rho} _0}\displaystyle{{\partial {\rm \varepsilon} \lpar x\comma \; t\rpar } \over {\partial t}} = \displaystyle{{\partial \lpar J + {q_e} + {q_i}\rpar } \over {\partial x}}\comma \; \quad {q_e} = {{{ \varkappa} }_e}\displaystyle{{\partial {T_e}} \over {\partial x}}\comma \; \quad {q_i} = {{{\varkappa} }_i}\displaystyle{{\partial {T_i}} \over {\partial x}}\comma \;$$ε = εe + εi. Integrating this equation with respect to time t from 0 to τpb and with respect to the spatial coordinate x from H − l p to H, where l p is the mean free-path length of protons, setting the heat fluxes q e = q i = 0 at the boundaries of the integration with respect to x and ignoring the initial energy of the fuel, obtain
$$ \bar {\rm \varepsilon} = \displaystyle{1 \over {{l_p}}} \int\limits_{H - {l_p}}^H \, {\rm \varepsilon} \lpar x\comma \; {{\rm \tau} _{\,pb}}\rpar {\rm d}x = \displaystyle{I \over {{{\rm \rho} _0}{l_p}}}\comma \; \quad I = \int\limits_0^{{{\rm \tau} _{\,pb}}} \, J\lpar t^{\prime}\rpar {\rm d}t^{\prime}\comma \;$$
where 
$\bar {\rm \varepsilon}$ is the average specific internal energy of plasma in the heating region at t = τpb, I is the beam energy per unit cross-sectional area. Since the dependence of the proton deceleration in plasma a on the density is close to the linear one, the free-path length of protons l p ~ ρ0−1. Therefore the averaged internal energy 
$\bar {\rm \varepsilon}$ is independent of ρ0 and is determined by only the beam energy I and by the initial proton energy (1 MeV) which determines the free-path length.
Epithermal protons give up the greater part of their energy to electrons, while ions are heated by both epithermal protons and more hot electrons. For the given beam energy I, the ion temperature at t = τpb increases both with increasing τpb and with increasing ρ0 because the relaxation time of the electron and ion temperatures τT ~ ρ−1. If the electron temperature at t = τpb is much greater than the ion temperature, the latter can be sufficient for the ignition after stopping the action of the proton beam at t > τpb.
The isochoric heating violation is connected with the rarefaction wave moving at the speed of sound from the free boundary into the target, and with the shock wave, arising due to the rapid growth of pressure in the heating region. As the condition of the isochoric heating violation, we take the inequality csτpb > l p, where c s is the speed of sound. Since l p ~ ρ0−1, with increasing ρ0 one should decrease the beam duration τpb to keep close to the isochoric heating, which, as it follows from Eq. (3.1), needs the appropriate increase in the beam intensity J 0 to conserve the average internal energy in the heating region.
The ion temperature and the mass velocity profiles by the Lagrangian variable with respect to the free boundary s − s b at the time t = τpb for the beams of the same energy and different durations τpb = 50 and 100 ps, as well as for the two values of the initial density ρ0 = 25ρs and 100ρs are presented in Figure 2. The parameter R α limiting the trajectory of α-particles is selected near the ignition boundary of the target by at least one of the beams (R α = 0.4 and 0.1 mm for ρ0 = 25ρs and 100ρs respectively).
Fig. 2. The ion temperature (a, b) and the mass velocity (c, d) as functions of the Lagrangian variable with respect to the free boundary s − s b at t = τpb: ρ0 = 25ρs, R α = 0.4 mm (a, c) and ρ0 = 100ρs, R α = 0.1 mm (b, d) for different beams of the same energy, J 0 = 1019 W/cm2, τpb = 50 ps (solid lines) and J 0 = 1018 W/cm2, τpb = 500 ps (dashed lines).
First, we examine the case of ρ0 = 25ρs (see Figs. 2a and 2c).
At τpb = 50 ps, as follows from the corresponding velocity profile (the solid line in Fig. 2c), the heating is almost isochoric: the formation of the shock wave has not started yet, and the rarefaction wave occupies only a small part of the heating region. Although the ion temperature (the solid line in Fig. 2a) is relatively small, its subsequent growth due to heating by more hot electrons leads to the target ignition.
With increasing the duration τpb up to 500 ps, the proton beam heating is no longer isochoric. As seen in the velocity profile (the dashed line in Fig. 2c), a continuous compression wave is formed inside the heating region. Then this wave should transform into a shock wave. A rarefaction wave adjoins the compression one, which in the absence of heat by α-particles should lead to a rapid drop in the amplitude of the shock wave. However, the ion temperature (dashed line in Fig. 2a) is sufficient to ignite the target. Analysis of the calculation results shows that a significant heating of the plasma by α-particles takes place in the field of the compression wave in Figure 2c. As a result, the compression wave is rapidly converted to the detonation wave of the well-known type (Landau & Lifshitz, Reference Landau and Lifshitz1987) with a rarefaction wave, which is adjacent to the detonation wave front.
We now turn to the case ρ0 = 100ρs (see Figs. 2b and 2d). When τpb = 50 ps, the ion temperature (the solid line in Fig. 2b) is much greater than in the case of the same beam and ρ0 = 25ρs due to the decrease of the temperature relaxation time τT ~ ρ−1. As can be seen from the velocity profile (the solid line in Fig. 2d), the shock wave begins to form, and the rarefaction wave takes until a small portion of the heating region. The ignition mechanism for these parameters will be discussed below.
For the beam with τpb = 500 ps, the target does not ignite. As can be seen from the velocity profile (dashed line in Fig. 2d), the shock wave at t = τpb is too far away from the heating region and no longer can transform into a detonation wave. After some time, the rarefaction wave overtakes the shock front and starts to reduce its amplitude. Calculations show that for the fuel density ρ0 = 103ρs, a similar flow pattern without the target ignition occurs for the beam with τpb = 50 ps.
The aforementioned lack of the ignition by the beam with the intensity of 1018 W/cm2 for the fuel density ρ0 = 100ρs ≈ 22 g/cm3, and the ignition by the beam with the intensity of 1019 W/cm2 for the same density, as well as the lack of the ignition by this beam for the fuel density ρ0 ≈ 220 g/cm3, correspond to the density dependence of the minimum beam intensity required for ignition (Atzeni, Reference Atzeni1999) based on results of two-dimensional calculations, which gives a value of approximately 6 × 1018 W/cm2 for ρ0 = 22 g/cm3 and about 5 × 1019 W/cm2 for ρ0 = 220 g/cm3.
Results of calculations for a target with ρ0 = 100ρs, R α = 0.1 mm, and for the beam with τpb = 50 ps, J 0 = 1019 W/cm2 are presented below. Consider first the case of the fuel layer half-width H = 0.5 mm, which corresponds to the value of the parameter Hρ0 ≈ 1 g/cm2.
Figure 3 shows the ignition process of the target. The functions T e(x), T i(x), ρ(x), u(x) and the burn-up rate χB(x), which determines the local burn-up factor B loc by Eq. (2.15), are given at the termination of the beam 50 ps and three subsequent times 100, 150 and 200 ps. A small discontinuity in the function χB(x) is caused by the discontinuity in the approximation formula for 〈σv〉DT(T i) from Brueckner and Jorna (Reference Brueckner and Jorna1973).
Fig. 3. The temperature (a, the solid lines correspond to ions, the dashed lines—to electrons), the burn-up rate (b), the density (c) and the mass velocity (d) spatial profiles in the formation of the detonation wave at t = 50 (1), 100 (2), 150 (3), and 200 ps (4) for the target with ρ0 = 100ρs, R α = 0.1 mm, H = 0.5 mm and for the beam with τpb = 50 ps.
At t = 50 ps, the maximal ion temperature is about 100 MK, which is sufficient for the ignition. The shock wave has not yet formed, and the burning wave, which can be identified with behavior of the function χB(x), is continuous. At t = 100 ps, the shock wave is formed, and the burning wave is a little behind it and still continuous. At t = 150 ps, the burning wave has caught up with the shock wave and transformed it into the detonation one. Comparing functions T i(x) and χB(x) at t = 150 and 200 ps, one can see that the temperature and the burn-up rate at the detonation wave front quickly grow with time. Before the detonation wave one can see a precursor caused by electron heat conduction and α-particles (near the wave front) and by the self-radiation of high-temperature plasma (at large distances from the front).
To conclude this section, we present the numerical results for the target with the initial density ρ0 = 100ρs at different values of R α, those limiting trajectory of α-particles, and for different models of α-particle transport: the Fokker–Plank equation and the track method. Table 1 shows the values of the burn-up factor for targets with two values of H = 0.5 (B 0) and 2.5 mm (B 1). These values correspond to a sufficiently late point in time, when their significant growth ceases. Zero values B 0 and B 1 means that their value is less than 0.01, and are interpreted as a lack of ignition.
Table 1. The burn-up factors for the targets with ρ0 = 100ρs, H = 0.5 (B 0) and 2.5 mm (B 1), for different values of the parameter R α and different models of α-particle transport: the Fokker–Plank equation (FP) and the track method (TM)
At R α = 0.15 mm, the factors B 0 and B 1 are the same for both models, and slightly differ from respective values at R α = 0.1 mm. At R α = 0.05 mm, the Fokker–Plank equation still leads to the ignition of the target with somewhat smaller values of B 0 and B 1 compared with the case of R α = 0.1 mm, and in the case of the track method B 1 = 0. This means that the detonation wave is not formed for H = 2.5 mm since the shock wave at large distances is extinguished by the rarefaction wave. If the symmetry plane is located close enough to the edge of the target (H = 0.5 mm), the shock wave at the moment of reflection is strong enough to ignite the target, creating a reflected detonation wave. At R α = 0.03 mm, the ignition is possible only if H = 0.5 mm with small values of B 0. At R α = 0.01 mm the ignition is absent. To estimate the energy of ignition, we set R α = 0.1 mm, wherein the reliable ignition occurs in both models.
4. REFLECTION OF DETONATION WAVE
In the case of H = 0.5 mm, reflection of the detonation wave from the symmetry plane x = 0 is shown in Figure 4. The functions T i(x), ρ(x), u(x) and the local burn-up factor B loc(x) are presented in the six time points. The first two points (200 and 250 ps) meet the detonation wave before its reflection from the symmetry plane. Note the relatively small values of the local burn-up factor (about 0.05). Pay attention to the velocity profile at t = 250 ps (curve 2 in Fig. 4d). One can see that the precursor of the detonation wave transforms into a flow with a linear velocity profile between the wave and the symmetry plane. Such flows, which are characterized by the dependence
Fig. 4. The spatial profiles of the ion temperature (a), the local burn-up factor (b), the density (c) and the mass velocity (d) in reflection of the detonation wave at points of time 200 (1), 250 (2), 300 (3), 350 (4), 400 (5), and 450 ps (6) for the parameters of the problem, as shown in Figure 3.
$$ u\lpar x\comma \; t\rpar = {\rm \varphi} \lpar t\rpar x\comma \;$$are well known in hydrodynamics (Sedov, Reference Sedov1972) and combined into a large family of solutions, each member of that is determined by an arbitrary function of one argument and three arbitrary constants. For all of these solutions, the following formula takes place
$$ {\rm \rho} \lpar s\comma \; t\rpar = {\rm \rho} \lpar s\comma \; {t_1}\rpar \exp\left({ - \int\limits_{{t_1}}^t \, {\rm \varphi} \lpar t^{\prime}\rpar {\rm d}t^{\prime}} \right)\comma \;$$which is obtained by integrating an ordinary differential equation arising under the substitution (4.1) to the equation of discontinuity (2.3). In the case under consideration, since ϕ(t) < 0, the density in the region of the precursor and, in particular, in the symmetry plane, grows with time.
The last four points of time in Figure 4 (300, 350, 400, and 450 ps) give the flow pattern after reflection of the detonation wave. One can see that the flow with a linear velocity profile arises again between the reflected detonation wave and the symmetry plane. Its special feature is the proximity of the thermodynamic functions to the constant upon x value, which depends on time. The local burn-up factor increases significantly after reflection of the detonation wave and continues to grow because the burn-up rate χB(x,t) decreases with time slowly enough, remaining almost constant function of x. Note that the spatial homogeneity of the thermodynamic functions after the reflected wave front is not typical for the solutions of the equations of hydrodynamics in the case of spherical or cylindrical geometry. For example, a well-known solution on the problem of a converging spherical shock wave (Guderley, Reference Guderley1942; Stanyukovich, Reference Stanyukovich1955), extended through time after the shock wave collapse, has zero density at the point of symmetry (Goldman, Reference Goldman1973).
Solution of the equations of hydrodynamics with a linear velocity profile and constant on x values of thermodynamic functions is contained in the above-mentioned family from Sedov (Reference Sedov1972). It also can be obtained directly from the equation of motion (2.4). Putting in Eq. (2.4) ∂p/∂x = 0 and substituting the velocity in the form (4.1), we obtain an ordinary differential equation 
$\dot{\rm \varphi} + {\rm \varphi}^2 = 0$, whose solution has the form
$$ {\rm \varphi} \lpar t\rpar = \displaystyle{1 \over {C + t}}\comma \;$$where C is an arbitrary constant. Substituting (4.3) to Eq. (4.2), we obtain
$$ {\rm \rho} \lpar t\rpar = \displaystyle{{{{\rm \rho} _0}C} \over {C + t}} = {{\rm \rho} _0}C{\rm \varphi} \lpar t\rpar \comma \;$$where constant ρ0 is chosen so that ρ(0) = ρ0.
Since the density ρ is independent of the spatial coordinate x, the velocity has a linear profile along the Lagrangian coordinate s = ρx, the slope of which, as follows from Eqs. (4.1), (4.3), and (4.4), u/s = (ρ0C)−1, is independent of time. Shown in Figure 5 velocity profiles along the Lagrangian coordinate indicate that this property is satisfied with good accuracy for the problem in the region between the symmetry plane and the reflected detonation wave.
Fig. 5. Mass velocity profiles along the Lagrangian coordinate after reflection of the detonation wave for the same parameters of the problem and the time points as in Figure 4.
Using the simulation results, one can determine the approximate value of constant C in Eq. (4.3). It is convenient to substitute time 
$\bar t = t - t_{\ast}$ for t in Eqs. (4.3) and (4.4), where 
$t_{\ast}$ is the time point of appearance of the reflected detonation wave. Note that the density on the symmetry plane as a function of time ρ(0, t) increases as the incident detonation wave approaches to the symmetry plane, and begins to decrease in accordance with Eq. (4.4) after appearance of the reflected wave. It is therefore natural to define the time of formation of the reflected wave 
$t_{\ast}$ and constant ρ0 in Eq. (4.4) by the condition maxt(ρ(0, t)) = ρ(0, 
$t_{\ast})={\rm \rho}_0$. The result is ρ0 ≈ 80 g/cm3, 
$t_{\ast}$ ≈ 270 ps.
The calculation results are compared with Eqs. (4.3) and (4.4) for the four time points shown in Figure 4 after reflection of the detonation wave, t = 300, 350, 400, and 450 ps. As the slope 
${\rm \varphi}$ and density ρ for each time, the values at the point x approximately two times less than the coordinate of the front of the reflected detonation wave are selected. We denote by 
$\bar{\rm \varphi}$ and 
$\bar {\rm \rho}$ the corresponding values included in Eqs. (4.3) and (4.4). At first, we put 
$\bar{\rm \varphi} = {\rm \varphi}$ and calculate C from Eq. (4.3), then 
$\bar {\rm \rho}$ by Eq. (4.4) and the relative error of calculating the density 
${\rm \delta \rho} = \vert 1 - \bar {\rm \rho}/{\rm \rho}\vert $. The results of calculations are shown in the first three columns of Table 2. For the first three values of t, constant C varies within 10% and the error δρ is within 6%. At the last time point t = 450 ps, constant C changes by about 10% more, and the error δρ grows to 14%. Nevertheless, it is possible to select a single value of constant C for all four points of time and get a relatively small magnitude errors 
${\rm \delta\varphi} = \vert 1 - \bar {\rm \varphi}/{\rm \varphi}\vert $ and δρ. Such a possibility is demonstrated in the last two columns of Table 2, where the values of these errors are presented for C = 40 ps. It can be seen that both errors are within 3%.
Table 2. Comparison of simulated values of φ and ρ with the analytical Eqs. (4.3) and (4.4): the value of C from Eq. (4.3) as well as the relative error δρ for the given value φ and the relative errors δρ40 and δφ40 for C = 40 ps at some points of time t
Let us consider now the target with the layer half-thickness H = 2.5 mm, which corresponds to the value of the parameter Hρ0 ≈ 5 g/cm2. Figure 6 shows the ion temperature, the density, the mass velocity, and the local burn-up factor as functions of x at eight points of time. The first four points give a picture of the flow before the detonation wave reflection from the symmetry plane. Detonation wave intensity increases with time, as can be seen from the profiles of the ion temperature and the mass velocity, reaching values of the order of 700 MK and 2 Mm/s at the instant of reflection. Such a high velocity in the incident detonation wave at the reflection leads to the ion temperature after reflection of about 1.3 GK (see curve 5 in Fig. 6a). Further, as in the above case H = 0.5 mm, the flow with a linear velocity profile and the weak dependence of the thermodynamic functions upon x arises between the reflected detonation wave and the symmetry plane.
Fig. 6. Spatial profiles of the ion temperature (a), the local burn-up factor (b), the density (c) an the mass velocity (d) at the time points 0.2 (1), 0.4 (2), 0.6 (3), 0.8 (4), 1 (5), 1.2 (6), 1.4 (7) and 1.8 ns (8), corresponding to the incident (1–4) and the reflected (5–8) detonation waves for H = 2.5 mm (another parameters of the problem are the same as in Fig. 3).
Note the rapid rise of the local burn-up factor after reflection of the detonation wave as close to the symmetry plane as in the periphery. At the last instant shown in Figure 6, the local burn-up factor in the region 0 ≤ x ≤ 3.5 mm, containing about 80% of the fuel mass, varies from about 0.7 to 0.6.
The burn-up factor B(t) is presented in Figure 7 for the two above values of H. The point of reflection of the detonation wave t = t * is indicated. One can see that the expansion stage t > t * gives the main contribution to the final value of the factor in the case of H = 0.5 mm. For H = 2.5 mm the burn-up factor increases from 0.4 to 0.7 at the expansion stage.
Fig. 7. The time dependence of the burn-up factor for H = 0.5 (1) and 2.5 mm (2) (another parameters of the problem are the same as in Fig. 3); asterisks denote the instant of the detonation wave reflection.
5. INTEGRATED CHARACTERISTICS
Generalized characteristics of the target for three values of the initial density ρ0 = 5ρs, 25ρs, and 100ρs are shown in Table 3. For ρ0 = 100ρs, the characteristics are obtained by the results of the above calculations with the parameter R α = 0.1 mm. For other values of ρ0, the parameter R α is chosen from the condition ρ0R α = const. For each value of ρ0 two values of the layer half-width H, corresponding to the two given values of the parameter Hρ0 ≈ 1 and 5 g/cm2, are chosen. Interpreting the parameter R α as the radius of the cylindrical target, shown in Figure 1, one can define the target ignition energy E ig = πR α2I(∞) and the mass M = 2HπR α2ρ0. These values are given in the table together with the values of the gain with respect to the energy of neutrons and the burn-up factor after the burning process.
Table 3. The ignition energy for one proton beam E ig, the target mass M, the burn-up factor B and the gain G for three values of ρ0, the cylinder radius R α~ρ0−1 and different values of the layer half-width H determined by two given values of the parameter Hρ0
As in the case of the well-known approximate formula for the expansion of a spherical target (Basko, Reference Basko2009), the burn-up factor B, as well as the gain G, are determined by the parameter Hρ0. If Hρ0 ≈ 1 g/cm2, B ≈ 0.3, G ≈ 200. If Hρ0 ≈ 5 g/cm2, the burn-up factor increases by more than a factor of two, which, together with the increase in the fuel mass five times, gives the gain G > 2000.
Since the parameter R α, determined by the condition R αρ0 = const, provides the target ignition, the ignition energy decreases with increasing ρ0 under the law E ig ~ ρ0−2 in accordance with the known theoretical estimate (Tabak et al., Reference Tabak, Hammer, Glinsky, Kruer, Wilks, Woodworth, Campbell, Perry and Mason1994). For ρ0 = 100ρs ≈ 22 g/cm3, the ignition energy E ig = 160 kJ.
6. CONCLUSIONS
Computational experiments with two proton beams of the same energy and different intensities show that the beam intensity required for the ignition increases with the initial fuel density in accordance with a known formula, obtained by results of two-dimensional calculations.
If the target is ignited, a detonation wave, which is adjacent to the front of the rarefaction wave, occurs. The formation of the detonation wave precedes a short time interval with a shock wave that goes ahead of a slow combustion wave. When the combustion front overtakes the shock wave, the latter transforms into the detonation wave.
There is no ignition if the distance between the shock front and the region of the fuel heating by protons is too large at the end of irradiation. Over time, the rarefaction wave overtakes the shock front and reduces the intensity of the shock wave.
For small values of R α, limiting trajectories of α-particles, and small target thickness (Hρ0 ≈ 1 g/cm2), the target ignition is possible after the shock wave reflection from the symmetry plane. With increasing H, starting from a certain value, the target ceases to ignite due to the aforementioned decrease in the intensity of the shock wave that is overtaken by the rarefaction wave.
Upon reflection of the detonation wave from the symmetry plane, the flow with the linear velocity profile along the spatial variable x and with a weak dependence of the thermodynamic functions upon x, close to one of the representatives of a certain family of solutions of hydrodynamic equations, occurs between this plane and the reflected detonation wave. Combustion efficiency at this stage is largely due to the spatial homogeneity of the burn-up rate (the product of the DT reaction rate and the triton number density) whose integral along the trajectory of a Lagrangian particle defines the local burn-up factor.
Solution of the equations of hydrodynamics with a linear velocity profile and constant in x values of thermodynamic functions is found analytically up to an arbitrary constant. The possibility to choose this constant so that the solution describes with good accuracy the change in the density with time after reflection of the detonation wave is shown.
For the parameter Hρ0 ≈ 1 g/cm2, the main contribution to the final value of the burn-up factor B ≈ 0.3 is given by the expansion stage after reflection of the detonation wave. For Hρ0 ≈ 5 g/cm2, the burn-up factor is increased by the expansion stage from 0.4 to 0.7. At Hρ0 ≈ 1 g/cm2, the gain with respect to the energy of neutrons G ≈ 200. At Hρ0 ≈ 5 g/cm2, the burn-up factor increases by more than a factor of two, which, together with the increase in the fuel mass five times, gives the gain G > 2000. The above values of the burn-up factor and the gain should be considered as the maximum possible values for the real cylindrical target of the same initial density and length.
The quasi-1D model, limiting trajectories of α-particles by a cylinder of a given radius, reproduces known theoretical dependence of the ignition energy E ig ~ ρ0−2. For ρ0 = 100ρs ≈ 22 g/cm3, the quasi-1D model gives E ig = 160 kJ.
ACKNOWLEDGMENTS
The work is supported by grants of the Russian Foundation for Basic Research (12-01-00130 and 14-08-00967) and the President of the Russian Federation (NSh-6614.2014.2), as well as by the programs of scientific investigations of the Russian Academy of Sciences (program 3 of the Division of Mathematics RAS and program 2P of the Presidium RAS).
Appendix: TRACK METHOD OF PLASMA HEATING BY CHARGED PRODUCTS OF A THERMONUCLEAR REACTION AND KINETIC FOKKER–PLANK EQUATION
We limit our consideration by plane 1D flows for which all of the plasma parameters depend upon only one space variable x. Without loss of generality, the domain 0 ≤ x ≤ 1 can be considered. The charged products of a thermonuclear reaction will be referred to as particles. The deceleration a(x, v) < 0, the number of the reaction events per unit time and unit volume F(x) and the velocity of particle thermalization v th(x) are supposed to be given.
To simplify calculations, at first we consider the problem without particles entering the domain. At the end of the section, we present the modification of the method for the condition of completely elastic reflection of particles from one of the boundaries, which simulates the symmetry condition at this boundary.
Consider the track method (Brueckner & Jorna, Reference Brueckner and Jorna1973; Duderstadt & Moses, Reference Duderstadt and Moses1982) based on simple physical reasons applied to a discrete media. F(x 0) particles with the same magnitude of velocity v 0 are produced in unit volume per unit time at the point x 0. The particles move in straight rays in all directions, which of them is given by the cosine μ of the angle between the ray direction and the axis x, −1 ≤ μ ≤ 1. Consider the particle velocity v(ξ, x 0, μ) determined by the Cauchy problem
$$ v\displaystyle{{\partial v} \over {\partial {\rm \xi} }} = a\lpar x\comma \; v\rpar \comma \; \quad x = {x_0} + {\rm \xi} {\rm \mu} \comma \; \quad v\lpar 0\comma \; {x_0}\comma \; {\rm \mu} \rpar = {v_0}\comma \; $$where ξ is the coordinate along the ray. Since a < 0, the particle velocity decreases. The maximal value of ξ for given x 0 and μ is determined by the conditions
$$ 0 \le {x_0} + {{\rm \xi} _{\max}}{\rm \mu} \le 1\comma \; \quad v\lpar {{\rm \xi} _{\max}}\comma \; {x_0}\comma \; {\rm \mu} \rpar = {v_{{\rm th}}}\lpar {x_0} + {{\rm \xi} _{\max}}{\rm \mu} \rpar .$$
Let the grid x l, l = 1, 2, …, N + 1, x 1 = 0, x N+1 = 1 divides the coordinate x into the cells of the size Δx l = x l+1 − x l, l = 1, 2, …, N. For definiteness, assume the plasma parameters are defined at the centers of the cells 
$\mathop {\bar x}\nolimits_l = 0.5\lpar {x_l} + {x_{l + 1}}\rpar $.
Number of particles produced in the cell j with unit cross-section area per unit time, 
${\rm \sigma}_j = F\lpar \mathop{\bar x}\nolimits_j\rpar \Delta x_j$. The particles are divided into M angular groups with the constant step of μ, Δμ = 2/M, basing on the isotropic angular distribution of particles (σjΔμ/2 particles in the group). All particles of one group m produced in the cell j move along the straight ray with μ = μm = −1 + (m − 0.5)Δμ and are decelerated in compliance with Eq. (A1) with 
$x_0 = \mathop{\bar x}\nolimits_j$. We restrict our consideration by the case of even values of M, for which μm ≠ 0 at any m.
Let the particles of a certain angular group produced in the cell j traverse the cell k of dimension Δx k. If k = j, the boundary condition for Eq. (A1) is given at the cell center, v = v 0. Otherwise, the boundary condition is given at one of the points x k or x k+1 with a smaller value of ξ, v = v b0, where the velocity v b0 has been obtained from computation of the previous cell. In the most simple variant of the method, which is considered in the present paper, the function a(x, v) within the cell supposes to be of constant value 
$a_k=a\lpar \mathop {\bar x}\nolimits_k\comma \; v_b\rpar $, where v b = v b0 at k ≠ j and v b = v b0 at k = j. The velocity v e at the point with the greater value of ξ is determined by the Cauchy problem (A1) and the thermalization condition 
${v_e} \ge v_{{\rm th}}^k = {v_{{\rm th}}}\lpar \mathop {\bar x}\nolimits_k \rpar $
$$ \mathop {\tilde v}\nolimits_e^2 = v_b^2 + 2a_k\Delta {\rm \xi} \comma \; \quad {v_e} = \max\lpar v_{{\rm th}}^k \comma \; \mathop {\tilde v}\nolimits_e \!\rpar \comma \;$$
where Δξ = Δξk = Δx k/|μm| at k ≠ j and Δξ = Δξk/2 at k = j. If particles are thermalized 
$\mathop {\tilde v}\nolimits_e \le v_{{\rm th}}^k $, computation of the given angular group is terminated.
The energy transferred by particles of the angular group m produced in the cell j to unit volume of the cell k in unit time
$$ {{\rm \psi} _{kjm}} = \displaystyle{{{m_p}{{\rm \sigma} _j}\Delta {\rm \mu} } \over {2\Delta {x_k}}}\displaystyle{{v_b^2 - v_e^2 } \over 2} =$$
$$ = \displaystyle{{{m_p}{{\rm \sigma} _j}\Delta {\rm \mu} } \over {4\vert {{\rm \mu} _m}\vert\Delta{\rm\xi}_k }}\left\{{\matrix{ {\min \lpar{-}2\Delta {{\rm \xi} _k}\comma \; v_{b0}^2 - {{\lpar v_{{\rm th}}^k \rpar }^2} \comma \; }&{k \neq j } \cr {\min\lpar {-}a_k\Delta {{\rm \xi} _k}\comma \; v_{0}^2 - {{\lpar v_{{\rm th}}^k \rpar }^2} \comma \; }&{k = j } \cr } } \right.$$where m p is the particle mass. The total energy transferred to unit volume of the cell k in unit time
$$ {w_k} = \sum\limits_{m = 1}^M \, \sum\limits_{\,j = k}^{k + \Delta k} \, {{\rm \psi} _{kjm}}\comma \;$$where k + Δk is the number of the cell in which the particles of group m, those are thermalized in the cell k, were produced. The sign of Δk depends of the sign of μm: Δk ≤ 0 at μm > 0 and Δk ≥ 0 at μm < 0. The evident restriction on Δk is 1 ≤ k + Δk ≤ N. To compute the sum (A4), the terms ψkjm are computed step-by-step for all of the cells 1 ≤ j ≤ N and the angular groups 1 ≤ m ≤ M.
To study connection between the track method and the kinetic steady-state equation in Fokker–Plank approximation (Gus'kov & Rozanov, Reference Gus'kov and Rozanov1982), it is necessary to formulate the track method in terms of the mathematical physics rather than algebraic equations, i.e. to find the limit of w k at h = max(Δx l), Δμ → 0.
Let, to begin with h → 0. Then, the number of terms in the internal sum of Eq. (A4) tends to infinity. Therefore, one can reject the first term with j = k and the last term with j = k + Δk. As a result, the internal sum in Eq. (A4) takes the form
$$ {\Psi _{km}} = - \displaystyle{{{m_p}\Delta {\rm \mu} } \over {2\vert {{\rm \mu} _m}\vert }}\sum\limits_{\,j = k + {\rm sign}\lpar \Delta k\rpar }^{k + \Delta k - {\rm sign}\lpar \Delta k\rpar } \, a\lpar \mathop {\bar x}\nolimits_k\!\comma \; {v_{b0}}\rpar F\lpar \mathop {\bar x}\nolimits_j \rpar \Delta {x_j}.$$ Denote 
$x = \lim_{h \to 0}\mathop {\bar x}\nolimits_k $, 
${x_0} = \lim_{h \to 0}\mathop {\bar x}\nolimits_j$. As long as limh→0x k= limh→0 x k+1 = x, limh→0v b0 = v(ξ, x 0, μm), ξ = (x − x 0)/μm, where v(ξ, x 0, μm) is the solution of the Cauchy problem (A1). As a result, obtain
$$ \mathop {\lim}\limits_{h \to 0} {\Psi _{km}} = \Psi \lpar x\comma \; {{\rm \mu} _m}\rpar =$$
$$ = \displaystyle{{{m_p}\Delta {\rm \mu} } \over {2{{\rm \mu} _m}}} \int\limits_x^{x + \Delta x} \, a\lpar v\lpar {\rm \xi} \comma \; {x_0}\comma \; {{\rm \mu} _m}\rpar \comma \; x\rpar F\lpar {x_0}\rpar {\rm d}{x_0}\comma \; $$where Δx is determined by thermalization condition at the point x for particles of group m produced at the point x + Δx:
$$ v\lpar \!- \Delta x/{{\rm \mu} _m}\comma \; x + \Delta x\comma \; {{\rm \mu} _m}\rpar = {v_{{\rm th}}}\lpar x\rpar .$$Disappearance of the modulus sign for μm under conversion from the sum (A5) to the integral (A6) is due to that the sum (A5) conserves its sign at changing the sign of Δk, while the integral (A6) changes its sign at changing the sign of Δx.
At μm → 0, the function Ψ(x, μm) remains finite though μm is in the denominator of the function. We limit our proof by the simplest case a = const. Then the condition (A7) and the Cauchy problem (A1) give Δx = μm(v 02 − v th2(x))/(2a), from which
$$ \Psi \lpar x\comma \; 0\rpar = \displaystyle{{{m_p}\Delta {\rm \mu} F\lpar x\rpar \lpar v_0^2 - v_{{\rm th}}^2 \lpar x\rpar \rpar } \over 4}.$$Finally, the required limit, which is the rate of energy transfer to unit volume of plasma at the point x, is as follows:
$$ W\lpar x\rpar = \mathop {\lim}\limits_{h\comma \Delta {\rm \mu} \to 0} {w_k} =$$
$$ = \displaystyle{{{m_p}} \over 2} \int\limits_{ - 1}^1 \, \int\limits_x^{x + \Delta x} \, \displaystyle{{a\lpar x\comma \; v\lpar {\rm \xi} \comma \; {x_0}\comma \; {\rm \mu} \rpar \rpar F\lpar {x_0}\rpar } \over {\rm \mu} }{\rm d}{x_0}{\rm d}{\rm \mu} \comma \;$$where v(ξ, x 0, μ) is the solution of the Cauchy problem (A1), Δx is determined by the condition (A7).
Now we consider the kinetic steady-state equation for the distribution function f(x, v, μ) in the Fokker–Plank approximation. If all of the produced particles have the same initial velocity v 0 and the isotropic angular distribution, as well as the diffusion of the function f(x, v, μ) in the velocity space can be ignored, the considered equation reduces to the following Cauchy problem for the kinetic homogeneous equation (Gus'kov & Rozanov, Reference Gus'kov and Rozanov1982)
$$ {\rm \mu} v\displaystyle{{\partial f} \over {\partial x}} + \displaystyle{{\partial af} \over {\partial v}} = 0\comma \; \quad {v_{{\rm th}}}\lpar x\rpar \le v \le {v_0}\comma \; \quad f\lpar x\comma \; {v_0}\comma \; {\rm \mu} \rpar = - \displaystyle{{F\lpar x\comma \; {\rm \mu} \rpar } \over {2a\lpar x\comma \; {v_0}\rpar }}.$$The rate of energy transfer to unit volume of plasma at the point x has the form
$$ {W_{{\rm FP}}}\lpar x\rpar = \int\limits_{ - 1}^1 \, {\Psi _{{\rm FP}}}\lpar x\comma \; {\rm \mu} \rpar {\rm d}{\rm \mu} \comma \; \quad {\Psi _{{\rm FP}}}\lpar x\comma \; {\rm \mu} \rpar = - {m_p} \int\limits_{{v_{{\rm th}}}\lpar x\rpar }^{{v_0}} \, f\lpar x\comma \; v\comma \; {\rm \mu} \rpar a\lpar x\comma \; v\rpar v{\rm d}v.$$Instead of the velocity v, we introduce the variable η = v 2/2, setting f = f(x, η, μ), a = a(x, η), η0 = v 02/2, ηth(x) = v th(x)2/2. Then the function ΨFP from Eq. (A10) takes the form
$$ {\Psi _{{\rm FP}}}\lpar x\comma \; {\rm \mu} \rpar = - {m_p} \int\limits_{{{\rm \eta} _{{\rm th}}}\lpar x\rpar }^{{{\rm \eta} _0}} \, f\lpar x\comma \; {\rm \eta} \comma \; {\rm \mu} \rpar a\lpar x\comma \; {\rm \eta} \rpar {\rm d}{\rm \eta} .$$Using the change of variable x = x 0 + ξμ, we introduce the new independent variable ξ and the parameter x 0, which sense will be specified below. Equation (A9) takes the form
$$ \displaystyle{{\partial f} \over {\partial {\rm \xi} }} + a\displaystyle{{\partial f} \over {\partial {\rm \eta} }} + f\,{a_{\rm \eta} } = 0\comma \; \quad f = f\lpar {x_0} + {\rm \xi} {\rm \mu} \comma \; {\rm \eta} \comma \; {\rm \mu} \rpar.$$Its characteristics are determined by the differential equation dη = adξ and have the form of the one-parameter family
$$ {\rm \eta} \lpar {\rm \xi} \comma \; {x_0}\comma \; {\rm \mu} \rpar = {{\rm \eta} _0} + \int\limits_0^{\rm \xi} \, a\lpar {x_0} + {\rm \xi} ^{\prime}{\rm \mu} \comma \; {\rm \eta} \lpar {\rm \xi} ^{\prime}\comma \; {x_0}\comma \; {\rm \mu} \rpar \rpar {\rm d}{\rm \xi} ^{\prime}.$$where the parameter x 0 defines the characteristic that pass through the point (x 0, η0). Note, that the above family rewritten in another notations coincides with the family of solutions of the Cauchy problem (A1).
Along the characteristics (A13), equation (A12) and the boundary condition from the Cauchy problem (A9) take the form
$$ \displaystyle{{{\rm d}f} \over {{\rm d}{\rm \xi} }} + f\,{a_{\rm \eta} } = 0\comma \; \quad f = f\lpar {x_0} + {\rm \xi} {\rm \mu} \comma \; {\rm \eta} \lpar {\rm \xi} \comma \; {x_0}\comma \; {\rm \mu} \rpar \comma \; {\rm \mu} \rpar \comma \; \quad f\lpar {x_0}\comma \; {{\rm \eta} _0}\comma \; {\rm \mu} \rpar = - \displaystyle{{F\lpar {x_0}\rpar } \over {2a\lpar {x_0}\comma \; {{\rm \eta} _0}\rpar }}\comma \;$$that gives
$$ f\lpar {x_0} + {\rm \xi} {\rm \mu} \comma \; {\rm \eta} \lpar {\rm \xi} \comma \; {x_0}\comma \; {\rm \mu} \rpar \comma \; {\rm \mu} \rpar = - \displaystyle{{F\lpar {x_0}\rpar {\rm \chi} \lpar {\rm \xi} \comma \; {x_0}\comma \; {\rm \mu} \rpar } \over {2a\lpar {x_0}\comma \; {{\rm \eta} _0}\rpar }}\comma \;$$
$$ {\rm \chi} \lpar {\rm \xi} \comma \; {x_0}\comma \; {\rm \mu} \rpar = \exp\left({ - \int\limits_0^{\rm \xi} \, {a_{\rm \eta} }\lpar {x_0} + {\rm \xi} ^{\prime}{\rm \mu} \comma \; {\rm \eta} \lpar {\rm \xi} ^{\prime}\comma \; {x_0}\comma \; {\rm \mu} \rpar \rpar {\rm d}{\rm \xi} ^{\prime}} \right).$$In Eq. (A11), the integration over η is replaced by the integration over x 0 connected with η by Eq. (A13) and by the relation x = x 0 + ξμ = const. To simplify calculations, the relation between the differentials dη and dx 0 can be found by going from the variable ξ to the variable x in Eq. (A13), taking the differential of Eq. (A13) and setting dx = 0. As a result, we obtain
$$ {\rm d}{\rm \eta} = - \displaystyle{{a\lpar {x_0}\comma \; {{\rm \eta} _0}\rpar } \over {\rm \mu} }{\rm d}{x_0}.$$The limits of integration η0 and ηth(x) are changed to x and x + Δx respectively, where Δx is determined by the condition η(−Δx/μ, x + Δx, μ) = ηth(x) similar to the condition (A7) in the track method. By replacing the limits in Eq. (A11) and using Eqs. (A14) and (A15), we obtain
$$ {\Psi _{{\rm FP}}}\lpar x\comma \; {\rm \mu} \rpar = \displaystyle{{{m_p}} \over 2} \int\limits_x^{x + \Delta x} \, \displaystyle{{F\lpar {x_0}\rpar {\rm \chi} \lpar {\rm \xi} \comma \; {x_0}\comma \; {\rm \mu} \rpar a\lpar x\comma \; {\rm \eta} \lpar {\rm \xi} \comma \; {x_0}\comma \; {\rm \mu} \rpar \rpar } \over {\rm \mu} }{\rm d}{x_0}.$$One can see that the inverse change of variable η by variable v in Eq. (A16) and the substitution (A16) into (A10) give the function W FP(x) differed from the function (A8) in the track method by only the multiplier χ(ξ, x 0, μ), which is defined by Eq. (A14). Therefore, the track method described at the beginning of the section for finite numbers of the grid cells on x and the angular groups on μ can be used as the numerical method for the problem (A9), (A10) with the additional computation of the function χ(ξ, x 0, μ) along the every trajectory of particles.
To compute the rate of the energy transfer to unit of plasma volume separately for electrons and ions, the decelerations of particles by plasma components a e(x, v) or a i(x, v) should be used in Eq. (A3) instead of a(x, v).
If one of the boundaries, x = 0, for example, is the symmetry plane, the condition of symmetrical reflection of the particle trajectory from the boundary at μ < 0
$$ x = \left\{{\matrix{ {{x_0} + {\rm \xi} {\rm \mu} \comma \; }&{{\rm \xi} \le {{\rm \xi} _ * } = - {x_0}/{\rm \mu} \comma \; } \cr { - \lpar {\rm \xi} - {{\rm \xi} _ * }\rpar {\rm \mu} \comma \; } &{{\rm \xi} \gt {{\rm \xi} _ * }\comma \; } \cr } } \right.$$is added to the problem (A1), and the condition x 0 + ξmaxμ ≥ 0, that determines maximal ξ value, is removed from the conditions (A2).
In the case of the quasi-1D model, the condition (2.17) is added to the conditions (A2).
