Auxiliary Physics-Informed Neural Networks for Forward, Inverse, and
  Coupled Radiative Transfer Problems
In this paper, we develop and employ auxiliary physics-informed neural
networks (APINNs) to solve forward, inverse, and coupled integro-differential
problems of radiative transfer theory (RTE). Specifically, by focusing on the
relevant slab geometry and scattering media described by different types of
phase functions, we show how the proposed APINN framework enables the efficient
solution of Boltzmann-type transport equations through multi-output neural
networks with multiple auxiliary variables associated to the Legendre expansion
terms of the considered phase functions. Furthermore, we demonstrate the
successful application of APINN to the coupled radiation-conduction problem of
a participating medium and find distinctive temperature profiles beyond the
Fourier thermal conduction limit. Finally, we solve the inverse problem for the
Schwarzschild-Milne integral equation and retrieve the single scattering albedo
based solely on the knowledge of boundary data, similar to what is often
available in experimental settings. The present work significantly expands the
current capabilities of physics-informed neural networks for radiative transfer
problems that are relevant to the design and understanding of complex
scattering media and photonic structures with applications to metamaterials,
biomedical imaging, thermal transport, and semiconductor device modeling.
Roberto Riganti, Luca Dal Negro
Auxiliary Physics-Informed Neural Networks for Forward, Inverse, and Coupled Radiative Transfer Problems

###### Abstract

In this paper, we develop and employ auxiliary physics-informed neural networks (APINNs) to solve forward, inverse, and coupled integro-differential problems of radiative transfer theory (RTE). Specifically, by focusing on the relevant slab geometry and scattering media described by different types of phase functions, we show how the proposed APINN framework enables the efficient solution of Boltzmann-type transport equations through multi-output neural networks with multiple auxiliary variables associated to the Legendre expansion terms of the considered phase functions. Furthermore, we demonstrate the successful application of APINN to the coupled radiation-conduction problem of a participating medium and find distinctive temperature profiles beyond the Fourier thermal conduction limit. Finally, we solve the inverse problem for the Schwarzschild-Milne integral equation and retrieve the single scattering albedo based solely on the knowledge of boundary data, similar to what is often available in experimental settings. The present work significantly expands the current capabilities of physics-informed neural networks for radiative transfer problems that are relevant to the design and understanding of complex scattering media and photonic structures with applications to metamaterials, biomedical imaging, thermal transport, and semiconductor device modeling.

## I Introduction

Over the past few years, there has been a growing interest in developing deep learning (DL) and artificial intelligence (AI) algorithms for electromagnetic wave engineering, metamaterials design, and radiative transport problems [1; 2; 3]. Rapidly emerging approaches include training artificial neural networks (ANNs) to solve complex inverse problems, parameter estimation in structured photonic environments, and in strongly scattering media [4; 5; 6; 7; 8]. Although successfully demonstrated with respect to several inverse design problems, traditional methods remain essentially data-driven techniques and require time-consuming training steps and massive datasets [9; 10]. In order to improve on purely data-driven methods, it is essential to constrain and regularize them by leveraging the underlying physics of the investigated problems, thus relaxing the burden on training and data acquisition. Building on the firm foundation of the universal approximation theorem for multi-layer ANNs [11; 12], physics-informed neural networks (PINNs) have recently emerged as a powerful framework for the efficient solution of both forward and inverse problems mathematically described by partial differential equations (PDEs) of integer or fractional orders [13; 14; 15; 16]. The approach of PINNs has been successfully applied to a number of differential problems in engineering ranging from Navier-Stokes fluid dynamics, solid mechanics, and thermal transport [17; 18; 19]. Moreover, PINNs have shown remarkable results and noise robustness in the solution of electromagnetic inverse problems for metamaterials design, radiative transfer, imaging, and in the parameter retrieval of resonant photonic nanostructures [20; 21; 22; 23]. However, the solution of Boltzmann-type, integro-differential transport equations using PINNs still poses significant challenges due to the need to resort to numerical quadrature methods such as Gauss-Legendre or Gauss-Chebyshev for the approximation of the integral terms [24]. Such methods add computational complexity and inevitably introduce quadrature errors in the numerical solutions [25; 26; 27].

In order to eliminate such problems, a new PINN framework called auxiliary physics-informed neural networks (APINNs) was recently introduced by Yuan et al. [28]. This approach allows one to recast integro-differential equations into equivalent differential ones through the introduction of a network architecture containing additional auxiliary variables at its output, each corresponding to an integral term in the original, constrained by suitable relations. Therefore, the APINN formulation avoids the numerical approximation of integrals that are instead directly "guessed" by the network at a minimal cost, significantly improving both the numerical accuracy and computational efficiency compared to traditional PINNs.

In this paper, we develop a general APINN framework for solving relevant forward and inverse integro-differential transport equations of radiative transfer theory, which is a domain of vital importance in science and engineering with applications to complex photonic devices, medical imaging, metamaterials, thermal transport, as well as astrophysics, climate dynamics, and nuclear engineering [29; 30; 31]. In particular, we address and demonstrate APINN formulations for the accurate solution of forward, inverse, and coupled radiation-conduction problems of radiative transport in the relevant slab geometry for different choices of scattering phase functions.

Our paper is organized as follows: in Section II, we will provide a brief introduction to the radiative transfer equation (RTE), along with a description of the general APINN employed throughout this paper. In Section III.1, we dis
cuss forward problems for different phase functions governing the scattering processes. Specifically, we present benchmarked solutions for isotropic, Rayleigh, and Henyey-Greenstein scattering phase functions that are often utilized in engineering applications [29; 30; 32]. In Section III.2, we discuss the APINN solution of a coupled radiation-conduction problem, enabling the accurate description of radiation transfer in a participating medium. Lastly, in Section III.3, we show the solution of a canonical inverse problem described by the Schwarzschild-Milne integral equation, and we show that the radiative intensity solution and the single scattering parameters are accurately retrieved solely based on intensity data at the boundaries of the slab.

Our work shows that APINNs possess the flexibility, accuracy, and robustness required to become a powerful tool for inverse scattering and thermal transport modeling beyond the limitations of Fourier theory. Therefore, this work expands significantly upon the current capabilities and range of applications of PINNs methods and paves the way to the study of higher-dimensional transport problems in strongly scattering media with applications to nanophotonics, metamaterials, biomedical imaging, and optoelectronic device modeling.

## II Apinns for radiative transfer problems

The framework of radiative transfer theory for the study of complex scattering media was originally developed in astrophysics as a way to quantitatively describe the radiative equilibrium in interstellar clouds, planetary and stellar atmospheres [33]. Radiative transfer theory has found a very wide range of applications beyond astrophysics, including biomedical optics [32], atmospheric science [31], radiation hydrodynamics [34; 35] and remote sensing [36; 37]. For example, the propagation of light through fogs and clouds, white pains or paper, milky and turbid liquids, human tissue, and the brain can be adequately described by the classical theory of radiation transfer that we discuss in this paper using APINNs. The radiation transfer theory is founded upon the RTE, which is a Boltzmann-type integro-differential equation expressing the detailed energy balance for the propagation of directed energy flow, or radiance, through a multiply scattering discrete random medium. For scalar waves in three spatial dimensions the RTE can be written as follows:

\[\begin{split}\frac{1}{c}\frac{\partial I(\mathbf{r},\hat{s},t)}{ \partial t}&=-\hat{s}\cdot\nabla I(\mathbf{r},\hat{s},t)-(\kappa+ \sigma)I(\mathbf{r},\hat{s},t)+\\ \sigma\int_{4\pi}I(\mathbf{r},\hat{s}^{\prime},t)p(\hat{s}^{\prime}, \hat{s})d\Omega^{\prime}+S(\mathbf{r},\hat{s},t)\end{split} \tag{1}\]

where \(\kappa\) and \(\sigma\) are the absorption and scattering coefficients, respectively. Here \(S(\mathbf{r},\hat{s},t)\) denotes a generic source term and \(p(\hat{s}^{\prime},\hat{s})\) is the phase function describing the angular distribution of the scattering process. Alternatively, after introducing the optical thickness \(\tau\) and the single scattering albedo \(\omega\) as:

\[\tau(S)=\int_{S^{\prime}=0}^{S}\beta(S^{\prime})dS^{\prime} =\int_{S^{\prime}=0}^{S}[\kappa(S^{\prime})+\sigma(S^{\prime})]dS^{\prime} \tag{2}\] \[\omega=\frac{\sigma}{\beta}=\frac{\sigma}{\kappa+\sigma} \tag{3}\]

one can rewrite Eq. 1 in the alternative form:

\[\begin{split}\frac{1}{\beta c}\frac{\partial I(\mathbf{\tau},\hat{ s},t)}{\partial t}&=-\hat{s}\cdot\nabla_{\tau}I(\mathbf{\tau},\hat{s},t)-I( \mathbf{\tau},\hat{s},t)+\\ \omega\int_{4\pi}I(\mathbf{\tau},\hat{s}^{\prime},t)p(\hat{s}^{ \prime},\hat{s})d\Omega^{\prime}+S(\mathbf{\tau},\hat{s},t)\end{split} \tag{4}\]

which is the RTE in its standard form. For a detailed discussion and derivation of the RTE, we refer the reader to references [29; 30; 33]. In essence, the RTE states that a directed beam of light in a uniform random medium loses energy through divergence and extinction, including both absorption and scattering away from the beam (i.e., out-scattering contributions), and it gains energy from radiation sources, fluorescence or scattering events that redirect it towards the beam (i.e., in-scattering contributions). In the standard formulation, wave interference effects, polarization and non-linearity in the medium are neglected. Radiative transfer theories for vector waves have also been developed but are outside the scope of this work and more details on these subjects can be found in references [31; 37]. Even for the relevant slab geometry, the RTE introduced above is generally very difficult to solve [38]. Analytic solutions only exist for very simple cases while in many realistic situations, numerical methods such as Monte Carlo transport simulations are usually employed [39]. For this reason, the RTE is often approximated, under suitable conditions, by the simpler but less accurate diffusion equation [32].

In our paper, we developed APINNs to obtain the forward and inverse solution of the scalar RTE in the steady-state and for different choices of phase functions. However, the developed framework can be naturally extended to time-dependent and vector RTE problems, anisotropic phase functions, and arbitrary nonlinear responses. All the implementations of the APINN algorithms developed in this paper are obtained in the powerful TensorFlow environment [40].

The general APINN network utilized to solve forward and inverse RTE problems in the slab geometry is illustrated in Fig. 1. We considered a fully connected neural network

Figure 1: Schematics of the APINN solving the RTE problem in a slab. The FCNN has \(N+1\) outputs where \(N\) is the number of auxiliary variables \(v_{i}\) required in the Legendre polynomial expansion of the scattering phase function. The outputs of the network are then used to satisfy the PDE, initial conditions (ICs), and boundary conditions (BCs) of the differential equation, which are then combined into the loss function \(\mathcal{L}\). During the training process, the loss function is minimized until it passes a threshold \(\sigma\).
(FCNN) with input vector \(\mathbf{x}=(\mathbf{\tau},\mathbf{\mu})\) with randomly distributed values of the optical thickness \(\tau\) and \(\mu=\cos\theta\) over a two-dimensional spatial-angular domain \(\Omega\) and output that is the predicted surrogate \(\hat{I}(\mathbf{\mu},\mathbf{\tau};\mathbf{\tilde{\theta}})\) of the RTE solution \(I(\mathbf{\mu},\mathbf{\tau};\mathbf{\theta})\). Here, \(\mathbf{\theta}\) denotes the angle of the directed energy flow with the axis \(z\) perpendicular slab's surface and \(\mathbf{\tilde{\theta}}\) is the vector of weights and biases of our FCNN. In addition, the FCNN outputs \(n\) auxiliary variables \(v_{i}(\mathbf{\mu},\mathbf{\tau};\mathbf{\tilde{\theta}})\), each corresponding to an integral expansion term in the RTE. The outputs of the APINN are then used to compute, by means of automatic differentiation (AD), the derivatives of \(\hat{I}(\mathbf{\mu},\mathbf{\tau};\mathbf{\tilde{\theta}})\) and \(v_{i}(\mathbf{\mu},\mathbf{\tau};\mathbf{\tilde{\theta}})\), along with the PDE, initial conditions, and boundary conditions, depending on the nature of the problem. Each calculated value is then combined into a term of the loss function \(\mathcal{L}(\mathbf{\tilde{\theta}})\) defined as:

\[\begin{split}\mathcal{L}(\mathbf{\tilde{\theta}})=& \mathcal{L}_{int}(\mathbf{\tilde{\theta}};\mathcal{N}_{int})+\mathcal{L}_{b}( \mathbf{\tilde{\theta}};\mathcal{N}_{b})\\ &+\mathcal{L}_{aux}(\mathbf{\tilde{\theta}};\mathcal{N}_{aux})+ \lambda\sum_{i}\mathbf{\tilde{\theta}}_{i}^{2}\end{split} \tag{5}\]

In the expression above,

\[\mathcal{L}_{int}(\mathbf{\tilde{\theta}};\mathcal{N}_{int})=\frac{1}{|\mathcal{N }_{int}|}\sum_{\mathbf{x}\in\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{ \mathcal{\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{ \mathcal{\mathcal{\mathcal{ \mathcal{ \mathcal{       \mathcal{    }}}}}}}}}}}}}} \left|\int\left(\mathbf{x};\mathbf{\hat{I}}} \frac{\partial\mathbf{\hat{I}}{\partial\mathbf{\tau}}}{v_{0},\ldots,v_{n}}\right) \right|\right|^{2} \tag{6}\]

denotes the loss term calculated in the interior domain \(\Omega\) and

\[\mathcal{L}_{b}(\mathbf{\tilde{\theta}};\mathcal{N}_{b})=\frac{1}{|\mathcal{N}_ {b}|}\sum_{\mathbf{x}\in\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{ \mathcal{\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{\mathcal{ \mathcal{\mathcal{\mathcal{\mathcal{ \mathcal{       \mathcal{    }}}}}}}}}}}}}}} \left|\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{ \mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\
which is constrained by the following system:

\[\mu\frac{dI}{d\tau}+I-\frac{\omega}{2}v(1)=0 \tag{11a}\] \[v(\mu;\tau)=\int_{-1}^{\mu}I(\mu^{\prime};\tau)d\mu^{\prime}\] (11b) \[v(-1;\tau)=0,\quad\frac{dv}{d\mu}(\mu;\tau)=I(\tau,\mu)\]

We then train the APINN to solve the problem for different values of the albedo \(\omega\) varying from 0.2 to 1.0. Table 2 shows the speed and accuracy of our APINN implementation in solving the Milne problem. In the large scattering limit of \(\omega\geq 0.9\), APINN minimized the loss function with values that are two orders of magnitude lower and for a fraction of the time than for the equivalent geometry displayed in Ref. [27], where a quadrature method was employed. Two representative APINN solutions for the spatial-angular distributions of the radiation intensity for \(\tau_{\max}=1.0\) are displayed in Fig. 2 (a) and (b). To benchmark our solutions using the tables calculated by Van de Hulst's in Ref. [42], we computed the zeroth moment or point-direction gain \(G(\tau)\) of the radiative intensity, which is defined as [29; 30; 33; 42]:

\[G(\tau)=\int_{-1}^{1}I(\tau,\mu)d\mu \tag{12}\]

Fig. 2 (c) displays the validation data of \(G(\tau)\) calculated by Van de Hulst and the solution from our network, showing an excellent agreement achieved by the APINN framework. This is further confirmed by the average relative error between the two solutions displayed in the last column of Table 2. Fig. 2 (d) shows a comparison between the APINN and the standard PINN quadrature loss function to solve the same problem, as implemented in Ref. [27]. In this figure, we display the loss function versus the number of epochs for the three largest scattering values of \(\omega\). We can immediately notice that the quadrature solution is heavily affected in its performance by the scattering strength, and the L-BFGS-B solver terminates the training early because the loss function has already saturated to its minimum value and is not decreasing further. In contrast, the APINN's loss function monotonically decreases independently of \(\omega\). This result confirms the robustness, flexibility, and accuracy of the APINN framework in solving transport problems for strongly scattering media. In a variety of engineering applications, however, the material's response is not isotropic. Therefore, in Section III.1.2, we employ the APINN framework to solve the RTE in a slab with an anisotropic Rayleigh scattering phase function.

#### ii.1.2 The Rayleigh scattering phase function

The Rayleigh phase function is employed to study anisotropic light scattering processes in various fields, from

Figure 2: (a) and (b) APINN solutions for the Milne equation in a slab, with the single scattering albedo set to \(\omega=1\) and \(\omega=0.2\), respectively. \(\omega\) sets the strength of the scattering term in the PDE. (c) Validation of the APINN solutions using Van de Hulst’s data results for the point direction gain \(G(\tau)\). (d) Comparison of the APINN training performance in minimizing the objective loss function with the quadrature solution: the APINN method consistently trains independently of scattering strength, whereas the quadrature method’s performance worsens as the scattering strength, represented by \(\omega\), increases.
optics to astronomy [33]. The phase function reads:

\[p(\cos\theta)=\frac{3}{4}(1+\cos\theta^{2}) \tag{13}\]

and because the scattering from spherically symmetric particles is cylindrically symmetric with respect to the incoming direction, this symmetry holds after averaging over all possible orientations. Therefore, in these situations, the phase function depends on \(\phi-\phi^{\prime}\) and one can compute this average resulting in the projected phase function [41]:

\[p_{0}(\mu,\mu^{\prime})=\int\frac{d\phi}{2\pi}\frac{d\phi^{\prime}}{2\pi}p(\mu, \phi;\mu^{\prime},\phi^{\prime}) \tag{14}\]

Using the equality \(\mu=\cos\Theta=\mathbf{n}\cdot\mathbf{n}^{\prime}=\sin\!\theta\!\sin\!\theta^{\prime} \!\cos(\phi-\phi^{\prime})+\!\cos\!\theta\!\cos\!\theta^{\prime}\) one obtains:

\[p_{0}(\mu,\mu^{\prime})=\frac{3}{8}(3-\mu^{2}-\mu^{\prime 2}+3\mu^{2}\mu^{ \prime 2}) \tag{15}\]

To facilitate the calculations and the auxiliary variable formulation of the APINN framework, one typically considers the expansion of the scattering phase function in Legendre polynomials:

\[\Phi(\mu,\mu^{\prime})=\sum_{\ell=0}^{\infty}w_{\ell}P_{\ell}(\mu)P_{\ell}( \mu^{\prime}) \tag{16}\]

Note that, for the Rayleigh phase function, the only nonzero \(w_{\ell}\) terms are \(w_{0}=1.0\) and \(w_{2}=0.1\). Therefore, Eq. 9 in a slab with Rayleigh scattering becomes

\[\mu\frac{dI(\tau,\mu)}{d\tau}+I(\tau,\mu)=\frac{\omega}{2}\int_{-1}^{1}I(\tau,\mu^{\prime})\sum_{\ell=0}^{\infty}w_{\ell}P_{\ell}(\mu)P_{\ell}(\mu^{\prime} )d\mu^{\prime} \tag{17}\]

and after rearranging terms and truncating the series expansion at \(\ell=2\) we get:

\[\mu\frac{dI(\tau,\mu)}{d\tau}+I(\tau,\mu)=\frac{\omega}{2}\Big{[} w_{0}P_{0}(\mu)\int_{-1}^{1}I(\tau,\mu^{\prime})P_{0}(\mu^{\prime})d \mu^{\prime}\] \[+w_{2}P_{2}(\mu)\int_{-1}^{1}I(\tau,\mu^{\prime})P_{2}(\mu^{\prime })d\mu^{\prime}\Big{]} \tag{18}\]

Finally, we recast the problem by adding two auxiliary variables to the network with their respective constraints as follows:

\[\mu\frac{dI(\tau,\mu)}{d\tau}+I(\tau,\mu)=\frac{\omega}{2}\left[w_{0}P_{0}( \mu)v_{0}(1)+w_{2}P_{2}(\mu)v_{2}(1)\right] \tag{19a}\] \[v_{0}(\mu;\tau)=\int_{-1}^{\mu}I(\tau,\mu^{\prime})P_{0}(\mu^{ \prime})d\mu^{\prime}\] (19b) \[v_{0}(-1;\tau)=0,\quad\frac{dv_{0}}{d\mu}(\mu;\tau)=I(\tau,\mu)P_{0}(\mu)\] \[v_{2}(\mu;\tau)=\int_{-1}^{\mu}I(\tau,\mu^{\prime})P_{2}(\mu^{ \prime})d\mu^{\prime}\] (19c) \[v_{2}(-1;\tau)=0,\quad\frac{dv_{2}}{d\mu}(\mu;\tau)=I(\tau,\mu)P_{2}(\mu)\]

Due to the lack of benchmark solutions for Rayleigh scattering in a slab, we decided to consider a physical system similar to the one studied by Mishra and Molinaro in Ref. [27], namely the case where the single scattering albedo depends on the optical thickness \(\tau\) of the material. In this case, Eq. 19a becomes

\[\frac{dI(\tau,\mu)}{d\tau}+I(\tau,\mu)=\frac{\omega(\tau)}{2} \Big{[} w_{0}P_{0}(\mu)v_{0}(1;\tau) \tag{20}\] \[+w_{2}P_{2}(\mu)v_{2}(1;\tau)\Big{]}\]

To solve this problem, we train APINN with the parameters specified in Table 1, using 40 neurons per layer. The training for this solution took 12 minutes, and the final value of the loss function \(\mathcal{L}\) was \(10^{-6}\), demonstrating the adaptivity and flexibility of APINN in solving anisotropic scattering problems. Fig. 3(a) displays the APINN radiative intensity solution as a function of \(\mu\) and the optical thickness. This result highlights the flexibility of APINN in finding the solution to an analytically intractable problem [27]. In turn, this motivates us to study the RTE with strongly anisotropic scattering properties modeled by the Henyey-Greenstein (HG) phase function.

#### ii.2.3 The Henyey-Greenstein phase function

Here we consider the forward RTE problem in the slab with the Henyey-Greenstein (HG) phase function governing the scattering processes. The HG phase function finds applications in astrophysics, atmospheric optics, and biomedical imaging, and it depends on both the cosine of the incident angle and the asymmetry factor \(g\in[0,1]\) that appears in the equation below [30; 43; 44; 29; 45]:

\[p(\mu,g)=\frac{1-g^{2}}{(1-2g\mu+g^{2})^{3/2}} \tag{21}\]

where \(\mu=\cos\!\theta\). In the limit of \(g\to 0\), the HG phase function reduces to isotropic scattering, while in the limit of \(g\to 1\), HG describes strongly anisotropic scattering events.

As for the Rayleigh phase function, the HG phase function can be rewritten using the Legendre polynomials expansion in Eq. (16). However, unlike the Rayleigh case, the Legendre

\begin{table}
\begin{tabular}{c c c c} \(\omega\) & Training time & Loss & Avg. Rel. Error \\ \hline
0.2 & 3 min & \(7\times 10^{-6}\) & \(7\times 10^{-2}\) \\
0.6 & 10 min & \(5\times 10^{-6}\) & \(8\times 10^{-2}\) \\
0.9 & 10 min & \(3\times 10^{-6}\) & \(7\times 10^{-2}\) \\
0.99 & 10 min & \(3\times 10^{-6}\) & \(7\times 10^{-2}\) \\
1.0 & 10 min & \(3\times 10^{-6}\) & \(7\times 10^{-2}\) \\ \end{tabular}
\end{table}
Table 2: APINN training information for the Milne equation in a slab. The last column shows the average relative error between the network solution of \(G(\tau)\) and Van de Hulst’s results.
expansion converges more slowly, and additional terms need to be included to achieve accurate numerical results:

\[\begin{split}\mu\frac{dI(\tau,\mu)}{d\tau}+I(\tau,\mu)=\frac{\omega}{ 2}\Big{[}& w_{0}P_{0}(\mu)\int_{-1}^{1}I(\tau,\mu^{\prime})P_{0}(\mu^{ \prime})d\mu^{\prime}\\ &+w_{1}P_{1}(\mu)\int_{-1}^{1}I(\tau,\mu^{\prime})P_{1}(\mu^{ \prime})d\mu^{\prime}\\ &+\ldots\\ &+w_{n}P_{n}(\mu)\int_{-1}^{1}I(\tau,\mu^{\prime})P_{n}(\mu^{ \prime})d\mu^{\prime}\Big{]}\end{split} \tag{22}\]

where:

\[w_{\ell}=(2n+1)g^{n} \tag{23}\]

In our numerical studies, we chose to benchmark the RTE with HG phase function, \(g=0.5\), which allowed us to utilize Van de Hulst's tables as validation data [46]. The polynomial expansion of the phase function was truncated after ten terms, introducing ten auxiliary variables and their corresponding constraint conditions in the simulation:

\[\begin{split}&\mu\frac{dI(\tau,\mu)}{d\tau}+I(\tau,\mu)=\frac{ \omega}{2}\Big{[}w_{0}P_{0}(\mu)v_{0}(1;\tau)\\ &+\cdots+w_{10}P_{10}(\mu)v_{10}(1;\tau)\Big{]}\\ & v_{0}(\mu;\tau)=\int_{-1}^{\mu}I(\tau,\mu^{\prime})P_{0}(\mu^{ \prime})d\mu^{\prime}\\ & v_{0}(-1;\tau)=0,\quad\frac{dv_{0}}{d\mu}(\mu;\tau)=I(\tau,\mu) P_{0}(\mu)\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \qquad\ldots\\ & v_{10}(\mu;\tau)=\int_{-1}^{\mu}I(\tau,\mu^{\prime})P_{10}(\mu^{ \prime})d\mu^{\prime}\\ & v_{10}(-1;\tau)=0,\quad\frac{dv_{10}}{d\mu}(\mu;\tau)=I(\tau, \mu)P_{10}(\mu)\end{split} \tag{24c}\]

Table 3 provides a summary of the APINN training for this problem. Considering the larger number of auxiliary variables, we trained using 80 neurons per layer instead of 40. Similarly to the isotropic and Rayleigh cases, the loss function is minimized to extremely low values with a minor trade-off in speed due to the larger number of auxiliary variables in the system, as the second and third columns of Table 3 demonstrate. The accuracy of these results, displayed in the last column of Table 3, confirms the versatility of the APINN framework, which excels in solving even strong anisotropic scattering problems. Fig. 3 (c) shows a representative solution of the radiation intensity when \(\omega=1.0\), and Fig. 3 (b) displays the benchmarked solutions for this problem by comparing the integrated radiative intensity \(G(\tau)\) calculated from the APINN network with the Van de Hulst's data. These results open the doors for multiple biomedical, metamaterials, and nano-optics applications where the HG phase function is often utilized to model realistic scattering processes [43, 44, 45].

### The coupled radiation-conduction problem of a participating medium

We now apply our APINN method to the solution of a coupled problem in radiative transfer theory. Specifically, we con

Figure 3: (a) APINN solution for the RTE with Rayleigh scattering phase function, where the single scattering albedo \(\omega(\tau)\) depends on the optical thickness. (b) Validation of the APINN solutions for the HG phase function with \(g=0.5\) using Van de Hulst’s data results for the point direction gain \(G(\tau)\). (c) Representative APINN solution for the RTE in a slab with HG scattering phase function, \(g=0.5\), \(\omega=1.0\).

\begin{table}
\begin{tabular}{c c c c} \(\omega\) & Training time & Loss & Avg. Rel. Error \\ \hline
0.2 & 49 min & \(8\times 10^{-6}\) & \(3\times 10^{-2}\) \\
0.6 & 50 min & \(2\times 10^{-5}\) & \(2\times 10^{-2}\) \\
0.9 & 65 min & \(3\times 10^{-5}\) & \(6\times 10^{-3}\) \\
0.99 & 65 min & \(6\times 10^{-5}\) & \(2\times 10^{-3}\) \\
1.0 & 65 min & \(5\times 10^{-5}\) & \(1\times 10^{-3}\) \\ \end{tabular}
\end{table}
Table 3: APINN training information for the HG phase function in a slab, \(g=0.5\). The last column shows the average relative error between the network solution of \(G(\tau)\) and Van de Hulst’s table.
sider a conducting and participating slab that couples to the radiation hitting the boundary in the steady-state. Such problems have been analyzed extensively in the literature [47; 48; 49; 50; 51; 52; 53; 54; 55; 56], but, to our knowledge, have never been solved using physics-informed neural networks. Here, we use the APINN framework to analyze this problem, where the slab's temperature profile is governed by a Poisson-like equation with a coupling term to the RTE [51]. We will further analyze how the conduction-radiation parameter \(N_{CR}\) affects the traditional Fourier temperature solution in the steady-state when significant temperature differences are imposed at the boundaries of the slab. The conduction-radiation parameter \(N_{CR}\) measures the ratio of conductive to radiative heat contributions in a given medium, and it is defined as [29]:

\[N_{CR}=\frac{k\beta}{4k_{B}T^{3}}=\frac{k(\kappa+\sigma)}{4k_{B}T^{3}} \tag{25}\]

For the simulations that follow, we chose to study coupled systems where \(N_{CR}\) varies from 10 (for \(N\rightarrow\infty\), we get the Fourier limit) to 0.001 (for \(N\to 0\), radiative processes dominate).

We consider the heat transfer problem due to conduction and radiation in a participating medium presented by Ref. [51] governed by the two following coupled integro-differential equations:

\[\frac{d^{2}\Theta}{d\tau^{2}}-\frac{(1-\omega)}{N_{CR}}\left[ \Theta^{4}(\tau)-\frac{1}{2}G(\tau)\right]=0 \tag{26}\] \[\mu\frac{dI(\tau,\mu)}{d\tau}+I(\tau,\mu)=H[\Theta(\tau)]+\frac{ \omega}{2}\int_{-1}^{1}I(\tau,\mu^{\prime})\Phi(\mu,\mu^{\prime})d\mu^{\prime}\] (27) \[\text{for}\quad 0<\tau<1,\quad-1\leq\mu\leq 1,\quad\Phi(\mu,\mu^{ \prime})=1,\quad\omega=0.9\]

where the temperature is being modeled by the normalized adimensional quantity \(\Theta=T/T_{1}\). The coupling terms are:

\[G(\tau)=\int_{-1}^{1}I(\tau,\mu)d\mu,\quad H[\Theta(\tau)]=(1-\omega)\Theta^{4} \tag{28}\]

and \(G(\tau)\) is the zeroth moment of the intensity \(I(\tau,\mu)\). The boundary conditions are:

\[I(0,\mu)=1,\mu\in(0,1],\quad\text{and}\quad I(1,\mu)=0,\mu\in[-1,0) \tag{29}\] \[\Theta(0)=1\quad\text{and}\quad\Theta(1)=T_{2}/T_{1} \tag{30}\]

Since the problem involves two undetermined coupled functions, we modified the architecture of the APINN framework. The changes are illustrated in Fig. 4: the input parameters are passed to the radiative intensity network \(\hat{I}(\tau,\mu)\) with auxiliary variables as for the uncoupled cases discussed so far, but the spatial variable \(\tau\) is also used to train simultaneously the adimensional temperature network \(\hat{\Theta}(\tau)\). The coupled problem recasted in the APINN formalism reads:

\[\frac{d^{2}\Theta}{d\tau^{2}}-\frac{(1-\omega)}{N_{CR}}\left[ \Theta^{4}(\tau)-\frac{1}{2}v(1;\tau)\right]=0 \tag{31}\] \[\mu\frac{dI(\tau,\mu)}{d\tau}+I(\tau,\mu)=H[\Theta(\tau)]+\frac{ \omega}{2}v(1;\tau) \tag{32}\]

where we introduced the auxiliary variable \(v(\mu;\tau)\) and its corresponding conditions like in Eq. (11):

\[v(\mu;\tau)=\int_{-1}^{\mu}I(\mu^{\prime};\tau)d\mu^{\prime}, \tag{33a}\] \[v(-1;\tau)=0,\quad\frac{dv}{d\mu}(\mu;\tau)=I(\tau,\mu) \tag{33b}\]

By means of automatic differentiation, the outputs of the two networks are then used to compute the required PDE conditions, initial conditions, and boundary conditions, which are then incorporated into the coupled loss function:

\[\mathcal{L}=\mathcal{L}_{I(\tau,\mu)}+\mathcal{L}_{\hat{\Theta}(\tau)} \tag{34}\]

To solve this problem, we coupled the APINN network for the radiative intensity with a PINN estimating the dimensionless

Figure 4: Schematics of the APINN solving the coupled radiation-conduction problem in a slab. Here, two FCNNs representing the intensity of radiation \(\hat{I}(\mu,\tau)\) and temperature \(\hat{\Theta}(\tau)\) are combined in a loss function that minimizes the PDE, initial conditions (ICs), and boundary conditions (BCs) of two coupled partial differential equations. The training process minimizes the loss function \(\mathcal{L}\) until it passes a threshold \(\sigma\).

\begin{table}
\begin{tabular}{c c c c} \(\Delta\Theta\) & \(N_{CR}\) & Training time & \(\mathcal{L}_{I(\tau,\mu)}\) & \(\mathcal{L}_{\hat{\Theta}(\tau)}\) \\ \hline \multirow{4}{*}{150K} & 10 & 18 min & \(6\times 10^{-4}\) & \(6\times 10^{-6}\) \\  & 0.1 & 18 min & \(5\times 10^{-4}\) & \(5\times 10^{-7}\) \\  & 0.01 & 18 min & \(5\times 10^{-4}\) & \(7\times 10^{-6}\) \\  & 0.001 & 18 min & \(5\times 10^{-4}\) & \(8\times 10^{-5}\) \\ \hline \multirow{4}{*}{270K} & 10 & 33 min & \(5\times 10^{-4}\) & \(2\times 10^{-5}\) \\  & 0.1 & 30 min & \(5\times 10^{-4}\) & \(2\times 10^{-6}\) \\ \cline{1-1}  & 0.01 & 31 min & \(5\times 10^{-4}\) & \(6\times 10^{-6}\) \\ \cline{1-1}  & 0.001 & 30 min & \(5\times 10^{-4}\) & \(4\times 10^{-5}\) \\ \end{tabular}
\end{table}
Table 4: APINN Training information for radiative-conductive coupled problem in a slab.
temperature \(\Theta(\tau)\), with parameters according to Table 1. Fig. 5 shows the solutions for the coupled problem when two different temperature jumps are imposed at the rightmost boundary. Fig. 5(a) displays a \(\Delta\Theta=150K\), whereas Fig. 5(b) a \(\Delta\Theta=270K\). Moreover, we analyze the dimensionless temperature behavior when the conduction-radiation parameter \(N_{CR}\) decreases, as previously investigated in Refs. [29; 30]. It is important to realize that both panels in Fig. 5 display a beyond-Fourier behavior as \(N_{CR}\) decreases, demonstrating that the temperature profile is significantly affected by radiative scattering phenomena. Lastly, Table 4 presents some relevant information regarding the APINN training. We note that, even for the coupled case, the APINN successfully minimizes both the temperature loss function \(\mathcal{L}_{\hat{\Theta}(\tau)}\) and the radiative intensity loss function \(\mathcal{L}_{\hat{\Theta}(\tau)}\) independently of the parameter \(N_{CR}\).

### Inverse problem: retrieval of the albedo from the boundary data

Finally, we present here the solution of an inverse problem of radiative transfer theory where we employ APINN to retrieve simultaneously the forward solution of the intensity \(I(\tau,\mu)\) and the single scattering albedo \(\omega\). We do not, however, introduce synthetic data everywhere in the domain, as it has been done previously in the literature [27], but we limit ourselves to introducing two data points representing the integrated intensity \(G(\tau)\) at the edges of the slab, simulating a lab environment with two detectors capturing integrated radiation entering and exiting the slab, respectively. The reason to present an inverse problem in such a fashion is to demonstrate the full potential and capabilities of physics-informed neural networks that, with no additional overhead and computing power, can solve a forward and parameter retrieval problem simultaneously. We thus modify the Schwarzschild-Milne equation for a slab discussed in an earlier section. In particular, Eq. (11) is changed to include the unknown albedo parameter \(\omega_{\theta}\):

\[\mu\frac{dI}{d\tau}+I-\frac{\omega_{\theta}}{2}v(1)=0 \tag{35}\]

and the loss function in Eq. (5) is modified to include the two synthetic detector data points at the boundaries of the slab:

\[\begin{split}\mathcal{L}(\mathbf{\theta},\mathbf{\omega}_{\theta})& =\mathcal{L}_{int}(\mathbf{\theta},\mathbf{\omega}_{\theta};\mathcal{N}_{ int})+\mathcal{L}_{b}(\mathbf{\theta},\mathbf{\omega}_{\theta};\mathcal{N}_{b})\\ &+\mathcal{L}_{aux}(\mathbf{\theta},\mathbf{\omega}_{\theta};\mathcal{ N}_{aux})+\mathcal{L}_{inv}(\mathbf{\theta},\mathbf{\omega}_{\theta};\mathcal{N}_{im}) \end{split} \tag{36}\]

where

\[\begin{split}&\mathcal{L}_{inv}(\mathbf{\theta},\mathbf{\omega}_{\theta };\mathcal{N}_{inv})=\\ &\frac{1}{|\mathcal{N}_{im}|}\sum_{(\tau,\mu)\in\mathcal{N}_{im} }\left|\left|\int_{-1}^{1}\hat{I}(\tau,\mu)d\mu-G(\tau)\right|\right|^{2}=\\ &\frac{1}{2}\left(\left|\left|\int_{-1}^{1}\hat{I}(0,\mu)d\mu-G (0)\right|\right|^{2}+\left|\left|\int_{-1}^{1}\hat{I}(1,\mu)d\mu-G(1)\right| \right|^{2}\right)\end{split} \tag{37}\]

Fig. 6 displays the fast convergence of the retrieved APINN parameter \(\omega_{\theta}\) to the actual value \(\omega\). Each line corresponds to a different APINN training procedure during which the only data points added were \(G(0)\) and \(G(1)\) obtained from the Van de Hulst's tables [42], and used to minimize \(\mathcal{L}_{inv}\) during the training process. Despite the loss term with two data points was not weighted differently from the interior or boundary ones, APINN achieved a precise inversion of the parameter of interest. In fact, as displayed in Table 5, the loss function converges independently of the albedo \(\omega\) and with great precision, as displayed by the relative error between the known

Figure 5: (a) and (b) APINN results displaying the temperature profile in the conductive participating medium of Eq. 33 for a temperature jump of 150K and 270K, respectively. The simulations have been conducted for four sets of values of the conduction-radiation parameter \(N_{CR}\): as \(N_{CR}\) decreases, the solution changes from the well-known Fourier temperature profile to a beyond-Fourier solution characterized by a dominant scattering contribution.

Figure 6: Results showing the convergence of each value of different trainable variables \(\omega_{\theta}\) for separate APINN simulations, each of which where fed different boundary values \(G(\tau)\) for the inverse problem, while all the other parameters and PDE conditions stayed the same. As the picture shows, the neural network provided very accurate results across all simulations, training \(\omega_{\theta}\) to the true albedo value \(\omega\), especially for high scattering values.
albedo and the predicted APINN albedo \(\omega_{\theta}\) shown in the last column of Table 5. Therefore, APINN retrieved the correct parameter of interest \(\omega_{\theta}\) when only two points were added during the training process.

## IV Conclusions

Throughout this paper, we have described different applications of APINN for solving the radiative transfer equation, which is a Boltzmann-type transport equation. We successfully solved forward problems in a slab with both isotropic and anisotropic scattering phase functions and irrespective of the albedo. The results presented improved upon previous attempts to use physics-informed neural networks for solving the RTE in both accuracy and speed [27]. Furthermore, we presented the solution of the first coupled radiation-conduction problem in a participating medium using the APINN framework and we showed that the loss functions of coupled neural networks quickly converged to a low minimum value below \(10^{-5}\). Our findings open the possibility to utilize APINN to analyze higher dimensional systems and discover more interesting physics with applications to metamaterials and semiconductor device modeling. Finally, we solved an inverse problem following a setup that replicates an experimental setting with data points at the boundary of the system. It will be interesting in future studies to build on the APINN platform to address higher dimensional coupled, inverse coupled, and strongly scattering forward systems with applications to biomedical imaging, nanophotonics, metamaterials, and thermal modeling of semiconductor devices.

###### Acknowledgements.

 We acknowledge the support from the U.S. Army Research Office, RF-Center managed by Dr. J. Qiu (Grant #W911NF-22-2-0158). We thank professors Mike Kirby, Akil Narayan, and Shandian Zhe for useful discussions on this topic.
