This paper is devoted to the solving of a number of important problems of neutron kinetics. The most significance here is that problems connecting with questions of neutron multiplication in the nuclear fission process are presented. Special interest induces the kinetic analysis of delayed neutrons in the general neutrons balance. In the work, determining influence of delayed neutrons on the average life time of one neutron generation in considered nuclear device (nuclegen) is noted. Author’s works  ,  ,  ,  also are devoted to the theoretical grounds of the nuclegen creation.
With the implementation of the project on the creation of toroidal nuclear electric generator the important is the question about the diffusion of neutrons and the question about spreading of diffusing the flow (flex) of neutrons in this nuclear device  ,  ,  .
The main attention in Sections 2-5 is given to the solving of neutron kinetics problems in the general view. With that purpose the conception of the matrix of complete neutron kinetics is introduced. Given matrix allows finding an overt form for solutions of kinetic equations by means of employment: 1) of the separation method for kinetic equations and 2) of the reducing method of integrodifferential equations equivalent for kinetic equations to corresponding integral equations and its solutions.
In Sections 6-8 the problem of neutron diffusion in a nuclegen is investigated. It is succeeded in solving of diffusion equation with the aid of Fourier’s separation method of variables. Obtained solutions allow confirming about the practically complete absence of some penetrating diffusion.
2. Problem of Neutron Kinetics. Basic Concepts and Preliminary Clarifications
In the process of nuclear fission two chief neutron groups are formed. In the first one of these groups neutrons are emitted almost instantaneously (within 10−14 s from the very process of fission) and because of that are called prompt neutrons. Prompt neutrons make about 99.36% of the total quantity of emitted neutrons. The second group is formed of the delayed neutrons, which come from the splinters of fission after their radioactive β-decay after some time, measured by fractions and tens of seconds.
10 known precursors of delayed neutrons can be pointed out  ,  and besides each one of them produces neutrons with a corresponding period of half-life and energy, which is lower than that of prompt neutrons. Of these precursors six  ,  ,  , having the largest output of delayed neutrons, are usually selected to be taken into account due to difficulty of analysis.
Under a stationary operating condition of a nuclear device (reactor, generator-nuclegen) a critical state of a system, when the multiplication factor of neutrons , corresponds with a self-sustaining fission reaction. Understandably during operation of the device within the limits of this condition the factor k is to be changing within the neighbourhood of one. Because of this, a new coefficient is introduced to describe the operation state of a nuclear device
and it is called excess multiplication factor or excess reactivity. It is obvious that the device remains in the critical state if , in supercritical if and in subcritical if .
Apart from the values k and in order to characterize state of the active zone a value
is often used, which is called reactivity and signifies the fraction of the change of neutron density in a new generation.
The average neutron generation time I is the time from the moment of the neutrons formation during fission to the moment of their absorption (with the challenge of a new act of fission or ultimately output from the reaction). The rate of increase or decrease of neutron density n (number of neutrons N in a unity of volume) in the active zone depends on the average generation time I of neutrons. Naturally, the less is the generation time I until capture with fission, the more is the rate of increase of density n, that is
This expression can be given a more precise form. Indeed, in the next generation at the presence of excess reactivity the neutron density is to be equal . Because of this, the rate of change of neutron density without taking delayed neutrons into account obeys the equation
from which , in which is neutron density at the initial time.
In kinetic analysis of nuclear devices the concept of period T of a given device is introduced, so that is the time, in the duration of which the neutron density n changes e times, that is
and, consequently, .
If delayed neutrons are to be taken into account in the general kinetics equation, then the following neutron balance is achieved. Specifically, for groups of delayed neutrons the rate of change of neutron density equals
in which the first summand in the right-hand side of the Equation (1) signifies the rate of formation of all neutrons by fission, the second one signifies the rate of neutron loss, the third one signifies the rate of formation of precursors, the fourth one signifies the rate of formation of delayed neutrons out of these precursors, the fifth one signifies the intensity S of an external source of neutrons.
Here β is the portion of delayed neutrons, is the portion of prompt neutrons, is a constant of decay of fission splinters (precursor nuclei), which emit delayed neutrons of group number i, is the lifetime of these precursor nuclei, is the number (concentration) of precursor nuclei of delayed neutrons of group number i.
The quantity satisfies the equation
in which is the portion of delayed neutrons of each group. In the Equation (2) the first summand in the right-hand side gives the decrease of the number of precursor nuclei due to decay, and the second one gives their formation. Also notable is that for the substantiation of the system of Equations (1) and (2) it is supposed that the lifetime of daughter delayed neutrons and their precursor nuclei coincides.
Obviously, Equation (1) can be written dawn in the form
if it is considered that . The system of Equations (2) and (3) describes the behavior of neutron density in the time and because of that it may be admitted as the most important in the task of kinetic neutron research. From that system it is clearly seen that the number (concentration) of precursor nuclei , forms feedback in regards to neutron density.
For simplification and greater visualization of the calculations an averaged approximation of one group of delayed neutrons with an averaged lifetime is typically used for kinetic analysis, in which
Then, obviously, the average lifetime of a generation of neutrons can be found as a sum of average lifetimes of delayed and prompt neutrons taking into account their portions in the total composition of neutrons:
In the expression (4) there is an estimation
from which follows an important ascertaining: an average lifetime of one neutron generation in a nuclear device is practically defined by the lifetime of delayed neutrons.
3. Matrix of Full Neutron Kinetics
Write a system of differential Equations (2) and (3) in vector-matrix form
where designations are introduced
For a constant square matrix A of dimension the Equation (5) has a single solution . However that solution cannot be used to write down the explicit form of all the solution components x because of the presence of a matrix exponent.
Matrix A can be called a matrix of full neutron kinetics. The calculation algorithm of the determinant of this matrix is adduced.
Theorem 1. Determinant of matrix A is calculated by the formula
Proof. With the goal of deriving the correlation (6) we shall utilize the expansion method of the determinant on columns with highest number of zeroes and reduction to calculation of lesser order determinants of triangular form. Then, starting with the second column, on the first step we get
Following this, in the equality (7) on 2-nd step the is again expanded along the elements of the second column:
This means that
Substituting the expression (8) into the formula (7) we get
Continuing this process and successively reducing it to calculation of and etc we come to the equality (6) on the final m-th step:
The Theorem 1 is proven.
Consequence. The expressions for are substituted into formula (6). Then, obviously, we get
From here an important conclusion ensues that the matrix A is degenerate for , that is, for the critical state of the nuclear device. In that case there is a stationary operating regime (rate of change of neutron density and their precursors equals zero). Then Equations (2) and (3) take on the form ( ) respectively:
in which . From the first equation . If this value of is substituted into the second equation of the system (9), then we come to an identity, which proves the rightness of the substantiation of initial conditions data of the process of nuclear (neutron) kinetics.
Considered further is a task of finding the explicit form of the solution of the Equation (5) while accounting for the results of Theorem 1. Apparently all elements of the matrix exponent need to be found for that. As it is known (refer, for instance, to the works  ,  ), this task is solved in a general case by means of finding eigenvalues and eigenvectors of the matrix A.
Various eigenvalues of matrix A are designated by . The left-hand side of the characteristic equation , where I is the identity matrix, for determination of takes the form of a polynomial by of degree , the coefficients of which (and thus also the roots) are rather difficult to calculate. As a vivid example polynomials for various groups of delayed neutrons are written down.
1) One-group kinetic model:
2) Two-groups kinetic model:
3) Three-groups kinetic model:
From here the search for an analytical solution with the help of the algorithm of finding roots of the equation in a general case is revealed as prospectless. In the meantime, the situation is not as hope less and changes in a radical way if the system Equation (2) and scalar Equation (3) are separated.
4. Separation Method. Integrodifferential Task
System (2) can be rewritten in the form
As a reminder, the initial conditions are defined by the expressions (9).
Integrating the Equation (10), we get
Because in the solution (11) the value is scalar, and the vector is constant, we have
In the expression (12) the matrix exponent is very simply found in the form
and that means, that the solution (12) can be written down element wise:
Then the obtained values (13) are substituted in the Equation (3), and we get
Obviously, this equation is representable in the form of a linear nonhomogeneous integrodifferential Volterra equation of 1-st order
where for the known functions of time and designations are introduced ( ):
Integrodifferential equations of the form (14) apparently have first appeared in the tasks of mathematical biology (20-s and 30-s of 20-th century: Verhulst’s task about populations and competition within a species, Lotka-Volterra task about “predator-prey” type associations and etc.  ). A feature of these equations is that all of them are caused by various effects of delay. Thus it is not surprising that the presence of delayed neutrons in the process of nuclear fission also caused a corresponding integrodifferential equation in its wake, in which the effect of the delay is given with the help of the function .
For the solution of the Equation (14) one of the approximate methods can be utilized, for instance, one of quadratures, one of successive approximations and etc  . The method of reduction of the initial equation to an integral one with consequent solution through a resolvent kernel is among the more redical analytical methods of solution of integrodifferential equations.
5. Method of Reduction to an Integral Task
However in our case the reduction of the initial integrodifferential equation to
an integral equation leads to big difficulties due to the presence of three constituents of various differential dimensions at once: and .
Nonetheless, let’s try to utilize one of the variants of the method of reduction  with the help of a fitting change of variables.
in which is the new sought function, is a function of a homogenous kernel, which is hereinafter selected from certain conditions of solvability of the equations. Substituting all of these expressions into the Equation (14) we get
Let’s further designate
and write the Equation (15) down in the form of a linear Volterra integral equation of the second kind:
in which hereinafter the kernel is to be selected appropriately.
Provided identity of the Equations (15) and (16) there is a need for the equality
to hold. This correlation holds, in turn, if the following equality of integrands holds:
that is, if there might be the equality
Obviously, in the integral condition (17) the first fraction should not depend on t. This happens, for example, if the function selected as the function is such
in which is a certain known function of time s. Thus, we come to a restriction
from which we get an equation of a linear holonomic differential constraint of 1-st order:
This constraint (18) may also be called a regulated one, because in it the functional coefficient may be selected by some purposeful way.
The simplest and most natural case of the solution of the Equation (18) is considered, in which . We have , in which . In connection with the presence of a series of new introduced functions (which are supposed to be continuously everywhere differentiable with respect to all arguments) and of the corresponding equations an important question of their solvability arises. The following statement in the form of the Theorem of matching provides an answer to it.
Theorem 2. For the selection of a function of time
in which is a certain given function is a known solution of the integral Equation (16), are known functions of time, the system of integral dependences (16)-(18) is a concordant one, that is consistent and possessing an unambiguous solution for given initial conditions.
Proof. Let’s proceed primarily from the presence of the Equations (16) and (18). On one hand, the function being given fully defines the solution in the integral Equation (16). On the other hand, the solution in the Equation (18) defines the nature of the dependence
if the function of time is known.
This is handled in the following way. Let . Then in the correlation (18) we have
Due to the definition of the function we need to write down
Functions and are assumed to be known functions of time. Let’s assign the function . Then from the correlation (20) we get
in which is an already known function of time. What is more, we have
in which is a known solution of the integral Equation (16). Substituting the expression for from the correlation (22) into the dependence (21), we find
from which the condition of matching (19) follows, quod erat demonstrandum.
In a general summary of the task of solving kinetic equations, it can be concluded that the integral scheme of solution consists of several stages: firstly, we define the function and the number , then we find the functions and , by them we determine functions and . The found values allow determining the solution of the task, specifically the neutron flux density and concentration of precursor nuclei .
6. Diffusion Problem of Neutrons and Its Solution by Furier’s Method
It is well known the process of formation of free neutrons inside a fission medium (thereby of formation of charged fission particles) for realization and maintenance of nuclear chain reaction is a decisive significance. Therefore, the question about diffusion of neutrons and charged particles flex over the whole volume of external magnetic reflecting field which created by toroidal nuclear electric generator (nuclegen) is presented sufficiently important.
In critical stationary generator the coefficient of neutron multiplication k from generation to generation doesn’t change . Hence it follows that for any point of active zone we can write the equation of neutron balance
where F is the “formation”, L is the “leakage”, A is the “absorption” and R is the “reflection”. Here through R in relation (23) the external magnetic pressure is indicated which is making in respect to the neutron flex the reflected and restored function and reflected effect of itself toroidal surface also.
The basic equation for diffusion effects (in this case neutron diffusion) is the Fick’s law, which establishes that resulting current (flex) of diffusive substance is proportional to the density gradient of this substance and is directed to the side of region with the lowest density
where is the operator of spatial gradient, J is the resultant flex of neutrons, is the scalar neutron flex, D is the coefficient of diffusion is determined by relation (24). Here n is the number of neutrons in the unit of volume (the spatial density), v is the average neutron density.
We denote the leakage L in Equation (23) through L too. The value L is express itself in the equality
where is the Laplace operator. If to count up that in Equation (23) “formation”, “absorption”, “reflection” each separately are proportional to the flex then taking into account relations (24), (25) the Equation (23) can be written for some proportionality coefficient R in the form:
where is the boundary condition on the from some region concerning the function . We solve Equation (26) with a given boundary condition (from region ) into supposition that the flex is twice continuously differentiable on spatial variables. With that end in view we introduce cylindrical coordinates . We have for the arbitrary function :
Owing to the symmetry of toroid the flex doesn’t depend on the angular coordinate and the diffusion problem (26) becomes two-dimensional: , i.e.
where . We use the Furier’s method variables separation for the solving of Equation (27) supposed that the flex is presented in the form of two functions product as
We substitute expression (28) in the equation of diffusion (27). After termwise division on we receive the equation
The left part of Equation (29) depends only from r and the right―only from z. The function (28) can be the solution of Equation (29) as soon as in the case if both parts of this equation are constant by the arbitrary changing of variables r and z, i.e. must be carried out the following relations
where are positive numbers connected by equality owing to Equation (29).
7. Solution of Equation on R(r)
Consider Equation (30). At first we multiply it by :
and then introduce the new variable . As a result we obtain the Bessel’s equation of zero order
The general solution of Equation (32) is presented in the form
where are the arbitrary constants, are the Bessel’s functions of zero order of first and second kinds determined with the help of power series
Here , is the Bessel’s function of n-th order with the constant , i.e.
where is the Euler gamma-function. It can be shown, that for small :
where is the Euler’s constant. Sometimes the function is called the Neumann’s (or Weber’s) function of zero order. If come back to the variable r then we can write as a result
8. Solution of Equation on Z(z) and the General Solution
Consider now Equation (31). Its solution is determined easily
where are the arbitrary constants. We take from conditions of symmetry
Then it is obvious that in solution (34) by conditions (35) we have and consequently
Multiply these two solutions (33) and (36). As a result we obtain the expression for flex (28):
Solution (37) should satisfy the basic boundary condition:
is the δ-toroid, and axial conditions of symmetry
It is obviously that conditions (38), (39) correspond to the neutron flex for which
The initial diffusion Equation (26) with boundary (16) and axial (17) conditions corresponds the flex , with components
The flex (40) corresponds to the absence of leakage in a stationary nuclear electric generator, to the presence of the zero parameter in Equation (26), and to satisfying of the neutron balance equation:
The provision of condition (41) means an important feature and advantage of toroidal vacuum-nuclear system―the almost complete circular motion of the streams within the toroidal active zone of the generator.
In condition (38) there is a parameter δ called the extrapolation length. We shall say a few words about his physical sense. It turns out that the neutron flex does not vanish on the boundary or the edge of the medium in which the reaction takes place. In a toroidal nuclear device of finite dimensions there is always, though insignificant, a leakage through this external boundary. This means that in the immediate vicinity of the boundary inside of active zone there should be neutrons that leave its limits. At the boundary, some of the neutrons leave outside and the other part diffuses inside and increases the value of the neutron flex. As a result of this, the neutron flex on the boundary doesn’t equal zero but is supposed equal zero from the outside on some little distance δ.
Let’s consider this case in more detail. Suppose that the boundary condition (38) is satisfied, and on the toroid itself T there is a small neutron leakage, i.e. , but
where is a small numerical parameter.
Then, obviously, the small external transverse diffusion in Equation (26) corresponds that , where is some small parameter. In fact, this means that
where are the coefficients of flex proportionality in the “formation”, “absorption”, “reflection” of neutrons, respectively. Hence we have are small parameters.
Further we can write for the solution (37):
Since for the Bessel and Neumann functions of zero order we have
then it is naturally we can choose: .
Thus, the case of small transverse neutron diffusion with the indicated characteristics in the toroid is quite probable and can be considered as a kind of nuclear analogy of Arnold diffusion  across the stochastic layer.
In conclusion, we investigate the nonstationary case by limiting to a simple situation in which there are no delayed neutrons and . The change in the neutron density with time is described by the relation , where is the neutron density in the initial moment of time, is the period. We assume from this that the neutron flex can be represented in the form
where is a continuously differentiable function of time. Then the general diffusion equation can be written as follows
where the rate of change of the nonstationary neutron flex is proportional to the total “penetrating diffusion” over the lifetime of one generation T. Here is the effective coefficient of reproduction.
It remains to substitute into Equation (42) and take into account that . We have
Note that due to the assignment of the flex takes place . However, it is not difficult to show that on the basis of the previous analysis, that in the case where w is the spatial variable, we have for the toroid . Thus, we can conclude from this that the cases of neutron leakage from the active toroidal zone are excluded. The same reasoning can also be given for an ionite gas. Other coefficients and notations appear, but the essence, must be assumed, remains the same: the ionites in the toroid will practically not be subject to diffusion.
First of all, it should be noted once more that in this paper the problem of neutron multiplication in the nuclear fission process is solved in general form. It may be achieved owing to: 1) the introduction of the matrix of full neutron kinetics and the investigation of its algebraic properties; 2) the reducing of integrodifferential task to integral task for the effective solution of kinetic neutron equations.
Flex (40) corresponds to the absence of the leakage in stationary nuclear electric generator, the availability of zero parameter in Equation (26) and the carrying-out of neutron balance equation: . The fulfillment of this condition means the important peculiarity and advantage of toroidal vacuum-nuclear systems, expressed into practically complete (without losses) circular motion of nuclear flexes inside of the toroidal active zone of nuclegen.
 Tertychny-Dauri, V.Yu. (2017) Diffusion of Neutrons in the Toroidal Nuclear Electrogenerator. VI International Conference on Mathematical Modeling in Physical Sciences, Pafos, Cyprus, 28-31 August 2017, IOP Journal of Physics: Conference Series, 936, Article ID: 012007.
 Tertychny-Dauri, V.Yu. (2018) Nuclear Electromagnetic Generator: Introduction in Charge Algebra and Elements of Charge Kinetics. Open Access Library Journal, 5, e4614.