SELECTED
TOPICS in REACTOR PHYSICS II.
From Budapest University of Technology and Economics
The adjoint function and its applications
Neutron field and detector field
The adjoint transport equation
The adjoint diffusion equation
2. Time dependence of the adjoint function
The Lagrangian of transport theory
Perturbation of the reactivity
The characteristic of a control rod
5. The point kinetic equation and the effective delayed
neutron fraction
Neither the transport equation nor the diffusion equation is selfadjoint – in contrast, for example, to the Schrödinger equation. It would not be difficult to formulate the equations adjoint to the transport or diffusion equations (at least formally). However, we are interested in a more fundamental question: do the solutions of the adjoint equations have some physical meaning? We shall see that the answer is: yes! As a matter of fact, a proper understanding of the transport equation is not possible without understanding the adjoint equation. The physical meaning can be best understood in the framework of the simplest model e.g. the onegroup diffusion equation. This is the only approximation of the transport equation which selfadjoint.
Consider a bare critical reactor in which the neutron flux is constant in time. This flux is proportional to F_{1}(r) that is the eigenfunction of the Helmholtz equation:
_{} 
(1) 
belonging to the smallest eigenvalue _{}. Assume that a neutron appears at the point r_{0} in the moment t = 0. After this, the flux will become time dependent. The initial condition for it can be written in the form
_{} 
(2) 
where c is some constant. For later moments t > 0, the time dependence of the flux is determined by the timedependent onegroup diffusion equation:
_{}. 
(3) 
We write its solution in the form of the series
_{} 
(4) 
where the functions F_{n}(r) are the eigenfunctions of Equation 1 (n = 1, 2, ...). The coefficients f_{n} are the expansion coefficients of the initial distribution given by Equation 2:
_{}. 

For a critical reactor, w_{1} = 0 and w_{n} < 0 for n > 1. Thus, for sufficiently large times t, the flux goes over in a constant value again:
_{}. 
(5) 
We have found that the appearance of a single neutron increases the steadystate flux by a value proportional to F_{1}(r_{0}). We can put it also in an other way: for the chain reaction, it is the worth of a neutron injected at point r_{0}. This is the idea of neutron importance. Equation 5 shows that, in onegroup diffusion theory, the importance is equal to the flux.
In the next section, we shall derive the equation satisfied by the neutron importance. We shall find that this equation is the adjoint of the transport equation. In the present section, our starting point was the onegroup diffusion equation (see Equation 3). All operators in the right hand side are selfadjoint. Thus, it was not by mere chance that the neutron importance was found to be equal to the neutron flux.
There are problems in transport theory which, principally, could be solved within the framework of transport theory, but, practically, they cannot be solved without the adjoint equation. An important set of such problems is the treatment of small perturbations: we change some characteristic of the reactor by some small amount (for example, we move a control rod a little), and we want to determine the corresponding changes of the reactivity and the neutron flux. Since we want to predict the response of the reactor to small changes, such problems can not be treated by the usual numerical methods. Let us consider the example of the reactivity change. Denote by _{} and _{} the values of the multiplication factor corresponding to the initial and final positions of the control rod, respectively. Then the reactivity changes by
_{}. 

The usual order of magnitude of Dr is 10^{–5} and we want to get it as the difference of two numbers which are both close to 1. Any numerical solution of the transport equation is only an approximation. We would get a rather poor accuracy by calculating the small difference of two large quantities. In such cases, it is always better to look for such methods which, although approximate as well, give the difference directly. This is the basic idea of perturbation theory. It was first elaborated for the Schrödinger equation. Later, it was developed for the transport equation as well. However, there is a difference with respect to the Schrödinger equation: it is necessary to solve both the transport equation and its adjoint.
Although the flux F is calculated for most practical applications, the solutions of the adjoint equation can be interpreted consequently only if we start out from the transport equation formulated in terms of the neutron density n. Onegroup theory is an exception: the neutron speed v is not a variable of the equations, thus, the flux and density are proportional to each other.
Since – as far as possible – we do not wish to repeat our considerations for all approximations of the transport equation, we put our starting equations in an operator form:
_{} 
(1.1) 
where
_{}. 
(1.2a) 
The operators are defined for the transport equation as follows:
_{} 

_{} 
(1.2b) 
is the destruction operator describing the disappearance of the neutrons while
_{} 
(1.2c) 
is the production operator describing the neutrons produced by fission. In Equation 1.2b, the kernel S(r,E˘®E,_{}) is the sum of elastic and inelastic scattering kernels. These operators have the following form in diffusion theory:
_{} 

_{} 
(1.3a) 
and
_{} 
(1.3b) 
where – as before – the kernel S(r,E˘®E) is the sum of elastic and inelastic scattering kernels. The solution of Equation 1.1 has to satisfy the boundary conditions belonging to the model chosen. The initial condition is written in the form
_{} 
(1.4) 
where the function n_{i} has to satisfy the same boundary conditions as the neutron density n itself.
There is no instrument which would allow the measurement either of the neutron density or of the neutron flux. All what can be done consists in putting some neutron detector in the reactor and in counting the nuclear reactions induced in it by the neutrons. The result of the counting will be called meter reading in the following. If the cross section of the detector is S_{d}(r,E) for neutrons of energy E at point r of the reactor, the number of nuclear reactions in unit of time, in volume element dV at r is given by the integral
_{}. 

The total meter reading at time t is then
_{} 

_{} 
(1.5a) 
where V is the volume of the reactor. In diffusion theory, the same thing can be written as
_{}. 
(1.5b) 
We cite two examples for explaining what we mean by “detector”. If it is a BF_{3} counter, S_{d} is proportional to the absorption cross section of ^{10}B, and the integrals for r in Equations 1.5 are extended only to the detector’s volume since S_{d} vanishes outside the detector’s volume. Each neutron absorbed in the detector results in an electric pulse, thus, the meter reading is obtained by counting the pulses. Another example is the reactor power: S_{d} is then proportional to the fission cross section of the reactor fuel, thus, the integrals in Equations 1.5 are extended to the whole volume of the reactor core.
It follows from these considerations that the observable quantities cannot be described by the neutron field alone but we need another field, too, that is the detector field. Any function belonging the detector field will be identified by a superscript “+” such as y^{+}. The variables of the detector functions are the same as those of the model actually used. This means that we have y^{+}(r,E,_{},t) and y^{+}(r,E,t) in transport theory and diffusion theory, respectively. The physical meaning of the former is the contribution of a neutron characterized by (r,E,_{}) at time t to the meter reading measured at time t.[1] The meaning of y^{+}(r,E,t) is analogous. According to this definition, we have in Equations 1.5
_{}. 

The integrals given in Equations 1.5 are inner products of the detector function and the neutron density. We use the same notation for them as in functional analysis:
_{}. 
(1.6) 
Variables r, E, and _{} are omitted in formulae like this but the time variable t is given when, in the actual context, it is necessary to specify for which moment of time the inner product is defined.
The detector functions y^{+} may be arbitrary in principle but we suppose them to satisfy some boundary conditions which follow from their physical meaning. Since a neutron flying outwards at the outer boundary of the reactor gives no contribution to any meter reading, we require:
_{} 
(1.7a) 
where N_{s} is the normal vector of the boundary pointing outwards (Figure 1.1). Note that this is just complementary to the condition set for the neutron density:
_{}. 
(1.7b) 
In diffusion theory, for the same reason, we require that y^{+}(r,E,t) vanished at the extrapolated reactor boundary. Note that the same boundary condition is set for the neutron density.
Figure 1.1. Illustration of the boundary condition set for the detector functions in transport theory
Assume that a measurement will be performed at t = t_{f} by using a detector with detector function _{}.[2] The meter reading will be by definition
_{}. 
(1.8) 
In order to calculate n(t_{f}), we need to know the initial neutron density n_{i} for some time t_{i} < t_{f}. Equation 1.1 determines how the neutron density develops between the two moments of time. Now we formulate two definitions:
· For t < t˘, consider a neutron in both of these moments of time. The latter is called a descendant of the former if its birth can be traced back to the former through translation, scattering and fission.
· We call n^{+}(r,E,_{},t) neutron importance if it gives the total contribution of all descendants of a neutron characterized by (r,E,_{}) at moment t to the meter reading measured in the moment t_{f} (t < t_{f}).
It follows from the definition that the importance has to satisfy boundary condition 1.7 as well as the following “final condition”:
_{}. 
(1.9) 
The neutron importance n^{+}(r,E,t) of diffusion theory can be defined in an analogous way.
The neutron importance defined in the introduction is a special case of this general definition. As the relation between the two definitions is not trivial, we give the following explanation. Firstly, we used onegroup theory there, thus, variables E and _{} are omitted. Secondly, we considered the solution of Equation 4 in the limit of t_{f} ® Ą. Thirdly, we considered the neutron flux at point r. Equation 1.8 gives it for
_{}. 

In order to get an equation for the neutron importance, we apply the same reasoning which leads to the differential form of the transport equation. Our equation will express the following obvious statement: the importance of a neutron is equal to the sum of the importances of its descendants. We can call it the principle of importance conservation.
Let a neutron be at point r at time t and having energy E, flying in direction _{}. After a time lap of dt, we find it at point (r + v_{}dt), unless it induces some nuclear reaction. The probability of the former eventuality is (1 – vS_{t}(r,E)dt). According to the definition of the kernel _{} (see Equation 1.2), the probability that, following an (elastic or inelastic) scattering, the energy of the neutron changed to a value in the interval (E˘, E˘+dE˘), its flight was deflected to a direction in the solid angle d_{} around _{}, is given by


Finally, the total number of fission neutrons of energy in the interval (E˘,E˘+dE˘), of flight direction in the solid angle d_{} around _{}, is
_{}. 

It follows from the conservation principle of the importance that
_{} 

_{} 

_{}. 

We get in the limit of dt ® 0:
_{} 


_{}. 



If the total time derivative is expressed by partial derivatives, we get an equation similar to Equation 1.1:
_{} 
(1.10) 
where
_{} 

_{} 
(1.11a) 
and
_{}. 
(1.11b) 
We have got an equation for n^{+}, in which time t passes “backwards”: if we want to compute n^{+}, we start from the final condition specified by Equation 1.9 for the moment t_{f}, and we solve Equation 1.10 backwards in time.
We will show that the operators _{}, _{}, and _{} defined by Equations 1.11 are adjoint to operators _{}, _{}, and _{} of Equations 1.2. This means that we have to prove the validity of the following equations:
_{}, 

_{} 

where n and y^{+} are arbitrary functions of the neutron and detector fields, respectively, which satisfy the respective boundary conditions. If these equations hold, it follows from this that
_{}. 

Let us look at the operators _{} and _{} terms by term:
· The kernels of the integral operators are the same except that variables E and E˘ exchanged: we have _{} in _{} (see Equation 1.2b) while the corresponding term of _{} is _{}. Since _{} and _{} appear symmetrically, this means that the integral operators are adjoint to each other. As an exercise, we advise to check the same thing for operators _{} and _{}.
· The term containing S_{t}(r,E) is trivially selfadjoint.
· The leakage terms are identical except the sign. We prove our statement by taking into account the following identity:
grad(gh)
= g×grad
h + h×grad g 

which holds for any differentiable functions g(r) and h(r). We may now write
_{} 

_{}. 

We apply the Gauss theorem in the first term of the right hand side:
_{}. 

This surface integral vanishes for every _{}: if _{} points outwards (see Figure 1.1), we have y^{+} = 0 by Equation 1.7a, if _{} points inwards, we have n = 0 by Equation 1.7b. This means that
_{}. 

This completes the proof of that the adjoint operator of _{}grad is –_{}grad.
We can now state the following conclusion: the neutron importance satisfies the adjoint transport equation. For this reason, the quantity n^{+}(r,E,_{},t) is sometimes called the adjoint function.
The formal relation between the transport equation and its adjoint leads certain authors to call the function n^{+} adjoint flux. We do not subscribe to that practice. It follows from our derivation and the physical meaning of n^{+} that there are no particles whose transport would be described by the adjoint transport equation. In addition to this, we have a simple formal objection, too. We have seen that the functions of the detector field are the form vS i.e. their dimension is 1/s. At the same time, the dimension of the flux is 1/(cm^{2}×s). This formal reason alone is sufficient to ban the expression “adjoint flux”.
Their are three ways to derive the adjoint diffusion equation: firstly, we can give some physical meaning to the adjoint function n^{+}(r,E,t) and apply the same derivation as above; secondly, we can start from the adjoint transport equation and derive the adjoint diffusion equation just as the diffusion equation can be derived from the transport equation; thirdly, we can formally find the operators adjoint to those defined in Equations 1.3. We choose the third way which is the simplest.
We shall show that operators adjoint to Equations 1.3 are the following:
_{} 

_{} 
(1.13a) 
and
_{}. 
(1.13b) 
The statement is trivial except for the leakage term which is the first term of _{}. The proof will be based on the identity
_{} 

which holds for any (differentiable) function h(r) and vector g(r). Let n and y^{+} be arbitrary functions of the neutron and detector fields, respectively, which satisfy the diffusion boundary conditions. Then we may write:
_{} 


_{}. 



We have applied the Gauss theorem for the first term of the right hand side. The surface integral is defined for the outer (extrapolated) surface of the reactor where y^{+} vanishes. The value of the surface integral is thus 0. In the second term of the right hand side, we apply again the above identity and the Gauss theorem:
_{} 

_{} 

_{}. 

n satisfies the diffusion boundary condition, thus, the first term on the right hand side vanishes. We conclude now that
_{}. 

This means that the operator div[vDgrad] is selfadjoint as stated in Equation 1.13a.
As a conclusion of the present section, we give some hints concerning the derivation of the adjoint diffusion equation starting from the adjoint transport equation. When the diffusion equation is derived from the transport equation, the first step is the series expansion of the neutron density in terms of the spherical functions:
_{} 

Here, the diffusiontheoretical density _{} is the integral of the transporttheoretical density _{} with respect to _{}. Starting from this series development, we get Fick’s law:
_{}. 

In case of the adjoint function, we could formally generalize this as follows:
_{} 

If we put the series development in this form, _{} would be the integral of the transporttheoretical adjoint function _{} with respect to _{}. Unfortunately, this is not acceptable for two reasons. Firstly, in contrast to the neutron density _{}, the neutron importance _{} is not a density with respect to _{}: it is the importance of a neutron flying in some selected direction _{}. Thus, we cannot assign any physical meaning to the integral of _{}. Secondly, the transporttheoretical final condition 1.9 can be transformed into a diffusiontheoretical final condition if _{} is the average of the transporttheoretical adjoint function _{} with respect to _{}:
_{}. 
(1.15) 
In order to show this, assume that the detector function appearing in the final condition is isotropic. (This is the usual case.) Then Equation 1.9 has the form
_{}. 
(1.16a) 
Its diffusiontheoretical equivalent should trivially read as
_{}. 
(1.16b) 
Equations 1.16a and 1.16b can hold simultaneously if we use the definition 1.15.
These considerations lead to the series development
_{} 

If we apply it in Equation 1.10 and 1.11, we derive the adjoint Fick’s law:
_{} 
(1.17) 
and the adjoint diffusion equation corresponding Equations 1.13. When comparing Equation 1.17 with the Fick’s law related to the neutron density, we find a remarkable difference: the negative sign is absent in Equation 1.17. This can be interpreted as follows: vector J^{+} of the adjoint current points in the direction of the increase of neutron importance.
Two kinds of eigenvalue equations can be derived from the timedependent equation 1.1: kinetic (or dynamic) and the static eigenvalue equations. Both are obtained from Equation 1.1 by setting S = 0. This corresponds to a reactor left alone without external neutron source.
The kinetic eigenvalue equation describes a reactor which develops in a persistent mode i.e. in which the time dependence of the neutron density is separable from the dependence of the other variables:
_{}. 

By substituting this in Equation 1.1, we can easily show that the time dependence can be only an exponential function:
_{} 

with some suitable w satisfying the kinetic eigenvalue equation:
_{}. 
(2.1a) 
The static eigenvalue equation is the well known
_{} 
(2.2a) 
equation which is solved by many computer codes.
It is simple to write the analogous adjoint equations:
_{} 
(2.1b) 
and
_{}. 
(2.2b) 
The negative sign in the kinetic eigenvalue equation is due to the negative sign of the time derivative in Equation 1.10. If necessary, we apply a subscript “s” in order to identify quantities belonging to the static eigenvalue equation.
In all what follows, we try to keep general notations by using operators because this allows to make statements which are independent of the model actually used. If it is necessary to specify the other variables in addition to the time variable t, we shall use vector x for that purpose. Its meaning depends on the actual model. For example:
· transport theory: x = (r,E,_{});
· onegroup transport theory: x = (r,_{});
·
diffusion
theory: x = (r,E);
· onegroup diffusion theory: x = (r);
· multigroup diffusion theory: x = (r, g) where g is the group subscript,
and so on. When calculating inner products, we have to integrate for all components of vector x. For arbitrary functions y(x) and y^{+}(x), we calculate the integral
_{}. 

For multigroup models, the integral with respect to energy E is substituted by a sum with respect to the group subscript g. For example, we have in multigroupdiffusion theory
_{}. 

The eigenfunctions belonging to different eigenvalues are biorthogonal. In order to show this, we write the kinetic eigenvalue equations:
_{} and _{}. 

We take the inner products of these equations with the eigenfunction of the other equations, and we subtract the resulting equations:
_{}. 

If _{}, the inner product _{} vanishes. Under quite general conditions, it is possible to choose the subscripts such that, for j = l, the equality _{} holds. We can show the biorthogonality for the static eigenfunctions too:
_{} for _{}. 

Note that the orthogonality holds only when the production operator is included in the inner product.
Consider a subcritical reactor, which, as a special case, can be also a nonmultiplying medium. For such a system, all kinetic eigenvalues are negative: w_{l} < 0, for l = 1, 2, ... As a consequence of this, the neutron density tends to zero for t ® Ą independently of the initial density. A steady state neutron density can prevail only in the presence of an external source:
_{} 
(2.3) 
where S(x) is a steady state external neutron source.
The neutron importance has an analogous time behavior. Let us develop the importance in a series of the eigenfunctions:
_{} for t Ł t_{f.} 

By virtue of the final condition (see Equation 1.9), we may write
_{}_{.} 

By the orthogonality proven above, we can obtain the coefficients of this expansion:
_{} 

leading to the result
_{}. 
(2.4) 
Since all eigenvalues are negative, the limit of this importance is zero for t ® –Ą independently of the final condition. This behavior allows to get an equation for the neutron importance similar to Equation 2.3.
Assume that the external neutron source injects neutrons into the system since t = –Ą. We ask now the following: what will be the total contribution of the descendants of these neutron to the meter reading at time t_{f}. The neutron importance goes to zero for t ® –Ą, thus, the contribution of the descendants of the neutrons injected sufficiently early with respect to t_{f} will be negligible. Let the time interval (t, t+dt) be close to the moment t_{f}. The contribution of the descendants of the neutrons injected during this interval is by definition
_{}. 

If this is integrated for all moments, the result is the total contribution:
_{} 
(2.5) 
where
_{}. 
(2.6) 
It follows from Equation 2.4 that this time integral is convergent. The quantity defined by Equation 2.7 is called source importance.
Equation 1.10 allows the derivation of an equation for the source importance. If we integrate Equation 1.10 for t in the interval (–Ą, t_{f}], we get
_{}. 

By taking into account Equations 2.6 and 1.9, this may be rewritten as
_{}. 
(2.8) 
This equation is similar to Equation 2.3. We have thus found that there is an adjoint equation with source in which the “source” is the detector function.
We give an example of applying our last result. Assume that a spatially extended neutron source injects neutrons into the system, and we are interested in the radiation dose in some dose point outside the system. In terms of our former notations, S(x) ą 0 in a large volume while _{} = 0 if x is other than the dose point. Such problems of biological shielding are often treated by Monte Carlo programs. Unfortunately, the calculation will not converge since the probability of that a neutron hits the dose point is zero if the neutrons are started from the neutron source S(x). The problem is lessened a little by not considering the detector a geometrical point but extending it to a small volume around the dose point. Then the probability of hitting the detector is no more zero but it is still small. The real solution is the adjoint Monte Carlo method which allows to solve Equation 2.8 by Monte Carlo simulation. Although we emphasized in Section 1 that there are no real particles the flux of which could be the adjoint function, we can consider the source importance, at least formally, the flux of some imaginary particles whose transport is described by the operator _{}. Now we start the “adjoint particles” from _{} as the “adjoint source”, then the payoff is computed according to Equation 2.5 when the “diffusing adjoint particles” cross the real neutron source. Since the real neutron source is extended in space, the adjoint particles cross it with a finite probability. Of course, even the adjoint Monte Carlo method is helpless if the neutron source S(x) is also concentrated in a point.
In practical applications, it is the static eigenfunction that is computed most frequently. This means the solution of the static eigenvalue problem defined by Equations 2.2a and 2.2b. The static eigenvalues constitute a monotonously decreasing sequence, the first one being the effective multiplication factor (k_{eff}). The static eigenvalue equations corresponding to it are
_{} 
(2.9a) 
and
_{}. 
(2.9b) 
The subscript “s” means “static” equations. There are several computer programs that solve Equations 2.9 in various models (i.e. finite difference diffusion, finite element diffusion). Consequently, we may assume in the following that n_{s} and _{} are known in every case considered.
In addition to what was discussed above, the adjoint function plays a further role, too: it allows to formulate the transport equation (and all of its approximate versions) in form of variational principles. We apply variational principles frequently in theoretical physics. Their role is twofold: for one hand, they contain the basic equations in an concise form, and, for the other hand, they allow to compute certain quantities accurately (up to first order) by means of approximate solutions of the basic equations. Since the transport equation is not selfadjoint, the variational principles necessarily contain booth the neutron density and the adjoint function.
In transport theory, the Lagrangian is the meter reading defined by Equation 1.8. Let n^{+} be a (for the time being arbitrary) detector function which satisfies the boundary condition formulated in Equations 1.7. Rewrite R(t_{f}) as follows:
_{} 

where t_{i} is the moment to which initial condition 1.4 applies. R(t) for t_{i} < t < t_{f} is given by the formula
_{}. 

If Equation 1.1 is also taken into account, the Lagrangian becomes
_{} 

_{}. 
(3.1) 
According to the variational principle, we ask the question: how to choose n^{+} in order to make R(t_{f}) stationary against small variations of the neutron density? In order to find the answer, we change n by dn. The neutron density is given for t = t_{i} (see Equation 1.4), consequently, dn(t_{i}) = 0, leading to
_{}_{.} 
(3.2) 
This vanishes for all dn if
_{}. 

In other words: the Lagrangian is stationary if n^{+} satisfies Equation 1.10 i.e. the adjoint equation.
We have found that the variational principle based on the Lagrangian 3.1 contains the adjoint equation. We can simply see that this principle contains the transport equation 1.1, too. We integrate in Equation 3.1 by parts:
_{}. 
(3.3) 
Now we take the variation of this with respect to n^{+}:
_{}_{.} 
(3.4) 
This is stationary for all dn^{+} if n satisfies Equation 1.1. Furthermore,
_{} 

needs to hold as well, i.e. n^{+} has to satisfy the final condition 1.9.
As a conclusion, we can state the following: the Lagrangian 3.1 is stationary with respect to both the neutron density and the adjoint function, and both equations can be obtained from the variational principle based on it.
The most important application of the variational principle formulated above is perturbation theory which allows to compute small reactivity changes with high accuracy. Perturbation theory can be applied to the computation of small changes in the neutron flux, too, but we will not go into the details of this. Before presenting perturbation theory, we present a simple example illustrating the power of the method.
Equation 1 (i.e. the Helmholtz equation) is selfadjoint. Consequently, we can show that the formula
_{} 
(3.5) 
is stationary against the variations of F . Assume that some function F(r) satisfies Equation 1. Then, its small variations dF(r) do not change B^{2} in first order:
_{} 

_{} 

_{}_{.} 

Let us take the example of a sphere with radius R = 1. The Laplace operator has the following form in spherical geometry:
_{}. 

Equation 3.5 has the following form in this case:
_{}. 
(3.6) 
The buckling is known for this case: B^{2} = p^{2} = 9,87. In order to illustrate the method, we take some trial functions satisfying the boundary condition F(1) = 0 and certain trivial requirements. One of the simplest functions which are differentiable twice is
_{}. 
(3.7a) 
Formula 3.6 yields the buckling _{} for it. Its error is 40%. This trial function has the drawback that it has no extreme for r = 0 its first derivative being equal to –2. We modify the trial function:
_{}. 
(3.7b) 
Its derivative vanishes for r = 0. Formula 3.6 leads to _{} which is a much better approximation, the error being only 6.5%.
We could imagine further trial functions. Now the question arises: which one is the best. We know from functional analysis that formulae like 3.5 have a minimum property: from among the trial functions invented, that one leads closest to the exact solution for which the buckling is the smallest. This observation allows to refine the trial functions further. Let us consider for example a linear combination of the function 3.7b and of its square:
_{} 

where c is not specified for the time being. Now we proceed as follows: substitute F_{3} in Equation 3.6, then find the minimum of the buckling in terms of c. We get two values minimizing B^{2}:
c_{1} = 1.098 and c_{2} = –1.41. 

The corresponding bucklings are
for _{} and for _{}. 

The error of the former value is only 0.1%. The other buckling belongs to the second harmonic which is an approximation for 4p^{2} = 39.5. Such methods are known as RayleighRitz methods in mathematics.
Perturbation theory is applied when we have to compute the effects of small changes directly. Most frequently, we are looking for the change in reactivity. There are cases when we are interested in the perturbation of the neutron flux but we discuss only reactivity perturbations in the following. We call perturbation every small change in the operators describing the reactor. We give two examples for illustrating what can be considered small:
(1) When the perturbation consists in moving a control rod a little, the cross sections change only in a small volume of the reactor. The changes in the cross sections can be locally large.
(2) When the concentration of boric acid dissolved in the moderator changes a little, the absorption cross section is perturbed. Here the perturbation concerns the whole volume of the reactor but the change is small in every point.
The reactivity is obtained by solving the static eigenvalue equation 2.9a. The first idea is to express the multiplication factor from the equation itself:
_{} 

leading to the reactivity:
_{} 

where y^{+}(x) is an arbitrary weight function (belonging to the detector field). The perturbation of the reactor means that the operators change in Equation 2.9a:
_{} changes to _{}, 

_{} changes to _{}. 

As a result of this, both the reactivity and the static neutron density change:
r changes to r + dr, 

n_{s}(x) changes to n_{s}(x) + dn_{s}(x). 

The reactivity perturbation is then obtained as the difference
_{}. 
(4.1) 
This formula yields the required results but it is not satisfactory for two reasons: firstly, its application requires that the density perturbation dn_{s}(x) is computed, secondly, there is now warranty that the effects of the approximations applied in solving Equation 2.9a were minimal. Both difficulties can be eliminated if the weight function y^{+} is chosen adequately. The variational principles discussed in the previous section suggest that the best choice is the static adjoint eigenfunction obtained as a solution of Equation 2.9b:
_{} or _{}. 
(4.2) 
We show that this is stationary against small variations of n_{s}(x):
_{} 

_{}. 

In the last step, we took into account Equation 4.2 and the fact that the adjoint function satisfies Equation 2.9b. It follows from this that this choice of the weight function eliminates – at least up to first order – the eventual approximations applied in computing the neutron density. Unfortunately, also the adjoint function can be computed only approximately. Thus, we could think that our method has the drawback that we have to know the adjoint function accurately. Things are not as bad as that: we could show by a similar derivation that Equation 4.2 is stationary against small variations of the adjoint function, too. We may thus draw the conclusion that the error of Equation 4.2 is of the order of
_{} 

i.e. the error is only of second order.
It follows from our considerations that the effects of small perturbations of the neutron density in Equation 4.1 are negligible, and it is sufficient to take into account the perturbations of the operators:
_{}. 
(4.3a) 
This leads to the usual perturbation formula for the reactivity:
_{} 
(4.3b) 
where
_{}. 
(4.4) 
Equations 4.3 are the fundamental formulae of perturbation theory. They demonstrate how we could realize our goal formulated in the introduction: we can determine the reactivity perturbation by solving Equations 2.9 for the unperturbed state of the reactor. In terms of these solutions, Equations 4.3 yield the reactivity change related to an arbitrary perturbation. It is an important aspect of the method that the eigenvalue equations 2.9 need not be solved for the perturbed reactor state at all.
The characteristic of a control rod is a wellknown application of perturbation theory. Assume that the rod moves in an empty tube parallel to the zaxis, and its horizontal coordinates are (x_{0}, y_{0}). Let the reactor be critical when the rod is its upper position. Since the rod moves in an empty tube, the cross section are (practically) zero in points inside the tube when the rod is its upper position. We solve the eigenvalue problem 2.9 for this reactor state. We choose the onegroup diffusion model according to which the solution is
_{} 
(4.5) 
where H is the reactor height and A is a coefficient depending on the horizontal rod position.
When the rod is inserted to some depth, we can say (approximately) that only the absorption cross section changed in the volume occupied by the control rod since the scattering cross of the rod’s material is generally small. When the rod is inserted to the height z, the perturbation of the operator is
_{} 
(4.6) 
where S_{a} is the absorption cross section of the control rod. This means that, in Equation 4.3b, we have to integrate the numerator only for the inserted volume of the control rod while the integral in the denominator is extended to the whole volume of the reactor. Thus, the dependence of the reactivity on z is determined by the numerator. The denominator, S_{a}, and the other quantities independent of z can be condensed in a single constant C. With this, we get
_{}. 
(4.7) 
If we plot this as a function of z, we get a curve like that presented in Figure 4.1. It has a typical shape recalling letter “S” for control rods for which the absorption cross section is independent of z. Such Scurves are regularly taken by the reactor operators. The curve has a relatively long section near z = H/2 which is linear in a good approximation called the linear section of the control rod. The slope of this section is
_{}. 

This defines the physical meaning of coefficient C: it is the slope of the rod’s characteristic on its linear section.
Figure 4.1. Characteristic
of a control rod according to Equation 4.7
(C = 0,05 $/cm;
H = 100 cm)
The point kinetic equation is a simple tool for discussing and interpreting time dependent phenomena. Its derivation is based on the assumption that the neutron density is separable:
_{} 
(5.1a) 
where y(x,t) is a function to be specified later. The point kinetic equation is related to the factor j(t). Such a separation is obviously arbitrary. It has some use only if the first factor is somehow restricted. Generally, the following is required:
_{} 
(5.1b) 
where n^{+} is the solution of Equation (2.9b). If we keep the time dependence of y(x,t), the parameters of the point kinetic equation (such as L and b_{eff}) will be also time dependent. In order to avoid this, the best choice is
_{} 
(5.2) 
where n is the solution of Equation (2.9a). With this, condition (5.1b) is automatically satisfied.
We omit the derivations because they are lengthy and not too interesting. We restrict ourselves to the final result concerning b_{eff}. As a matter of fact, the main interest of the derivations consists in obtaining this formula, which is fundamental for reactor kinetics. We shall need some notations.
Introduce subscript m which runs through the fissile and fertile isotopes (m = ^{235}U, ^{238}U, ^{239}Pu etc.). f_{i}(E) is the spectrum of delayed neutron group i, f_{p}(E) is the spectrum of the prompt neutrons. All these spectra may be assumed as independent of m. All others, in turn, do depend on m: b_{im}, nS_{fm}. By means of these notations, we define the production operator describing the prompt neutrons emitted by isotope m:
_{} 
(5.3) 
where b_{m} is the nuclearphysical delayed neutron fraction for isotope m:
_{}. 
(5.4) 
In an analogous way, we define the production operators for the delayed neutrons (i = 1, 2, ..., 6):
_{}. 
(5.5) 
The total production operator is thus given as the sum of all these operators:
_{}. 
(5.6) 
The static fission spectrum appearing in the operators (1.2c) and (1.3b) is a linear combination of these spectra. The latter should be averaged for all isotopes involved. The averaging factors are the integrals
_{} 
(5.7a) 
where
_{}. 
(5.7b) 
Strictly speaking, both F and F_{m} depend on r but this is not shown in the formulae. Now the average static fission spectrum is
_{}. 
(5.7c) 
After some lengthy derivations, we get the following formula for the effective delayed neutron fraction:
_{} 
(5.8) 
where
_{}. 
(5.9) 
In order to apply these formulae, we have to solve the transport equation (or one of its approximations) and the corresponding adjoint equation. This is usually done for the critical state of the reactor. Formula (5.9) involves integration for the spatial variable r as well as summation (or integration) for the energy variable. Thus b_{eff} depends on the spatial and energy dependence of the neutron density and the adjoint function. The main effect is the latter i.e. the energy dependence: b_{eff} needs to be computed with a large number of energy groups. We show an example illustrating this statement. If the core composition of the considered reactor is simple like in a critical assembly, the spatial effects are much less important than the spectral ones. Therefore, the best is to treat them by asymptotic reactor theory. However, formula (5.9) should be applied explicitly for energetic reactors in which the fuel composition can change from fuel pin to fuel pin.
Let us write Equation (5.9) explicitly for the energydependent diffusion model:
_{} 
(5.10a) 
where V is the volume of the reactor. Usually, such calculations are performed in multigroup approximation. We rewrite this formula in a multigroup picture for completeness:
_{}. 
(5.10b) 
G is the total number of energy groups, v_{g} is the characteristic speed of energy group g. The asymptotic approximation means the following:
_{} and _{}. 

If this is substituted in Equation (5.10a), the integrals for the variable r are the same in the numerator and the denominator leading to
_{}. 

This can be further simplified by virtue Equations (5.7a) and (5.7b):
_{}. 
(5.11) 
In these inner products, the integration is done only for variable E.
The spectra y(E) and y^{+}(E) appearing in Equation (5.11) are computed by means of some asymptotic slowing down program. y(E) is the solution of an equation whose operator form is (cf. Equation (2.9a)):
_{} 
(5.12a) 
with
_{}. 
(5.12b) 
The form of the destruction operator depends on the slowing down model actually used and on the way of treating the leakage. In order to have something in hand, we give its form for the simplest case i.e. for the case when the leakage is treated in diffusion theory:
_{} 
(5.12c) 
where B^{2} is the buckling of the fundamental mode. Equations (5.12) constitute an eigenvalue equation for y(E) which has a nontrivial solution only when k_{eff} and B^{2} take some special values. If, for example, we set k_{eff} = 1, B^{2} should be equal to the material buckling of the core composition. For any selected value of k_{eff}, the slowing down program computes the corresponding value of B^{2} and vice versa. The eigenfunction can be normalized in an arbitrary way. The most comfortable normalization is
_{} 
(5.13) 
since it results in the following slowing down equation:
_{}. 
(5.14) 
The advantage of this “trick” (which, by the way, is possible only in asymptotic theory) consists in that the slowing down equation has become an inhomogeneous integral equation which can be solved for any value of the buckling B^{2}. If, for some given value of B^{2}, we compute the solution y(E), Equation (5.13) yields the multiplication factor k_{eff} as a function of B^{2}. This has a further byproduct: the solution allows to compute the fission fractions F_{m} on the basis of Equation (5.7b) which reads as
_{} 
(5.15) 
in the actual model. According to Equation (5.7c), the spectrum f(E) depends on these fractions. Thus, strictly speaking, we ought to perform an iteration for this dependence. This is rarely done since its effect is rather small in most cases.
Equation (5.11) can be applied only if the adjoint spectrum y^{+}(E) is known. This can be eventually a severe problem because adjoint slowing down programs are not common place. Fortunately, this difficulty can be easily overcome in asymptotic theory. Let us write the adjoint slowing down equation (cf. Equation (2.9b)):
_{} 
with
_{}. 

The adjoint function can be normalized similarly to Equation (5.13):
_{} 
(5.15) 
leading to the adjoint equation
_{}. 
(5.17) 
We see from this that the adjoint function does not depend on the spectrum f(E). This allows to compute the quantity
_{} 

rather easily. Equation (5.17) remains the same when we put the following factor in its right hand side since it is equal to 1:
_{}. 

This equation, in turn, is adjoint to the
_{} 
(5.18) 
equation. If its solution is normalized as
_{} 
(5.19a) 
(cf. Equation (5.13)), Equation (5.18) becomes
_{} 
(5.19b) 
(cf. Equation (5.14)).
As a conclusion of our considerations, we formulate the following algorithm of computing the effective delayed neutron fraction:
1. We first compute _{} i.e. that value of B^{2} for which k_{eff} = 1. The spectrum y(E) obtained in this step is used to compute F_{m} and F according to Equations (5.7a) and (5.7b), see also Equation (5.15).
2. Keeping the value of B^{2} constant, Equation (5.19b) is solved for i = 1, 2, ..., 6, and the k_{i} factors are computed by means of Equation (5.19a).
3. The effective delayed neutron fraction for group i and isotope m is computed by using equation
_{} 
(5.20) 
(cf. Equation (5.11)).
4. The effective delayed neutron fraction for group i is
_{} 

while the total b_{eff} is their sum:
_{}. 

It is noted finally that Step 1 can be generalized: we first find that value of B^{2} for which k_{eff} is an arbitrary csosen positive number (not necessarily equal to 1).
We give a numerical example obtained by a 40group B_{1} slowing down program using the GruelingGoertzel approximation. The example will illustrate how important the number of energy groups can be. b_{eff} was computed in four approximations:
Approximation 1: asymptotic slowing down in 38 groups out of which one is the thermal group;
Approximation 2: asymptotic slowing down in 6 groups out of which one is the thermal group;
Approximation 3: diffusion calculation of a bare reactor in 6 groups, using finite differences; the geometric buckling of the bare reactor was equal to the B^{2} obtained in Approximation 1;
Approximation 4: reflected reactor calculated in 6group finite difference diffusion approximation.
The considered reactor was critical in all cases. For Approximations 1 and 2, the method described was applied while Equation (5.10b) was used for Approximations 3 and 4. Only ^{235}U was considered as fissile isotope. The nuclearphysical value is b = 0.0065 for it. The obtained b_{eff} values are summarized in the following table.
Approximation 1 yielded a large effect since the ratio of the effective and nuclearphysical delayed neutron fractions is
_{} 

i.e. 11%. The difference between the values obtained by Approximations 1 and 2 is also large (2,9%). This means that we loose almost one third (» 2.9%/11%) of the total effect (11%) when the number of energy groups is reduced from 38 to 6. The differences between Approximations 2, 3 and 4 are much less ( 0,6% at maximum). It is remarkable that the effect of the reflector is practically negligible (Approximation 3 vs. Approximation 4).
Approximation 
b_{eff} 
1 
0.007217 
2 
0.007016 
3 
0.006973 
4 
0.006984 
[1] It is important to emphasize that the characteristics (r,E,) of the neutron are specified in the same moment t as the measurement is done. In case of the neutron importance to be defined later, these two moments will be different.
[2] Subscript “f” recalls final – just as subscript “i” recalls initial in Equation 1.4.