Traveling Wave Reactor and Condition of Existence of Nuclear Burning Soliton-Like Wave in Neutron-Multiplying Media

Physical fundamentals of traveling wave reactor are considered. We show that the condition of existence of nuclear burning soliton-like wave in a neutron-multiplying medium is determined in general by two conditions. The first condition (necessary) is determined by relationship between the equilibrium concentration and critical concentration of active (fissionable) isotope that is a consequence of the Bohr–Sommerfeld quantization condition. The second condition (sufficient) is set by the so-called Wigner quantum statistics, or more accurately, by a statistics of the Gaussian simplectic ensembles with respect to the parameter that describes the squared width of burning wave front of nuclear fuel active component.


Introduction
In spite of obvious efficiency of the nuclear power engineering of new generation, the main difficulties of its acceptance are predetermined by non-trivial properties which a future ideal nuclear reactor must possess.At first, the natural, i.e., unenriched uranium or thorium must be used as a nuclear fuel.Secondly, the traditional control rods must be completely absent in the control system of nuclear reactor.Thirdly, in spite of the absence of control rods a reactor must possess the property of so-called internal safety.It means that a reactor core must always be in the critical state, i.e., the normal operation of reactor is automatically provided not as a result of personnel activity, but by virtue of underlying physical laws, which prevent the explosive development of chain reaction in a natural way.Figuratively speaking, the reactor with internal safety it is "the nuclear installation which never explodes" [1].It seems to be strange, but the reactors satisfying such unusual requirements are really possible.For the first time the idea of such self-regulating fast reactor in general terms (so-called breed-and-burn mode) was proposed by Russian physicists Feinberg and Kunegin in 1958 at II Genevan conference [2] and relatively recently was "resuscitated" as an idea of the self-regulating fast reactor in traveling-wave mode by Russian physicist L. Feoktistov [3] and independently by American physicists Teller, Ishikawa and Wood [4].
The main idea of reactor with internal safety consists in the selection of fuel composition so that, at first, the characteristic time τ β of the nuclear burning of fuel active (fissionable) component is considerably greater than the characteristic time of delayed neutrons production and, secondly, necessary self-regulation conditions are maintained during the reactor operation (that always take place, when the equilibrium concentration ñfiss of fuel active component is greater than its critical concentration [5] n crit [3]).These very important conditions can practically always to be attained, if in the reactor the chain of nuclear transformations of the Feoktistov uranium-plutonium cycle type [3] 238 U (n, γ) → 239 U or the Teller-Ishikawa-Wood thorium-uranium cycle type [4] 232 T h(n, γ) → 233 T h will be present at sufficient level.
In both cases the produced fissionable isotopes 239 P u in (1) or 233 U in (2) are the active components of nuclear fuel.The characteristic time of such a reaction corresponds to the proper β-decays time and is approximately equal to τ β = 2.3/ ln 2 3.3 days for the reaction (1) and τ β 39.5 days for the reaction (2).This is several orders greater than the time of delayed neutron production.
The self-regulation of nuclear burning process consists in the fact that such a system left by itself can not pass from a critical state to reactor runaway mode, because a critical concentration is bounded above by the finite equilibrium concentration of nuclear fuel fissionable component (plutonium in (1) or uranium in (2)), i.e., ñfiss > n crit (the condition for existence of Feoktistovs wave mode [3]).At phenomenological level the self-regulation of nuclear burning is manifested as follows.The increase of neutron flux for some reason or other leads to the rapid burn-up of nuclear fuel fissionable component (plutonium in (1) or uranium in (2)), i.e., its concentration and, as a result, the neutron flux will decrease, while the new nuclei of corresponding fissionable component of nuclear fuel are produced with the Energies 2011, 4 same rate of production during time τ β .And vice versa, when the neutron flux is sharply decreased due to external action, the fuel burn-up rate decreases too, and the accumulation rate of fuel fissionable component will increase as well as the number of neutron production in a reactor during the same time τ β .
However, as is known [3], the Feoktistov condition for existence of wave mode is only necessary but not sufficient condition.Therefore complete generalization of the condition for existence of wave mode for critical waves of nuclear burning in neutron-multiplying media is the purpose of this article.

The Features of Condition for Existence of Wave Mode of Nuclear Burning Critical Wave
According to Feoktistov As is known, to verify any physical hypothesis one can use simplified models.One of the possible simplifications consists in a separate consideration of neutron-nuclear processes and heat removal process.Such a simplification is especially justified at a long control time.At the same time, neutron processes can be studied in the one-dimensional geometry rather than in the three-dimensional geometry and also can be considered in the diffusion approximation and one-group approximation.This means that neutron spectral characteristics are averaged in an appropriate manner over the fixed neutron energy interval, and the problem is solved at the fixed neutron energy.
Following References [3,6], let us consider the kinetics of self-regulating fast uranium-plutonium reactor, where the Feoktistov self-propagating neutron-fission wave of nuclear burning is generated.Such a wave mode meets all requirements, which are appropriate to the nuclear reactor with internal safety.
Main transmutation chain corresponds to the uranium-plutonium fuel cycle (1).We consider the 238 U -filled half-space, which is irradiated by neutrons from the open surface.We also assume for simplicity that the neutron spectrum and fission spectrum are identical.The typical neutron energy in the medium strongly depends on moderator properties.Following Feoktistov [3,7], we will consider the case, when the moderator is absent or is present in small amounts and the neutron spectrum practically coincides with initial one.
The main goal of such a simplified model consists in finding of an autowave solution for the transmutation chain (1) under the indispensable condition n crit < ñPu .
The balance equation for plutonium concentration looks like where N P u is the 239 P u concentration, N 8 is the 238 U concentration, n is the neutron concentration, v is the neutron velocity in the one-group approximation, σ a and σ f are the neutron-capture cross-section and fission cross-section, respectively.From here it follows that the plutonium relative equilibrium concentration, when the derivative becomes zero, has the form ñPu = ÑPu Here N P u | t=0 = 0 and the current concentration N P u (t) cannot exceed ÑPu .
Recall that the value of constant ñPu strongly varies depending on the neutron energy, for example, for a thermal reactor ñPu = 0.25%, while for a fast reactor ñPu ≈ 10%.

1340
The other characteristic of the uranium-plutonium medium is the plutonium critical concentration n crit .At n P u > n crit the system becomes supercritical and capable of self-multiplication, and, conversely, at n P u < n crit the system becomes subcritical and the neutron flux density attenuates with time.
The value n crit can be obtained from the neutron balance where n i are the relative concentrations of the elements taking part in the reaction, ν is the average number of prompt neutrons per 239 P u fission, i σ a is the absorption cross-section for i-th element including U and Pu (neutron escape can be also included).
The magnitude determines the plutonium concentration when the multiplying medium is in critical state.The magnitude n crit is also a function of neutron energy.Since the two dimensionless numbers ñPu and n crit are composed from the different combinations of constants, the following variants ñPu > n crit , ñPu < n crit and ñPu = n crit are possible.It appeared that for thermal neutrons ñPu < n crit and for fast neutrons ñPu > n crit .In the first case the system is viable only in the presence of neutron source.If this external source will be switched off, the reaction immediately stops.In other case, which we will consider below, the asymptotical and independent of initial conditions solution in the form of stationary traveling wave is possible.This is not paradoxical, if to take into account the foregoing explanation of physical sense of Equation (1).Thus, in case of 238 U slow burning induced by fast neutrons it is impossible to overstep the criticality threshold.On formal level this makes possible to neglect the partial time derivative in the neutron transport equation.Further developments can be presented in the following way.Neutrons from the external source are absorbed in the nuclear fuel layer, whose thickness does not exceed the neutron free path, and due to fission the uranium is transmuted into the plutonium.With the plutonium accumulation the fission process intensifies, the neutron number grows and fission reaction begins to develop farther.In time the center of energy release shifts away from the neutron source, whose influence decreases, and the system comes to the stationary traveling-wave mode.In other words, all functions become depending on the argument z = x + ut (the wave is directed from right to left and u is its velocity).At the same time the wave velocity u is of order L/τ , where L ≈ 5 cm is the neutron diffusion length and τ = 2.3/ ln 2 = 3.3 days is the 239 P u-production time from 238 U .
To obtain the system of kinetic equations for neutrons and nuclei taking part in the transmutation chain (1) with respect to the autowave variable z = x + ut, we will write down at first this system of equations in the traditional coordinates {x; t}.
So, the kinetic equation for neutrons has the form dñ(x, t) dt = D∆ñ(x, t) + q(x, t) Energies 2011, 4 1341 where ñ(x, t) is a neutron density; D = v/3Σ S is the neutron diffusion coefficient, cm 2 /s; v is the neutron velocity in one-group approximation, cm/s; Σ S is the neutron macroscopic scattering cross-section, cm −1 ; q is the neutron source strength, cm −3 s −1 .
It is obvious that within the framework of our problem the expression for the neutron source strength looks like where N i is the concentration of i-th isotope in the reaction chain (2), P u σ f is the microscopic fission cross-section for 239 94 P u, i σ a is the neutron-capture cross-section by i-th isotope, ν is the average number of prompt neutrons per 239 94 P u fission.Using Equation (6) for the plutonium relative critical concentration and normalizing the neutron concentration and isotope concentration to the 238 U concentration, we obtain the diffusion Equation (7) in the following form where or, in other words, n(x, t) and n i (x, t) are the relative concentration of neutrons and isotopes, respectively, n P u and n P u crit are the plutonium relative equilibrium and critical concentration, respectively.We have mentioned above, that the inequality ñPu > n crit formally predetermines the stationary form of the kinetic equation.So long as for qualitative analysis of Equation (7) it is sufficient to find an approximate solution of steady-state equation in the region in front of wave [3,7], i.e., in the asymptotical region x −→ −∞ we can neglect the summands n 9 and n P u , which are finding together with n 8 ∼ = 1 .This means that the stationary form of the kinetic Equation (9) in the simplified form looks like where the constant C ∞ is equal to τ is the neutron lifetime, D = L 2 /τ is the neutron diffusion coefficient, cm 2 s −1 and concentration N 8 (x, t) is equal to the initial concentration N 8 (x, 0).Following Refs.[3,7], we will seek the approximate solution of the steady-state Equation (11) in the region in front of wave (x −→ −∞).
Energies 2011, 4 1342 Now we write down the kinetic equations for each of isotopes taking part in the reaction chain (1).First we write down an equation for the 238 92 U isotope In References [3,7] the second summand on the right-hand-side of Equation ( 13) is introduced under the assumptions that the 238 U -neutron and 239 P u-neutron capture leads to 240 P u production, whose properties are assumed to be identical to the properties of 238 U initial nucleus.Obviously, this is done to simplify the problem or, in other words, to satisfactorily close the systems of kinetic equations for neutrons and nuclei taking part in the reaction chain (1).
Assuming that 8 σ f 8 σ a [3,7], the kinetic equation for 238 U takes the following form where The kinetic equation for 239 U looks like where τ β = τ 9 β + τ N p β is the 239 U lifetime which is equal by assumption [3,7] to the sum of half-lives of 239 U and 239 N p β-radioactive nuclei.
Taking into account that 9 σ f 9 σ a and normalizing to initial concentration N 8 (x, 0)we obtain Equation (16) in the following form Finally we write down the kinetic equation for 239 P u as Taking into consideration the expression for the plutonium equilibrium concentration (4) and normalizing to the initial concentration N 8 (x, 0), we have Equation (18) in the form Now we are ready to write down the system of kinetic equations for neutrons and nuclei that take part in the transmutation chain (2) with respect to the dimensionless autowave variable z = ξ/L = (x + ut)/L, where u is the velocity of steady-state traveling wave going from left to right (as in References [3,7]) and L is the neutron average diffusion length.
Energies 2011, 4 1343 For this we use the following operators Following [3,7], we assume, without loss of generality, that 8 σ a ∼ = 9 σ a ∼ = P u σ a .From here it follows that c 8 = c 9 = c P u .Then, introducing the dimensionless constant Λ = uτ β /L and the variable n * (z) = c 1 τ β n(z), and simultaneously performing the transformation of coordinates in Equations ( 11), ( 14), ( 17) and ( 19) we obtain the system of kinetic equations for neutrons and nuclei taking part in the reaction chain (1) with respect to the dimensionless autowave variable z where u is the phase velocity of stationary traveling wave, L is the neutron average diffusion length, v is the neutron velocity in the one-group approximation, cm • s −1 ; Σ s is the the neutron macroscopic scattering cross-section, cm −1 ; τ = 1/v σ i a N i is the neutron lifetime in medium, N i is the nucleus concentration of 8( 238 U ), 9( 239 N p) and 239 P u; Λ = uτ β /L is the dimensionless constant, is the plutonium relative critical concentration, n i is the nucleus concentration of 8( 238 U ) and 9( 239 N p) normalized to 238 U initial concentration, i.e., to N 8 (−∞), σ a and σ f are the microscopic neutron capture cross-section and fission cross-section, respectively, ν is the average number of prompt neutrons produced per plutonium nucleus fission.Solving these equations Feoktistov has used the analogy between the diffusion equation and the Schrödinger steady-state equation in the quasi-classical approximation [3].Naturally, in this case (see Equation ( 22)) the solution stationarity condition is satisfied integrally, because there are points where n P u > n crit , and there are points where n P u < n crit .In this sense, the region n P u > n crit corresponds to kind of allowed region, whereas the region n P u < n crit corresponds to subbarrier region.In other words, the inverted profile of plutonium concentration in the 238 U medium plays the role of potential well (Figure 1 [8]).For the region in front of wave (z = −∞) the approximate solution looks like It will be recalled that searching this solution, we have neglected summands n 9 and n P u whose values are determined by edge condition n 8 ∼ = 1.Then, assuming that the subbarrier region ends at z = 0, we have n P u = n crit at this point.This makes it possible to determinate the constant C. According to the Bohr-Sommerfeld quantization condition, we have the following equality at the point z = a a 0 where the integral is taken over the supercritical region (n P u > n crit ).
At the same time the condition (30) plays also the role of condition for finding the point a at n P u = n crit , i.e., when the transition into subbarrier region happens due to burn-up (see Figure 1a [8] and Figure 2) [9].
Figure 2. The schematic view of the permitted and subbarier (gray colored) region corresponding to the conditions n P u > n crit and n P u < n crit , respectively.The delineated by square region is considered more particularly in Figure 3. n Pu Performing the ordinary for quasi-classical approximation matching with the supercriticality region As a critical state is automatically maintained at n P u > n crit [3] (that is the direct consequence of the Bohr-Sommerfeld quantization condition), we can use this fact for generalization of the following inequality: ñPu > n P u > n P u crit (31) where ñPu is the plutonium equilibrium concentration (see Figure 1).Thus, Feoktistov have shown for the first time that the soliton-like propagation of neutron-fission wave of nuclear burning is possible in 238 U medium only at definite ratio between the equilibrium and critical plutonium concentrations (ñ P u > n crit ), which is imposed by the Bohr-Sommerfeld quantization Energies 2011, 4 1346 condition.In other words, only in this case the critical (quasi-stationary) state of system (reactor core) can automatically maintained without any external intervention, and, consequently, only in this case the reactor completely and unambiguously possesses the property of internal safety.
It is appropriate here to pay an attention to very important Feoktistovs parameter, which, as shown below, is basis for the existence of soliton-like wave of nuclear burning: where a is the width of allowed range of integration in the Bohr-Sommerfeld condition (30), where the inequality n P u > n crit (Figure 2) and ñPu > n crit , respectively, are satisfied; Λ(a) is the dimensionless coefficient, which appears within the framework of simplified diffusion model of the Feoktistov reactor ( 22)- (25).
It is obvious that Equation (30) due to its physical meaning is a key factor, which predetermines the phase velocity of soliton-like burning wave.Therefore, this equation exists regardless of level of idealization of reactor core model and should appear in the explicit or implicit form in any model in which the system of kinetics equations for neutrons and nuclei has soliton-like solutions for neutrons.At the same time, as the average width of soliton wave is about 2L, the maximum values for the dimensionless coefficient Λ(a) and wave velocity u are determined by the following approximate equality where the coefficient b is about 2, although a final estimate will be given below.From analysis of Equation (33) it follows that the velocity of soliton-like wave propagation does not always equal to the diffusion speed u = L/τ β .It can be considerably slower or faster due to very strong domination of the nonlinearity parameter or, conversely, the variance parameter, which reflect the peculiarities of nuclear transformation kinetics (for example, in the chain (1) and/or in (2)).In practice, they manifest itself as different degree of fuel burn-up.
In other words, when the wave velocity and consequently the degree of fuel burn-up are low, the wave stops due to the following reasons.Neutrons from an external source, which take place in the initial stage of wave initiation, burn out the plutonium on the medium boundary and simultaneously transmute the uranium into 239 N p.Neptunium with time begins to produce the plutonium, but it can not create the required high concentration, while the 239 P u accumulation decreases due to the uranium burn-up.More and more thick layer without both 238 U and 239 P u grows near to the medium boundary.The neutron diffusion through this layer does not provide the increase of plutonium concentration in next layers, and the wave does not arise even at n P u (x, 0) = n crit .And, vice versa, when the wave velocity and consequently degree of fuel burn-up are high, the wave stops also because of the scarce (or more exactly, delayed) plutonium production which takes place due to another reason.Figuratively speaking, the situation resembles the forest fire by strong wind, when only crown posts burn.When the wind velocity increases, it could extinguish the fire at all.We have the similar situation, when there is a velocity, at which in the early stage (when x ≈ 0) the neutron soliton wave front outruns the plutonium (sigmoid) wave front, and this advance exceeds the neutron diffusion length.This leads, in fact, to transformation of fast wave into slow wave or to its full stop.It is interesting to note that this case is not Energies 2011, 4 1347 discussed in the literature (with the exception of [8,10]), but it is possible to assume that it corresponds to some hypothetical situation, when the nuclear burning wave forms in highly enriched fuel which has the ultra-low critical concentration of fuel fissionable component.
Thus, when the lag (Figure 1b [8]) or advance of neutron wave front relative to the plutonium wave front is considerably greater than the neutron diffusion length, these waves stop and totally degrade.This means that the degradation of waves with very low or very high initial phase velocity will exhibit in the fact that Equation (32) tends to zero at some low or very high values of a. Therefore taking into account Equation (33), we can conclude that Equation (32) is true in the range 0 (1/b)Λ(a) 1.Based on this generalization, we can make an important assumption that the expression (1/b)Λ(a) means the certain probability density distribution p(a) with respect to a: Let us consider and substantiate the type and main properties of such a statistics and also show the results of its verification based on the known computational experiments on simulation of nuclear burning wave in the U-Pu (1) and Th-U (2) fuel cycles.

Chaos and Integrability in Nonlinear Dynamics of Reactor Core
In order to solve the assigned task we use the known analogy between the neutron diffusion equation and the Schrödinger steady-state equation in quasi-classical approximation.We would remind that earlier we have used this analogy to search the solution of the system of kinetics equation for neutrons and nuclei ( 22)- (25) in the reaction chain (1) of the U-Pu fuel cycle.Since the system of equations for neutrons and nuclei in the Th-U fuel cycle (2) is structurally identical to the system equation for the U-Pu fuel cycle (1), the computed "quantum mechanical" solution, which describes the statistics (34), will be general for both fuel cycles, except for a few inessential details.
So, due to mentioned analogy we use the Bohr-Sommerfeld quantization condition, which in the case of the one-dimensional systems determines the energy eigenvalues E n in the explicit form where m and p(x) are the particle mass and particle momentum in the field of some smooth potential V (x).
For the Feoktistov nearly integrable system of the equations ( 22)-( 25) or for the similar Teller system of equations, for which it is assumed that m = 1/2, = 1, V (x) = 1 and n = 0, this condition is applied in the following form where index f is signifies the fissionable isotope, for example, the 239 P u in the Feoktistov U-Pu fuel cycle (1) or the 233 U in the Teller Th-U fuel cycle (2).However, in describing the real evolution of fast reactor core, the corresponding systems of equations for neutrons and nuclei are non-integrable almost without exception.This, in its turn, means that Energies 2011, 4 1348 according to the Kolmogorov-Arnold-Moser theorem [11,12] quasi-classical quantization formulas are inapplicable for the system, where the motion in phase space is not restricted by multidimensional tori.This is caused by the fact that in the Hamiltonian non-integrable systems the more and more number of tori is collapsed in phase space with perturbation (non-integrability) growth.As a result, the trajectories of majority of coupled states get entangled in phase space, the motion becomes mainly chaotic, and coupled states themselves and their energies can not be described by the rules of quasi-classical quantization, for example, such as the Einstein-Brillouin-Keller (EBK) quantization rule for multidimensional case [12,13], which generalizes the Bohr-Sommerfeld quantization rule.Note that now a notion "quantum chaos" joins the range of problems related to quantum-mechanical description of systems chaotic in a classic limit [14,15].
Since the results of random matrices theory will be used for research of chaotic properties of the statistics (30), we give in advance an overview of the main concepts of this theory.
First, following [14,15], let us shortly consider a nature of so-called universality classes and kinds of the Gaussian ensembles.It is known, if the Hamilton operator matrix has any symmetry, it can be reduced to the block-diagonal form.At the same time, the matrix elements in each block are specified by a certain quantum number set.For simplicity, we assume below that the Schrödinger equation i (∂ψ/∂t) = Ĥψ describes the states belonging only to the one block.At the same time, the size of the operator matrix Ĥ is finite and equals to some integer.
As shown in [14,15], these universality classes divide physical systems into groups depending on their relation to orthogonal, unitary or symplectic transformations, after which the matrix Ĥ remains invariant.In other words, as it ascertained in [14]: • the Hamiltonian of spinless system, which is symmetrical with respect to time inversion, is invariant under orthogonal transformations and can be represented by a real matrix; • the Hamiltonian of spinless system, which is not symmetrical with respect to time inversion, is invariant under unitary transformations and can be represented by the Hermitian matrix; • the Hamiltonian of the system with the spin of 1/2, which is symmetrical with respect to time inversion, is invariant under symplectic transformations and can be represented by a quaternion real matrix.Now let us talk about the Gaussian ensembles.If the matrix element distribution function is invariant under one of mentioned transformations, this means that the sets of all matrices with elements, which are described by these distribution functions, form the Gaussian orthogonal ensemble (GOE), the Gaussian unitary ensemble (GUE) and the Gaussian symplectic ensemble (GSE), respectively.
At the same time it should be noted the one very important detail.The matrix element distribution function of the Gaussian ensembles can not be directly measured, because the experiment can give us information about the energy levels of investigated quantum-mechanical system only.In other words, just the distribution function of energy eigenvalues is of greater interest from the practical point of view.
Energies 2011, 4 1349 Derivation of corresponding equations for the considered types of the Gaussian ensembles can be found in [15].At the same time, the correlated distribution function of energy eigenvalues for all ensemble types it is possible to write down in the sufficiently universal form: where ν is the universality index, which takes on the value of 1, 2 and 4 for GOE, GUE and GSE statistics, respectively.At ν = 0 the energy eigenvalues are not correlated.In this case, the energy level spacing distribution function is described by the Poisson statistics, and the matrix ensemble itself is called the ensemble.
So long as the energy level spacing distribution function is the most studied characteristic of chaotic systems, following [14], we give a calculation only for relatively simple case of the Gaussian ensemble of the matrices 2 × 2 in size.We calculate the energy level spacing distribution function p W (s) substituting the function P (E1, E2) into (37): Constants A and C are determined by the two normalization conditions: The first condition is the total probability normalization to unit, and the second condition is the average energy level spacing normalization to unit.Integration of (38) gives the so-called Wigner energy level spacing distribution functions, which correspond to the different Gaussian ensembles: Despite the fact that these functions were obtained for the Gaussian ensembles of matrices 2 × 2 in size, they describe with good accuracy the spectra of arbitrary size matrices [14].
Note that random matrix theory at first was developed to find some regularities for heavy nucleus energy spectra [15,16], but it attracted keen interest after the Bohigas, Giannoni and Schmit conclusion [17] that this theory is applicable to any chaotic system.
Energies 2011, 4 1350 Returning to our problem, we will attempt to use the considered statistical properties of the Gaussian ensembles to determine the Equation (34) statistics type.

The Wigner Quantum Statistics and Generalized Condition for Existence of Traveling Wave Mode of Nuclear Burning
To apply the results of previous section, in the framework of nearly integrable system, to which the system of equations describing the kinetics of nuclear burning in the Feoktistov U-Pu fuel cycle (1) or the Taylor Th-U fuel cycle (2) belongs, we use some formal quantum analogy of this system [18].
It is obvious that to research the Wigner statistics type a two-level quantum system is needed, at least.Therefore we assume that a quantum energy system equivalent to the analyzed area of nuclear burning is two-level (Figure 2) and formally introduce the "energy" eigenvalue of stationary ground state as (n f iss ) 0 /n crit = E 0 and "energy" eigenvalue of quasistationary state as (n f iss ) quasi /n crit = E quasi (Figure 3).At the same time E 0 > E quasi and (n f iss ) 0 is some specified value of fissionable isotope concentration limited from above by the fissionable isotope equilibrium concentration, i.e., (n f iss ) 0 < ñfiss .We assume also that the quasistationary level is situated near the potential well bottom, i.e., E quasi → 1 (see Figure 3), and is strongly unstable [19].Therefore, the nuclear burning mode "lives" most of the time in the ground energy state, i.e., on the level E 0 .
So, to describe the wave mode of nuclear burning we use below the quantum-mechanical analogy, in the framework of which the "energy" spectrum of nuclear burning in allowed region is described by some quasiequivalent two-level scheme (Figure 3).
Then, for the nearly integrable system which describes the kinetics of nuclear transformations in the Feoktistov (1) or the Teller (2) fuel cycle in general case we can use the Bohr-Sommerfeld approximate condition in the following form From here the obvious and important assertion follows: by virtue of the Bohr-Sommerfeld condition (42) the type of the Wigner energy level spacing statistics unambiguously predetermines the analogous statistics type of parameter, which characterizes the squared width (a 2 ) of concentration wave front of active (fissionable) material.
Note that we have not any information about the value of "energy" E 0 before the experiment, whereas it is possible to consider that E quasi = 1.If, additionally, in the steady-state mode all wave kinetic parameters are predetermined by the equilibrium ñfiss and critical n f iss crit concentration of active (fissionable) isotope (whose values are known before experiment), the physical meaning and the necessity of following change It is obvious that the conditions ( 42) and ( 43) make it possible to obtain the expression for the parameter a * : The next step for determining the statistics p(a * ) of Equation ( 34) type consists in the experimental verification of proposed hypothesis.With that end in view we have compared the Gaussian ensemble statistics (41) with well-known computational experimental data [8,[23][24][25][26][27] and have obtained a good accordance of calculation data with the theoretical dependence, which is described by the Gaussian symplectic ensemble statistics (see Table 1 and Figure 4).Thus, we can conclude that the wave velocity u (see Equation ( 34)) is predetermined by the following approximate equality where coefficient b = 2 (see Equation ( 34)); τ β is the delay time caused by active (fissionable) isotope production, which is equal to the β-decay period of compound nuclei in the Feoktistov (1) or the Teller (2) fuel cycle; p s W (a * ) is the Wigner symplectic statistics.Thus, based on the verification results of Equation we can make a conclusion, which generalizes the physical conditions for existence of Feoktistovs wave mode.The velocity of soliton-like wave propagation in neutron-multiplying medium must be determined in general case by the two conditions.The first condition (necessary) is predetermined by relationship between the equilibrium concentration and critical concentration of active (fissionable) isotope (ñ P u /n crit > 1) or, more exactly, by the corresponding Bohr-Sommerfeld quantization condition.The second condition (sufficient) is set by statistics of the Gaussian symplectic ensembles with respect to the parameter a, which describes the burning wave thickness of active (fissionable) component of nuclear fuel.

Computation 3D-Experiment and Verification of the Wigner Quantum Statistics
We consider here the simplified diffusion model for neutrons and nuclei kinetics in the chain (1) in the one-group approximation (neutrons energy is ∼1 M eV ) and cylindrical geometry.The corresponding system of differential equations, which describes the kinetics of Feoktistovs U-Pu fuel cycle with consideration of delayed neutrons, i.e., the kinetics of initiation and propagation of neutron-fission wave n(x, t), is as follows [23] ∂n(x, t) ∂t = D∆n(x, t) + q(x, t) where Ñi ln 2 Energies 2011, 4

1353
To specify the last term q(x, t) on the right-hand-side of Equation ( 46), we use the approach of effective additional neutron absorber: Taking into account the fact that fission with two fragment formation is most probable, the kinetic equation for N (x, t) becomes Ñi ln 2 n(x, t) is the neutron density; D is the neutron diffusion constant; v n is the neutron velocity (E n = 1 M eV , the one-group approximation); Ñi are the concentrations of neutron-rich fission fragments of 239 P u nuclei; N 8 , N 9 , N u are the concentrations of 238 U , 239 U , 239 P u, respectively; Ni are the concentrations of remaining fission fragments of 239 P u nuclei; σ a is the neutron-capture microcross-section; σ f is the fission microcross-section; τ β is the nucleus life time with respect to the β-decay; p i are the parameters characterizing delayed neutron groups for main fuel fissionable nuclides [27].
The boundary conditions for the system of differential Equations ( 46)-(50) are where Φ 0 is the neutron density of plane diffusion neutron source, which is located on the boundary x = 0; l is the uranium bar length.An estimation of the neutron flux density from the inner source on the boundary Φ 0 can be obtained from an estimation of the P u critical concentration which is of order of 10%: and therefore we have Note that Equation ( 55) is only an estimation of Φ 0 .The results of the computational experiment show that it can be substantially smaller in reality.
In general, different boundary conditions can be used, depending on the physical conditions under which nuclear burning is initiated by the neutron source, for example, the Dirichlet condition of (29) type, the Neumann condition or a so-called third-kind boundary condition, which generalizes first two conditions.The use of the third-kind boundary condition is recommended in neutron transport theory [27].In the simple case this condition (which is known as the Milne problem) is a linear combination of the neutron concentration n(x, t) and its spatial derivative ∂n/∂x(x, t) on the boundary n(0, t) − 0.7104λn (1,0) where λ is the neutron free path and n (1,0) (0, t) ≡ ∂n/∂x(0, t).
Energies 2011, 4 1354 Although the "neutron source-nuclear fuel" system behavior depends on the boundary conditions near the boundary, computational experiments show that inside the reactor core, i.e., far from the boundary, the system behavior is asymptotically invariant.This confirms the independence of wave propagation in reactor volume on the boundary conditions and the way in which the nuclear burning is initiated.In this sense, the problem of determining the optimum parameters of nuclear fuel "ignition" in the "neutron source-nuclear fuel" system is a nontrivial and extraordinarily vital issue, which requires a separate examination.
The initial conditions for the system of differential Equations ( 46)-(50) are where ρ 8 is the density, g•cm −3 ; µ 8 is the gram-molecular weight, g•mole −1 ; N A is the Avogadro number.
The following values of constants were used for the simulation: .38 • 10 −26 cm 2 ; σ 9 a = σ P u a = 2.12 • 10 −26 cm 2 (61) The system of Equations ( 46)-(51) with the boundary conditions (53)-(56), initial conditions (57)-( 59) and the values of constants (60)-( 62) is solved numerically using the software package Fortran Power Station 4.0.At the same time we use the DMOLCH subprogram from the IMSL Fortran Library.The DMOLCH subprogram solves a system of partial differential equations of the form u t = f (x, t, u x , u xx ) by the method of straight lines [23,28].The solutions for diffusion model of neutrons and nuclei kinetics in the chain (1) in the one-group approximation and cylindrical geometry are presented in Figure 5.
Verification of the Wigner symplectic statistics consists in comparison of the experimental velocity of nuclear burning wave obtained by a computational 3D-experiment with its theoretical value obtained by Equation (45).For this purpose we at first find the plutonium critical concentration n crit from the profile of its experimental concentration distribution (Figure 5).It is obvious that the absolute value of critical concentration approximately equals to N P u crit ∼ = 8 • 10 20 cm −3 (see Figure 6b).It follows from here that the plutonium normalized critical concentration is where by virtue of Equation (58) the initial uranium concentration is N 8 (x, 0) = 4.79 • 10 22 cm −3 and the value of a * is equal to 0.704 by virtue of Equation (44).In other words, the important case when a * < 1 takes a place (see Figure 4).y the method of straight lines [19,24].The solutions for diffusion model of neutrons and nuclei kinetics the chain (1) in the one-group approximation and cylindrical geometry are presented in Fig. 5. Verification of the Wigner symplectic statistics consists in comparison of the experimental velocity f nuclear burning wave obtained by a computational 3D-experiment with its theoretical value obtained y Eq. (45).For this purpose we at first find the plutonium critical concentration n crit from the profile f its experimental concentration distribution (Fig. 5).It is obvious, that the absolute value of critica oncentration approximately equals to N P u crit ∼ = 8 • 10 20 cm −3 (see Fig. 6b).It follows from here that the lutonium normalized critical concentration is here by virtue of Eq.(58) the initial uranium concentration is N 8 (x, 0) = 4.79 • 10 22 cm −3 and the value f a * is equal to 0.704 by virtue of Eq. (44).In other words, the important case when a * < 1 takes a lace (see Fig. 4).
Taking into account the plutonium normalized equilibrium concentration ñfiss = 0.1, by virtue of Eq 5) we obtain the Wigner theoretical symplectic probability: Taking into account the plutonium normalized equilibrium concentration ñfiss = 0.1, by virtue of Equation ( 45) we obtain the Wigner theoretical symplectic probability: which corresponds to the nuclear burning wave velocity u theor = 2.82 cm/day at given parameters L = 5 cm and τ β = 3.3 days.Now it is easy to determine the experimental value of nuclear burning wave velocity and, accordingly, the Wigner symplectic probability.In Figure 6a the profile of experimental concentration distribution of neutrons is shown.We can see that the wave crest has covered the distance of 600 cm during t = 217 days.Thus, the velocity of nuclear burning neutron wave is u simul = 600/217 2.77 cm/day (65) This, in its turn, corresponds to the value of 1/2Λ(a * ) = p s W (a * ) = 0.9141.Thus, the approximate equality of the experimental and theoretical velocity of nuclear burning wave (u theor ∼ = u simul ) makes it possible to conclude that the Wigner quantum (symplectic) statistics verified by computing 3D-experiment (see Figure 4) satisfactorily describes experimental data characterized by the parameter Λ(a * ).Here we note that computing experiments show that the conditions of wave blocking, which describe the degradation and subsequent stop of wave, are predetermined by the degree of burn-up of the main nonfissionable ( 238 U ) and fissionable ( 239 P u) components of nuclear fuel in front of the wave by neutrons from an external source in the initial stage of wave "ignition".This process is very important, since when the degree of fuel burn-up in front of the wave is high, the wave can not overcome this "scorched" region just like a fire in the steppe can not overcome the plowed stripe of the land in front of combustion wave.It is obvious that in the initial stage of wave initiation the degree of fuel burn-up is determined, first of all, by the energy and intensity of neutrons from the external source and by the properties of nuclear fuel.The most important of these properties is the delay time τ β of active (fissionable) isotope production, which is equal to the effective period of β-decay of compound nuclei in the Feoktistov U-Pu fuel cycle (1) or the Teller Th-U fuel cycle (2).

1357
In spite of the general understanding of physics of nuclear burning wave blocking, it is obvious that indicated above difficulties in the describing this process testify to nontriviality of given problem.Unfortunately, the solving of this problem exceeds the scope of this work, but it will be a subject of future research.

On the Depth of Analogy between Diffusion Equation for Neutrons and the Schrödinger Equation in the WKB Approximation
It is known [29], that the solution of differential equation of the form in the Liouville-Green approximation looks like where |f −1/4 | is the sufficiently small and slowly varying function, A and B are arbitrary constants.
From Equation (67) the famous expression for the phase with consideration of turning point follows Note that physicists often name Equation (67) the WKB approximation after Wentzel, Kramers and Brillouin, who had developed a method for finding approximate solutions to linear partial differential equations with spatially varying coefficients.However, as is known [29], their contribution consists not in construction of the Liouville-Green approximation which had been already obtained by that time, but in establishing of the equations which connect exponential and oscillatory solutions in turning points on real axis.Thus, it is obvious that to obtain the necessary condition for existence of the nuclear burning wave (31) the use of the Bohr-Sommerfeld quantization condition (30) or, more precisely, the analogy between the diffusion equation and the Schrödinger equation in the WKB approximation was not quite obligatory.The use of Equation ( 68) is quite enough in this case, because, as is well-known [29], this equation is the independent result of the asymptotic decomposition theory and exists regardless of goals and ideology of the mathematical apparatus of quantum mechanics.
However the purpose of our research consists in finding of the generalized condition for existence of nuclear burning wave including both the necessary condition and sufficient condition.In this sense a question arises: "How will the nuclear burning wave velocity change, when the values of current concentration n f iss and critical concentration n f iss crit of nuclear fuel active components in Equation ( 22) change in a random manner?"It is obvious that in this case formulation of the problem (with consideration of Equations ( 34) and ( 68)) comes to finding the probability density distribution function p(a * ) of the random value of effective thickness of a nuclear burning concentration wave.This function, in its turn, is predetermined by a "random" nature of the equilibrium ñfiss and critical n f iss crit concentration of nuclear fuel active component, and thereby sets a random nature of wave velocity (see Equation (34)): And here is the key moment of the article.To avoid the use of random matrix theory by "brute force" in researching the random process in (22), we have applied its results in the form of quantum analogy stated above.In this sense there are not serious reasons to search any hidden physical sense of the found analogy.In other words, interpreting the results of present paper, it is not advisable to talk about some real quantum effects and, all the more, about observation of the so-called quantum chaos mode in our case.Indeed we show that the probabilistic nature of nuclear burning wave velocity (see Equation (45) and Figure 4) is immediately predetermined by random nature of the initial values of the equilibrium ñfiss and critical n f iss crit concentration of nuclear fuel active component.In this connection, a natural question arises, how in such a neutron-multiplying system, on the one hand, dynamic and chaotic modes coexist and, on the other hand, the chaotic mode in the form of the Wigner statistics is observed, although we have considered, at first sight, absolutely different fuel media, which in addition have the different initial fuel composition (see Figure 4).
Indeed it is known that the dynamic and chaotic modes can be observed under certain conditions in the same nonlinear dynamic system.This is easy to show on an example of the evolution of motion phase trajectories of the following hypothetical system.
Let us imagine some cylinder volume which is 238 U filled.We mentally divide it into identical disks.Each of these disks is described by the eigenvalue of equilibrium concentration ñfiss and critical concentration n f iss crit of 239 P u under the obligatory condition ñfiss > n f iss crit .More specifically, each of the disks is described by the eigenvalue a * n = ∆ • n (where n = 0, 1, 2, . . ., n max = 3/∆ (see Figure 4) and ∆ is the subinterval of the segment [0, a * max = 3]).At the same time, the disks with the different values ałn are equally spaced in a cylinder.Then for the traveling wave of nuclear burning the dynamic mode takes place in each of disks, but when the wave passes from one disk to another, the change of the dynamic modes which are characterized by the different pair of equilibrium ñfiss and critical n f iss crit concentration of 239 P u takes place due to the chaotic mode.In this sense, such an example visually demonstrates the random nature of changes of the different modes of motion phase trajectories, when the chaotic mode (strange as it may seem) plays the role of constructive chaos [30], which "instantly" transfers the system from one dynamic mode to another.
Another question related to the different fuel media and their initial composition can be reformulated in the following way: "What does the nuclear burning wave velocity mainly depend on in the uranium-thorium cycle (1) and thorium-uranium cycle (2)?"An answer is obvious and rather simple.In both cycles the nuclear burning wave velocity (far from a source which initiates the firing process) is fully characterized by the equilibrium ñfiss and critical n f iss crit concentration of nuclear fuel active component.First of all, this is explained by the fact that the equilibrium ñfiss and critical n f iss crit concentration of nuclear fuel active component completely identify the neutron-multiplying properties of fuel medium, because they are the conjugate pair of integral parameters, which due to their physical content completely and sufficiently characterize all physics of nuclear transformations predetermined by initial fuel composition.This also follows from the simple analysis of the solutions of the system of kinetic equations for neutrons and nuclei ( 22)- (25).From this it follows that regardless of the type of nuclear cycle and initial fuel composition the nuclear burning wave velocity will be determined by the values of the equilibrium ñfiss and critical n f iss crit concentration of nuclear fuel active component and, consequently, as calculating experiments show (Figure 4), will obey the Wigner statistics.

Conclusions
The solutions of the system of diffusion type equations for neutrons and concomitant kinetic equations for nuclei obtained by numerical 3D-simulation persistently point to the regions, where the stable soliton-like solutions for neutrons and solitary wave solutions for nuclei are existed.This is no wonder for nearly intergrable systems, to which the investigated system of equations for neutrons and nuclei belongs, whereas the existence of stable soliton-like solutions in three spatial dimensions causes a surprise for the following reason.
As is known, the derivation and solution of integrable nonlinear evolution partial differential equations in three spatial dimensions has been the Holy Grail in the field of integrability since the late 1970s.The celebrated Korteveg-de Vries and nonlinear Schrödinger equations, as well as Kadomtsev-Petviashvili and Davey-Stewertson equations, are prototypical examples of integrable evolution equations in the one and two spatial dimensions, respectively.Do there integrable analogs of these equations in three spatial dimensions exist?
As it has turned out, quite recently, in 2006, the method for finding of analytical solutions of indicated above partial differential equations in three spatial dimensions was developed [31].Therefore, the natural question arises: "To which from this equations does the diffusion equation for neutrons correspond, or, maybe, his is perfectly a new type of soliton partial differential equations in three spatial dimensions?"

Figure 1 .
Figure 1.Time dependence of neutron concentration.Propagating wave (a) and locked wave; (b) a segment of the curve n P u (z) above the n cr line is the reactor core; the scales of n cr and n P u are given with a × 10 magnification [8]; t is expressed in arbitrary units.

Figure 3 .
Figure 3. Schematic description of the permitted and forbidden region boundaries of nuclear burning according to the Borh-Sommerfeld condition (a) and the corresponding quasi-equivalent two-level scheme (b).

Figure 4 .
Figure 4.The theoretical (solid line) and calculational (points) dependence of Λ(a * ) on the parameter a * .

Figure 5 . 25 Figure 5 .
Figure 5. Concentration kinetics of neutrons, 238 U , 239 U and 239 P u in the cylindrical reactor core with radius of 125 cm and 1000 cm long over the time 240 days.Here r is transverse spatial coordinate axis (cylinder radius), z is longitudinal spatial coordinate axis (cylinder length).

Figure 6 .
Figure 6.The concentration distribution of neutrons at the velocity of wave propagation u simul ≈ 2.77 cm/day (a) and 239 P u at ñPu = 0.1 and n P u crit = 0.0167 (b) along the axis of a cylinder over the period t = 217 days.

Table 1 .
The parameters of nuclear burning wave.Forecast for the Th-U fuel cycle in infinite medium at the 10% enrichment of 233 U . *