Probability of Initiation in Neutron Transport

Abstract We discuss the numerical solution of the nonlinear integro-differential equation for the probability of a divergent neutron chain in a stationary system (i.e., the probability of initiation (POI)). We follow the development described in Bell’s classic paper on the stochastic theory of neutron transport. As noted by Bell, the linearized form of this equation resembles the linear adjoint neutron transport equation. A matrix formalism for the discretized steady state (or forward) neutron equation in slab geometry is first developed, and is then used to derive the discrete adjoint equation. A main advantage of this discrete development is that the resulting discrete adjoint equation does not depend upon how the multigroup cross sections for the forward problem are obtained. That is, we derive the discrete adjoint directly from the discrete forward equations rather than discretizing directly the adjoint equation. This also guarantees that the discrete adjoint operator is consistent with the inner product used to define the adjoint operator. We discuss three approaches for the numerical solution of the POI equations, and present numerical results on several test problems. The three solution methods are a simple fixed point iteration, a second approach that is akin to a nonlinear Power iteration, and a third approach which uses a Newton-Krylov nonlinear solver. We also give sufficient conditions to guarantee the existence and uniqueness of nontrivial solutions to our discrete POI equations when the discrete system is supercritical, and that only the trivial solution exists when the discrete system is subcritical. Our approach is modeled after the analysis presented for the continuous POI equations by Mokhtar-Kharroubi and Jarmouni-Idrissi, and by Pazy and Rabinowitz.

probability of initiation (POI)). We follow the development described in Bell's classic paper on the stochastic theory of neutron transport (Bell 1965). As noted by Bell, the linearized form of this equation resembles the linear adjoint neutron transport equation. A matrix formalism for the discretized steady state (or forward) neutron equation in slab geometry is first developed, and is then used to derive the discrete adjoint equation. A main advantage of this discrete development is that the resulting discrete adjoint equation does not depend upon how the multigroup cross sections for the forward problem are obtained. That is, we derive the discrete adjoint directly from the discrete forward equations rather than discretizing directly the adjoint equation. This also guarantees that the discrete adjoint operator is consistent with the inner product used to define the adjoint operator (see Equation (30)).
We discuss three approaches for the numerical solution of the POI equations, and present numerical results on several test problems. The three solution methods are a simple fixed point iteration, a second approach that is akin to a nonlinear Power iteration, and a third approach which uses a Newton-Krylov nonlinear solver. The nonlinear Power iteration has also been presented by Baker (2005), but was first developed by Bell and coworkers. The nonlinear Power and the Newton-Krylov methods demonstrate the most robustness. Our main purpose in writing this article is the consistent development of the discrete adjoint (and hence POI) equations. Hence, we are most interested in robustness and exactness of the numerical methods, and have not been concerned with developing the most efficient implementations of the methods. We also give sufficient conditions to guarantee the existence and uniqueness of nontrivial solutions to our discrete POI equations when the discrete system is supercritical, and that only the trivial solution exists when the discrete system is subcritical. Our approach is modeled after the analysis presented for the continuous POI equations by Mokhtar-Kharroubi and Jarmouni-Idrissi (1998), and by Pazy and Rabinowitz (1969).
We review some of the notation used throughout the article. For matrices A 2 R mÂn and B 2 R kÂl , the Kronecker (or tensor) product of A and B is the mk Â nl matrix denoted by where A ¼ ða ij Þ: Kronecker products have many interesting properties. We list here the ones relevant to our discussion: If A and B are nonsingular, then A B is nonsingular with ðA BÞ À1 ¼ A À1 B À1 , ðA BÞ T ¼ A T B T , Given matrices A, B, C, and D, ðA BÞ Á ðC DÞ ¼ AC BD, as long as both sides of the equation make sense, ðA þ BÞ C ¼ A C þ B C, and A ðB þ CÞ ¼ A B þ A C: For any matrix A ¼ ða ij Þ and real numbers b cðb < cÞ, by b A cðb < A < cÞ we mean b a ij cðb < a ij < cÞ for all i and j. The rest of the article is organized as follows. In §2, we consider the time independent neutron transport equation in slab geometry and review the standard multigroup semi-discretization. In §3, we review the Petrov-Galerkin spatial discretization discussed in Brown (1995) and Step discretization (Lewis and Miller 1993), and then develop a matrix form of the fully discretized multigroup, discrete ordinates equations. In §4, we derive the discrete adjoint problem to the discretized multigroup equations in §3, and then in §5 we introduce the time independent probability of initiation (POI) equations, derive a discretized POI system along with a matrix form, and discuss its numerical solution. We also discuss existence and uniqueness of nontrivial solutions to the discrete POI in this section. We present some numerical results in §6, and then in §7 make some concluding remarks.

The neutron transport equation
We begin with the linear time-independent neutron transport equation in one-dimensional slab geometry with P N scattering (Lewis and Miller 1993). The spatial domain is the interval ½a, b in x, l is an angle cosine in ½À1, 1, the energy variable is E 2 ð0, 1Þ, and the equation in the flux wðx, l, EÞ is given by l @w @x ðx, l, EÞ þ rðx, EÞwðx, l, EÞ ¼ ð 1 0 X N s n¼0 2n þ 1 2 r s, n ðx, E 0 ! EÞP n ðlÞ ð 1 À1 wðx, l 0 , E 0 ÞP n ðl 0 Þdl 0 dE 0 þ qðx, l, EÞ: (1) Here, the P n are the (unnormalized) Legendre polynomials (P 0 ¼ 1, P 1 ¼ l, etc.), N s is the scattering order, and q is a source term. The r s, n are the Legendre moments of the scattering kernel r s ðx, l, E 0 ! EÞ, given by r s, n ðx, E 0 ! EÞ 1 2 ð 1 À1 r s ðx, l 0 , E 0 ! EÞP n ðl 0 Þdl 0 , where l 0 is the cosine of the scattering angle. The total cross section is r.

Discrete formulation of the 1-D slab problem
In this subsection we describe in some detail the Petrov Galerkin discretization for (1)-(2), as well as the Step discretization (cf. Lewis and Miller 1993), and also present matrix formulations of the resulting discretized system. Our 1D Petrov Galerkin method is equivalent to the well-known diamond difference discretization scheme, which was discussed in Ashby et al. (1995). We begin by introducing a spatial grid and let Dx i ¼ x i À x iÀ1 : The x i are referred to as nodes, and function values at these points are called nodal values. Assume that r g and the r s, g, g 0 , n have constant values on each zone x iÀ1 < x < x i , denoted by r g, i and r s, g, g 0 , n, i respectively. Function values that are constant on zones will be referred to as zone-centered values. We use a discrete ordinates collocation of (1) at (an even number of) Gauss points l ' with and l Lþ1À' ¼ Àl ' : The integral in (5) is then approximated by 1 2 for each x 2 ½a, b, where the w ' are the corresponding Gaussian weights. We use w g, ', i to denote the approximation to w g ðx i , l ' Þ, the true solution at ðx i , l ' Þ: Note that P L '¼1 w ' ¼ 1: The particular symmetry assumption we will need is symmetry through the origin. Namely, if l ' is a quadrature point with corresponding weight, w ' , then Àl ' is also a quadrature point, and letting ' À denote its index (i.e., l ' À ¼ Àl ' , it is also true that w ' À ¼ w ' : We also assume that the quadrature rule is of high enough order to integrate exactly the Legendre polynomials of order N s and below. (See the Discrete Ordinates section below.)

The Petrov Galerkin method for slab geometry
Using piecewise linear basis functions, and testing against constant functions, integrating (ignoring the group index g for the moment) l @w @x ðx, lÞ þ rðxÞwðx, lÞ ¼ f ðx, lÞ over the ith zone yields where Using (8) in the full system (1), and collocating at l ¼ l ' , we obtain The discretized boundary conditions are given by It then follows that (10)-(12) give a system of GLðM þ 1Þ equations in GLðM þ 1Þ unknowns.
2.3. The step method for slab geometry Starting from the system (10), we modify the removal and source terms to create the Step method (cf. Lewis and Miller 1993). For l ' < 0, we use and for l ' > 0, for The discretized boundary conditions are the same as for the Petrov Galerkin method, namely (11) and (12). It then follows that (13)- (14) and (11)-(12) give a system of GLðM þ 1Þ equations in GLðM þ 1Þ unknowns.

A matrix formulation for slab geometry
We next wish to write system (10)-(12) in matrix notation. To accomplish this, we first consider the matrix form of the discretized operator l ' @=@x þ r for a single group g. Define the vector Next, for the Petrov Galerkin method, define the block diagonal matrix The matrices S þ and S À define the zone-centered fluxes for the Step method as either the upwind or downwind nodal flux. Next, define matrices To isolate the boundary values, define row vectors of length M þ 1 by where indices on the standard basis vectors e ' run from 0 to M. Define the matrices Z and Z b by is the discrete diamond difference representation of l ' @=@x þ r g , with both x and l discretized. Note that H g, ' operates on nodal vectors.
One can easily show that H g is nonsingular for either the Petrov Galerkin or Step methods. Indeed, from its block-triangular form, H g will be nonsingular if for each ' ¼ 1, Á Á Á , L, the matrix ! is nonsingular. For ' L=2, we have l ' < 0 and B ' ¼ e T M : For this case, the matrix H g, ' is an upper triangular matrix with positive diagonal entries, and is therefore nonsingular. For ' > L=2, we have l ' > 0 and B ' ¼ e T 0 : Moving the last row of H g, ' up to the first gives a lower triangular matrix with positive diagonal entries. Thus, H g is nonsingular.

The discrete ordinates method
Continuing the matrix development of the overall discretization of (1), we begin by defining discretized representations of the operations of taking moments of the flux. As operators on zone-centered vectors, these are easily seen to be given by the M Â LM matrices L n ðl n WÞ I M where l n ðP n ðl 1 Þ, P n ðl 2 Þ, Á Á Á , P n ðl L ÞÞ and W diagðw 1 , Á Á Á , w L Þ: If the vector W g approximates w g ðx, lÞ, then L n W g will approximate the n th moment of w g ðx, lÞ, namely / g, n ðxÞ: Similarly, we define the LM Â M matrices L þ n ð2n þ 1Þl T n I M : For a vector U approximating /ðxÞ, L þ n U will approximate P n ðlÞ/ðxÞ: We also will find it useful to define the grouped matrices Given an N ¼ N s , the number of terms in the scattering kernel, we have assumed that the quadrature rule is symmetric through the origin (see remarks above) and such that the Legendre polynomials of order N s and less satisfy for all 0 n, n 0 N s : This can be written more compactly as To represent the source term, define the vectors Next, define matrices of scattering coefficients, R s, g, g 0 , n diagðr s, g, g 0 , n, 1 , Á Á Á , r s, g, g 0 , n, M Þ 2 R MÂM , and the matrices The matrix Z injects zone-centered vectors into the nodal vector space, and the matrix S averages nodal vectors to obtain zone-centered ones. Note that Z T Z ¼ I and Z T Z B ¼ 0: Using the above matrices, the matrix H g can be written as The matrices Z and Z B are needed since H g operates on nodal vectors, while the scattering matrix operates on zone-centered vectors. (Recall that f was assumed to be zone-centered in the development of the Petrov-Galerkin method discussed earlier.) The complete discretization of (17)-(18) can be written in the compact form with F g ZQ g : where R N s s, gg 0 diagðR s, g, g 0 , 0 , Á Á Á , R s, g, g 0 , N s Þ, and define S I G S, Z I G Z, H diagðH 1 , H 2 , Á Á Á , H G Þ, L þ I G L N s , þ , L I G L N s , then (27) can be written as with F ¼ ZQ:

Discrete form for adjoint problem
In this section we derive the discrete adjoint problem to the system (12). We define the matrix T by and so Equation (12) can be written as TW ¼ F: Given a discrete inner product that approximates the integral over all of space x, direction l, and energy E, the discrete adjoint system can be derived in a way analogous to that for the continuous system (cf. Lewis and Miller 1993). For vectors X, X Ã 2 R MLG , we define the inner product with DE diagðDE 1 , Á Á Á , DE G Þ, and DE g ¼ E gÀ1 À E g : Note that the ðDEÞ À1 factor is needed since the multigroup formulation solves for and hence a factor DE À1 g is needed to approximate the integral over energy of the product of two such terms. (We note that one could alternatively solve forw g, ', i w g, ', i =DE g , in which case the above inner product would then have a factor of DE instead of DE À1 : However, the cross section matrices would then need to be scaled by DE À1 on the left and DE on the right. We use the familiar form of the multigroup equations here.) Additionally, the inner product is defined for zone-centered vectors, whereas the unknown W in (29) is a nodal vector. To remedy this situation we rewrite (29) using a zone-centered unknown. This is possible since for adjoint problems we only consider vacuum boundary conditions on outgoing faces. Assuming this has been done, and we denote W z as the zonecentered unknown with T z the zone-centered version of T, then the adjoint T Ã z is defined as the matrix satisfying for all W z and W Ã z : We now begin with the discussion of how to rewrite (29) using a zonecentered unknown. We start with a lemma.
Proof. For l ' < 0Z l ' < 0, the matrix ZS þ Z b B ' is upper triangular with positive diagonal entries, and hence is nonsingular. For l ' > 0, reordering the last row first the matrix is lower triangular, again with positive diagonal entries, and is nonsingular. Hence, the matrix is nonsingular for all ' ¼ 1, Á Á Á , L: Thus, Since Z T Z ¼ I M and Z T Z b ¼ 0, multiplying (32) on the left by Z T gives and multiplying on the right by Z gives the second assertion in (31). Similarly, multiplying (32) on the left by Z T b gives and then multiplying on the right by Z b gives the first assertion in (31). w Next, given a zone-centered vector W z , the nodal vector W defined by For the above defined W, we have from Lemma 3.1 that W ¼ ðZS þ Z B BÞ À1 ZW z : Thus, substituting into (29) and multiplying on the left by Z T gives (recall that Z T Z ¼ I) which is (29) reformulated in terms of a zone-centered unknown. We have the following relationship between the inverses of H and H z : Lemma 3.2. The matrices H and H z are nonsingular, and we have Proof. From the definition of C z , we have Since H g, ' is nonsingular for all g and ', it follows that the matrix H is nonsingular. Therefore, multiplying this last equation on the left by SH À1 gives SH À1 ZH z ¼ I, which proves that H z is nonsingular and gives (34). w Returning now to the adjoint Equation (30), we have In defining the adjoint matrix T Ã z , we start by considering the case in which all of the cross sections are zero. For this case, From the definition of the matrices in this last equation, note that the matrix C z ¼ CðZS þ Z B BÞ À1 Z is block-diagonal for each g and ' with size M Â M. In fact, there is no g dependence on the blocks, and for each ' we have a block of the form If we choose the vectors W z and W Ã z so that they have nonzero components for only one g and ', then Also for this case, we will have We next need some additional lemmas.
Lemma 3.3. Assume the 1D Petrov Galerkin discretization is being used. Then for all ' ¼ 1, Á Á Á , L, the following relationship holds Proof. As B ' equals e T M or e T 0 , we consider only the case where B ' ¼ e T M , with the other case just a permutation of this one. For this case, note that B ' À ¼ e T 0 : We also take M ¼ 4 to illustrate the proof for the general case. Then, it is easily seen that Hence, it also follows easily that Comparing to a similar calculation of the right side of (36) with B ' À ¼ e T 0 gives the result.
w Lemma 3.4. Assume the 1D Step discretization is being used. Then for all ' ¼ 1, Á Á Á , L, the following relationship holds Proof. For the Step method, the proof is even simpler. An easy calculation gives that with the other case just a permutation of this one. For this case, note that B ' À ¼ e T 0 : We also take M ¼ 4 to illustrate the general case. Then, it is easily seen that Comparing to a similar calculation of the right side of (36) with B ' À ¼ e T 0 gives the result. w Using Lemmas 3.3 and 3.4, we can now define T Ã z for this particular case as and we have Returning now to the other terms in (35), we have To simplify the discussion, we assume isotropic scattering. The general case follows similarly. In this case, Thus, In general, defining DE N s DE I N s þ1 I M and we can write the full zone-centered adjoint matrix as with C Ã z ÀCðZS þ Z B B À Þ À1 Z: The final step in defining the adjoint to the original nodal matrix T involves using the relationship Thus, we can define the nodal adjoint matrix via where H Ã ZðÀC þ RSÞ þ Z B B À : Using this definition we have the two identities We also note that the analogue of Lemma 3.2 holds for the adjoint matrices, namely ðH Ã z Þ À1 ¼ SðH Ã Þ À1 Z: Given a code that solves the forward problem (29), to modify it to solve the adjoint system (assuming that the quadrature rule is symmetric through the origin, i.e., if l is a point in the quadrature rule, then so is Àl), no changes are required to invert the ZðÀC þ RSÞ þ Z B B À matrix since it is just a permutation of the directions in the forward matrix ZðC þ RSÞ þ Z B B, and the scattering matrix ZL þ R Ã s LS is invariant under this same permutation that replaces an l ' by its corresponding l ' À : The only change has to do with using the transpose of the scattering matrix, R s , and scaling on the left and right, respectively, by DE and DE À1 : It does not matter how the scattering matrix R s was generated.

Probability of initiation problem
In this section we describe the form of the probability of initiation problem that is to be solved, and then discuss its numerical solution. First, let pðx, l, EÞ be the probability that a neutron at position r, traveling in direction l, with energy E, will initiate a divergent chain. From Bell (1965) where the superscript i is an index over isotopes, and where r i , r i s , and r i f are the total, scattering, and fission, microscopic cross sections for isotope i: As the fission emission spectrum, v i ðE 0 Þ, is different for each isotope, the nonlinear terms are a weighted sum over i of nonlinear terms, each similar in form. Note that where i ðEÞ is the mean number of neutrons produced in a fission of isotope i caused by a neutron with energy E. We also note this formulation assumes that fission is the only collision process which produces more than one neutron per collision. From Terrell (1957),Ĉ i n ðEÞ depends upon E through i ðEÞ through the formula where r p ¼ 1:08 and b is a small parameter (b < 10 À2 ) which Terrell recommends setting to b ¼ 0:5 À 0:5f ðð i ðEÞ þ 0:5Þ=r p Þ, where f ðxÞ ð2pÞ À1 Ð x À1 exp ðÀt 2 =2Þdt: Equation (43) is nonlinear, but if M f ¼ 1, then it reduces to the linear adjoint equation This equation (with vacuum boundary conditions on outgoing faces) has a nontrivial solution if it is exactly critical. Otherwise, a 1=k Ã factor is introduced in front of the fission term, and one solves a criticality eigenvalue problem instead, with the system being subcritical for k Ã < 1 and supercritical for k Ã > 1: Note that if k is the criticality eigenvalue for the forward criticality problem, then k ¼ k Ã : See Lewis and Miller (1993) for more details about criticality and adjoint problems.

Discrete POI problem
Our solution procedure proceeds as for the forward problem, using a multigroup discretization in energy, discrete ordinates in direction, and a Petrov-Galerkin or Step discretization in space. As Bell (1965) has noted, when the physical system is only slightly supercritical, the nonlinear terms in (43) are very small. Hence, the solution pðx, l, EÞ should be wellapproximated by the solution w Ã ðx, l, EÞ to the adjoint problem (45), with an appropriate scaling found by integrating (43) over all of phase space with p ¼ w Ã : We describe the energy discretization in some detail, as the discrete form of the nonlinear term must be performed in such a way that the discrete POI problem is consistent with the discrete adjoint equation derived above. Discretizing in energy leads to the multigroup equations and this is exactly the same relationship found in the discussion of the discrete adjoint problem in (40). Continuing with the discrete form, we now have for g ¼ 1, Á Á Á , G, where p g 0 , n ðxÞ 1 2 ð 1 À1 p g 0 ðx, lÞP n ðlÞdl is the n th moment of p g 0 : Next, from the forward problem, we would have Note that P G g¼1 v g ¼ 1 in the multigroup formulation since Ð 1 0 vðEÞdE ¼ 1: Next, the coefficients A i m ðEÞ depend upon the energy E only through ðEÞ, which follows from the definition of theĈ i n ðEÞ in (44). DefiningĈ i n, g We note that from (44), As p g ðx, lÞ is the unknown being solved for, we can approximate ð g 1 2 as P 0 ðlÞ 1: Thus, the multigroup Equation (48) can be written as for g ¼ 1, Á Á Á , G: If M f ¼ 1, then this last equation reduces to the multigroup discretization of the adjoint problem (45), for g ¼ 1, Á Á Á , G: The spatial and directional discretizations proceed as before, and thus we can write the fully discretized systems in matrix form. For the discrete adjoint criticality problem, we have where , and R f is the discretized fission matrix for the forward criticality problem. In order to describe the nonlinear term in (50), we must describe the matrix R f in more detail. Defining we can write R f as Thus, and for any vector U, define diagðUÞ to be a square diagonal matrix with the number of rows and columns equal to the number of entries in U: Next, for any vector U 2 R MG , definê for each isotope i, and let P P 1 . . .
where p g, ', i % p g ðx i , l ' Þ: With the notation defined above, the discrete POI problem can be written in the form Á T for any vector U: We note that since the discrete adjoint matrices R Ã s and R Ã f do not depend upon how the forward multigroup cross sections are defined, the discrete POI system also does not depend upon how the multigroup cross sections for the forward problem are defined. We have also implicitly assumed that the functions n i ðxÞ are piecewise constant on the spatial zones, and when evaluating the spatial integral of the nonlinear terms over each zone Z i to derive the finite element equations, we use the cell-centered (i.e., nodal average) value of p g, 0 ðxÞ, (i.e., 1 2 ðP g, 0, i þ P g, 0, iþ1 Þ, or L 0 SP).

Solution of the discrete POI problem
In this subsection, we discuss three different approaches to solving (54). The first is a simple fixed point iteration, the second can best be described as a nonlinear power iteration, and the third approach uses a Newton-Krylov nonlinear solver.

Algorithm: POI Fixed Point Iteration
(1) Choose an initial guess, P 0 : (2) For each l ¼ 0, 1, Á Á Á , until convergence, solve Given a vector F Ã , we solve the system T Ã W Ã ¼ F Ã by the Block-Jacobi iteration Algorithm: Block-Jacobi Iteration (1) Choose an initial guess, W Ã, 0 : (2) For each k ¼ 0, 1, Á Á Á , until convergence, solve Note that inverting H Ã involves only sweeps. In the context of the POI fixed point iteration, Ã f ðL 0 SP l ÞL 0 SP l , and we use the initial guess W Ã, 0 ¼ P l : This becomes a better and better initial guess as the POI fixed point iteration gets close to convergence. Since the iteration (55) is nonlinear, some care must be taken with the initial guess P 0 : Otherwise, the iteration will diverge. We use the solution of the adjoint criticality problem, W Ã , scaled so that it satisfies nonlinear system (54) in an integral sense. Let the vector v ¼ ðDE I L I M Þ1, where 1 2 R MLG is a vector with all 1's. Then for a zone-centered vector W z approximating the function wðx, l, EÞ, we have hv, W z i % Given the adjoint criticality vector, W Ã , and eigenvalue, k Ã , that solves (52), we set P 0 ¼ k 0 W Ã , where k 0 is chosen so that P 0 satisfies (54) in an integral sense via From the polynomial form ofF i f ðUÞ, Equation (57) is a polynomial in k 0 of degree ðM f À 1Þ: For this initial guess, it is sufficient to take M f ¼ 4, and then set k 0 equal to the real root closest to unity. As remarked above, for systems that are just slightly supercritical, this P 0 should be close to the solution P: However, the convergence rate of this fixed point iteration is extremely slow for systems just slightly supercritical.

Nonlinear power iteration
The second solution approach is akin to a nonlinear power method iteration. Essentially this same approach is discussed in Baker (2005). We present it here for completeness.

Algorithm: POI Nonlinear Power Iteration
(1) Solve the adjoint criticality problem (52) for W Ã and k Ã : (2) Find k 0 , the root closest to unity of the polynomial (57), using M f ¼ 4: As the iteration convergences, we will have k l ! 1: The roots of the above polynomials are obtained by computing the eigenvalues of the related companion matrix, and then choosing the root closest to unity. This iteration works extremely well, even for systems only slightly supercritical.

Newton-Krylov nonlinear solver
This last solution approach uses a Newton-Krylov nonlinear iteration. The nonlinear problem to solve has the form We use the KINSOL Newton-Krylov solver from the SUNDIALS suite (Hindmarsh et al. 2005). For a preconditioner, we must approximate the Jacobian matrix F 0 ðPÞ: From the polynomial form of the nonlinear term, we have that The preconditioner linear system F 0 ðPÞW ¼ w is approximately solved using 5 iterations of the Block-Jacobi iteration: for l ¼ 0, 1, Á Á Á , 5, For an initial guess, we take P 0 ¼ W Ã , where W Ã is the adjoint criticality eigenvector. Since physically, pðx, l, EÞ 0 for a critical or subcritical system, we only need to solve the nonlinear problem for supercritical systems.

Convergence, existence, and uniqueness results
In this subsection, we give sufficient conditions that guarantee convergence of the above Block-Jacobi and Fixed Point iterations. Furthermore, we demonstrate that the POI problem (54) has a unique nontrivial solution when the system is supercritical, i.e., when k Ã > 1, and has only the trivial solution when k Ã 1: In our existence and uniqueness proofs for the POI problem (54), we employ techniques that closely follow those used in Mokhtar-Kharroubi and Jarmouni-Idrissi (1998) and Pazy and Rabinowitz (1969) for the study of the continuous probability of initiation equations. Sufficient conditions are given therein for the existence and uniqueness of nontrivial solutions to the continuous POI equations. We show similar results for the discrete POI problem (54). For technical reasons, we will work with an equivalent zone-centered version of (54), which we derive first. Consider Equation (54) written in terms of the zonal unknown P z ¼ SP, Equation (58) can be rewritten as Now, ðDE À1 I LM ÞT Ã z ðDE I LM Þ ¼ H Ã z À L þ R T s L, and hence, (60) can be written as Multiplying on the left by ðH Ã z Þ À1 gives We give conditions below that will imply the spectral radius qððH Ã z Þ À1 L þ R T s LÞ < 1: Thus, the matrix we have (defining R 0 diagðR 1 , Á Á Á , R G Þ and noting that R À1 L þ which is a discrete integral equation for the scalar version ofP z , namelỹ R z : IfR z is a solution of (63), then the solutionP z to (61) can be found by solving the system It is for Equation (63) that we will prove our results. Note thatR z has pointwise values in energy, and hence as a probability should have values between 0 and 1.
As the criticality eigenvalue problem will be important later, we write it in the form Specifically, we will need that the eigenvectorŨ Ã z, 0 > 0: We give sufficient conditions on the cross sections, quadrature rule, and spatial grid to guar-  Greenbaum (1997)) will imply that k Ã > 0, is a simple eigenvalue, and has corresponding eigenvectorŨ Ã z, 0 > 0: From results in Mokhtar-Kharroubi and Jarmouni-Idrissi (1998) and Pazy and Rabinowitz (1969), it is clear that when the discrete problem is supercritical (k Ã > 1), then Equation (63) should have a nontrivial solution between 0 and 1, and when critical or subcritical (k Ã 1) should have only the identically 0 solution. We will need some additional assumptions. We assume that a multigroup approximation in energy has been made leading to a system of the form (50), and that a quadrature rule has been chosen that is symmetric through the origin (as discussed above). We also assume the total cross sections r g, i > 0 for all groups g and zones Z i : We make the following additional assumptions: z, g Þ À1 Q g 1 for all g and any Q g 2 R M with 0 Q g 1, and for each g the corresponding diagonal block of L 0 ðH Ã z Þ À1 L þ 0 is a positive matrix. H2 The vectors A i 1 , F i f , X i > 0, and N i ! 0, with P i n i ðx i Þ > 0 for every zone Z i (i.e., each zone has fissionable material), and define the diagonal matrix recalling that R¼diagð R 1 ,ÁÁÁ, R G Þ, with R g ¼I L R g , and R 0 diagðR 1 ,ÁÁÁ,R G Þ: From the definitions of the various vectors in H2, the diagonal matrix ðR 0 Þ À1R T f has main diagonal entries P i n i i r i f, g =r g, i for all zones. To satisfy L þ R T s L ! 0 in H2 when there is anisotropic scattering, it may be necessary to solve a least squares problem as discussed in Dahl, Ganapol, and Morel (1999) and Moskalev (1993). We also note the condition, 0 ðH Ã z, g Þ À1 Q g 1 for all g and any Q g 2 R M with 0 Q g 1, in Assumption H1 implies that kðH Ã z Þ À1 k 1 1: To see this, let A ¼ ða ij Þ be a matrix satisfying 0 AQ 1 for any 0 Q 1: If a kl > 1, let Q ¼ e l , the l-th unit vector. Then AQ equals the l-th column of A whose k-th entry is greater than 1, a contradiction. Thus, a ij 1 for all i and j. Next, if a kl < 0, again let Q ¼ e l : Then AQ again equals the l-th column of A whose k-th entry is negative, again a contradiction. Thus, a ij ! 0 for all i and j. Finally, taking Q to be a vector of all 1's, AQ has entries the row sums of A, all of which are 1: As the maximum of these equals kAk 1 , we have kAk 1 1: We verify assumption H1 for the Petrov Galerkin (under an additional assumption) and the Step discretizations.
Lemma 4.1. Suppose we are given functions r g ðxÞ such that a r g ðxÞ b on a x b for all g ¼ 1, Á Á Á , G for some constants 0 < a < b < 1, and an L > 0, the order of a symmetric Gaussian quadrature on ½À1, 1. Then it is possible to choose the spatial grid so that for each i ¼ 1, Á Á Á , M, with Furthermore, we then have that Assumption H1 holds for the 1D Petrov Galerkin discretization.
Proof. Since the functions r g ðxÞ are bounded, and since for our Gaussian quadrature none of the quadrature points are zero, it follows that the spatial grid can be chosen to satisfy (67). We next show that for each g and ', ZQ 1 for any 0 Q 1: We suppress for the moment the g and ' dependence, and assume that l ' > 0, for ease of exposition. A close look at the slab Petrov Galerkin equations gives for each zone with w 0 ¼ 0 and r i ¼ r g ðx iÀ1=2 Þ: Thus, Therefore, if 0 q i , w iÀ1 1, then 0 w iÀ1=2 1: By (67), l=ðr i Dx i Þ > 1=2, and so we also have that 0 w i 1: By assumption 0 q i 1 for all i. For i ¼ 1, w 0 ¼ 0, and so 0 w 1=2 , w 1 1: Continuing gives 0 w iÀ1=2 , w i 1 for all i ¼ 1, Á Á Á , M: Furthermore, note that if q i 0 > 0 for some i 0 , then w iÀ1=2 > 0 for all i ¼ i 0 , Á Á Á , M: A similar result holds for l < 0, and if q i 0 > 0 for some i 0 , then w iÀ1=2 > 0 for all i ¼ 1, Á Á Á , i 0 : Hence, if 0 Q 1 and ZQ, then 0 W z, g, ' 1: Next, for each g, the corresponding diagonal block of Thus, if 0 Q 1, then Furthermore, if q i > 0 for some i, then from the remarks above it is clear thatQ > 0, and hence the main diagonal blocks of L 0 ðH Ã z Þ À1 L þ 0 are positive matrices. Therefore, Assumption H1 holds. w We omit the proof for the 1D Step method as it is similar to the above proof. However, we note that in the case of the Step method, Assumption H1 holds given any discretization of the POI problem. That is, the assumption holds for any spatial mesh.
Lemma 4.2. Suppose we are given a discretization of (63) such that Assumption H2 holds. Then (62). Furthermore, if Assumption H1 also holds, then K 0 ! 0, and the M Â M main diagonal blocks of K 0 are positive.
Proof. From H1 and H2, as noted in the proof of Lemma 4.2 we have that K 0 ðR 0 Þ À1 is a nonnegative matrix with strictly positive main diagonal blocks. Furthermore, by Assumption H2, P i n i ðx i Þ > 0 for every zone Z i , and hence the block matrix (with block size M Â M) is such that each column of each block A g, g 0 has at least one positive entry. To see this, note that by reordering the unknowns so that the group index is fastest, there exists a permutation matrix O such that where each V i 2 R GÂG and has the form by H2, and noting that OðX i I M Þ ¼ I M X i : Therefore, letting V i ¼ ðv i, g, g 0 Þ, we have , and each v i, g, g 0 > 0: LetK ¼ K 0 ðR 0 Þ À1 ¼ ðK g, g 0 Þ: ThenK g, g > 0 and K g, g 0 ! 0 by Lemma 4.2. Therefore,K g, g A g, g 0 > 0 for all g, g 0 ¼ 1, Á Á Á , G, and hence the matrix g, g 00 A g 00 , g 0 @ 1 A > 0: w Lemma 4.4. Assume H1 and H2 hold for the discrete POI system (63) and that (49)  Proof. For each g and isotope, define From the definitions of the A i m, g , it follows that f g, i ðxÞ can also be written as 1 for all g and isotope by (49). Therefore, 0 f g,i ðxÞ 1 for all 0 x 1: Also, f 0 g,i ðxÞ ¼ P M f n¼1 nĈ i n,g ð1 À xÞ nÀ1 , and so f 0 g,i ð0Þ > 0 and f 0 g,i ðxÞ > 0 for 0 x < 1: Therefore, f g,i is monotone on 0 x 1: From the definition of X i , we have 0 ðX i I M Þ T U 1 for all 0 U 1: Thus, it follows that for all 0 U 1, where 1 2 R GM is a vector of all 1's. Hence, using (59) for all 0 U 1: Next, we need to show that K 0 ðR 0 Þ À1R T f 1 1: This will be true if kK 0 ðR 0 Þ À1R T f k 1 < 1, which holds by condition (66) in H2: Since K 0 ! 0 and ðR 0 Þ À1 ! 0, it follows that 0 N ðUÞ ¼ K 0 ðR 0 Þ À1R T f ðUÞU K 0 ðR 0 Þ À1R T f 1 1, for all 0 U 1, which proves the first part of the lemma.
For the monotonicity of N , note that for each isotope the monotonicity of f g, i ðxÞ for all g implies 0 for all 0 U U 0 1: Thus, it follows that 0 R T f ðUÞU R T f ðU 0 ÞU 0 1: The monotonicity of N now follows from the nonnegativity of the matrices ðR 0 Þ À1 and K 0 : w Theorem 4.5. Assume H1 and H2 hold for the discrete POI system (63). Then (63) posseses a maximal solution.
From the monotonicity of N , one obtains inductively Thus, the vectors fU l g form a bounded monotonically nonincreasing sequence. Hence, it has a subsequence that converges to U, and the monotonicity of the sequence fU l g implies that the entire sequence converges to U: By (69), U satisfies N ð UÞ ¼ U: To show that U is maximal, let U be any solution of (63). Then 0 U 1 ¼ U 0 : The monotonicity of N implies that U ¼ N ðUÞ N ðU 0 Þ ¼ U 1 , and so inductively U U l for all l ! 0: This implies U U, and hence U is maximal.
w Without further assumptions on the cross section matrices, (63) may have only the trivial solution U ¼ 0: We give sufficient conditions for the existence and uniqueness of a nontrivial solution. Consider the adjoint criticality problem (64), and let k Ã andŨ Ã z, 0 be the maximal eigenvalue-eigenvector pair. Suppose that k Ã > 1: Then with assumptions H1 and H2, Equation (63) will have a nontrivial solution.
Theorem 4.6. Assume H1 and H2 hold for the discrete POI system (63). Assume also that the system is supercritical, i.e., k Ã > 1. Then (63) possesses a nontrivial solution. In particular, U is nontrivial.
Proof. By Lemma 4.3, K 0 ðR 0 Þ À1 R T f is a positive matrix. Therefore, Perron's Theorem (Theorem 10.2.3 in Greenbaum (1997) f Þ is a simple eigenvalue with corresponding eigenvector U k > 0, and we assume that it has been scaled so that U k 1: Next, let V i and the permutation matrix O be as defined in the proof of Lemma 4.3. Since by assumption k Ã > 1, for all 0 U 0 , using an argument similar to that used in Lemma 4.3. For any 0 < < 0 , we have Thus, let U 0 ¼ U k and U lþ1 ¼ N ðU l Þ: Then we have U 1 ¼ N ðU 0 Þ ! U 0 , and the monotonicity of N implies U l N ðU l Þ ¼ U lþ1 1 for all l ! 0: As above, the sequence fU l g converges monotonically to U, with U ¼ N ðUÞ and 0 < U k U U 1: 1 for every isotope, and this ensures that the nonlinear terms are strictly concave. As in Pazy and Rabinowitz (1969), a function f(x) for z 2 ½0, 1, is said to be strictly concave if for all 0 < t < 1 and x 1 , x 2 2 ½0, 1 f ðtz 1 þ ð1 À tÞx 2 Þ > tf ðx 1 Þ þ ð1 À tÞf ðx 2 Þ: If f is differentiable, a sufficient condition for f to be strictly concave is if f 00 ðxÞ < 0: Now for each g and isotope i, as done earlier, let f g,i ðxÞ ¼ n,g ð1 À ð1 À xÞ n Þ: Then f g,i 0 ðxÞ ¼ P M f n¼1 nĈ i n,g ð1 À xÞ nÀ1 and f g,i 00 ðxÞ ¼ À P M f n¼2 nðn À 1ÞĈ i n,g ð1 À xÞ nÀ2 : IfĈ i 0,g þĈ i 1,g < 1, then at least one ofĈ i n,g for n>1 must be nonzero, and hence f g,i 00 ðxÞ < 0, which implies that f g,i is strictly concave. Additionally, we have f g,i 0 ðx 1 Þ ! f g,i 0 ðx 2 Þ if x 1 x 2 , with equality only if x 1 ¼ x 2 : Therefore, f g,i ðxÞ < f g,i 0 ð0Þx for 0 < x 1, and note that f g,i 0 ð0Þ ¼ A i 1 : Theorem 4.7. Assume H1 and H2 hold for the discrete POI system (63). Assume also that the system is supercritical, i.e., k Ã > 1, and for all g ¼ 1, Á Á Á , G and every isotope (49) holds withĈ i 0, g þĈ i 1, g < 1. Then (63) possesses a unique nontrivial solution.
Proof. By Theorem 4.5, (63) has a maximal solution U > 0: Now, suppose there exists another nontrivial solution, denoted by U: We show that U > 0: Since U is nontrivial, there exists a g and zone Z i with U g, i > 0: Define Ç ¼ OU, where O is the permutation matrix defined in the proof of Theorem 4.5. Since X i > 0 for all isotopes by H2, we have 0 < ðX i Þ T Ç i for all isotopes. Hence, using theV i defined in the proof of Theorem 4.5, we have > 0, using the strict positivity of the main diagonal blocks of K 0 : Next, since U is maximal, we must have U U: If U 6 ¼ U, then there exists a g and zone Z i such that U g, i < U g, i : Since X i > 0 for all isotopes by H2, we have In a way similar to that done above, we have U ¼ N ðUÞ < N ð UÞ ¼ U: Next, for any 0 < t < 1, define UðtÞ tU þ ð1 À tÞ U, and let Ç ðtÞ ¼ OUðtÞ: Using the strict concavity of the functions f g, i ðxÞ, for each zone Z i AÛ >Û > 0; where A K 0 ðR 0 Þ À1R T f ðt 0Û Þ: Therefore, A kÛ >Û for all k ¼ 1, 2, Á Á Á : Hence, for all k kA k k 1 Á kÛk 1 ! kA kÛ k 1 > kÛk 1 > 0, or that kA k k 1 > 1 for all k. Thus, qðAÞ ¼ lim k!1 kA k k 1=k 1 ! 1, which is a contradiction since qðAÞ < 1 by (70). Thus,Û ¼ 0: w Theorem 4.9. Assume H1 and H2 hold for the discrete POI system (63). Assume also that the system is critical, i.e., k Ã ¼ 1, and for all g ¼ 1, Á Á Á , G and every isotope (49) holds withĈ i 0, g þĈ i 1, g < 1. Then (63) possesses only the trivial solution U ¼ 0: Proof. We argue again by contradiction. Suppose thatÛ is a nontrivial solution as in the previous proof. Then we must haveÛ > 0: From the strict concavity of the functions f g, i , we have N ðtÛÞ > tN ðÛÞ ¼ tÛ > 0 for all 0 < t < 1: Since 0 < f g, i ðxÞ < f g, i 0 ð0Þx for 0 < x 1, and f g, for all 0 < U 1: Since ðR 0 Þ À1 ! 0 and K 0 ! 0 with positive main diagonal blocks, it follows that N ðUÞ K 0 ðR 0 Þ À1 R T f U for all 0 < U 1: Next, let U Ã > 0 be the eigenvector corresponding to k Ã ¼ 1 given by Perron's theorem for the matrix ðK 0 ðR 0 Þ À1 R T f Þ T > 0: Then, for any 0 < t < 1 0 < tðU Ã Þ TÛ < ðU Ã Þ T N ðtÛÞ ðU Ã Þ T ðK 0 ðR 0 Þ À1 R T f ÞðtÛÞ ¼ tðU Ã Þ TÛ , which is a contradiction toÛ being nontrivial. Therefore,Û ¼ 0: w From the results presented in this subsection, it is clear that the Fixed Point POI iteration discussed above will converge if the Block-Jacobi iteration is performed using a tight enough tolerance, and the initial guess is chosen properly. From the proof of Theorem 4.6, the choice of 2 ð0, 0 Þ is akin to the procedure outlined after the presentation of the Fixed Point iteration. Basically, 2 ð0, 0 Þ in the theorem should be chosen so that We have not been able to verify that the procedure discussed after the Fixed Point iteration will always give an acceptable initial guess, but it seems to work well in practice.

Numerical tests
In this section we present results for several tests. We first verify our adjoint problem formulation, and then consider some POI test problems.
As the Newton-Krylov solver approach is quite different from the other two solvers, having all three solvers give the same solution on the same test problem is a good verification of the correctness of the solver implementations. The slab test problems were solved using MATLAB. Both the Petrov-Galerkin and Step methods were used, but we only report the Petrov-Galerkin methods here.

Adjoint problem tests
We present some tests that verify our adjoint formulation. Consider solving the steady state forward and adjoint problems with the same source and vacuum boundary conditions, namely, TW ¼ ZQ and T Ã W Ã ¼ ZQ: For the solutions W and W Ã , let W z ¼ SW and W Ã z ¼ SW Ã : Then we have Note that (71) holds for any source Q: We consider a slab of aluminum of width 12 cm and density 2.7 g=cm 3 : A spatial mesh of 40 zones was used for the slab À6: x 6:, 8 directions (i.e., S 8 ), P 2 scattering, and 18 energy groups with group boundaries given in Table 1 The source is a volume source contained in the aluminum ball, uniform in energy, isotropic in direction, and with an intensity of 1 mole/sec. Cross sections for all tests were obtained using the LLNL ENDL2009 evaluated nuclear database and access routines described in Beck (2002, McNabb 2003. For this test, we have hW Ã z , Qi ¼ hW z , Qi ¼ 18:38836 moles=sec: We have also run this test with the source restricted to only one energy group and obtained similar results, as well as with higher numbers of energy groups and various flux weightings when generating the cross sections.

POI solution implementation details
For the fixed point and nonlinear Power iteration methods discussed above, we iterate until the Euclidean norm of the relative change in the directionally integrated POI solution, i.e., Ð S 2 p g ðr, XÞdX, from one iterate to the next is less than or equal to a stopping tolerance of 10 À6 : More specifically, for iterates P l and P lþ1 , we calculate the Euclidean norm of the vector E whose components are given by For the Block-Jacobi iterations, we stop iterating when kW Ã, lþ1 À W Ã, l k 1 10 À8 Á kW Ã, l k 1 : We found it necessary to use a tighter stopping tolerance for the inner iteration. Otherwise, for some problems the outer iteration failed to converge.
For the Newton-Krylov method, we iterate until the Euclidean norm of the nonlinear residual, F ðP l Þ, is less than or equal to 10 À8 : Note that because of the vacuum boundary equations the nonzero terms in the residual all have the same units. To be on the conservative side, we use a slightly tighter tolerance for the Newton-Krylov iteration. As noted above, we use the KINSOL Newton-Krylov solver from the SUNDIALS suite (Hindmarsh et al. 2005). Each preconditioner solve is accomplished by 5 Block-Jacobi iterations. We used the GMRES solver option in KINSOL, with a maximum subspace size of 10, and linesearch backtracking turned on. Some runs were also done using the BiCGSTAB solver, but the statistics did not differ appreciably from those with GMRES, so we have not reported them. Our main interest in this study is robustness, and so we also have not reported run times.
In the tables below, FP refers to the fixed point iteration, NP to the nonlinear Power iteration, and NK to the Newton-Krylov solution approach. We report the total number of Block-Jacobi iterations for all three solvers, where for the Newton-Krylov solver this is obtained by multiplying the number of preconditioner solves by 5. We also report the total number of GMRES iterations for the Newton-Krylov solver.

POI tests
For this test we consider a slab of width 5.8252 containing Pu 239 at a density of 14.6 g=cm 3 : We used the 18 energy groups as in Table 1, 40 spatial zones, 4 directions (i.e., S 4 ), and isotropic scattering. The system is supercritical with a criticality value of k ¼ 1:129752: For this problem, condition (67) holds, and kðH Ã z Þ À1 k 1 ¼ 0:9999987, kR À1 L þ R T s Lk 1 ¼ 0:8405877, and kK 0 ðR 0 Þ À1R T f k 1 ¼ 0:5007058: Hence, both conditions H1 and H2 hold. Table 2 shows the results. While the nonlinear Power method and Newton-Krylov method performed well, the fixed point iteration took roughly ten times more iterations to converge. All three methods produced exactly the same POI solution (to 5 significant digits) on all the problems (with the one exception below where the fixed point method failed to converge).
For a second test, we changed only the density to 12.5 g=cm 3 : The system is now only slightly supercritical with a criticality value of k ¼ 1.003148. For this problem, condition (67) holds, and kðH Ã z Þ À1 k 1 ¼ 0:9999905, kR À1 L þ R T s Lk 1 ¼ 0:8405877, and kK 0 ðR 0 Þ À1R T f k 1 ¼ 0:4798445: Hence, both conditions H1 and H2 hold. We used all three solution methods on this second problem. Table 3 shows the results for each. Note that the fixed point iteration failed to converge after 400 iterations, while the nonlinear power method and the Newton-Krylov solver had no problems.
Our third test again changed only the density to 20.6 g=cm 3 : The system is now very supercritical with a criticality value of k ¼ 1.438712. For this problem, condition (67) holds, and kðH Ã z Þ À1 k 1 ¼ 1:000000, kR À1 L þ R T s Lk 1 ¼ 0:8405877, and kK 0 ðR 0 Þ À1R T f k 1 ¼ 0:6220586: Hence, both conditions H1 and H2 hold. All three solution methods performed well on this problem, and Table 4 shows the results. On these first three problems, we also tried using the initial guess from the fixed point iteration for the Newton-Krylov solver. The results did not change significantly, except for the second problem where the Newton-Krylov solver converged to the identically zero solution with the fixed point initial guess. The adjoint criticality eigenvector appears to be the best initial guess for the Newton-Krylov solver.
For our final test, we consider the same slab containing Pu 239 at a density of 14.6 g=cm 3 : Here, we used 175 energy groups, 40 spatial zones, 20 directions (i.e., S 20 ), and P 2 scattering. The system is supercritical with a criticality value of k ¼ 1.183168. We only used the nonlinear Power iteration for this test, with an outer iteration count of 14, and a total of 248 Block-Jacobi iterations. Figures 1 and 2 show contour plots of the scalar POI solution, ð1=DE g Þ Ð S 2 p g ðx, XÞdX, as a pointwise function of energy and space. The first plot focuses on the low energies, while the second focuses on the high energy dependence. Figure 3 shows the lower energy spectrum at a specific point in space.

Discussion
We have presented a matrix framework for the development and solution of the forward and adjoint multigroup neutron transport equations in slab  geometry, and then used this framework to develop a matrix form and solution methods for the neutron probability of initiation (POI) equations. We also presented sufficient conditions guaranteeing the existence and uniqueness of a nontrivial solution to the discrete POI system when it is supercritical, and only trivial solutions when it is critical or subcritical. Several tests were performed to verify the consistency of the adjoint development, and numerical results were presented for several POI problems. To our knowledge, we have not found these results presented in the literature. While the nonlinear Power iteration appears to work extremely well  for all of the POI problems we've tested it on, there appears to be a lack of theory discussing it's convergence properties. While our development has been restricted to the Petrov-Galerkin (or equivalenty, Diamond Difference in 1D slab geometry) and Step spatial discretizations, the matrix formalism easily extends to other spatial discretization approaches in slab and to higher spatial dimensions. We have implemented 1D spherical and 2D r-z geometry adjoint solvers using the same principles developed for slab and geometry, and the 1D, 2D and 3D results agree well on problems with the appropriate symmetries. As the scattering and fission operators are all local operators, the matrix formalism applies directly to these discrete operators in higher dimensions. Additionally, much of our analysis of the discrete equations does not depend on the particular spatial discretization. In a related technical report (Brown 2008), many of the results here have been extended to 3D Cartesian geometry. Finally, we also note that the discrete adjoint systems do not depend upon how the forward problem's cross sections were obtained.