An Extended Closure Relation by LightGBM for Neutrino Radiation Transport in Core-collapse Supernovae
Shota Takahashi Akira Harada Shoichi Yamada




tkst1228@g.ecc.u-tokyo.ac.jp


Department of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan


0000-0003-1409-0695]Akira Harada
National Institute of Technology, Ibaraki College, Hitachinaka 312-8508, Japan
Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS), RIKEN, Wako, Saitama 351-0198, Japan


Advanced Research Institute for Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan

§ ABSTRACT 
We developed a machine learning model using LightGBM, one of the most popular gradient boosting decision tree methods these days, to predict the Eddington tensor, or the second-order angular moment, for neutrino radiation transport in core-collapse supernova simulations. We use not only zero-th and first moments as in ordinary closure relations but also information on the background matter configuration extensively. For training the model we utilize some post-bounce snapshots from one of our previous Boltzmann radiation-hydrodynamics simulations. The Eddington tensor as well as the zero-th and first angular moments are calculated from the neutrino distribution function obtained in the simulation. LightGBM is light indeed and its high efficiency in training enables us to feed a large number of features and figure out which features are more important than others. We report in this paper the results of the training and validation as well as the generalization of our model: it can reproduce the Eddington factor better in general than the M1-closure relation, one of the most commonly employed algebraic closure relations at present; the generalization performance is also much improved from our previous model based on the deep neural network.







§ INTRODUCTION

Core-collapse supernovae (CCSNe) are explosive demise of massive stars. Despite extensive investigations, the detailed mechanism driving this explosion remains to be established. The neutrino heating <cit.> is the mechanism that most researchers in the society think is most likely at present. It goes as follows. The gravitational collapse of a massive star is halted when the core reaches nuclear density at the center. Then a shock wave is formed and sweeps the accreting outer core, photo-dissociating heavy nuclei therein. The shock energy depleted through this endothermic reaction, leading eventually to the shock stall. In the meantime, a proto-neutron star (PNS) is formed at the center and starts to emit neutrinos copiously. These neutrinos heat up matter behind the shock and re-energize the stagnant shock wave eventually. It should be obvious that the good understanding of neutrino transport is imperative for firmly establishing the theory of the CCSN mechanism.

Numerical simulations have been extensively employed for this purpose <cit.>. They comprise the computations of hydrodynamics, gravity, and neutrino transport and reactions. The neutrino transport, the target of the this study, is described by the Boltzmann equation. Its faithful solution, although preferable, is numerically demanding because of the high dimension. The majority of CCSN simulations these days resorts to be established to some approximations. Among them, the truncated moment method is most widely used <cit.>. Instead of the distribution function in momentum space at each spatial point, it considers the angular moments thereof and integrates the Boltzmann equation with respect to the two angles that specify the neutrino momentum, angular moments of the Boltzmann equation, thereby reducing the dimension.

The resultant equations for the moments of different orders are coupled with one another and form an infinite hierarchy. In numerical simulations, it is truncated at some order, typically the first order. The efficacy of this truncated moment method hinges on the accuracy of the closure relation, which one imposes by hand to obtain unsolved moments from the solved ones. If one truncates the hierarchy at the first order, the Eddington tensor, or the second moment of the distribution function, is normally given by the closure relation as a function of the zero-th and first order moments. Among various closure relations proposed in the literature to date, the M1-closure relation is the most popular these days <cit.>. According to the comparison of the Eddington tensors obtained by the directly integrating the distribution functions obtained in the CCSN simulations that solved the Boltzmann equation faithfully, which we refer to as the Boltzmann simulations hereafter, and the Eddington tensors given by the M1-closure relation for the same data, the diagonal and off-diagonal components of the Eddington tensor may develop discrepancies of several tens of percent and a factor of 2, respectively <cit.>.

Some efforts are still being made to enhance the accuracy of closure relation. <cit.> suggested a fitting formula for the Eddington factor, the largest eigenvalue of the Eddington tensor, based on the spherically symmetric Boltzmann simulations. Though it offers a reasonably faithful formula, more improvement is required to apply for the multi-dimensional simulations. In a separate endeavor, <cit.> developed a tensor basis neural network (TBNN), a deep neural network that provides the Eddington tensor from the neutrino energy density and flux, i.e., the zeroth and first order moments and the local matter velocity. The network was trained by using the data of a Boltzmann simulation. Although the TBNN-closure outperforms the M1-closure in terms of accuracy, there remains room for improvement both in the accuracy and in the generalization performance.

This paper aims at improving the closure relation by machine learning techniques. From our experience in the previous attempt with the TBNN <cit.>, we think that a substantial improvement will be achieved only by providing more information other than the two lowest order moments to the machine learning model. In fact, <cit.> found that the information on local matter motions was useful as an input to the TBNN. One particularly powerful technique to identify relevant features is the decision tree for regression problems (see Section 2.2 for details). We adopt Light Gradient Boosting Machine (LightGBM) in this paper. It is known that LightGBM can handle a large number of features very efficiently. As one of the decision tree algorithms, on the other hand, it may lack a sufficient smoothness in expressing functions. We attach more importance to  the former capability in this study.

This paper is organized as follows. In section 2, we briefly review the M1-closure relation and LightGBM, the machine learning model employed in this paper. Section 3 provides the details of the training: data, feature engineering, validation and hyper parameters. Then, in section 4, we present the main results on the performance of our novel models making comparisons with the Eddington tensors calculated directly from the original data obtained by our Boltzmann simulation as well as with the Eddington tensors derived from the M1-closure relation applied to the zero-th and first angular moments calculated from the same original data. Finally, section 5 is the conclusion. Throughout this paper, we employ the units of c=1, where c is the light speed. The metric signature is (-+++). Greek and Latin indices run over 0–3 (spacetime) and 1–3 (spatial).



§ OVERVIEWS OF THE M1-CLOSURE RELATION AND LIGHTGBM

The purpose of this paper is to improve by machine learning techniques the estimation of the Eddington tensor in the truncated moment method for the neutrino transport in CCSN simulations. This section is devoted to a brief review of the M1-closure relation in section <ref> and LightGBM in section <ref>.



 §.§ M1-closure relation


The neutrino population in CCSNe is described by the distribution function f(x,p) in phase space, in which x and p are the spacetime coordinate and momentum of neutrinos. It is governed by the Boltzmann equation. The high dimension and stiff source term of this equation poses a numerical challenge under the current computational resources. To the numerical cost, most researchers are opting the truncated moment method nowadays.

Since the formulation of the moment method is described in <cit.> in detail, we give only its outline in the following. We first define the unprojected second moment of the distribution function as[We follow the definition in <cit.>, which is slightly different from that in <cit.>.]

    M^αβ(ε, x ) = ∫ f(x, p^') δ(ε^3/3-ε^' 3/3) p^'α p^'β d V_p^',

where ε^'=-u_μ p^'μ is the neutrino energy measured in the fluid-rest (FR) frame with u being the fluid four-velocity. By applying appropriate projections to this moment, we can derive the energy densities E, energy flux 
F^i, and the stress tensors P^ij in both the FR and the laboratory (LB) frames:

    E_ FR = M^αβ u_α u_β,  E_ LB = M^αβ n_α n_β,


    F_ FR^i = -M^αβ u_α h_β^i,   F_ LB^i = -M^αβ n_αγ_β^i,

and

    P_ FR^ij = M^αβ h_α^i h_β^j,   P_ LB^ij = M^αβγ_α^i γ_β^j,

where n is the unit normal vector to the spatial hypersurface with t = const., γ_α^i = δ_α^i + n_α n^i is the projection onto that surface with δ_α^β being the Kronecker's delta, and h_α^i = δ_α^i + u_α u^i is the projection onto the plane perpendicular to u^α. In the truncated moment formalism, the time evolutions of E_ LB and F_ LB^i are solved normally by setting a closure relation to give P_ LB^ij as a function of E_ LB and F_ LB^i[In order to solve the spectral evolution of spectrum in general relativity, the third moment of the distribution function is also required. We ignore it in this paper for simplicity, though.]. These three quantities are essentially the lowest three angular moments. The exact relations are given, e.g., in <cit.>.

The closure relation is normally imposed to the Eddington tensor k^ij defined as the stress tensor divided by the energy density. A common way to construct it is to employ the interpolation
of the optically thick and thin limits:

    P_ LB^ij = 3p_ν-1/2 P_ thin^i j + 3(1-p_ν)/2 P_ thick^i j.

In the above expression, P^ij_ thin and P^ij_ thick are the stress tensors in the optically thick and thin limits, resepctively, and can be expressed with the energy density and energy flux density as given shortly; p_ν is the Eddington factor and should be also prescribed in the truncated moment method. Then the Eddington tensor is obtained as

    k^ij = P_ LB^ij/E_ LB.

The functional forms of the stress tensor in the two limits and of the Eddington factor is given as follows <cit.>. In the optically thin limit, neutrinos move freely in one direction specified by the flux vector. The stress tensor P_ thin^ij is then given in the LB frame as

    P_ thin^i j = E_ LBF_ LB^i F_ LB^j/|F_ LB|^2.

In the optically thick limit, on the other hand, neutrinos are in thermal equilibrium with matter, and are isotropic in the FR frame. In the LB frame, P_ thick^ij is given by the Lorentz transformation as

    P_ thick^i j   =   1/3 E_ FR(γ^ij + 4 v^i v^j) 
            + F_ FR^i v^j +  F_ FR^j v^i.

As for the Eddington factor, the M1 closure scheme adopts the following expression proposed by <cit.>:

    p_ν = 3 + 4 χ^2/5 + 2 √(4 - 3 χ^2),

where χ=|F_FR|/E_FR is referred to as the flux factor.



 §.§ Light Gradient Boosting Machine


The decision tree (DT) is a machine learning algorithm that performs classification or regression tasks by recursively partitioning a dataset based on certain criteria <cit.>. Suppose that we find a functional relation between an n-dimensional input vector x and an output real number[If the output is integer or discrete number, the task is classification.] y from a dataset of size N: {( x^i, y^i)}_i=1^N. In the DT method, each data ( x^j, y^j) is binary-classified according to x^j, and the the two destinations are called nodes. The data classified to a certain node is further classified to the daughter nodes recursively. This process hence forms a tree graph. In the DT regression, each node has the prediction value ŷ_k = ∑_j∈ D_k y^j/|D_k|, where D_k is a subset of the data classified to node k and |D_k| is its size. The classification criterion at each node is determined so that it would minimize the mean squared error of the prediction value ŷ_R(L)k for the data in the right (left) daughter node of node k. This binary classification is repeated until the error thus obtained gets small enough. Then the output of the DT is the prediction value ŷ^j of the node, to which the input x^j is finally classified.

The Gradient Boosting DT <cit.> is an improved version of DT. The output of single DT for regression has large errors normally. We can improve the prediction by training another DT for these errors. With recursive correctios with many DTs, the GBDT achieves an accurate prediction. Suppose that f_ℓ (x) is the output of the ℓ-th DT for the input vector x. In the least-square GBDT, each DT is trained so that the following functions

    ∑_j=1^N{(y^j-η∑_m=0^ℓ-1f_m(x^j))-f_ℓ(x^j)}^2

are employed from ℓ=1; in the above expression, 0< η < 1 is the learning rate employed to avoid overfitting and f_0(x) is given conventionally by ∑_j=1^N y^j/N. Then, the prediction of GBDT is given by ŷ^j = η∑_ℓ=0^M f_ℓ (x^j), where M is the total number of DTs.

LightGBM (Light Gradient Boosting Machine) is a
variant of GBDT introduced by <cit.>. LightGBM incorporates four techniques to enhance GBDT:  (1) leaf-wise tree growth that prioritizes creations of daughter and granddaughter nodes than those of cousin nodes, suppressing the growth of 
the nodes that focus on irrelervant variables for regression, as shown in figure <ref>, (2) a histogram-based algorithm, which reduces numerical cost by binning data, (3) the gradient-based one-sided sampling (GOSS) to reduce the training data by focusing on the critical elements, and (4) the exclusive feature bundling, which combines (EFB), which combines features that are not zero simultaneously. Thanks to these techniques, LightGBM excels in computational and memory efficiency compared to traditional models, hence it accelerates the search for the relevant variables of the closure relation.




§ METHOD

This section provides a detailed description of the model developed in this paper. Section 3.1 concerns the data used in this study, section 3.2 discusses the feature engineering, and section 3.3 explains the model.



 §.§ Data Description

The data set used in this study are five snapshots taken from a 2D axisymmetric CCSN simulations performed with the Boltzmann-radiative-hydrodynamics code <cit.>. This code solves the Boltzmann equation for neutrinos, the hydrodynamics equations for stellar gas and the Poisson equation for gravity. Simultaneously, we adopt Furusawa–Togashi EOS <cit.>, which is based on the variational method and is extended to handle a large number of nuclei in the nuclear statistical equilibrium. The progenitor is a non-rotating star with the ZAMS mass of 15 M_⊙ modeled by <cit.>. For the detailed description of this code we refer readers to <cit.>, <cit.>, <cit.> and <cit.>. Whereas we consider three neutrino species: electron-type neutrinos and anti-neutrino as well as heavy-lepton-type neutrino, we employ the results for the electron-type neutrino alone in this study. The grid sizes for radius N_r, zenith angle N_θ, energy N_ε, momentum angles N_θ_ν, N_ϕ_ν are (N_r, N_θ, N_ε, N_θ_ν, N_ϕ_ν) = (384, 128, 20, 10, 6). The computational zone covers the region up to 5000 km from the stellar center. In this study, on the other hand, we utilize the data within the ranges of 0<r(km)<200, 0<θ(rad)<π, and 0<ε(MeV)<300 because this is the region of the greatest relevance for the neutrino heating and shock revival. We extract the distribution function of electron-type neutrino in this region from the simulation data at the post-bounce time of 300  ms and use it as the training data. For validation, we employ four other distribution functions sampled at t=100 ms, t=150 ms, t=200 ms and t=250 ms after bounce. Note that the CCSN simulation we adopt in this paper is identical to what was used in <cit.>, but the radial extent and the number of snapshots are expanded.

We will construct a machine learning model to give the Eddington tensor as a function of some variables, or features, that include the lower moments. Those we refer to as the basic features are listed in table <ref>. As discussed later, additional features are generated from these basic ones through feature engineering. The basic and additional features are collectively treated as inputs to the machine.

The output of the machine is the Eddington tensor. In this paper, we do not project it onto the local orthonormal vectors: e_r, e_θ and e_ϕ, of the spherical coordinates in the LB frame. Instead we employ the orthonormal vectors, e_ FP, i, one of which is aligned with the neutrino flux and the other two are orthogonal to it and also to each other. We refer to the 3-dimenstional frame defined by these orthonormal vectors as the flux-projected (FP) frame, and the Eddington tensor represented in this frame is denoted as k^ij_ FP hereafter. 
The FP frame is constructed for each neutrino energy as follows: e_ FP, 1 is the unit vector directed along the neutrino flux. The choice of the other two unit vectors is rather arbitrary. In this paper we obtain e_ FP, 2 from e_θ by subtracting the component parallel to e_ FP, 1, i.e., via the Gram-Schmidt orthogonalization. As a result, e_ FP, 2 so obtained lies in the plane spanned by the flux and e_θ. The last vector is derived so that e_ FP,1, e_ FP, 2 and e_ FP,3 should form a right-handed orthonormal system. The Eddington tensor in the FP frame is then related with that in the LB frame as k_ FP^ij e_ FP,i^l e_ FP,j^m = k_ LB^lm, where the superscript of e_ FR specifies the components and the repeated indices are summed over.

The reason why we employ the FP frame is the following. The Eddington tensor is not “aligned” in general with the coordinates we choose rather arbitrarily. Since the Eddington tensor is a symmetric second-rank tensor, it becomes diagonal in the appropriately chosen orthonormal frame, which is not equal to the LB frame in general. Suppose that the neutrino distribution in momentum space is axisymmetric with respective to a certain direction misaligned with any one of the orthonormal vectors in the LB frame. Instead the symmetry axis is parallel to the flux vector in this case. In fact, the Eddington tensor is diagonal in the FP frame. It should be stressed that the neutrino distribution in momentum space is not axisymmetric in general, except for the notable case with spherical symmetry. For example, the neutrino flux from a rotating oblate PNS is non-radial and has even nonvanishing ϕ-component, since neutrinos carry non-zero angular momenta. Then the Eddington tensor in the LB frame has substantial off-diagonal components <cit.>. We expect, on the other hand, that the off-diagonal components will be smaller in the FP frame even in this case.





 §.§ Feature Engineering

The closure relation is originally meant to give  high-order moments such as the Eddington tensor in terms of the lower-order moments like the energy density and flux. In our previous attempt to build a deep neural network to do this task, we realized that the lower-order moments are just not sufficient and added some other local features such as the matter velocity and its shear to the input. In this paper, we extend this idea. Thanks to the efficiency of LightGBM, we are able to incorporate much more information, some of them nonlocal. It is true that the more model-specific information we employ, the less general the machine becomes, but the purpose of this paper is to construct a machine learning model that we can apply to the CCSN simulation with the truncated moment method. 

In providing  a large number of features to the machine, we find feature engineering very important. By feature engineering, we mean finding better combinations of the original features, normalizing features and, if necessary, removing the basic features after using them for feature engineering. In principle, this process does not increase the information included in the original features. If the output depends on a function of rather involved combinations of these features, the machine may need a long time or even fail to train itself even if those features contain sufficient information. On the other hand, if the human knows a priori those combinations of the original features and teaches them to the machine, we can expect the training efficiency to be improved substantially. Unfortunately, we do not know what should be the best combinations. Thanks to the efficiency of LightGBM, we can try various combinations. The engineered features given below are indeed obtained that way.

As mentioned above, we will incorporate nonlocal quantities to the input features. This is probably understandable if one recalls the fact that the radiative transfer is most difficult to treat in the semi-transparent region, where the mean free path of neutrino becomes comparable to the scale height of matter and the neutrino distribution is not determined locally indeed. It is noted that the basic features in Table 1 include the optical depth, which is nonlocal information. As explained below, we add more.

The engineered features thus produced are classified into three types: advanced features, spatial-shift features, and spatial-difference features. The advanced features, listed in Table <ref>, are those features given as functions of the basic features. Particularly important for the improvement of the model are the flux factor in the LB frame χ_ LB = |F_ LB^i|/E_ LB, the colinearity parameter F_LB^j v_j/E_LB |v^i| that measures the alignment of the neutrino flux and matter velocity, and the parameter we call the“neutrinosphere indicator,” which is a binary indicator of whether the point is well inside the neutrinosphere or not, i.e., it is 1 if the optical depth τ is greater than 5 and 0 otherwise.



The spatial-shift features are the features concerning the neighboring grids. The closure relation is normally local, i.e., it is a relation among physical quantities at the same location. It turns out, however, that the inclusion of the neighboring physical quantities as input improves the accuracy of the inference. More concretely, we incorporate the flux factor χ_ LB=|F_ LB^i|/E_ LB and the colinearity parameter F_LB^jv_j/E_LB|v^i| not only on the grid point of concern but also on the neighboring grid points. For instance, let Q_k,l represent one of the aforementioned quantities at the k-th radial grid point and the l-th angular grid point; we then incorporate Q_k-5i, l (where i=-1, 1, 2, …, 6), Q_k+1, l, Q_k, l+1, Q_k-1, l and Q_k, l-1 as components of the input feature vector for Q_k,l.[If it goes beyond the spatial boundaries, those points are simply ignored.] In particular, it was found that incorporating information from the inner regions of the star, such as Q_k-5i, l (where i=1, 2, …, 6) and Q_k-1, l, significantly improves the prediction accuracy. The reasoning here is that since neutrinos predominantly move outward, upstream information from radially inner points is particularly valuable.

We also find it useful to employ the difference between the feature on the grid of concern and that on the inner adjacent grid:

    Δ Q_k,ℓ = Q_k,ℓ-Q_k-1,ℓ.

to which we refer as the spatial-difference feature. Although it may seem redundant, we include both the spatial-shift and difference features in the input as we do not know which feature is more effective. Note that irrelevant features are automatically ignored in LightGBM.



 §.§ Models

The pipeline of the model employed in this study is shown in Figure <ref>. After creating new features from the basic features through feature engineering, we conduct training with the 8-fold cross validation (CV), in which the data set is randomly divided into 8 subsets to build 8 models. One of them serves as the validation set and the remaining 7 as the training data. They make predictions individually and the arithmetic mean of those predictions is adopted as the final prediction of this model. Such an approach is known to effectively prevent the model from over-fitting a particular dataset, thereby improving generalization.



The main hyperparameters and their values in the LightGBM model are listed in Table <ref>. We employ the early stopping training, i.e., the training is terminated once the mean squared error is no longer lowered remarkably. It is also useful to avoid overfitting. As a consequence, we can set a sufficiently large value to , the number of decision trees to be built. Though it is not our main goal in this paper to fine-tune the hyperparameters, we set different values to some hyperparameters for the diagonal components of the Eddington tensor from those for the off-diagonal components as given in the table.





§ RESULTS


In this paper, the data extracted from the snapshot at 300 ms after the bounce are used as the training set. The target variable is the Eddington tensor in the FP frame k_FP^ij. We focus on the estimation of its 11-, 22-, 33- and 12-components, since they are the diagonal and main off-diagonal elements, respectively. The test data were collected from the snapshots at 100 ms, 150 ms, 200 ms and 250 ms after the bounce. In the following, the diagonal and off-diagonal components are considered separately in turn.




 §.§ Diagonal Components

The diagonal components of the Eddington tensor play the main role in neutrino transport, and hence their accurate estimation is crucially important. The conventional M1-closure relation does not reproduce it perfectly <cit.>, and the machine-learning closure is expected to provide better predictions. 

Figure <ref> shows the learning curve in the training of the 11-component of the Eddington tensor. The training is conducted to minimize the Mean Squared Error (MSE):

    L_ε^(2) = 1/|D_ε|∑_D_ε|k_ boltz, FP^ij-k_ infer, FP^ij|^2,

where the sum is taken over the subset D_ε of the data corresponding to a particular neutrino energy and |D_ε| is the number of its elements; k_ boltz, FP^ij is the Eddington tensor obtained by the Boltzmann-radiation-hydrodynamics simulation whereas k_ infer, FP^ij is the one obtained with our LightGBM model. As the solid lines show, the training is successful with the values of MSE decreasing monotonically down to <10^-5 for all 8 models. These lines are terminated at different numbers of estimator because of the early stopping. The evaluation of the validation data also presented in the figure as dashed lines indicates that the overfitting is avoided reasonably well indeed. Similar results are obtained for other components of the Eddington tensor.



We move on to the predictions of these models. In Figures 4, 5, 6 we plot the mean absolute error (MAE) for the diagonal components of the Eddington tensor for different post-bounce times: 250, 200, 150, and 100 ms. MAE is defined as

    L_ε^(1) = 1/|D_ε|∑_D_ε |k_ boltz, FP^ij - k_ infer, FP^ij|,

where k_ infer, FP^ij is the ij-component of the Eddington tensor obtained either with our Light-GBM model or via the M1-closure relation. It is apparent that our LightGBM closure gives more accurate predictions than the M1 closure for almost all energies at these times. The behavior of MAE is essentially the same among the three diagonal components. As the time goes back to earlier times, the accuracy is degraded somewhat as expected. It is encouraging, however, that the results are still better than those for the M1 closure except at ∼ 30 MeV, where the M1 closure gives exceptionally good results.







Although MAE is useful to see the overall convergence of the model to the target, what is more important is to what extent the individual components of the Eddington tensor are reproduced at different points in space. We hence compare next the radial profiles of the components of the Eddington tensors along the radial line at θ=π / 2 in figures <ref>–<ref>. The results are not much different for other angles.


























We can confirm again that the LightGBM closure outperforms the M1 closure in general. Of particular note is that our LightGBM model can reproduce the values of the 11-component of the Eddington tensor smaller than 1/3 and the 22- and 33-components greater than 1/3 at rather small radii. Those values are obtained when the opacity decreases rapidly with the radius in the vicinity of the neutrinosphere; then the neutrino distribution is almost isotropic in the outward hemisphere whereas it  is heavily depleted in the opposite hemisphere <cit.>. By its construction the M1 closure cannot treat such angular distributions properly. It is noted that our LightGBM fairs even better than the TBNN model in <cit.>.

It is also worth mentioning that the LightGBM model can handle the shock wave with reasonable accuracies. The shock wave is located at 170  km, 140  km, 130  km, 120  km, and 110  km at 100  ms, 150  ms, 200  ms, 250  ms, and 300  ms post bounce, respectively. The Eddington tensor is discontinuous at the shock wave. The jump gets larger as the energy increases. The radial profiles of the Eddington tensor for the energy of 54  MeV presented in figures <ref>–<ref> demonstrate indeed that the predictions of the LightGBM model are closer to the results of the Boltzmann simulations than the M1-closure relation in general.

However, the Eddington tensor predicted by LightGBM has several issues. Firstly, the predicted values are not so smooth as those obtained with the M1 closure. This tendency becomes more pronounced as we go back in time and see the predictions themselves degraded. In fact, the Eddington tensor at 250 ms post bounce plotted in Figures <ref>, <ref>, <ref> is relatively smooth, whereas the lines are more jagged at 100 ms after bounce as shown in Figures <ref>, <ref>, <ref>. We think that this phenomenon occurs because the current predictions are made for each grid point individually and that the issue could be resolved by adopting a machine learning model capable of incorporating non-local information into the predictions. Machine learning models such as Convolutional Neural Networks (CNNs) <cit.> or Recurrent Neural Networks (RNNs) <cit.> may have robustness to local noises and outliers in the input data, thereby giving smoother predictions. 

It should be also mentioned that the current model does not take into account physical requirements. For instance, the trace of the Eddington tensor is unity by definition. There is no guarantee, however, that the predicted results satisfy this constraint although the violation is quite minor in the current models, as can be seen in Figure <ref>. This issue may be addressed with physics-informed neural networks <cit.>: for example, one may incorporate the violations of constraints into the loss function so that they could be minimized. Although these possibilities are interesting on their own rights and warrant further investigations, they are much beyond the scope of this paper and will be future works.





 §.§ Off-Diagonal Component

Although the diagonal components, particularly the 11-component, is dominant in neutrino transport, the off-diagonal components should not be forgotten <cit.>. Designed to reproduce the dominant component, the M1-closure is known to give large errors to the off-diagonal components sometimes <cit.>. It is hence one of our goals in the machine learning modeling of the Eddington tensor to better reproduce the off-diagonal components. It is admittedly true in this respect that the previous TBNN closure was not so good, either <cit.>. In the following, we examine the accuracy of the off-diagonal components in our LightGBM model.
 
We focus on the 12-component, which is actually the most important off-diagonal component. Figure <ref> shows the MAE for this component as a function of the neutrino energy at different postbounce times. We find again that our LightGBM model achieves higher accuracies in general compared to the M1 closure. The radial profiles of the 12-component of the Eddington tensor are presented in figures <ref>–<ref>. They are obtained for the radial ray at θ=π/2. At 250 ms post bounce, our LightGBM model fairs better than the M1 closure except at ε_ν=54.0MeV, where the two models give similar errors. As we go back in time, the deviation from the results of the Boltzmann simulation gets larger for the LightGBM model and its advantage becomes not so remarkable.











A potential factor that contributes to the compromized success is the variation in the statistical properties among the data at different time steps. Tables <ref> and 
<ref> give the the mean and standard deviation of k_ boltz, FP^11 and k_ boltz, FP^12 at each time step. It is evident that the statistical properties of k_ boltz, FP^12 are significantly different from time to time compared to those of k_ boltz, FP^11. It may be that the generalization should not have been expected for the  off-diagonal component in the first place. Standardisation of the dataset across all time steps may be a solution but again is out of the scope of this paper.

It is worth noting that different features were utilized for the diagonal and off-diagonal components in our LightGBM model. Specifically for the off-diagonal components, it was unclear which features held the most significance. In predicting the diagonal components, we utilized 35 features, whereas for the off-diagonal components, 317 features were employed. For example, features such as fluid velocity, which were removed when predicting diagonal components, are included. When predicting diagonal components, the prediction accuracy significantly decreased; however, in predicting off-diagonal components, the accuracy was not greatly compromised. Off-diagonal components are more significantly affected by factors such as fluid motion compared to diagonal components, which suggests that background information may be more important for their prediction.







 §.§ Importance of feature engineering

Feature engineering plays a crucial role in predicting the Eddington tensor accurately. In this work we incorporated many advanced features indeed as listed in table <ref>. It is possible because LightGBM is very efficient both in training and prediction, ignoring irrelevant features automatically. As explained in section 3.2, we find that the spatial-shift and difference features are two most important features for the improvement of accuracy, which we demonstrate here. 

In Figure <ref>, we show how MAE gets better from the one for the model with only the basic features included. The MAE for the model with the basic features alone is actually larger than that for the M1 closure method for most of the neutrino energy. The MAE becomes smaller, however, once the advanced features are incorporated, particularly for energies close to the mean energy. It is further improved if the spatial-shift and difference features are taken into account in addition. The same trend is observed for the Eddington tensor k_FP^11 itself as a function of radius as demonstrated in Figure <ref>.

What we learned from these experiments on feature engineering may be three-fold. Firstly, the local low-order moments, the inputs for the ordinary closure methods, are not sufficient to reproduce the Eddington tensor accurately. We find that the information on the background matter such as the baryon density or the optical depth is particularly important. Secondly, the nonlocal information on the neutrino distribution is probably the most crucial. Although their incorporation in LightGBM is rather straightforward, it will not be so easy for the analytic closure relation to take them into account. Thirdly, finding nice combinations of the basic features is also an important thing. Our LightGBM model is useful in this respect, since we can try many possible combinations in relatively short times. The new features so obtained may be employed in other machine learning models.







§ CONCLUSION

In this paper we develop a machine learning model using LightGBM to predict the Eddington tensor, or the second moment, of the angular distribution of neutrinos often employed in the core-collapse supernova simulation. Unlike the ordinary closure relations, we employ not only low-order moments of the neutrino angular distribution in momentum space but exploit also the information on the matter distribution as well as non-local features of the neutrino distribution. For the training and validation of the machine, we utilize the numerical result of the Boltzmann-radiation-hydrodynamics simulations of the core-collapse of a non-rotating 15M_⊙ progenitor model. The training data are a snapshot at 300 ms after core bounce, while the validation data are snapshots taken at 100 ms, 150 ms, 200 ms and 250 ms post-bounce. For all these times, our machine learning model shows better accuracy than the M1 closure in predicting the Eddington tensor, highlighting the potential of machine learning to be an alternative closure relation.

The key factor in the model building is the feature engineering process. It is true that finding features of the greatest relevance in improving the prediction is the most critical but it is equally important to produce nice combinations of the basic features that the machine can easily understand. Those features should have some physical meanings so that they could be applied to different evolutionary stages. The LightGBM, a variant of the Gradient Boosting Decision Tree, is efficient in the management of memory and time, allowing us to handle a large number of features at a time quickly and, as a result, identify the features that improve generalization ability substantially. In fact, LightGBM is able to reproduce better both diagonal and off-diagonal components of the Eddington tensor in most cases than the M1 closure although the accuracy degrades somewhat as we go back to earlier post-bounce times. In particular, the dip in the 11-component near the neutrino sphere is represented quite well, which not only the M1 closure but our previous tensor-basis neural network model is unable to do. The hump in the 12-component in the same region is represented reasonably as well.

Although these successes are encouraging, there are admittedly many unsatisfactory things. For instance, the resultant Eddington factor is not very smooth in radius, which may cause a trouble when we substitute it for the M1 closure relation in the truncated moment method and perform core-collapse simulations. Although it is much improved from our previous neural network model, the generalization performance is not good enough. This may suggest that we should use all data collected from different times for the training. We are currently trying such possibilities. The non-local features investigated in this paper are based on physical intuitions and admittedly very primitive, and there is a much room for improvement. The ultimate goal of this project is to create a versatile model capable of predicting the Eddington tensor (and even higher moments as well) for stars with different masses and/or degrees of rotation. It may turn out that the machine based on the decision tree is not the best one. Even in that case, the features we identify in this work may be used as the inputs to other machine learning models. The current model is useful to find other features to improve the prediction further.



§ ACKNOWLEDGMENTS


We acknowledge Hiroki Nagakura and Keiya Hirashima for fruitful discussions. 
This work was supported by KAKENHI Grant Numbers JP21K13913, 21H01083.
This work was also supported by MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (Toward a unified view of the universe: from large scale structures to planets).
S. T. is supported by International Graduate Program of Innovation for Intelligent World.
S. Y. is supported by Institute for Advanced Theoretical and Experimental Physics, Waseda University and the Waseda University Grant for Special Research Projects (project number: 2023-C141, 2024-C56, 2024-Q014).
We acknowledge the high-performance computing resources of the K-computer / the supercomputer Fugaku provided by RIKEN, the FX10 provided by Tokyo University, the FX100 provided by Nagoya University, the Grand Chariot provided by Hokkaido University, and Oakforest-PACS / Wisteria Odyssey provided by JCAHPC through the HPCI System Research Project (Project ID: hp130025, 140211, 150225, 150262, 160071, 160211, 170031, 170230, 170304, 180111, 180179, 180239, 190100, 190160, 200102, 200124, 210050, 220047, 220223, 230056, 230270, 240041) for producing and processing the supervisor data.


LightGBM (<cit.>)


aasjournal

