Efficient reduced order model based on the proper orthogonal decomposition for time-dependent MOC calculations

ABSTRACT An efficient reduced order model (ROM) for time-dependent transport calculations using the method of characteristics (MOC) is proposed. In the present ROM, the flux distributions and the net neutron currents between the adjacent unstructured mesh regions are taken from the MOC solution. Then, the coefficient matrices for the MOC-equivalent diffusion equation are reconstructed from them. The proper orthogonal decomposition (POD) is applied for the MOC-equivalent diffusion equation to reduce the degree of freedoms (DOFs) using the orthogonal bases obtained by the singular value decomposition (SVD) for the sampled MOC solution. The accuracy and computation time of the present ROM are verified in the C5G7-TD 2D benchmark problem. The calculation results show that the present ROM enables us approximately 5000–6000 times faster computation than the full order model (FOM) for kinetic calculations itself in the present calculation condition. The present method can be substituted as real-time simulations without the spatial homogenization when typical flux distributions and the net neutron currents of a target problem can be precalculated.


Introduction
The pin-resolved transport calculation in large heterogeneous geometry is one of the state-ofthe-art simulations in reactor analysis. The method of characteristics (MOC) [1] is widely used for such simulations because of its accuracy and fidelity. In the field of reactor kinetics, various efficient numerical algorithms are developed and are successful to achieve reasonable computation time for the time-dependent MOC calculations [2][3][4]. However, it is still difficult to apply such simulations to real-time applications due to a large number of degrees of freedom (DOF) in the spatial domain. Thus, lower-order efficient numerical models that reproduce the original solution of the full order models (FOMs) are desirable. Dimensionality reduction (DR) techniques can address this issue and the proper orthogonal basis decomposition (POD) [5][6][7][8] is one of the most successful methods to construct a reduced order model (ROM).
The basic concept of the POD is to expand the flux distributions in the target problem with a suitable orthogonal basis and to convert the neutron balance equation into the equation for the expansion coefficients. The DOFs to represent flux distributions in the target problem are dramatically reduced through the expansion. Such a basis is obtained by the singular value decomposition (SVD) for the flux vectors taken from the FOM solution of the target problem. In past studies, the POD has been successfully applied to diffusion and S N transport calculations and verified in simplified or homogeneous geometries [6][7][8]. However, there is a constraint because the neutron balance equation must be described in a matrix form to apply the POD. Since the neutron balance equation in the MOC is not generally described in a matrix form, it is difficult to directly apply the POD for MOC calculations.
Thus, in the present study, we aim to construct a ROM that reproduces the time-dependent MOC solution in a heterogeneous geometry using the POD. In the present study, the coefficient matrices which reproduce a MOC solution are reconstructed using diffusion calculation with a correction term for neutron current, and then the POD is applied to it. The accuracy and the computation time of the present ROM are verified with the C5G7-TD 2D benchmark problem [9].
The remaining part of this paper is organized as follows. The theoretical basis of the ROM for the time-dependent MOC calculations using the POD is described in section 2. The numerical results and discussions are shown in section 3. Finally, concluding remarks are summarized in section 4.

Dimensionality reduction for time-dependent MOC using POD
In the present study, the coefficient matrices that reproduce the MOC solution are reconstructed based on the diffusion calculation with a correction term. Then, the MOC-equivalent discretized diffusion equation is dimensionally reduced using the POD. In this section, the theoretical bases for the MOC-equivalent diffusion equation and the DR using the POD are described in sections 2.1 and 2.2, respectively.

Reconstruction of coefficient matrices for MOC
The time-dependent transport equation along a neutron flight path is written as: 1 v g;r @ψ g;m;k;i @t ¼ À @ @s þ Σ t;g;r � � ψ g;m;k;i þ 1 4π Q g;r ; (1) Q g;r ¼ X j λ j χ d g;j;r C j;r þ X g 0 Σ s;g 0 !g;r ϕ g 0 ;r þ 1 À β r À � χ p g;r X g 0 νΣ f ;g 0 ;r ϕ g 0 ;r ; (2) @C j;r @t ¼ β j;r X g 0 νΣ f ;g 0 ϕ g 0 ;r À λ j C j;r ; where s: coordinate along the neutron flight direction, v g;r : averaged neutron velocity, ϕ g;r : scalar flux, ψ g;m;k;i : angular flux, Σ t;g;r : total cross section, νΣ f ;g 0 : production cross section, Σ s;g 0 !g : scattering cross section, χ p g;r : prompt fission spectrum, χ d g;j;r : delayed fission spectrum, C j;r : precursor density, β j;r : delayed neutron fraction, λ j : decay constant, where, g, j,m, k, i and r are the subscripts for energy group, delayed precursor family, neutron flight direction, sequential number of the neutron flight path, segment, and flux region, respectively. Figure 1 shows their definitions.
The leakage term in Equation (1) is the matrix-free operator, which makes it difficult to directly apply the POD to the MOC. Integrating both sides of Equation (1) for the whole solid angle and spatially averaging it within each flux region, where J MOC g;r!r 0 : net neutron current between the flux regions, S r!r 0 : surface area between the adjacent flux regions, V r : volume of the flux region, The net neutron current, J MOC g;r!r 0 , is calculated as a summation of the incoming and outgoing angular fluxes between the adjacent regions as: ΔJ out g;m;k;i ¼ ΔA m;k ω m ψ out g;m;k;i sin θ m ; ΔJ in g;m;k;i ¼ ΔA m;k ω m ψ in g;m;k;i sin θ m ; where ΔA m;k : width of the neutron flight path, θ m : polar angle, ω m : solid angle weight, ψ out g;m;k;i : outgoing angular flux from segment i, ψ in g;m;k;i : incoming angular flux for segment i. Note that Equations (6)-(8) can be calculated using the ray-tracing information and the present approach can be applied to any geometries in which the neutron flight paths can be defined. By describing the net neutron current with the finite difference method with a correction term as it is done in the coarse mesh finite difference (CMFD) method as [10]: where D FD g;r!r 0 is a coupling coefficient and D COR g;r!r 0 is the correction factor to reproduce the net neutron current, J MOC g;r!r 0 , between the adjacent flux regions, r and r 0 . In the present study, the coupling coefficient is calculated as: where D g;r and S r are the diffusion constant and the cross-sectional area of the flux region in the radial direction, respectively. Note that S r in Equation (10) is different from the definition in the ordinary CMFD method to simplify the definition of the coupling coefficient for the unstructured mesh but the approximation error for it is corrected by D COR g;r!r 0 . As for the boundary, the coefficient for neutron current is described as where J MOC g;r!BC and S r!BC are the net neutron current and surface area at the boundary, respectively. The coefficient for the neutron current, D BC g;r , is calculated so that it preserves the relationship in Equation (11). Substituting Equations (9), (10) and (11) into Equation (5), the MOC-equivalent neutron balance equation for each flux region is derived as: where A g;r 0 !r is the annihilation operator described as: A g;r 0 !r ¼ D COR g;r!r 0 À D FD g;r!r 0 r�r 0 ð Þ: where D BC g;r ¼ 0 for the flux region that is not facing the boundary.
Employing the fully implicit method for time integration [11], Equation (12) is described as: where n is the superscript for a time step and δ r 0 r is the Kronecker delta. μ j , η j , and � j are the coefficients for temporal integration of precursors, which are derived assuming the linear transition of the fission source between the successive time steps as [11]: Finally, the neutron balance equation for each flux region, which reproduces the neutron balance calculated with the MOC, is written as a matrix form as: where I is the identity matrix and the other matrices and vectors are described as follows: where R is the number of the flux regions in the calculation geometry. The size of the identity matrix in Equation (26) is also R � R.

Dimensionality reduction using POD
Once the MOC-equivalent diffusion equation is derived in a matrix form, its DOFs can be reduced by the POD. In the POD approach, the flux and precursor density distributions are expanded by orthogonal bases suitable for the target problem, which are numerically obtained with the SVD for the distributions taken from the detailed FOM solutions. By substituting them for the neutron balance equation, the equation for the expansion coefficients, which represents the target problem with fewer DOFs is obtained. In the previous study, an orthogonal basis that represents the flux distributions for all energy groups was employed to apply the POD [12]. However, in the present study, the orthogonal bases are respectively constructed for the scalar flux and precursor distributions at each energy group/precursor family as proposed in [7] to obtain more suitable bases that represent the parameter distributions with fewer DOFs. In this section, the equations for the expansion coefficients are derived by applying the POD to the MOC-equivalent diffusion equation. The solution space of the scalar flux distributions, R F;g , is defined as: where N is the number of the scalar flux distributions taken from the FOM solution and N < R is assumed in the following derivation. Applying the SVD for R F;g , an orthogonal basis for the target problem is obtained.
In the present study, the thin SVD [13] is employed to reduce the computational cost of the decomposition as: where � and K denote the singular value and the rank of the solution space R F;g , respectively. The matrix U F contains the first K left singular vectors of group g and its size is R � K. V F consists of the first K right singular vectors of group g and its size is N � K. They are orthogonal bases for the column or row spaces of the solution space R F;g so that where I in Equations (42) and (43) is K � K identity matrix. Since the scalar flux distributions are represented in column space of R F;g , the flux distribution is expanded with U F;g as: where φ nþ1 g is the column vector of the expansion coefficients. Substituting Equation (44) into Equation (24) and multiplying the transposed orthogonal basis, the neutron balance equation for each flux region is converted into the equation for the expansion coefficients as: Note that the size of the coefficient matrices on the left side of Equations (47) and (48) is reduced from R � R to K � K. The source term, Q nþ1 g , is described as: where the left side of Equations (52) and (53) are also K � K matrices. In Equation (52), the production cross sections are compressed together with the fission spectrums. However, χ d g;j and β j νΣ f ;g must be compressed separately to transform Equations (27) and (28) with the POD. Thus, in the present study, the matrices for the production cross section and the fission spectrum for the prompt fission are also separately compressed for the symmetry of the derivation. Then, Equation (52) is redefined from them. They are derived in the following procedure: (1) The prompt and delayed fission rate distribution vectors, P nþ1 PF and P nþ1 DF;j , are employed to describe the fission source distributions (Equations (54) and (55)).
(2) Employ the thin SVD for the solution spaces for the prompt and delayed fission rate distributions and obtain the orthogonal bases for them, U PF and U DF;j , respectively. (Equations (57)-(60)) (3) Expand the prompt and delayed fission rate distributions, P nþ1 PF and P nþ1 DF;j , with the orthogonal bases, U PF and U DF;j . (Equations (61) and (62) The vectors for the prompt and delayed fission rate distributions, P nþ1 PF and P nþ1 DF;j , are defined as: Using P nþ1 PF and P nþ1 DF;j , the fission source term is described as: Note that the matrices for the fission spectrum can be separately compressed with the production cross section if the fission rate distributions are expanded with their orthogonal bases. Thus, the thin SVD is also employed to obtain the orthogonal bases for the solution spaces for the prompt and delayed fission rate distributions, R PF and R DF;j , as: where the matrices U PF , S PF , V PF , U DF;j , S DF;j , and V DF;j denote the R � K left singular vectors, K � K diagonal singular values and the N � K right singular vectors for the solution space of the prompt and delayed fission rate distributions, respectively. The left and right singular vectors contain the first K vectors for each solution space. The rank of the solution space of the fission rate distributions is assumed as the same as that of flux distributions to simplify the notation, where it is not a necessary constraint. Using the orthogonal bases, the prompt and delayed fission rate distributions are expanded as: where and φ nþ1 PF and φ nþ1 DF;j are the expansion coefficients for the prompt and delayed fission rate distributions, respectively. Substituting Equations (61) and (62) into Equation (56), the fission source is described as: where b χ p g and b χ d g;j are the compressed matrices for the prompt and delayed fission spectrum described as: Substituting Equations (54) and (55) into Equations (61) and (62) respectively, the matrices for the production cross section are also compressed as: where b Σ nþ1 PF;g 0 and b Σ nþ1 DF;g 0 ;j are the compressed matrices for the prompt and delayed production cross sections described as: Note that both the production cross sections and fission spectrum for the prompt and delayed fission are described as K � K matrices for each energy group. The coefficient matrix for the fission source can be reconstructed from them by substituting Equations (66) and (67) into Equation (63) as: As for the pseudo-source term, substituting Equation (27) into Equation (50) and using Equations (65) and (69), Since the DOFs of the precursor density distribution, C n j , is still that of the FOM, assuming the precursor density can be expanded with the orthogonal basis for the delayed fission rate distribution, U DF;j , as: where ς n j denotes the expansion coefficients for the precursor density distribution. Substituting Equations (72) into Equation (71) and using Equation (65), the pseudo source term is described as: Finally, multiplying the transposed orthogonal basis for delayed fission rate distribution, U T DF;j , for Equation (28) and using Equations (44), (69) and (72) for it, the temporal integration for the expansion coefficient of the precursor density is described as: Note that Equations (46), (51), (73), and (74) are the equations for expansion coefficients and they can be solved with fewer DOFs than the balance equation in the FOM. However, since the coefficient matrices before the DR are taken from the FOM solution, the compressed coefficient matrices are only available at the sampled time steps. Thus, the compressed coefficient matrices are interpolated between the successive sampled time steps. The accuracy of the interpolation is verified in section 3.1.

Numerical results and discussion
In the present study, the C5G7-TD 2D benchmark problem is used to verify the accuracy and performance of the present ROM. Figures S1 and S2 in the supplemental online material accompanying this paper show the top view of the benchmark geometry. The cross sections for 7 energy groups and the delayed neutron parameters for 8 precursor families for each material specified in the benchmark problem are used in the present verification. The following 3 exercises provided for the 2D benchmark calculation are used in the present verification. Each exercise includes multiple test problems with different perturbation conditions, i.e. TD1-1 -TD3-4. The detail of them is summarized in Table S1, Figures S3 and S4 in the supplemental online material accompanying this paper. In all test problems, the cross sections linearly varies for t = 0.0 sec -1.0 sec and 1.0 sec -2.0 sec, respectively. The present verifications consist of 2 parts. Section 3.1 shows verification for the time dependence on the compressed coefficient matrices. Section 3.2 shows the verification of the computational efficiency of the ROM construction. In the present verification, all FOM calculations are carried out with the in-house 2D MOC code which has been verified in a previous study [4].

Verification for time dependence on coefficient matrices
In this section, the time dependence on the compressed coefficient matrices is verified. Table 1 and Figure 2 show the calculation condition of the reference FOM solution and the spatial mesh structure, respectively. The reference solutions are calculated with the same calculation conditions for each test problem in the present verification.
The multi-grid amplitude function (MAF) method shown in Table 1 is one of the factorization methods which treats the spatial dependence of the amplitude function to reduce the temporal discretization error [4,15]. The linear source approximation [16,17] is also employed to reduce the spatial discretization error in the FOM calculation. Figure 3 shows the core power transitions of the reference solutions.
Since the C5G7-TD benchmark problem is the transient calculations for 10 sec and the time step size for the shape function is 500 msec in the present calculation condition, N ¼ 21 flux distributions and the coefficient matrices are taken from the reference solutions for each test problem, i.e. TD1-1, TD1-2, etc. Thus, the solution space for each test problem is described as a 31,501 � 21 matrix, where 31,501 is the number of the flux regions in the reference solution. In the present verification, the thin SVD is applied for each solution space and the orthogonal bases are independently constructed, where the size of the orthogonal bases is also 31,501 � 21 for each test problem. Thus, note that the present orthogonal   bases focus to represent the transitions of the parameter distributions in each test problem with minimum DOFs. Figure 4 shows the singular values of solution space for the flux, prompt and delayed fission rate distributions in the test problem, TD3-4. As shown in Figure 4, the singular values, which indicate the contribution of each column vector in the orthogonal basis to the solution space, are rapidly decreasing. Since the column vectors corresponding to small singular values have less impact on the representation of the solution space, it can be seen that the solution space is accurately represented by a small number of dimensions. The major trend of the singular values is almost the same in the other test problems.
In the present verification, 3 ROMs are constructed based on different calculation conditions. Table 2 shows the case matrix for the ROM construction. Note that the low-rank approximation for the orthogonal bases is not employed in the present verification, i.e. the size of the orthogonal bases used for the DR is also 31,501 � 21 and that of compressed coefficient matrices b A nþ1 g is 21 � 21, respectively. Since the size of the coefficient matrices before the DR with the POD, A nþ1 g , is 31,501 � 31501, the DOF the target problem is dramatically reduced with the POD.
As shown in Table 2, case 1 utilizes the coefficient matrices generated by 1.0 sec interval samplings and the matrices are linearly interpolated within the sampled time steps. Thus, it includes a larger interpolation error for the compressed coefficient matrices than the other cases. Case 2 is constructed with the finer sampling intervals for the coefficient matrices, 0.5 sec. However, if the coefficients in the compressed matrices vary non-linearly in time, the interpolation error for the coefficient matrices degrade the accuracy of the present ROM. Case 3 is constructed with the same sampling intervals as case 2 but the coefficient matrices are quadratically interpolated. Figure 5 shows the typical transition of a coefficient in the compressed matrix b A nþ1 g for the test problem, TD3-4.    As shown in Figure 5, the coefficient in the compressed matrix non-linearly varies during perturbation, where the discontinuity of the coefficient at t = 1.0 sec comes from the perturbation condition that the moderator density reaches a minimum at that point. Since the cross sections linearly vary in the C5G7-TD benchmark problem, this non-linearity comes from the correction factor for neutron current, D COR g;r!r 0 . In case 3, the interpolation intervals are selected as follows to avoid the overshoot of the interpolation: The time step size for the present ROMs is 10 msec for all test cases, which is the same as that of the amplitude function in the reference solution. Figures 6 and 7 show the relative difference for the core power and the maximum relative difference for the flux distributions with respect to the reference FOM solution, where the relative differences are calculated as: max g2Gφ n;ROM g Àφ n;FOM g ϕ n;FOM where CP denotes the core power. The superscripts, FOM and ROM, denote the FOM and the ROM solutions, respectively. As shown in Figures 6 and 7, all cases in TD1 are in good agreement with the reference solution since the reactivity insertion and the variations of the coefficient matrices are small. However, there is a noticeable trend for TD2 and TD3 among the test cases especially from t = 0 sec to 2 sec. These results show the impact of the non-linearity of the coefficient matrices. Since all cross sections are linearly perturbed in the exercises TD1-TD3 and the spatial homogenization across the material boundary is not employed in the present study, the nonlinearity appears only in the correction term for the neutron current in the present calculation condition. Thus, case 3 which employs non-linear interpolation for the coefficient matrices shows a better agreement with the reference FOM solution. These results suggest the effectiveness of the non-linear interpolation for the compressed coefficient matrices.

Verification for computational efficiency of ROM construction
In the previous section, the present ROM is constructed from the reference solution to see the impact on the temporal discretization error for the coefficient matrices. However, the FOM calculation used for the ROM construction should be carried out more efficiently in practical applications. Thus, an efficient way to construct the present ROM is investigated in this section. The previous studies [18,19] reveal that coarse step calculations are an effective way to construct an orthogonal basis when the relative shape of the solution space, i.e. flux distributions, are comparable between the reference and coarse step calculations. However, the temporal discretization error will also affect the reconstructed coefficient matrices for the leakage term in the pre- sent ROM. Thus, the accuracy of the present ROM constructed from coarse time step calculation is compared with the reference FOM solution in the present verification. Table 3 shows the calculation conditions for the coarse step FOM calculation.
As shown in Table 3, the fully implicit method and the analytical solution assuming the linear transition of the fission sources between the successive time steps [11] are employed for the time integration for fluxes and precursors to reduce the computation time, respectively. The flux distributions are available only in 8 steps from the coarse time step FOM solution. The other calculation conditions are the same as the reference solution. Figure 8 shows the transition of the maximum relative difference of the flux distributions with respect to the reference solution, which is calculated as follows: max g2Gφ n 0 ;CTS g Àφ n;ref where the superscript CTS and ref denote the coarse time step calculation and the reference solution, respectively. n 0 and n denote the time step in the coarse time step calculation and the reference solution, where n is selected so that t n 0 ;CTS ¼ t n;ref . Note that this comparison focuses on the applicability of the coarse time step calculation for the orthogonal basis construction and both calculations are performed with the FOM. As shown in Figure 8, the flux distributions in the coarse time step solution include large temporal discretization error. However, the normalized flux distributions, which are employed for the orthogonal basis construction in the present study, are in good agreement with the reference solution. These results suggest that the coarse time step calculation can be an alternative to the reference solution for the orthogonal basis construction. Table 4 shows the case matrix for the ROM construction.
As shown in Table 4, there are 3 test cases. Case 4 is the ROM constructed from the coarse time step solution. In case 5, the ROM is constructed from the reference solution, but with only 8 flux distributions at the same times used for the coarse step calculation (t = 0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 6.5, 10.0 sec). The orthogonal basis for case 6 is constructed with all flux distributions in the reference solution. For all test cases, the coefficient matrices are evaluated at the same times used for the coarse step calculation. The compressed coefficient matrices are quadratically interpolated in all test cases as follows:  Table 5 shows the size of the solution space, the orthogonal bases and the coefficient matrices in the present calculation conditions. The low-rank approximation for the orthogonal bases is not applied in all cases. Figures 9 and 10 show the relative difference for the core power and the maximum relative difference for the flux distributions with respect to the reference FOM solution, where the relative differences are calculated as Equations (75) and (76).
Since the calculation conditions of case 3 in section 3.1 and case 6 are consistent except for the interpolation intervals for the coefficient matrices, the calculation results of cases 3 and 6 show that the interpolation error for the coefficient matrices has a negligible impact on the accuracy in the present calculation conditions. Comparing the case 5 and 6, the small degradation of the accuracy for the core power in TD3 indicates the degradation of the orthogonal basis due to fewer sampling of the flux distributions. A similar trend also appears in case 4 but the degradation of the flux distributions in Figure 10(d) and (g) shows the impact of the temporal discretization error for the orthogonal basis and the compressed coefficient matrices, which is reasonably small in the present calculation conditions. The maximum relative differences of the core power and flux transition between case 4 and the reference solution in the present verification are about 0.15% and 0.06%, respectively. These results suggest that the coarse step solution also can be an alternative to the reference solution in the present ROM by sampling a sufficient number of the typical flux distributions and coefficient matrices.  Table 6 shows the computation time of the reference solution and case 4. All calculations are carried out on Intel(R) Core TM i9-9900 K (3.6-5.0 GHz) with 16 GB memory. Parallel calculation based on ray-trace-wise decomposition is carried out using 16 threads for the FOM calculation. The ROM calculation is not parallelized in the present verification.
As shown in Table 6, the present ROM enables us much faster computation than the FOM, i.e. about 5000-6000 times faster computation is achieved for the ROM calculation itself. It is still about 2 times faster than the FOM if we consider the computation time for the coarse step calculations, SVD and DR to obtain the almost same accuracy. These results suggest the effectiveness of the present approach to construct a ROM for the time-dependent MOC calculations. As for the application of the present ROM, the online core monitoring systems could be one of the candidates because the calculation geometry is fixed during the operations and the flux distributions close to the current core status are periodically calculated. Using those distributions effectively, there is a possibility to construct a ROM with the accuracy of high fidelity simulations.

Conclusion
In the present study, an efficient ROM which reproduces the time-dependent MOC solution using the POD is proposed. In the present approach, the MOC-equivalent diffusion equation for each flux region is derived using the net neutron current calculated with MOC, and then the POD is applied to it. The flux and prompt/delayed fission rate distributions taken from a FOM solution are employed to make the orthogonal bases, which are separately defined for the scalar flux at each energy group and the precursor density at    Figure 9. Accuracy of core power transition (case 4-6).
each precursor family in the present study. The accuracy of the present ROM is verified in the C5G7-TD 2D benchmark problem. The calculation results show that the non-linear time dependence of the reconstructed coefficient matrices degrades the accuracy in the present approach. However, it can be reasonably mitigated with the quadratic interpolation for them and the present ROM accurately reproduces the reference MOC solution with high accuracy in the present verification. The maximum relative differences for the core power and flux distributions in the present verification are about 0.15% and 0.06%, respectively. The computation speed of the present ROM is much faster than the MOC calculations. The maximum speed-up ratio in the present verification is about 5000-6000 for the ROM calculation itself and the computation time for the ROM construction is almost half of the reference solution. The present results also indicate that coarse step FOM calculations are effective to construct the present ROM and the temporal discretization error included in the reconstructed coefficient matrices has less impact on the accuracy. These results suggest the effectiveness of the present approach  for the time-dependent MOC calculations. In the future study, the present approach is planned to be applied to the 3-dimensional MOC calculations and the treatment of the feedback effect and the partial control rod insertion during transient will be investigated. That could enable us real-time applications such as online core monitoring with the high accuracy of the pin-resolved transport calculations.

Disclosure statement
No potential conflict of interest was reported by the author(s).