Symmetry and separability of the neutron diffusion equation

Separation of variables is one of the oldest techniques for solving certain classes of partial differential equations (PDEs). As is the case with many other solution techniques for differential equations, separation of variables may be codified within the broader framework of symmetry analysis. Though the separation of variables technique is frequently used in the nuclear engineering context with various equations describing neutron transport, its connection to the symmetries of those equations has not yet been thoroughly established. It is thus the purpose of this work to establish that connection using neutron diffusion as both an initial step toward analysis of more generally applicable equations, and as a connection to previous results in related problems. Using Lie group analysis, it is found that the traditional space-time separable solution of the neutron diffusion equation (featuring a single α-eigenvalue) corresponds to time translation and flux scaling symmetries. Additional solutions of this equation are also constructed using its broader symmetry set.

Glasstone [15], Weinberg and Wigner [16], Lewins [17], and briefly by Keepin [18]. The same concepts in the context of neutron diffusion equations are motivated by Duderstadt and Hamilton [19], Hetrick [20], and many other texts that include a derivation of the point kinetics equations from diffusion theory. A time-separable assumption is key in these developments, as among other concepts it may be used to construct definitions of criticality. It can also be used in various schemes for the numerical solution of the neutron transport equation or its surrogates.
Symmetry analysis of the neutron transport equation or its surrogates has been conducted by, for example, Tsyfra and Czyzycki [21] and Giron et al [22]. Otherwise, the existing literature on the subject appears to be sparse, and aside from some commentary by Case and Zweifel [23] no specific connection appears to have been established between these results and those obtainable via the classical separation of variables technique. This outcome is the basis of the primary motivations of this work, which are twofold: (1) The explicit connection between the separation of variables ansatz commonly employed in the analysis of the neutron transport equation (or its surrogates) and the outcome of invariance under symmetry groups, (2) As a corollary to the above, an explicit characterization of the transformation properties of separable solutions of the neutron transport equation (or its surrogates), and the associated physical implications.
Along these lines, exhaustive symmetry analyses have been conducted on the linear heat equation (see, for example, Bluman and Anco [2] or Olver [5]), and categorization as outcomes of symmetry properties of its separable solutions in a variety of coordinate systems accomplished by Miller [8][9][10]. Given this existing body of work, a neutron diffusion equation-which is reducible to the linear heat equation by a point transformationis chosen as a natural starting point for the broader investigation. Given the nature of this equation, emphasis will be placed throughout on space-time (multiplicative) separability.
The results of this work are novel and expected to be useful in that (1) They will provide an explicit understanding of the Lie symmetries that result in the usual space-time separable solution of the neutron diffusion equation or, conversely, the symmetries that solution should be expected to possess (e.g., whether or not they are scale-invariant, for the motivation of scaled experiments or the assessment of counterpart computational results), (2) Scenarios will also be identified where the usual space-time separable solutions do not exist, but others may be available based on the symmetry properties of the underlying PDEs. While a variety of such solutions have been demonstrated to exist in the context of the linear heat equation, they have not been widely disseminated in the context of neutron diffusion. 'Unusual' solutions may thus prove useful as test or benchmark problems for the verification of codes used to simulate the underlying equations, (3) Finally, this work is intended to lay the groundwork for future studies of more complicated (and realistic) mathematical models, which may manifest similar properties to those considered here (given that neutron diffusion may be viewed as a limit of more complicated neutron transport processes).
Though these results are rightfully expected to bear much resemblance to similar outcomes derived in the context of the linear heat equation, they represent a new, site-specific application of general ideas in the spirit of the studies performed by Chou and Qu [12], Estevez et al [13], Polyanin and Zhurov [14], and many other researchers who have investigated the fit of separation of variables techniques within the broader symmetry analysis framework.
In support of the above, a brief review of the relevant mathematical model is provided in section 2 (included a self-contained primer in differential geometry, which facilitates some of the analysis which follows; this is also an original contribution of this work). Section 3 includes the necessary background and outcomes of symmetry analysis, further substantiated by the examples provided in section 4. Conclusions and recommendations for future study are provided in section 5. Some expanded solutions of the underlying mathematical model are provided in two appendices.

Mathematical model
Of principal interest to this work is the one-dimensional (1D) time-dependent neutron diffusion equation as discussed by Duderstadt and Hamilton [19], Hetrick [20], and many other authors, where r=generalized 1D space coordinate, t=time coordinate, j(r, t)=neutron scalar flux, n=0, 1, 2 for 1D planar, cylindrical, or spherical symmetries, l ∞ =infinite medium neutron lifetime, L 2 =diffusion area, Δk ∞ =excess infinite medium multiplication factor. The parameters l ∞ and L 2 will be taken as constants throughout this work. In section 4, two cases will be considered where Δk ∞ is allowed to vary as a function of time (consistent with similar work done in the point kinetics context by Hetrick [20], Keepin [18], and others).
For a well-posed problem, equation (1) must be supplemented by one initial condition and two boundary conditions. In general, these may be written as, respectively, where t 0 is an arbitrary reference time where j 0 (r) specifies the initial r-dependence of the flux, and r 1 and r 2 are arbitrary reference positions where j 1 (t) and j 2 (t) specify the t-dependence of linear combinations of the fluxes and their first derivatives, respectively (i.e., c 1 -c 4 are arbitrary constants). As written, equations (1)-(3) represent a 'monoenergetic' model that is regarded as the purely isotropic scattering and flux limit of the more complicated neutron transport equation. The material properties appearing in equation (1) are implicitly taken as averages with respect to a prescribed neutron energy spectrum, and the energy dependence of j is not explicitly represented; the scalar flux is assumed to be evaluated at either a discrete or averaged value of the neutron energy. While equation (1) may be of limited practical utility, the forthcoming analysis conducted in its context will reveal many concepts that may eventually be extended to more complicated and realistic models (e.g., multi-group P N or S N approximations).
In 1971, using the language of differential forms, Harrison and Estabrook [24] formulated an alternative approach to traditional symmetry analysis methods. In this so-called 'isovector' context, the usual infinitesimal Lie symmetry generators are interpreted as Lie derivative or dragging operators in a manifold spanned by all independent and dependent variables native to a problem under investigation. Any differential equations being analyzed are re-cast as exterior differential systems (EDS), thus allowing for the leveraging of the commutability between the Lie and exterior derivative operations (among other advantages).
The isovector method has since been developed in rigorous detail and applied to multitudinous problems in science and engineering by authors including Harrison [25,26], Edelen [27], and Suhubi [28]. It has been found that the principal advantages of the isovector method include (but are not necessarily limited to) a natural geometric setting from which to initiate symmetry analysis, and enabling a more elegant concurrent evaluation of extended infinitesimal symmetry generators in the analysis of differential structures.
In the spirit of this previous work: to write equation (1) in the language of differential forms, it must first be reduced to a first order system of PDEs. With Equations (4) and (5) 5 may then be multiplied by the differential volume element  dr dt to yield The procedure of recasting equation (1) as equations (4) and (5) is also widely used in the identification of potential symmetries of differential equations (see, for example, Bluman and Kumei [29]). However, the focus of the current work will be on the determination of point symmetries. Equations (6) and (7) are referred to as a system of 2-forms, where the operator  used to multiply differentials is known as the wedge product. Its relevant properties are for any variables x and y. As written, equations (6) and (7) represent a system of differential objects on a higherorder manifold where r, t, j, and w are regarded as entirely independent. To recover equations (4) and (5) from equations (6) and (7), let j=j(r, t) and w=w(r, t) (i.e., restricting to the sub-manifold or graph space in which equations (4) and (5) are defined). Expanding the total differentials and setting μ 1 =0 and μ 2 =0, equations (6) and (7) become  (8) and (9), the nontrivial solution of which is given by equations (4) and (5). Thus, equations (6) and (7) will be used for the symmetry analysis studies to follow.

Symmetry analysis
Symmetry analysis techniques as applied to the study of differential equations are in a sense a unification theory for disparate, otherwise ad hoc methods for the solution of those equations. The key point surrounding these techniques is that if a differential equation possesses symmetries, they enable a change of coordinates through which the original equation may either be reduced to a simpler structure [e.g., from a PDE to an ordinary differential equation (ODE)] or solved outright. As discussed in section 1, of principal interest to this work is the correspondence between symmetry reduction and separation of variables (which effectively reduces PDEs to systems of ODEs), in the context of equation (1). This exercise represents an application of the general ideas laid out by Miller [8][9][10] to a specific mathematical model, in the spirit of previous studies performed by Chou and Qu [12], Estevez et al [13], and Polyanin and Zhurov [14] for a variety of physically-motivated PDEs. Put simply, a differential equation (or corresponding collection of objects such as an EDS) has a symmetry when there exists a coordinate transformation that leaves it invariant, or unchanged in structure when written in terms of the new variables. For example, in the context of section 2, for a generalized set of point transformations parameterized in terms of a parameter a given by  (14)- (17) when  (14)-(17) (e.g., translations or scaling), equation (19) often proves straightforward to evaluate directly. This will obviously not be the case when equations (14)-(17) assume either nonlinear or generalized forms. In these cases, it often proves advantageous to formulate an infinitesimal representation of equation (19). If the transformations given by equations (14)- (17) possess an identity element a 0 such that f 1 (r, t, j; a 0 )=r, f 2 (r, t, j; a 0 )=t, etc, then the left side of equation (19) may be expanded in a Taylor series about that element: for all i, and F new taken as shorthand for the left side of equation (19). Using the chain rule, and the functions ξ i are referred to as coordinate functions, defined by With equation (21), the invariance condition given by equation (19) becomes (27) and (28) show that the linear infinitesimal representation of invariance is sufficient for the determination of the global invariance properties of the functions F i . If it is further required that F i =0 for all i when equation (28) is satisfied, the symmetries arising from the evaluation of equation (28) will correspond to those possessed by the differential equation structure or equivalent EDS under investigation.
The operator X (pr) is variously referred to as the prolonged vector field generated by the group of transformations given by equations (14)-(17) (or, more compactly, the 'prolonged group generator'). The terminology 'prolonged' refers to the presence of the basis element in the w-derivative; as w is interpreted as a derivative of j via equation (4), its presence in equation (22) represents the extension of a 'kernel' group generator to a higher-order manifold containing derivative information.
As discussed by Olver [5], Harrison and Estabrook [24], Harrison [25,26], Edelen [27], and Suhubi [28], given the suggestive vector basis notation in which the group generators X and X (pr) are written, they may also be interpreted as Lie derivative or dragging operators. This interpretation is convenient when coupled to an EDS structure, as exterior derivatives d and Lie derivatives X (prolonged or otherwise) commute: where F is any smooth function. As a result, the evaluation of structures such as equations (6) and (7) for invariance becomes straightforward: all original independent and dependent variables appearing in these relations are regarded as entirely independent of one another in the appropriate manifold, and any differential terms appearing within these expressions are treated using equation (30).

Similarity variables
Given a prolonged group generator of the form given by equation (22), a change of variables may be constructed from its kernel given by equation (29). The invariants of the group generated by equation (29) may be regarded as similarity variables in terms of which the governing differential equation structure may be reformulated. Assuming the group generated by equation (29) is nontrivial (i.e., at least two of the ξ i are nonzero), the change of variables arising from it leads to a reduction in order of the original differential equation. The invariants of equation (29) are determined by evaluating the condition where F is an arbitrary smooth function of the indicated arguments. As equation (31) is a quasi-linear first-order PDE, it may be solved using the Method of Characteristics. The characteristic equations associated with equation (31) are For physical problems it is often the case that ξ 1 =ξ 1 (r, t), ξ 2 =ξ 2 (r, t), and ξ 3 =ξ 3 (r, t, j). It is then convenient to solve the first two members of equation (32), and either the first and third or second and third members as a simultaneous system. Assuming this system can be solved in closed form (which depends on the exact form of the coordinate functions), that solution can be represented as const. and , , const., 33 1 2 where either the r or t dependence of the second member of equation (33) can be eliminated using the first member. The two constants of integration appearing in equation (33) are the invariants of the group generated by equation (29), and the functions g 1 and g 2 are in turn interpreted as a change of variables in terms of which the original differential equation may be reformulated.
In the current example, the first member of equation (33) leads to a new independent variable that is written in terms of the two original independent variables r and t. The second member of equation (33) defines a new dependent variable that is written in terms of the original dependent variable and one of the original independent variables.

Separation of variables
The classical (multiplicative; assumed throughout the remainder of this work) space-time separation of variables technique as applied to equation (1) assumes For certain forms of the material parameters l ∞ , L 2 , and Δk ∞ , the assumption of equation (34) allows for the decoupling of equation (1) into two ODEs: one each for R(r) and T(t). The reduction of the PDE in this manner is similar in outcome to the symmetry analysis procedure discussed in sections 3 and 3.1.
The group generator corresponding to the general group of point transformations admitted by equation (1) is given by equation (29). The goal is to determine the specific forms of the coordinate functions ξ 1 , ξ 2 , and ξ 3 (and their global/physical interpretations) such that a structure equivalent to equation (34) is recovered as the corresponding change of variables. In this case, it may then be shown that equation (34) follows from a subset of the full range of symmetries represented by equation (29).
To begin, equation (34) indicates that combinations of r and t (or functions of those combinations; e.g., r-ct, where c is a constant) do not appear in space-time separable solutions. However, since the invariant g 1 appearing in equation (33) must be a function of at least one original independent variable, it then follows that either g 1 =r or g 1 =t. With equation (32), the only possible cases giving rise to these solutions are Following the discussion surrounding equation (33), the first member of equation (35) leaves r as the new independent similarity variable, while the second member correspondingly leaves t. As a result the function g 1 appearing in equation (35) is trivial; the only achievable reduction using the symmetry analysis formalism is through the definition of a new dependent variable arising from g 2 (obtained using the nontrivial solutions for j of equation (35)).
From equation (35), the corresponding nontrivial ODEs for the new dependent similarity variables are where τ and ρ are arbitrary functions of the indicated arguments. The representation of equation (29) associated with equation (38) is not unique, but two physically intuitive possibilities are Equation (

Examples
Following the developments of section 3, a general prolonged group generator from which space-time separable solutions of equation (1) may be constructed is provided that, in the course of the analysis, one and only one of ρ(r) and τ(t) vanish. Utilizing the differential form representation, the condition that equations (6) and (7)  Assuming Δk ∞ =Δk ∞ (t) and using equation (30), equations (41) and (42) become, respectively,  Given that ρ is only a function of r and τ is only a function of t, equation (49) is satisfied only when The solution of equation (56) is given by where a 4 is an arbitrary constant of integration. Equation (57) represents the functional form in time the excess reactivity must take so as to enable either space/time scaling (a 1 ) or time translation (a 3 ) symmetries. If Δk ∞ (t) assumes any other functional form, it can still be made to satisfy the constraint given by equation (56) if a 1 =a 3 =0. In general, with the above equation (40) becomes

Constant material property general solution
For a constant excess multiplication factor D = ¹ ¥ ( ) ( ) k t const. 0, 59 a 1 =0 as appearing in equation (57), and Δk ∞ =a 4 /a 3 . The group generator X becomes indicating the presence of space and time translations (a 2 and a 3 , respectively) in addition to flux scaling. Equation (60) holds for planar symmetry (n=0); for cylindrical or spherical symmetries (n=1 or 2, respectively) a 2 =0 from equation (55). However, following from the developments of section 3.2, in order for equation (60) to produce similarity variables that result in a space-time separable solution, one of a 2 or a 3 must be zero (see appendix A for a brief discussion of the non-separable case where both a 2 and a 3 are non-zero). For example, let a 2 =0 (generalizing the cylindrical/spherical cases to planar geometry). Following the developments of section 3.1, the similarity variables associated with equation (60)

Linear transient general solution A linear excess multiplication transient is given by
for b=const.>0. Given that this form of Δk ∞ (t) does not correspond to equation (57), its consistency with equation (56) requires a 1 =a 3 =0. The associated group generator X becomes where equation (4) has also been used to express the r-derivative of j as w. For both members of equation (75) to collapse to ODEs for the functions j 1 (t) and j 2 (t), it is immediately evident that either a 1 =0, or c 2 =0 and c 4 =0. In either case, the solutions of equation (75) are then given by where c is an arbitrary constant. The boundary time-dependence of j must be of the form indicated by equation (76) for the space-time separable solution generated by equation (60) to exist.
As noted in section 4, for a space-time separable solution to exist, either ρ(r) or τ(t) as appearing in equation (40) must vanish. As a result, either a 1 r/2+a 2 =0 or a 1 t+a 3 =0 in equations (73)-(76), depending on both the underlying constraints of the specific problem being investigated, or the character of the desired space-time separable solution. These conditions more strongly constrain the allowable forms of j 0 , j 1 , and j 2 appearing in equations (74) and (76).
Appendices A and B include examples of application of these boundary conditions to some problems of interest.

Conclusions
The space-time separable solution of the neutron diffusion equation (1) is pervasive in nuclear engineering texts, but its connection to the broader outcomes of symmetry analysis of that equation is heretofore unexplored directly, well-established connections to the linear heat equation notwithstanding. In the current work, spacetime separable solutions of equation (1) are found to correspond to an infinitesimal symmetry generator of the form given by equation (40). The symmetries that generate space-time separable solutions include scaling in the dependent flux variable plus (1) A transformation in the independent space variable (r) that depends only on r, but is otherwise arbitrary, or, (2) A transformation in the independent time variable (t) that depends only on t, but is otherwise arbitrary.
The correspondences noted above are constructed by viewing the usual space-time separable ansatz of equation (1) as a change of variables that reduces the original PDE to two ODEs. This outcome is similar to that of symmetry analysis, which also results in the reduction of more 'complicated' structures to ostensibly simpler forms. By interpreting the separability assumption given by equation (34) as a similarity transformation, it may be viewed as originating from a group generator structure that is itself a representation of a symmetry.
Two examples of this phenomenon are provided in sections 4.1 and 4.2, indicating that the usual space-time separable solution of equation (1) follows from Case (2) above, with the transformation in time corresponding to translation. However, neither this solution nor its feature of space-time separability is unique. Though this outcome has not been widely disseminated, it is not entirely unexpected in light of the multitudinous solutions of the linear heat equation as derived, for example, by Bluman and Anco [2] and Olver [5]. The practical relevance or computational utility of these additional solutions remains to be investigated.

Recommendations for future work
The practical utility of the mathematical model represented by equation (1) is limited, as discussed in section 2. However, the results of this work may likely be extended to a variety of scenarios, including those most accurately modeled by the linear neutron transport equation: are often employed in the analytical treatment or numerical solution of equation (77). As seen in section 4.1, the existence of a time-separable solution of this form is connected with the invariance of equation (77) under time translation and flux scaling symmetries: If the material parameters appearing in equation (77) (v, χ, ν, and the various Σ's) do not depend on t, it is clear by inspection that equation (77) admits these symmetries (equation (77) also admits space translation symmetries if its material properties do not depend on r). A more rigorous investigation of this outcome in the context of the broader symmetry set admitted by equation (77) constitutes a logical course of future study.
Moreover, while spherical harmonics expansions given by å å y  (77), their explicit connection to a symmetry reduction has yet to be investigated. The same is true of S N methods or multi-group methods for handling the angular and energy-dependence of equation (77), respectively; given that these methods are approximate, their connection to exact symmetries may not be as straightforward to apprehend. The solution indicated by equation (A.10) is not space-time separable in general. It instead includes terms of the form exp(r±ct), where the propagation speed c of these traveling wave type solutions is comprised of the various material constant data appearing in the problem formulation. In equation (A.10), these traveling wave components are multiplied by the global scale factor exp(αt), which is typical of parabolic PDEs. These types of solutions have been shown to exist for diffusion-type equations by Olver [5], and their physical interpretation has been discussed by Moore [30]. In short, it appears solutions of this type feature not only the usual global time dependence, but some traveling wave-like structure.
The traditional separable solution of equation (1)  where i is the imaginary unit. Naturally, equation (A.10) may also be found using classical methods for solving linear PDEs. Perhaps the only advantage to constructing such solutions using symmetry analysis methods lies in the ability to a priori build solutions with desired properties (in the current case, the translation and scaling symmetries from which the solution follows).
Numerical results from equation (A.11) are readily obtainable after application of the correct boundary conditions. Specifically, consider the alternative representation of equation (A.11):  factor of the order exp(t 2 ) is present. This 'faster' rate of rise in j owes to the time dependence imposed on Δk ∞ . Since the excess multiplication is increasing linearly in time, the associated flux rises at a concordantly greater rate.
As was the case with the solution derived in appendix A, equation (B.5) may also be found using classical methods for solving linear PDEs. The potential advantage in using symmetry analysis techniques is the same as discussed previously.
Finally, as demonstrated in appendix A, the invariant solution given by equation (B.5) can be expressed numerically given the correct boundary conditions. For example, for an insertion rate given by b=0.0001 s −1 , ν=2.35 neutrons per fission, l ∞ =10 −6 s, diffusion coefficient D=1.1 cm, Σ f =0.073 cm −1 , and β=1 at thermal neutron energies, it follows that Σ a =1/(vl ∞ ), and L 2 =D/Σ a so that