Application of Finite Symmetry Groups to Reactor Calculations

Group theory is a vast mathematical discipline that has found applications in most of physical science and particularly in physics and chemistry. We introduce a few of the basic concepts and tools that have been found to be useful in some nuclear engineering problems. In particular those problems that exhibit some symmetry in the form of material distribution and boundaries. We present the material on a very elementary level; an undergraduate student well versed in harmonic analysis of boundary value problems should be able to easily grasp and appreciate the central concepts.


Introduction
Group theory is a vast mathematical discipline that has found applications in most of physical science and particularly in physics and chemistry. We introduce a few of the basic concepts and tools that have been found to be useful in some nuclear engineering problems. In particular those problems that exhibit some symmetry in the form of material distribution and boundaries. We present the material on a very elementary level; an undergraduate student well versed in harmonic analysis of boundary value problems should be able to easily grasp and appreciate the central concepts.
The application of group theory to the solution of physical problems has had a curious history. In the first half of the 20th century it has been called by some the "Gruppen Pest" , while others embraced it and went on to win Noble prizes. This dichotomy in attitudes to a formal method for the solution of physical problems is possible in light of the fact that the results obtained with the application of group theory can also be obtained by standard methods. In the second half of the 20th century, however, it has been shown that the formal application of symmetry and invariance through group theory leads in complicated problems not only to deeper physical insight but also is a powerful tool in simplifying some solution methods.
In this chapter we present the essential group theoretic elements in the context of crystallographic point groups. Furthermore we present only a very small subset of group theory that generally forms the first third of the texts on group theory and its physical applications. In this way we hope, in short order, to answer some of the basic questions the reader might have with regard to the mechanical aspects of the application of group theory, in particular to the solution of boundary value problems in nuclear engineering, and the benefits that can accrue through its formal application. This we hope will stimulate the reader to look more deeply into the subject is some of the myriad of available texts.
The main illustration of the application of group theory to Nuclear Engineering is presented in Section 4 of this chapter through the development of an algorithm for the solution of the neutron diffusion equation. This problem has been central to Nuclear Engineering from the very beginning, and is thereby a useful platform for demonstrating the mechanics of bringing group theoretic information to bear. The benefits of group theory in Nuclear Engineering are not restricted to solving the diffusion equation. We wish to also point the interested reader to other areas of Nuclear Engineering were group theory has proven useful.
An early application of group theory to Nuclear Engineering has been in the design of control systems for nuclear reactors (Nieva, 1997). Symmetry considerations allow the decoupling of the linear reactor model into decoupled models of lower order. Thereby, control systems can be developed for each submodel independently.
Similarly, group theoretic principles have been shown to allow the decomposition of solution algorithms of boundary value problems in Nuclear Engineering to be specified over decoupled symmetric domain. This decomposition makes the the problem amenable to implementation for parallel computation (Orechwa & Makai, 1997).
Group theory is applicable in the investigation of the homogenization problem. D. S. Selengut addressed the following problem (Selengut, 1960) in 1960. He formulated the following principle: If the response matrix of a homogeneous material distribution in a volume V can be substituted by the response matrix of a homogeneous material distribution in V, then there exists a homogeneous material with which one may replace V in the core. The validity of this principle is widely used in reactor physics, was investigated applying group theoretic principles (Makai, 1992), (Makai, 2010). It was shown that Selengut's principle is not exact; it is only a good approximation under specific circumstances. These are that the homogenization recipes preserve only specific reaction rates, but do not provide general equivalence.
Group theory has also been fruitfully applied to in-core signal processing (Makai & Orechwa, 2000). Core surveillance and monitoring are implemented in power reactors to detect any deviation from the nominal design state of the core. This state is defined by a field that is the solution of an equation that describes the physical system. Based on measurements of the field at limited positions the following issues can be addressed: 1. Determine whether the operating state is consistent with the design state. 2. Find out-of-calibration measurements. 3. Give an estimate of the values at non-metered locations. 4. Detect loss-of-margin as early as possible. 5. Obtain information as to the cause of a departure from the design state.
The solution to these problems requires a complex approach that incorporates numerical calculations incorporating group theoretic considerations and statistical analysis.
The benefits of group theory are not restricted to numerical problems. In 1985 Toshikazu Sunada (Sunada, 1985) made the following observation: If the operator of the equation over a volume V commutes with a symmetry group G, and the Green's function for the volume V is known and volume V can be tiled with copies tile t (subvolumes of V), then the Green's function of t can be obtained by a summation over the elements of the symmetry group G. Thus by means of group theory, one can separate the solution of a boundary value problem into a geometry dependent part, and a problem dependent part. The former one carries information on the structure of the volume in which the boundary value problem is studied, the latter on the physical processes taking place in the volume. That separation allows for extending the usage of the Green's function technique, as it is possible to derive Green's functions for a number of finite geometrical objects (square, rectangle, and regular triangle) as well as to relate Green's functions of finite objects, such as a disk, or disk sector, a regular hexagon and a trapezoid, etc. Such relations are needed in problems in heat conduction, diffusion, etc. as well.
For illustrative purposes let us consider a simple set of three elements {E, A, B}. A law of composition for these three elements can be expressed in the form of a multiplication table, see Table 1. In position i, j of Table 1. we find the product of element i and element j with the numbering 1 → E,2 → A,3 → B. From the table we can read out that B = AA because element 2, 2 is B and the third line contains the products AE, AA, AB. The table is symmetric therefore AB = BA. Such a group is formed for example by the even permutations of three objects: E =(a, b, c), A =(c, a, b), B =(b, c, a). The multiplication table reflects four necessary conditions that a set of elements must satisfy to form a group G. These four conditions are: 1. The product of any two elements of G is also an element of G. Such as for example AB = E. 2. Multiplication is associative. For example (AB)E = A(BE). 3. G contains a unique element E called the identity element, such that for example AE = EA = A and the same holds for every element of G. 4. For every element in G there exists another element in G, such that their product is the identity element. In our example AB = E therefore B is called the inverse of A and is denoted B = A −1 .
The application of group theory to physical problems arises from the fact that many characteristics of physical problems, in particular symmetries and invariance, conform to the definition of groups, and thereby allows us to bring to bear on the solution of physical problems the machinery of abstract group theory.
For example, if we consider a characteristic of an equilateral triangle we observe the following with regard to the counter clockwise rotations by 120 degrees. Let us give the operations the following symbols: E-no rotations, C 3 -rotation by 120 o , C 3 C 3 = C 2 3 -rotation by 240 o . The group operation is the sequential application of these operations, the leftmost operator should be applied first. The Table 4. Multiplication table of a permutation group group G = {E, C 3 , C 2 3 }. We see immediately that the multiplication table of the rotations of an equilateral triangle is identical to the multiplication table of the previous abstract group G = {E, C 3 , C 2 3 }. Thus there is a one-to-one correspondence (called isomorphism) between the abstract group of the previous example, and its rules.

Subgroups and classes
Groups can have more properties than just a multiplication table. The illustration in the previous subsection is not amenable to illustrating this; the groups are to small. However, if we again consider the equilateral triangle, we note that it has further symmetry operations. Namely, those associated with reflections through a vertical plane through each vertex. Let us give these reflection operations the symbols σ v , σ ′ v and σ" v , reflection through planes through vertex a, b, c, respectively. We may describe the operation by the transformations of the vertices a, b, c. For example σ v : (a, b, c) → (a, c, b). By adding these three reflections operations to the rotations, we form the larger group of symmetry operations of the equilateral The multiplication table of the new group is given in  Table 3. Another group with the same multiplication table as above can be constructed by considering the six permutations of the three letters a, b, c. This leads to the multiplication table 4, which is isomorphic to Table 3. The two multiplication tables illustrate the concept of a subgroup that is defined as: A set S of elements in group G is considered as a subgroup of G if: 1. all elements in S are also elements in G 2. for any two elements in S their product is in S 3. all elements in S satisfy the four group postulates.
From the presented multiplication tables we see that {E, C 3 , c 2 3 } and {E, A, B} are subgroups; they are the only subgroups in their respective groups. These subgroups are the rotations in the former example, and the even permutations in the latter. While the remaining elements are associated with reflections and odd permutations, respectively. Furthermore, we see that two reflections are equivalent to a rotation, and two odd permutations are equivalent to an even permutation. Thus the operations {C 3 , C 2 3 } and σ v , σ ′ v , σ" v and similarly {A, B} and {a ′ , b ′ , c ′ } belong in some sense to different sets. This property is illustrated by taking the transform we obtain the following table of the transforms T −1 QT: themselves and are thereby called classes. Classes play a leading role in the application of group theory to the solution of physical problems. In general physically significant properties can be associated with each class. In the solution of boundary value problems, different subspaces of the solution function space are assigned to each class.

Group representations
The application of the information in an abstract group to a physical problem, especially to the calculation of the solution of the boundary value problem that models the physical setting, requires a mathematical "connection" between the two. This connection originates with the transformations of coordinates that define the symmetry operations reflected in the actions of a point group.
As a simple illustration, let us again consider the abstract group G = {E, A, B, a ′ , b ′ , c ′ } in the form of its realization in forms of rotations and reflections of an equilateral triangle, namely the point group Let this group be consistent with a physical problem in terms of, for example, material distribution and the geometry of the boundary. Furthermore, let us consider a two-dimensional vector space with an orthonormal basis {e 1 , e 2 } relative to which the physical model is defined. Each operation by an element g of the group C 3v can be represented by its action on an arbitrary vector r in a two-dimensional vector space. In the usual symbolic form we have r ′ = D(g)r for all g ∈ C 3v and where r ′ = r 1 e 1 + r 2 e 2 is the transformed vector, and D(g) is the matrix operator associated with the action of group element g ∈ C 3v . It is well known from linear algebra that the matrix representation of operator D(g) for each g ∈ C 3v is obtained by its action on the basis vectors, and that the transpose of matrix D ji (g) gives the action of the group element g on the coordinates of the vector r as For the point group C 3v we obtain the following six matrix representations. To spare room we replace the matrices by permutations: These matrices satisfy the group multiplication table of C 3v , and therefore also the multiplication table of the abstract group G that is isomorphic to C 3v . We note that this is not the only matrix representation of C 3v . There are two one-dimensional representations, in particular that also satisfy the multiplication table of C 3v , and will be of interest later. These are where E 2 is the 2 × 2 identity matrix; and The role played by these representations will become clear in later discussions of irreducible representations of groups, and their actions on function spaces.

Generation of group representations
To this point we have constructed the matrix representations the group elements of point groups such as C 3v in the usual physical space (two dimensional in our case). These representations were based on the transformations of the coordinates of an arbitrary vector in a physical space due to physical operations on the vector. Mathematical solutions to physical problems, however, are represented by functions in function spaces whose dimensions are generally much greater than three. Thus to bring the group matrix representations that act on coordinates to bear on the solution of physical problems in terms of functions, we need one more "connection" between symmetry operators on coordinates and symmetry operators on functions. This connection is defined as follows.
Let f (r) be a function of a position vector r =(x, y) and D(g −1 ) be the matrix transformation associated with group element g ∈ G, such that (x, y) → (x ′ , y ′ ) through What we need is an algorithm that uses D −1 (g) to obtain a new function h(r) from f (r).T o this end we define an operator O g as That is, operator O g gives a new function h(r) from f (r) at r, while f is unchanged at r ′ . For example, let f (x, y)=ax + by and For two group elements g 1 and g 2 in G, we obtain Note: operator O g acts on the coordinates of function f and not on the argument of f . Therefore and thus we get in words: the consecutive application of O g1 and O g2 is the same as the application of the transformation O g2g 1 belonging to group element g 1 g 2 , and the operators O g , g ∈ G have the same multiplication table as G and any group isomorphic with G.

Invariant subspaces and regular representations
A common approach to the solution of physical problems is harmonic analysis, where a solution to the problem is sought in terms of functions that span the solution space. If the problem exhibits some symmetry, we would expect this symmetry to be reflected in the solution for this particular problem. Intuitively we would expect therefore the solution to belong to a subspace of the general solution space, and that the subspace be invariant under the symmetry operations exhibited by the problem.
As an illustration of this notion, we assume the problem has the symmetry of the cyclic permutation group C 3 = {E, C 3 , C 2 3 } that was discussed previously. Let f E (r) be an arbitrary function that allows the operation of the operators in the group C 3 as discussed above. The action of each operator on f E defines a new function that, is Based on this and the group multiplication table we get relations such as etc. These observations can be summarized in a table: From that table we can construct matrix (permutation) representations of the operators O E , O C 3 , O C 2 3 as for example D(C 3 )=(2, 3, 1). (2.8) This procedure gives the so-called regular representation for the group C 3 as (2.9) The matrices, in general, satisfy the group multiplication table, and are characterized by only the one integer one in each column, the rest zeros, and the dimension of the matrix equals to the number of elements in the group. The functions f E , f C 3 , f C 2 3 that generate the regular representation, span the invariant subspace. They are not necessarily linearly independent basis functions.

Complete sets of linearly independent basis functions and irreducible representations
As was mentioned at the outset, symmetry as exemplified through group theory brings added information to the solution of physical problems, especially in the application of harmonic analysis. The heart of this information is encapsulated in the so called irreducible representations of the group elements. It should be stated at the outset that the irreducible representations used in most applications are readily available in tabulated form. Yet much of mathematical group theory is devoted to the derivation and properties of irreducible representations. We do not minimize in any way the importance of that material; it is necessary for a clear understanding of the applicability of the mathematical machinery and its physical interpretation. Our objective here is only to touch on a few of the central results used in the applications. Perhaps this may motivate the reader to look further into the subject.
The key property for the application of point groups to physical problems is that for a finite group all representations may be "built up" from a finite number of "distinct" irreducible representations. The number of distinct irreducible representations is equal to the number of classes in the group. Furthermore, the regular representation contains each irregular representation a number of times equal to the number of dimensions of that irreducible representation. Thus, if ℓ α is the dimension of the α-th irreducible representation, where |G| is the order of the group G to be satisfied.
Let us illustrate this with the group C 3 that was discussed previously. To identify the classes in C 3 , as before, we compute a table of T −1 QT, see Table 5. The elements that transform into can only be satisfied by ℓ 1 = ℓ 2 = ℓ 3 = 1. Therefore, there are three distinct one-dimensional representations. These are the building blocks for decomposing the regular representation to irreducible representations, and can be found in tables: where ω = exp(2πi/3). The element in each of the three irreducible representation conform to the multiplication of point group C 3 .
These low dimension irreducible representations are used to build an irreducible representation from the regular representation of the operator O C 3 for example, as follows.
The regular representation has the form of a full matrix, The irreducible representation has the form of a diagonal (block diagonal in the general case) matrix, ⎛ The mathematical relationship is discussed at length in all texts on the subject, and will not be repeated here. We assume the irreducible representations are known. Of interest is the information for the solution of physical problem, that is associated with irreducible representations.
Recall that starting with an arbitrary function f (r) belonging to a function space L (a Hilbert space for example), we can generate a set of functions f 1 ,..., f |G| that span an invariant subspace L s ⊂ L. This process requires the matrices of coordinate transformations g 1 ,...,g |G| that form the symmetry group G of interest. The diagonal structure of the irreducible representations of G tells us that there exists a set of basis functions { f 1 , f 2 ,..., f n } that split the subspace L s further into subspaces invariant under the symmetry group G, and are associated with each irreducible representation D (1) (g), D (2) (g),...,D (n c ) (g) where n c is the number of classes in G. That is and thus an arbitrary function f (r) ∈ L s is expressible as a sum of functions that act as basis function in the invariant subspaces associated with each irreducible representation D (α) (g), α = 1,...,n c as If the decomposition of the regular representation contains irreducible representations of dimension greater than one, we have for each basis function that "belongs to the α-th irreducible representation" where ℓ α is the dimension of the α-th irreducible representation.
The question now remains how do we obtain f α (r), the basis function of each irreducible representation?
To this end we can apply a projection operator that resolves a given function f (r) into basis functions associated with each irreducible representation. This projection operator is defined as The information needed to construct this operator-the coordinate transformations, the irreducible representations-are known in the case of the point groups encountered in practice. So, for example, the i-th basis function of the α irreducible representation that is ℓ α dimensional for a symmetry group with |G| elements is constructed from an arbitrary function This decomposition creates a complete finite set of orthogonal basis functions.
In practice, a more simple projection operator is generally sufficient. This is due to the fact that the D α ii (g)'s (the diagonal elements of a multidimensional irreducible representation) are quantities that are intrinsic properties of the irreducible representation D α (g). That is they are invariant under the change of coordinates.
Furthermore, the sum of the diagonal elements, or trace, of the irreducible representation D α (g) is also invariant under a change of coordinates. In group theory this trace is denoted by the symbol χ α (g) and and referred to as the character of element g ∈ G in the α-th irreducible representation. There are tables of characters for all the point groups of physical interest.
The projection operator in terms of characters is given as so that the basis functions are and f (r) is decomposed into a complete finite set of orthogonal functions, with one for each irreducible representation irrespective of its dimension.

Symmetries of a boundary value problem
Let us consider the following boundary value problem: where A and B are linear operators. Group theory is not a panacea to the solution of boundary value problems; its application is limited. The main condition that must be met in nuclear engineering problems is that material distributions have symmetry. This is generally true in reactor cores, core cells and cell nodes.
In the following we give a heuristic outline of how the machinery presented above enters into the solution algorithm of a boundary value problem, and what benefits can be expected.
Symmetry is the key. If we have determined that the physical problem has symmetries these symmetries must form a group G. The symmetry operator O g must commute for all g ∈ G with the linear operators A and B for group theory to be applicable. That is must hold for all g ∈ G. If this condition is met, the boundary value problem can be written as We can now use the projection operator (2.20) to form a set of boundary value problems Since the projection operator creates linearly independent components, we have decomposed the boundary value problem into a number (equal to the number of irreducible components) of independent boundary value problems. These are whose solution ψ α (r) belongs to the α-th irreducible representation. From this complete set of linearly independent orthogonal functions we reconstruct the solution to the original problem as where n c is the number of classes in G.
Why is this better? Recall that we are applying harmonic analysis. The usual approach is to use some series that forms an incomplete set of expansion functions and results a coupled set of equations; one "large" matrix problem. With group theory, we find a relatively small set of complete basis functions that form the solution from symmetry considerations. These are found by solving a set of "small" boundary value problems. It is clear that the effectiveness of group theory is problem dependent. However, experience over the past half century has proven group theory's effectiveness in both nuclear engineering and other fields.
We present an especially simple example (Allgover et al., 1992) that demonstrates the advantages of symmetry considerations. The example is the solution of a linear system of equations with six unknowns: (3.11) The example has been constructed so that the basis of the reduction is the observation that the matrix is invariant under the following permutations: p 1 =( 1, 6)(2, 5)(3, 4) and p 2 = (1, 5, 3)(2, 6, 4).A sp 1 and p 2 generate a group D 6 of six element, the matrix commutes with the representation of group D 6 by matrices of order six. This suggests the application of group theory: decompose the matrix and the vector on the right hand side of the equation into irreducible components, and solve the resulting equations in the irreducible subspaces. The D 6 group is isomorphic to the symmetry group of the regular triangle discussed in Section 2.2.
The character table of the group D 6 can be found in tables (Atkins, 1970;Conway, 2003;Landau & Lifshitz, 1980), or, can be looked up in computer programs, or libraries (GAP, 2008).
Using the character table, and projector (2.17), one can carry out the following calculations. The observation that D 6 is isomorphic to the symmetry group of the equilateral triangle makes the problem easier. (Mackey, 1980) has made the observation: There is an analogy of the group characters and the Fourier transform. This allows the construction of irreducible vectors by the following ad hoc method. Form the following N-tuples (N = |G|): e 2k−1 =( cos(2π/N * (2k − 1) * 1), . . . , cos(2π/N * (2k − 1) * N), e 2k =( sin(2π/N * (2k) * 1), . . . , sin(2π/N * (2k) * N), k = 1, 2, . . . N. (3.12) These vectors are orthonormal and can serve as an irreducible basis. After normalization, one gets a set of irreducible vectors in the N copies of the fundamental domain. Here one may exploit the isomorphism with the symmetry group of an equilateral triangle with the points positioned as shown in Fig. 1. Applying the above recipe to the points in the triangle, we get the following irreducible basis: e 1 =(1, 1, 1, 1, 1, 1) e 2 =(2, −1, −1, 2, −1, −1) e 3 =(0, 1, −1, 0, 1, −1) (3.13) e 4 =(2, 1, −1, −2, −1, 1) e 5 =(0, 1, 1, 0, −1, −1)) e 6 =(1, −1, 1, −1, 1, −1). (3.14) We note that the points in the vectors e i do not follow the order shown in Fig. 1. Thus we need to renumber the points, and normalize the vectors. For ease of interpolation we also renumber the vectors given above. It is clear that the vectors formed from cos and sin transform together. Thus they form a two-dimensional representation. We bring forward the one-dimensional representations. The projection to the irreducible basis is through a 6 × 6 matrix that contains where the prime indicates rearranging in accordance with Fig. 1. Using the rearranging 2 1 00000 0 −1 0000 00−62 a 00 00−a −10 0 0000 −62 a 0000 −a −1 where a = √ 3. Compare the structure of the above matrix with that given in Section 3, where the similar form is achieved by geometrical similarity. In the present example there is no geometry, just a matrix invariant under a group of transformations.
In order to solve the resulting equations, we need the transformed right hand side of the equation: Finally, note that instead of solving one equation with six unknowns, we have four equations, two of them are solved by one division for each, and we have to solve two pairs of equations with two unknowns for each. At the end, we have to transform back from Ox to x.
The Reader may ask: What is the benefit of the reduction? In a problem which is at the verge of solvability, that kind of reduction may become important.
A more favorable situation is when there are geometric transformations leaving the equation and the volume under consideration, invariant. But before immersing into the symmetry hunting, we investigate the diffusion equation.

The multigroup diffusion equation
The diffusion equation is one of the most widely used reactor physical models. It describes the neutron balance in a volume V, the neutron energy may be continuous or discretized (multi group model). The multi group version is: where the processes leading to energy change are collected in T kk ′ : where subscripts k, k ′ label the energy groups, v k is the speed of neutrons in energy group k, Ψ k (r) is the space dependent neutron flux in group k, and k ef f = 1. In general, the cross-sections D k , Σ tk , Σ k ′ →k , Σ fk ′ are the space dependent diffusion constant, the total cross-section, the scattering cross-section, and the fission cross-section. χ k is called the fission spectrum. Equation (4.1) is a set of partial differencial equations, to which the initial condition Ψ k (r,0), r ∈ V and a suitable boundary condition, e.g. Ψ k (r, t), r ∈ ∂V are given for every energy group k and every time t. The boundary conditions used in diffusion problems are of the type (∇n)Ψ k (r)+b k (r)Ψ k (r)=h k (r) k = 1,...,G. for r ∈ ∂V. Here b k (r) depends on the boundary condition and may contain material properties, for example albedo.
The diffusion equation is a relationship between the cross-sections in V and the neutron flux Ψ k (r, t). The equation is linear in Ψ k (r, t). The main variants of equation (4.1) that are of interest in reactor physics are: 1. Static eigenvalue problem: When the flux does not depend on t, the left hand side is zero, and (4.1) has a nontrivial solution only if the cross-sections are interrelated. To this end, we free k ef f and the static diffusion equation is put in the form of an eigenvalue problem: where the eigenvalue k ef f introduced as a parameter in T kk ′ thus allowing for a non-trivial solution Ψ k (r). That usage is typical in core design calculations. 2. Time dependent solution allowing time dependence in some cross-sections. A typical application is transient analysis. 3. Equation (4.1) is homogeneous but it is possible to add an external source and to seek the response of V to the source.
The structure of the diffusion equation is simple. Mathematical operations, like summation and differentiation, and multiplication by material parameters (cross-sections) are applied to the neutron flux. In such equations the symmetries are mostly determined by the space dependence of the material properties. In the next subsection we investigate the possible symmetries of equation (4.4) and the exploitation of those symmetries.
When the solutions Ψ k (r), k = 1, . . . , G are known, not only the reaction rates, and netand partial currents can be determined, but also matrices can be created to transform these quantities into each other. From diffusion theory it is known that the solution is determined by specifying the entering current along the boundary ∂V. Thus the boundary flux is also determined. But the given boundary flux also determines the solution everywhere in V. The solution is given formally by a Green's function as follows: Here G k0,k (r 0 → r) is the Green's function, it gives the neutron flux created at point r in energy group k by one neutron entering V at r 0 in energy group k 0 ; and f k 0 (r 0 ) is the given flux in energy group k 0 at boundary point r 0 . Similarly the net current is obtained as where the ∇ operator acts on variable r.

Symmetries of the diffusion equation
First, the symmetry properties of the solution do not change in time because (4.1) is linear. This is not true for nonlinear equations. Secondly, the equations in (3.3) need to be satisfied. That is, the operations of the equation (4.4) and the boundary conditions must commute with the symmetry group elements. The symmetries of equation (4.4) are determined by the operators, the material parameters (cross-sections) and the geometry of V. The first term involves derivatives: Here the first term contains a dot product which is invariant under rotations and reflections. The second term involves the laplace operator, which is also invariant under rotations and reflections. Thus, the major limiting symmetry factors are the material distributions, or the associated cross-sections as functions of space, and the shape of V. We assume the material distribution to be completely symmetric, thus for any cross-section Σ(r) we assume the transformation property Here O g is an operator applicable to the possible solutions. D(g) is a matrix representation of the symmetry group of the diffusion equation applicable to r. The following operators are encountered in diffusion theory. The general form of a reaction rate at point r ∈ V can be expressed as Here subscript 1 refers to the symmetric component. Since because the material distribution is assumed symmetric hence O g Σ(r)=Σ(r) for every symmetry g, the transformation properties of a reaction rate are completely determined by the transformation properties of the flux Ψ k1 (r). The normal component of the net current at r ∈ ∂V is J nk (r)=−D k (r)(n∇)Ψ k (r), (4.9) where n is the normal vector at r. We apply O g to J nk (r) to obtain: Thus, the transformation properties of the normal component of the net current agree with the transformation properties of the flux. In diffusion theory, the partial currents are defined as From (4.9) it follows that the transformation properties of the partial currents correspond to the transformation properties of the flux.
The boundary condition (4.3) commutes with rotations and reflections provided the material properties do. The same is true for the diffusion equation (4.1). Our first conclusion is that the material distribution may set a limit to the symmetry properties. As to the symmetries, the volume V under consideration may also be a limiting factor. Let O g be an operator that commutes with the operations of the diffusion equation (4.1) and (4.3). Furthermore, the representation D(g) maps V into itself. The set of operators form a group; the group operation is the repeated application. That group is called the symmetry group of the diffusion equation.
Example 4.1 (Symmetries in a homogeneous square). This symmetry group has eight elements, four rotations: E, C 4 , C 2 4 , C 3 4 and four reflections σ x , σ y , called of type σ v and σ d1 , σ d2 called of type σ d . Characters of a given class have identical values. This group is known as the symmetry group of the square and denoted as C 4v . The first column of a character table gives a mnemonic name to each representation, and a typical expression transforming according to the given representation. The first line is reserved for the most symmetric representation called unit representation. From the character table of the group C 4v , we learn that there are groups with the same character tables, there are five irreducible representations labeled A 1 , A 2 , B 1 , B 2 , E where As and Bs are one.dimensional and E is two-dimensional, it has two linearly independent components transforming as the x and y coordinates.
Example 4.2 (Symmetries in a homogeneous equilateral triangle). The group has six elements, three rotations: E, C 3 , C 2 3 , and three reflections through axis passing one edge: σ a , σ b , σ c called type σ v . The symmetry group is isomorphic to the C 3v group and its character table is the same as that of the group D 3 . The C 3v group is the symmetry group of the equilateral triangle, it has two one-dimensional and one two-dimensional representations.
The key observation concerning the applications of symmetry considerations in boundary value problems is as follows. For a homogeneous problem (4.4) where there is no external source, the boundary condition is homogeneous, and every macroscopic cross-section Σ(r), r ∈ V is such that O g Σ(r)=Σ(r) for all O g mapping V into itself. When the boundary conditions h k (r) in the expressions (4.3) transform according to an irreducible subspace f α (r) then the neutron flux Φ(r), the partial currents I(r), J(r), the reaction rate all transform under the automorphism group of V as do the boundary conditions h k (r).
The symmetry group of the volume V makes it possible to reduce the domain on which we have to determine the solution of the diffusion theory problem. Once we know the transformation rule of the flux, for example, it suffices to calculate the flux in a part of V and exploit the transformation rules. That observation is formulated in the following concise way. Let r ∈ V a point in V and let g · r be the image of r under g ∈ G. Then the set of points g · r, g ∈ G is called the orbit of r under the group G. If there is a set V 0 ∈ V such that the orbits of r 0 ∈ V 0 give every point 2 of V we call V 0 the fundamental domain of V. It is thus sufficient to solve the problem on the fundamental domain V 0 , and "continue" the solution to the whole volume V.
When the boundary condition is not homogeneous or there is an external source, we exploit the linearity of the diffusion equation. The general solution is the sum of two terms: one with external source but homogeneous boundary condition and one with no external source but with non-homogeneous boundary condition. In either case, it is the external term that determines the transformation properties of the respective solution component.

Selection of basis functions
The purely geometric symmetries of a suitable equation lead to a decomposition (2.16) of an arbitrary function in a function space, and thus the decomposition of the function space itself. The decomposed elements are linearly independent and can be arranged to form an orthonormal system. This can be exploited in the calculations.
In a homogeneous material one can readily construct trial functions that fulfill the diffusion equation at each point of V. For example consider The general solution to (4.12) takes the form of where the weight functions W k (e) are arbitrary suitable functions, i 2 = −1, and t k signify the eigenvectors of matrix A: At k = λ 2 k t k . (4.14) When e(θ)=( cos θ, sin θ), using (2.19), we build up a regular representation from (4.13) so that ψ 0 (r)= G ∑ k=1 t k 2π/|G| 0 e iλ k r·e(θ) W k (e(θ))dθ, (4.15) and the action of operators O g on ψ 0 (r) is defined as follows. O g acts on variable r, see (2.6), but in (4.13), r occurs only in the form of the dot product re(θ), therefore action of O g can be transferred to an action on θ. As a result, each O g acts as where O g maps the interval 0 ≤ θ ≤ 2π/|G| into the interval I g . In this manner we get the irreducible components of the solution as a linear combination of |G| exponential function, it is only the coefficients in the linear combination that determine the irreducible components. The weight function W k (θ) makes it possible to match the entering currents at given points of the boundary. Let θ = 0 correspond to the middle of a side. Then choosing we get by (4.15) the solution at face midpoints. The last step is the formation of the irreducible components. Observe that in projection (2.20) the solutions at different images of r are used in a linear combination, the coefficients of the linear combinations are the rows of the character table. But in the images (4.16), only the weight function changes. In each I g interval the image of W k (e) is involved, which is a Dirac-delta function, only the place of the singularity changes as the group elements map the place of the singularity. A symmetry of the square maps a face center into another face center thus there will be four distinct positions and the space dependent part of the irreducible component of ψ 0 will contain four exponentials: ± e iλ k x , ±e −iλ k x , ±e iλ k y , ±e −iλ k y . (4.18) From these expressions the following irreducible combinations can be formed: A 1 : cos λ k x + cos λ k y; A 2 : cos λ k x − cos λ k y; E 1 : sin λ k x; E 2 : sin λ k y. It is not surprising that when we represent a side by its midpoint the odd functions along the side are missing.
The above method may serve as a starting point for developing efficient numerical methods. The only approximation is in the continuity of the partial currents at the boundary of adjacent homogeneous nodes.
If elements of the function space are defined for all r ∈ V, and if f 1 , f 2 ∈ L V , then the following inner product is applicable: (4.20) Let f α ℓ (r), ℓ = 1,...,n α be a regular representation of group G. Then furthermore, for the reactions rates formed with the help of the cross-sections in (4.1), similar orthogonality relation holds. For the volume integrated reaction rates we have (4.22) in other words: solely the most symmetric, one dimensional representation contributes to the volume integrated reaction rates. Note that as a result of the decomposition of the solution or its approximation into irreducible components not only that irreducible components of a given physical quantity (like flux, reaction rate, net current) but also the given irreducible component of every physical quantity fall into the same linearly independent irreducible subspace. As a consequence, the operators (matrices) mapping the flux into net currents (or vice versa) fall into the same irreducible subspace, therefore the mapping matrix automatically becomes diagonal. subspaces α i , i < 5 are one-dimensional whereas the subspace α = 5 is two-dimensional, and in a two-dimensional representation there are two pairs of basis functions that are identical as to symmetry properties. Thus, we have altogether eight linearly independent basis functions.
The physical meaning of the irreducible components is that the flux distribution of a node is a combination of the flux distributions established by eight boundary condition types. The component A 1 represents a complete symmetry that is the same even distribution along each side. Component A 2 is also symmetric, but the boundary condition is an odd function on each side. Components B 1 and B 2 represent entering neutrons along one axes and exiting neutrons along the perpendicular axes, a realization of a second derivative with even functions over a face. B 2 is the same but with odd functions along a face. E 1 and E 4 represent streaming in the x and y directions with even distributions along a face, whereas E 2 and E 3 with odd distributions along a face.
The symmetry transformations of the square, map the functions given along the half faces into each other but they do not say anything about the function shape along a half face. Therefore, the functions in Fig. 2 serve only as patterns, the function shape is arbitrary along a half face. The corresponding mathematical term is the direct product; each function may be multiplied by a function f (ξ), −h/2 ≤ ξ ≤ +h/2. It is well known that a function along an interval can be approximated by a suitable polynomial (Weierstrass's theorem). We know from practice that in reactor calculations a second order polynomial suffices on a face for the precision needed in a power plant.
The invariant subspace means that the boundary flux, the net current, the partial currents must follow one of the patterns shown in Figure 2, the only difference may be in the shape function f (ξ), −h/2 ≤ ξ ≤ +h/2. This means the a constant flux may create a quadratic position dependent current, but the global structure of the flux, and current should belong to the same pattern of Figure 2.
Moreover, if we are interested in the solution inside the square, its pattern must also be the same although there the freedom allows a continuous function along 1/8-th of the square. These features are exploited in the calculation.

Iteration
It is known that the diffusion (as well as the transport) equation has a well defined solution in V provided the entering current is given along the boundary ∂V. From the Green's function and from the operators in (4.11) we set up the following iteration scheme. To formalize this, we write the solution as (4.23) Applying operator F that forms the exiting current from the flux, we obtain that can be put into the concise form where we have suppressed that the partial currents depend on position along the boundary and the response matrix R includes an integration over variable r ′ .
When volume V is large, we subdivide it into subvolumes (nodes) and determine the response matrices for each subvolume. At internal boundaries, the exiting current is the incident current of the adjacent subvolume. Thus in a composite volume the partial currents are connected by response matrices and adjacency. We collect the response matrices and adjacency into two big response matrices: (4.26) and because the adjacency is an invertible relationship, we multiply the first expression by H and get I = HRI. (4.27) Since there is a free parameter k ef f in matrix R, it makes the equation solvable. At external boundaries there is no adjacency, but the boundary condition there provides a rule to determine the entering current from the exiting current. With these supplements, the solution of equation (4.27) proceeds I (m+1) = HR (m) I (m) . (4.28) The iteration starts with m = 0 with an initial guess for the k ef f and the entering currents I. Let us assume that the needed matrices are available, their determinations are discussed in the subsequent Subsection. The iteration proceeds as follows. We sweep through the subvolumes in a given sequence and carry out the following actions (in node m): • collect the actual incoming currents of subvolume m.
• determine the actual response matrix to calculate the new exiting currents and contributions to volume integrals 3 . • determine the new exiting currents (J) from the entering currents and the response matrices using equation (4.26) and the contributions to the volume integrals.
After this, pass on to the next node. When the iteration reaches the last node, the sweep ends and the maximal difference is determined between the entering currents of the last two iterations. At the end of an iteration step, the parameter k ef f is re-evaluated from the condition that the largest eigenvalue of HR should equal one. If the difference of the last two estimates is greater than the given tolerance limit, a new iteration cycle is started, otherwise the iteration terminates. If we have a large number of nodes, the improvement after the calculations of a given node is small. This shows that the iteration process is rather slow, acceleration methods are required.
It has been proven (Mika, 1972) that the outlined iteration is convergent. The goal of the iteration is to determine the partial current vector. The length of vector I is N node × n F × G. From the point of view of mathematics, the iteration is a transformation of the following type: (4.29) where m is the number of the iteration, matrix A(k ef f ) makes the new entering current vector x m+1 from the old entering current vector x m . In the case of neutron diffusion or transport, operator A(k ef f ) maps positive vectors into positive vectors. In accordance with the Krein-Ruthman theorem, A(k ef f ) has a dominant eigenvalue and the associated eigenfunction 4 . When k ef f is a given value, the power method is a simple iteration technique to find a good estimate of x = lim i→∞ x i . Solution methods have been worked out for practical problems in nuclear reactor theory: for the solution of the diffusion and transport equations in the core of a power reactor. The original numerical method is described elsewhere, see Refs. (Weiss, 1977), (Hegedus, 1991).
Note that the iteration (4.29) is just an example of the maps transforming an element of the solution space into another element. Thus in principle one can observe chaotic behavior, divergence, strange attractors 5 etc. Therefore it is especially important to design carefully the iteration scheme. The iteration includes derived quantities of two types: volume integrated and surface integrated. When you work with an analytical solution, the two are derived from the same analytical solution. But when you are using approximations (such as polynomial approximation), it has to be checked if the polynomials used inside the node and at the surface of the node are consistent. In an eigenvalue problem, parameter k ef f in equation (4.29)should be determined from the condition that the dominant eigenvalue a in (4.29) should equal one. First we deal with the general features of the iteration.
As has been mentioned, one iteration step (4.29) sweeps through all the subvolumes. The number of subvolumes (Gadó et al., 1994) varies between 590 and 7980, the number of unknowns is 9440 and 111680. At the boundary of two adjacent subvolumes, continuity of Φ and D∂ n Φ (the normal current) is prescribed .
In node m in iteration i. In the derivation of the analytical solution we have assumed the node to be invariant under the group G V . Actually, not the material properties are stored in a program because the material properties depend on: • actual temperature of the node; • the initial composition of the fuel (e.g. enrichment); • the actual composition of the fuel as it may change with burn-up; • the void content of the moderator; • the power level; In the calculations, app. 50-60% of the time is spent on finding the actual response matrix elements, because those depend on a number of local material parameters (e.g. density, temperature, void content). We mention this datum to underline how important it is to reduce the parametrization work in a production code.

Exploiting symmetries
In a given node, the response matrices are determined based on the analytical solution (4.19). We need an efficient recipe for decomposing the entering currents into irreps and reconstructing the exiting currents on the faces. Since the only approximation in the procedure requires the continuity of the partial currents, we need to specify the representation of the partial currents and how to represent them. The simplest is a representation by discrete points along the boundary, the minimal number is four, the maximal number depends on the computer capacity. An alternative choice is to represent the partial currents by moments over the faces. Usually average, first and second moment suffice to get the accuracy needed by practice. The representation fixes the number of points we need on a side and the number of points (n) on the node boundary.
To project the irreps, we may use (Mackey, 1980) the cos((k − 1)2π/n), k = 1, . . . , n/2 and sin((k)2π/n), k = 1, . . . , n/2 vectors (after normalization). The following illustration shows the case with n = 4, i.e. one value per face. In a square node we need the following matrix (4.30) to project the irreducible components from the side-wise values. As (2.20) shows, irreducible components are linear combinations of the decomposable quantity 6 . The coefficients are given as rows in matrices Ω 4 .
In a regular n-gonal node the response matrix has 7 Ent [(n + 2)/2] free parameters. The response matrix also has to be decomposed into irreps, this is done by a basis change. Let the response matrix give J = RI Multiply this expression by Ω from the left: and we see that for irreducible representations the response matrix is given by ΩRΩ + .I na square node: (4.32) and the irreducible representation of R 4 is diagonal: where We summarize the following advantages of applying group theory: • Irreducible components of various items play a central role in the method. The irreducible representations often have a physical meaning and make the calculations more effective (e.g. matrices transforming one irreducible component into another are diagonal). • The irreducible representations of a given quantity are linearly independent and that is exploited in the analysis of convergence. • The usage of linearly independent irreducible components is rather useful in the analysis of the iteration of a numerical process. • In several problems of practical importance, the problem is almost symmetric, some perturbations occur. This makes the calculation more effective.
• It is more efficient to break up a problem into parts and solve each subproblem independently. Results have been reported for operational codes (Gadó et al., 1994).
The above considerations dealt with the local symmetries. However, if we decompose the partial currents into irreps, we get a decomposition of the global vector x in equation (4.28) as well. We exploit the linear independence of the irreducible components further on the global scale.
For most physical problems we have a priori knowledge about the solution to a given boundary value problem in the form of smoothness and boundedness. This is brought to bear through the choice of solution space. In the following, we introduce via group theoretical principles the additional information of the particular geometric symmetry of the node. This allows the decomposition of the solution space into irreducible subspaces, and leads, for a given geometry, not only to a rule for choosing the optimum combination of polynomial expansions on the surface and in the volume, but also elucidates the subtle effect that the geometry of the physical system can have on the algorithm for the solution of the associated mathematical boundary problem.
Consider the iteration (4.29) and decompose the iterated vector into irreducible component where because of the orthogonality of the irreducible components x β+ x α = 0 when α = β. The convergence of the iteration means that lim N→∞ x N+k1 − x N+k2 = 0 (4.35) for any k 1 , k 2 . But that entails that as the iteration proceeds, the difference between two iterated vector must tend to zero. In other words, the iteration must converge in every irreducible subspace. This observation may be violated when the iteration process has not been carefully designed.
Let us assume a method, see (Palmiotti, 1995), in which N basis functions are used to expand the solution along the boundary of a node and M basis functions to expand the solution inside the node. It is reasonable to use the approximation of same order along each face, hence, in a square node N is a multiple of four. For an Mth order approximation inside the node, the number of free coefficients is (M + 1)(M + 2)/2. It has been shown that an algorithm (Palmiotti, 1995) with a linear (N = 1) approximation along the four faces, with 8 free coefficients, of the boundary did not result in convergent algorithm unless M = 4 quartic polynomial, with 15 free coefficients, was used inside the node.
In such a code each node is considered to be homogeneous in composition. Central to the accuracy of the method are two approximations. In the first, we assume the solution on the boundary surface of the node to be expanded in a set of basis functions ( f i (ξ); i = 1,...,N). In the second, the solution inside the volume is expanded in another set of basis functions (F j (r); j = 1,...,M). Clearly the independent variable ξ is a limit of the independent variable r.
Any iteration procedure, in principle, connects neighboring nodes through continuity and smoothness conditions. For an efficient numerical algorithm it is therefore desirable to have i/Order 0 1 2 3 4 1 1 -(x 2 + y 2 ) -(x 2 y 2 ), ( Table 6. Irreducible components of at most fourth order polynomials under the symmetries of a square C 4v the same number of degrees of freedom (i.e. coefficients in the expansion) on the surface of the node as within the node. With the help of Table 6, for the case of a square node, we compare the required number of coefficients for different orders of polynomial expansion. A linear approximation along the four faces of the square has at least one component in each irreducible subspace. At the same time the first polynomial contributing to the second irrep is fourth order. Convergence requires the convergence in each subspace thus the approximation inside the square must be at least of fourth order. There is no linear polynomial approximation that would use the same number of coefficients on the surface as inside the volume. The appropriate choice of order of expansion is thus not straightforward but it is important to the accuracy of the solution, because a mismatch of degrees of freedom inside and on the surface of the node is likely to lead to a loss of information in the computational step that passes from one node to the next. A lack of convergence has been observed, see (Palmiotti, 1995), in the case of calculations with a square node when using first order polynomials on the surface. A convergent solution is obtained only with fourth or higher order polynomial interpolation inside the node. Similar relationships apply to nodes of other geometry. For a hexagonal node that there is no polynomial where the number of coefficients on the surface matches the number of coefficients inside the node.
In a hexagonal node in (Palmiotti, 1995), the first convergent solution with a linear approximation on the surface requires at least a sixth order polynomial expansion within the node. Thus, in the case of a linear approximation on the surface, in the case of a square node a third order polynomial within the node does not lead to a convergent solution, although the number of coefficients is greater than those on the surface. In the case of the regular hexagonal node, a convergent solution is obtained only for the sixth order polynomial expansion in the node, while both a fourth and a fifth order polynomial have a greater number of coefficients inside the node than on the surface. It appears that some terms of the polynomial expansion contain less information than others, and are thus superfluous in the computational algorithm. If these terms can be "filtered out", a more efficient and convergent solution should result. The explanation becomes immediately clear from the decomposition of the trial functions inside the volume and on the boundary. In both the square and hexagon nodes, the first order approximation on the boundary is sufficient to furnish all irreducible subspaces whereas this is true for the interpolating polynomials inside V for surprisingly high order polynomials.

Reactor physics
In analogy with the application of group theory in particle physics, where group theory leads to insights into the relationships between elementary particles, we present an application of group theory to the solution to a specific reactor physics problem. The question is whether it is possible to replace a part of a heterogeneous core by a homogeneous material so that the solution outside the homogeneous region remains the same? This old problem is known as homogenization (Selengut, 1960).
In particular, for non-uniform lattices, asymptotic theory has shown that a lattice composed of identical cells has a solution that is composed of a periodic microflux and a slowly varying macroflux. What happens if the cell geometry is the same but the material composition varies?
In reactor calculations, we solve an equation derived from neutron balance. In that equation, we encounter reaction rates, currents or partial currents. It is reasonable to derive all the quantities we need from one given basic quantity, say from the neutron flux at given points of the boundary. The archetype of such relation is the exiting current determined from the entering current by a response matrix. We show that by using irreducible components of the partial currents, the response matrix becomes diagonal.
The Selengut principle is formulated: if the response matrix of a given heterogeneous material in V can be substituted by the response matrix of a homogeneous material in V, there exist an equivalent homogeneous material with which one may replace V. This principle simplifies calculations considerably, and, therefore, has been widely used in reactor physics. We investigate the Selengut principle more closely (Makai, 2010), (Makai, 1992).
The analysis is based on the analytical solution of the diffusion equation derived in the previous Section. The problem is considered in a few energy groups, the boundary flux F is a vector, as well as the volume averaged fluxΦ. Using that solution, we are able to derive matrices mapping into each other the volume integrated fluxes, the surface integrated partial and net currents. The derivation of the corresponding matrices is as follows. Our basis is the boundary flux, that we derive for each irrep i from (4.13). The expression (4.13) has three components. The first one is vector t k which is independent of the position r and is multiplied by an exponential function with λ k r in the exponent. The third component is the weight W k which is independent of r but varies with subscript k. The product is summed for subscript k, that labels the eigenvalues of the cross-section matrix in (4.14). That expression can be put into the following concise form: where c i comprises the third component. Here < f i (r) > is a diagonal matrix. Note that position dependent quantities like reaction rates, follow that structure. The normal component of the net current is J net obtained from the flux by taking the derivative and is given in irrep i as J net,i = −DT < g i > c i .
We eliminate c i to get Here n is the outward normal to face F i , The volume integrated fluxΦ is obtained after integration from (4.13) as Φ = T < F A1 > c A1 ,F A1 = f A1 (r)d 3 r, (5.5) and the integration runs over volume V of the node. Note that only irrep A1 (i.e. complete symmetry) contributes to the average flux because of the orthogonality of the irreducible flux components. After eliminating c A1 from (5.1), we get the response matrix for determining the volume integrated fluxΦ from the face integrated flux F A1 : (5.6) This assures that V is completely described by matrix W and the diagonal matrices < F(r) > , < f (r) >, < g(r) > for each irrep. For example, we are able to reconstruct the cross-section matrix Σ from them. Note that WT = T < F/ f >, the eigenvectors of matrix W are the eigenvectors of A. Now we need only a numerical procedure to find the eigenvalues λ k from < F/ f >.
The question is, under what conditions are the above calculations feasible. We count the number of response matrices. The matrix elements we need to characterize V may be all different and the number of matrices depends on the shape of V, since the number of irreducible components of the involved matrices depends on the geometry. In a square shaped homogeneous V, we have four R i matrices and one W. Altogether we have to determine 5 * G * G elements. In an inhomogeneous hexagonal volume, there are 6 * G * G matrix elements, whereas the homogeneous material is described by G * (G + 1) parameters as in a homogeneous material there are altogether G * G cross-sections and G diffusion constants 8 . Therefore the Selengut principle is not exact it may only be a good approximation under specific circumstances. Homogenization recipes preserve only specific reaction rates, but they do not provide general equivalence.

Conclusions
The basic elements of the theory of finite symmetry groups has been introduced. In particular, the use of the machinery associated with the decomposition into irreducible representations, in analogy with harmonic analysis of functions in function space, in the analysis of Nuclear Engineering problems. The physical settings of many Nuclear Engineering problems exhibit symmetry, as for example in the solution of the multi-group neutron diffusion equation. This symmetry can be systematically exploited via group theory, and elicit information that leads to more efficient numerical algorithms and also to useful insights. This is a result due to the added information inherent in symmetry, and the ability of group theory to define the "rules" of the symmetry and allows one to exploit them.