MNRAS 000

, 1-10 (2024)

## Approximating Rayleigh Scattering in Exoplanetary Atmospheres using Physics-informed Neural Networks (PINNs)

David Dahlbüdding 1 2 3 , , ★ , Karan Molaverdikhani 1 2 , Barbara Ercolano 1 2 , , and Tommaso Grassi 2 3 ,

1 Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679 München, Germany

- 2 Exzellenzcluster 'Origins', Boltzmannstr. 2, D-85748 Garching, Germany

3 Centre for Astrochemical Studies, Max-Planck-Institut für extraterrestrische Physik, Gießenbachstr. 1, D-85749 Garching, Germany

Accepted 2024 July 30. Received 2024 July 8; in original form 2024 March 6

## ABSTRACT

This research introduces an innovative application of physics-informed neural networks (PINNs) to tackle the intricate challenges of radiative transfer (RT) modeling in exoplanetary atmospheres, with a special focus on efficiently handling scattering phenomena. Traditional RT models often simplify scattering as absorption, leading to inaccuracies. Our approach utilizes PINNs, noted for their ability to incorporate the governing differential equations of RT directly into their loss function, thus offering a more precise yet potentially fast modeling technique. The core of our method involves the development of a parameterized PINN tailored for a modified RT equation, enhancing its adaptability to various atmospheric scenarios. We focus on RT in transiting exoplanet atmospheres using a simplified 1D isothermal model with pressure-dependent coefficients for absorption and Rayleigh scattering. In scenarios of pure absorption, the PINN demonstrates its effectiveness in predicting transmission spectra for diverse absorption profiles. For Rayleigh scattering, the network successfully computes the RT equation, addressing both direct and diffuse stellar light components. While our preliminary results with simplified models are promising, indicating the potential of PINNs in improving RT calculations, we acknowledge the errors stemming from our approximations as well as the challenges in applying this technique to more complex atmospheric conditions. Specifically, extending our approach to atmospheres with intricate temperature-pressure profiles and varying scattering properties, such as those introduced by clouds and hazes, remains a significant area for future development.

Key words: radiative transfer - scattering - planets and satellites: atmospheres - methods: numerical

## 1 INTRODUCTION

ral posterior estimation (Vasist et al. 2023), or flow matching and importance sampling (Gebhard et al. 2023).

The use of Machine Learning (ML) in exoplanetary science has grown increasingly popular during past years (Nixon &amp; Madhusudhan 2020). Hereby, it is mostly used to accelerate the time-consuming process of finding the posterior distributions of atmospheric parameters given some observed spectrum. Usually, this inference is done using a Bayesian scheme, like nested sampling (e.g. Skilling 2004) where many so-called forward models simulate data which can then be compared to observations.

Pioneered by Waldmann (2016), who used a simple multi-layer perceptron (MLP) to predict the existence of molecules in exoplanetary atmospheres from their transit spectra, the approaches have diversified since then, from Random Forests (Márquez-Neila et al. 2018; Fisher et al. 2020; Nixon &amp; Madhusudhan 2020) to generative adversarial networks like ExoGAN (Zingales &amp; Waldmann 2018) and convolutional neural networks, often in combination with Monte Carlo dropout (Soboczenski et al. 2018; Cobb et al. 2019; Ardévol Martínez et al. 2022).

More recent approaches achieve generating more accurate and versatile posterior distributions by employing normalizing flows in combination with variational inference (Hou Yip et al. 2022), neu-

★ E-Mail: ddahlb@mpe.mpg.de

The approaches most similar to our own are by Himes et al. (2022) and Hendrix et al. (2023), who both try to speed up the forward model while leaving the inference scheme untouched. While Himes et al. (2022) use a neural network to emulate radiative transfer by predicting the spectrum given some atmospheric parameters, Hendrix et al. (2023) focus on accelerating computationally expensive disequilibrium chemistry calculations using an autoencoder and an LSTM-like neural network.

Most of these models rely on traditional numerical methods to generate the data that the ML model uses for training. Since all ML models cannot contain more information than their training data, they cannot generalize well beyond the used forward model. Including the time to train a model, the respective method only offers a speed-up if enough exoplanets can be analyzed using the same model trained on the same data. This can be particularly challenging for exoplanetary atmospheres, given the diversity of exoplanets that have been observed. Additionally, as soon as a forward model is updated to include some new insight, the ML model might be completely retrained. It can also prove difficult to use one model for observations with varying instruments, that, e.g., have a different number of data points or an incompatible spectral range.

In the future, as data become more abundant and precise and the

models more detailed, the computational burden of running a retrieval will further increase. Therefore, in this work, we will assess the feasibility of using pre-trained physics-informed neural networks (PINNs) to cheaply model complex phenomena in exoplanetary atmospheres. We focus on the simple example of Rayleigh scattering in the isothermal atmosphere of a transiting exoplanet.

Building upon the work of Mishra &amp; Molinaro (2021), who used PINNs to solve the radiative transfer (RT) equation, we first model RT considering only absorption. Since for transiting exoplanets, the traditional numerical method treats scattering the same as absorption in that the light scattered out of the beam is analogous to absorbed light (e.g. Mollière et al. 2019), we can compare our PINN to existing models to establish its speed and accuracy.

Secondly, we include Rayleigh scattering in the RT equation. As has been shown by Sengupta et al. (2020), scattered light can have a significant effect on the transit spectrum at short wavelengths. Using PINNs could help us study this effect in more detail and potentially improve our estimates of the hydrogen-to-helium ratio or the existence and properties of clouds and hazes.

In Section 2, after briefly introducing PINNs, we discuss how the respective problem is presented to both the absorption and scattering PINN, such that they are able to solve it efficiently. This involves generating and normalizing the inputs, the formulation of the loss function, and the exact training procedure. The results of both models are assessed in Section 3, while drawbacks are summarized in Section 4. Section 5 discusses key insights and possible enhancements.

## 2 METHODS

## 2.1 PINNs

While there are many different ways of incorporating physical laws into the training of a neural network, PINNs, as introduced by Raissi et al. (2019), achieve this through their loss function, which is directly based on a differential equation.

To see how a PINN is able to solve a nonlinear partial differential equation (PDE), as explained in Cuomo et al. (2022), let us first describe a PDE in its most general form defined on the domain Ω ∈ R d with its boundary 𝜕 Ω :

$$\mathcal { F } ( u ( z ) ; \gamma ) = f ( z ), \ \ z \in \Omega$$

$$\mathcal { B } ( u ( z ) ) & = g ( z ), \quad z \in \partial \Omega \,, & ( 2 ) \quad 2. 2 \, \text{Se}$$

where 𝑢 describes the unknown solution to the differential equation, dependent on the space-time coordinate vector 𝑧 : = [ 𝑥 1 , · · · , 𝑥 d -1 , 𝑡 ] . F refers to the nonlinear differential operator dependent on physics-related parameters 𝛾 , while 𝑓 is the function identifying the data of the problem (if available). Eq. (2) defines the boundary and initial conditions, where B is again an operator and 𝑔 is the boundary function.

While the PINN can also be used for an inverse problem of finding the parameters 𝛾 , we use it for approximating the solution 𝑢 for every 𝑧 , dependent on the (given) parameters 𝛾 . Hereby, the PINN itself functions as the approximation to this solution. The PINN therefore has inputs 𝑧 and output ˆ 𝑢 𝜃 ( 𝑧 ) ≈ 𝑢 𝑧 ( ) dependent on the network parameters 𝜃 . The corresponding loss function that should be minimized is made up of different components L 𝑖 with respective weights 𝜔 𝑖 :

$$\mathcal { L } & = \omega _ { \mathcal { F } } \mathcal { L } _ { \mathcal { F } } + \omega _ { \mathcal { B } } \mathcal { L } _ { \mathcal { B } } \,. & ( 3 ) \quad \Phi ( \Delta ) \, \colon$$

The first term L is the differential equation (1) itself, with all terms on one side of the equation such that it is equal to zero. Often, the following form analogous to the mean squared error is used:

$$o \, \underset { \text{model} } { \text{$\mathcal{L}_{\mathcal{F}=\sum_{i=1}^{n}(\mathcal{F}(\hat{\mathfrak{a}_{\theta}(z_{i})':\gamma) -f(z_{i}))\mathcal{$}^{2}$} } } \,, \quad z _ { i } \in \Omega \,, \quad \quad ( 4 )$$

where the square guarantees both positive values, such that the minimumliesatzero,anddifferentiability (compared to the mean absolute error).

Nowwecalculate the value of the differential equation F( 𝑢 𝑧 ( ) ; 𝛾 ) with automatic differentiation. We can not only use it to calculate the gradient of the loss function with respect to the weights of the network, but also calculate the gradient of the output ˆ 𝑢 𝜃 ( 𝑧 ) of the network with respect to the inputs 𝑧 , the space and time coordinates. Hence, we can insert all the values of the variables and derivatives into the differential equation to calculate the loss.

The second term can be calculated accordingly to make sure that the boundary conditions are fulfilled:

$$\underset { \text{iection} } { \mathrm e } \, \text{less} \quad \mathcal { L } _ { \mathcal { B } } = \sum _ { i = 1 } ^ { n } \left ( \mathcal { B } ( \hat { u } _ { \theta } ( z _ { i } ) ) - g ( z _ { i } ) \right ) ^ { 2 }, \quad z _ { i } \in \partial \Omega \,. \quad \quad \, ( 5 )$$

In this work, our approach to employing PINNs deviates from their conventional application. Typically, PINNs are tasked with solving specific instances of differential equations, where their advantage over traditional numerical methods is most pronounced in scenarios involving complex geometries or boundary conditions. However, we leverage PINNs in a novel way by training them on a parameterized form of the radiative transfer equation. This means that the solution approximation provided by the PINN, denoted as ˆ 𝑢 𝜃 ( 𝑧 ) , is influenced not only by the space-time coordinates 𝑧 but also by a set of physical parameters 𝛾 . This unique application allows the parameterized PINNto generate solutions, ˆ 𝑢 𝜃 ( 𝑧, 𝛾 ) , for varying physical conditions without the need for retraining, showcasing a flexible and powerful tool for modeling radiative transfer processes in exoplanetary atmospheres.

## 2.2 Setup

In this study, our objective extends beyond solving a PDE; we aim to tackle an integro-differential equation, specifically the stationary radiative transfer equation, which operates in two spatial dimensions:

$$\text{inc} \quad & \text{inc} ( \phi ) \frac { d u } { d x } + \sin ( \phi ) \frac { d u } { d y } + ( \alpha _ { a } + \alpha _ { s } ) \cdot u \\ \text{and} \ g \quad & - \frac { \alpha _ { s } } { 4 \pi } \cdot \int \Phi ( \theta, \phi, \theta ^ { \prime }, \phi ^ { \prime } ) u ( \theta ^ { \prime }, \phi ^ { \prime } ) d \Omega ^ { \prime } = 0 \,,$$

where the angle 𝜙 = 0 is defined as the positive x-direction and lies in the range 𝜙 ∈ [-𝜋, 𝜋 ] . 𝛼 𝑎 and 𝛼 𝑠 are the absorption and scattering coefficient, respectively and Φ the scattering phase function, depending on the incoming ( 𝜃 , 𝜙 ′ ′ ) and outgoing angles ( 𝜃, 𝜙 ) of the light. In the case of Rayleigh scattering, the scattering phase function depends only on the angular distance Δ :

$$\Phi ( \Delta ) = \frac { 3 } { 4 } \left ( 1 + \cos ^ { 2 } \Delta \right ) \,,$$

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00084v1-figure-ba0ebb61-a1d4-4e00-89be-94793094e573.jpg)

Figure 1. Left: The transit geometry. The x (and y) axis was scaled by a factor of 𝑙 𝑥 (and 𝑙 𝑦 ), such that both coordinates are normalized. Right: Example differential optical depth d 𝜏 d 𝑥 = 𝑙 𝑥 𝛼 𝑥, 𝑦 ( ) , where 𝛼 ∼ 𝑃 .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00084v1-figure-f01ef39c-c322-40d5-9745-713a2e49c95e.jpg)

where cos Δ = cos 𝜃 cos 𝜃 ′ + sin 𝜃 sin 𝜃 ′ cos ( 𝜙 -𝜙 ′ ) according to the law of cosines. For both the absorbed and scattered light component sensible symmetries for small angular extents of the star are exploited such that the integral can be evaluated numerically over only one angle (see Appendix B for the explicit integral).

The function of interest that we want to solve for is 𝑢 𝑥, 𝑦, 𝜙 ( ) , the transmission through the atmosphere. For this, we also need to define our geometry, as shown in Fig. 1. More specifically, we want to find the transmission 𝑢 𝑥 ( = 𝑥 max , 𝑦, 𝜙 = 0 , which we can then integrate ) over the 𝑦 or radial direction to find the effective area of the planet.

This depends on a number of other parameters. So, to reduce the dimension of the problem, we assume an isothermal atmosphere and constant conditions throughout the day-night-terminator region of the planet. The models are trained on an atmosphere spanning pressures from 10 -6 to 10 2 bar and consisting of 100 log-spaced layers, where the absorption and scattering coefficient is interpolated logarithmically between them.

Furthermore, for the PINN with scattering, we assume that the planet transits its star with an orbital inclination 𝑖 = 90 , simplifying the ◦ boundary condition. The scattering PINN also only calculates the solution for pure scattering without any absorption ( 𝛼 𝑎 = 0). The total transmission function, when simulating a spectrum, is then calculated by multiplying the solution for purely absorbing species with the solution for purely scattering species, which should hold as long as one of the two is dominating.

## 2.3 Data Generation

Before we are able to train a neural network on this problem, we have to generate and normalize the necessary input data.

For both the absorption and the scattering model, this includes the spacial coordinates 𝑥 and 𝑦 , which are normalized to be within the interval [-1 1 , ] by scaling them with factors 𝑙 𝑥 and 𝑙 𝑦 , respectively. These are defined such that the innermost atmospheric layer touches the lower boundary at ( 𝑥, 𝑦 ) = ( 0 , - ) 1 and the outermost layer touches the upper boundary at ( 𝑥, 𝑦 ) = ( 0 1 , ) as well as going through the lower corners ( 𝑥, 𝑦 ) = (-1 , - ) 1 and ( 𝑥, 𝑦 ) = ( 1 , - ) 1 :

$$l _ { x } & = \sqrt { r _ { 0 } ^ { 2 } - r _ { 9 9 } ^ { 2 } }, \quad l _ { y } & = \frac { r _ { 0 } - r _ { 9 9 } } { 2 } \,, & ( 8 ) & \quad \text{ et al. } \vec { 1 } \\ & \text{profile}$$

where 𝑟 0 refers to the outermost and 𝑟 99 to the innermost layer. Note that we have to account for these scaling factors inside of the differential equation by dividing the respective gradients by this factor.

The other two input parameters which are given to both models include the radius of the planet 𝑅 0 at the reference pressure 𝑃 0 = 10 -2 bar, spanning from 0 2 . 𝑅 J to 2 0 . 𝑅 J , normalized to also be within

[-1 1 . Secondly, the scale height , ] 𝐻 : = 𝑘 B 𝑇 /( 𝜇 𝑔 ) , where T refers to the temperature, 𝜇 to the mean molecular weight, and g is the gravity at the reference pressure 𝑃 0 . log 10 ( 𝐻 [ cm ]) is sampled from a normal distribution with a mean of 7.1 and a standard deviation of 0.35 and is normalized accordingly. Extreme values of the scale height can change the geometry of the problem significantly and hence, make it harder for the PINN to learn a solution. The normal distribution ensures that this does not hinder the training of the PINN, while maintaining a broad range of values in the training data, allowing the neural network to approximate regions of higher data density more quickly and accurately (e.g., Basri et al. 2020). Together, 𝐻 and 𝑅 0 allow us to calculate the radius profile 𝑟 ( 𝑃 ) assuming a constant temperature but variable gravity:

$$\underset { \text{local} } { \text{anning} } \text{ spaced} \quad r ( P ) = R _ { 0 } \cdot \left [ 1 + \frac { H } { R _ { 0 } } \ln \left ( \frac { P } { P _ { 0 } } \right ) \right ] ^ { - 1 } \,.$$

Although isothermality is often an inadequate approximation for exoplanetary atmospheres (e.g., Rocchetto et al. 2016), we use it to assess if the PINN can manage the simplest scenario. That is why we further assume constant mass fractions (or mixing ratios) of all species throughout the atmosphere.

For the pure-absorption PINN, we also need to generate the 𝛼 𝑎 ( 𝑃 ) -profiles, where the absorption coefficient 𝛼 𝑎 is extended to encompass both absorption and scattering sources. In this approach, scattering is effectively treated as a form of absorption, consistent with the traditional model. Ideally, these profiles should not just consist of a database of absorption coefficients from certain molecules at certain wavelengths, but should be more general by generating them randomly on the fly during training. For this, we first investigated example 𝛼 𝑎 ( 𝑃 ) -profiles using petitRADTRANS (Mollière et al. 2019), which calculates the absorption coefficient 𝛼 𝑎 = 𝜅𝑋𝜌 for every species at every layer as the product of opacity 𝜅 , mass fraction 𝑋 and density 𝜌 . Our analysis included the correlated-k opacities of H2 O and CO (Rothman et al. 2010), CH 4 and CO 2 (Chubb et al. 2021) and Na and K (Piskunov et al. 1995), as well as Rayleigh scattering from H 2 (Dalgarno &amp; Williams 1962) and He (Chan &amp; Dalgarno 1965) and collisionally-induced absorption from H2 -H 2 (Borysow et al. 2001; Borysow 2002) and H 2 -He (Borysow et al. 1988, 1989; Borysow &amp; Frommhold 1989). We found that most profiles follow either 𝛼 𝑎 ∼ 𝜌 ∼ 𝑃 or 𝛼 𝑎 ∼ 𝑃 2 (if 𝜅 ∼ 𝑃 ) and hence d ln ( 𝛼 𝑎 ) d ln ( 𝑃 ) either equates to 1 or 2. Therefore, a truncated exponential profile around 1 and 2 ( ≤ 2) was used to create the profiles from their logarithmic gradient. During training, we observed that the PINN improves when using more accurate but diverse profiles. Therefore, we include 6 different types of randomly generated profiles in total:

- · d ln ( 𝛼 𝑎 ) d ln ( 𝑃 ) ≈ 1

Figure 2. Testing the 𝛼 𝑎 ( 𝑃 ) -profile generation algorithm by comparing the resulting distribution of logarithmic gradients between every layer with the same distribution for all opacity sources of the example spectrum shown in Fig. 6. The portion of gradients of the example spectrum that lay outside of the shown interval accumulate to less than 0.1% for d ln ( 𝛼𝑎 ) d ln ( 𝑃 ) &lt; -0 5 and to . about 1.1% for d ln ( 𝛼𝑎 ) d ln ( 𝑃 ) &gt; 3 5. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00084v1-figure-903f269b-e4bb-4056-a945-ccff1ce8dfa3.jpg)

- · d ln ( 𝛼 𝑎 ) d ln ( 𝑃 ) ≈ 2
- · random transition from d ln ( 𝛼 𝑎 ) d ln ( 𝑃 ) ≈ 1 to 2
- · random transition from d ln ( 𝛼 𝑎 ) d ln ( 𝑃 ) ≈ 2 to 1
- · randomly switching between d ln ( 𝛼 𝑎 ) d ln ( 𝑃 ) ≈ 1 and d ln ( 𝛼 𝑎 ) d ln ( 𝑃 ) ≈ 2
- · random negative values ( ≥ -0 5) for . d ln ( 𝛼 𝑎 ) d ln ( 𝑃 ) in some profiles.

In Fig. 2, we compare the distribution of the generated logarithmic gradients with actual ones from all opacity sources and layers of the example atmosphere shown in Fig. 6. The more diverse set of gradient-profiles improved the error of a spectrum from initially around 3-5% to 1-2%.

The final 𝛼 𝑎 ( 𝑃 ) -profiles are then created by cumulatively summing up the gradients in both directions from a middle layer (with index 49), which gets assigned a log 10 ( 𝛼 𝑎 ) -value by a uniform distribution and hence acts as a random shift of the profiles.

To inform the PINN about the profile, we give it all the absorption coefficients 𝛼 , 𝑖 𝑖 ∈ [ 0 99 , ] from the 100 atmospheric layers, such that it has a total of 104 inputs.

In the case of Rayleigh scattering from hydrogen and helium, the scattering cross-section does not depend on the pressure. Hence, the scattering coefficient is simply proportional to the density 𝛼 ∼ 𝜌 ∼ 𝑃 , an example of which can be seen in Fig. 1, and can be parameterized with just one input variable.

The scattering PINN has additional inputs, namely the angle 𝜙 ∈ [-𝜋, 𝜋 ] as well as (half) the angular extent of the star as seen from the planet Δ ∗ ∈ [ 0 0 3 , . 𝜋 ] . Both are normalized to be within [-1 1 . , ] A summary of the parameter bounds and how they are sampled can be found in Table 1.

## 2.4 Training

Theconstruction of the loss function for the absorption PINN initially appears straightforward, but involves intricate considerations. The boundary condition is established by the incident stellar radiation, characterized by the function 𝑢 𝑥 ( = -1 , 𝑦 ) = 1 0, representing the . intensity of starlight entering the system from the left boundary. This boundary condition is coupled with the RT equation governing the system:

Table 1. Ranges of the input parameters to the absorption and scattering PINN. Most are sampled uniformly within their respective intervals, with the exceptions of the angle 𝜙 , half of which is sampled uniformly within [-Δ ∗ , Δ ∗ ] , the scale height 𝐻 , which is sampled via a normal distribution in log-space, and the absorption profile 𝛼 𝑎 (see the Section 2.3 for details on the generation of the 𝛼 𝑎 ( 𝑃 ) -profiles). The 49 refers to the absorption or scattering coefficient at the layer with index 49 (ranging from 0 to 99).

| Input                           | Range                   |
|---------------------------------|-------------------------|
| 𝑥                               | [- 𝑙 𝑥 , 𝑙 𝑥 ]          |
| 𝑦                               | [- 𝑙 𝑦 , 𝑙 𝑦 ]          |
| 𝜙 †                             | [- 𝜋, 𝜋 ]               |
| 𝑅 P [ 𝑅 J ]                     | [ 0 . 2 , 2 . 0 ]       |
| log 10 ( 𝐻 [ cm ] )             | N( 7 . 1 , 0 . 35 2 )   |
| log 10 ( 𝛼 𝑎, 49 [ cm - 1 ] ) ∗ | [- 20 . 2 , - 4 . 5 ]   |
| log 10 ( 𝛼 𝑠, 49 [ cm - 1 ] ) † | [- 18 . 74 , - 6 . 74 ] |
| Δ † ∗                           | [ 0 , 0 . 3 𝜋 ]         |

∗ : input only to absorption PINN

† : input only to scattering PINN

$$\frac { d u } { d x } + \alpha ^ { \prime } \cdot u = 0 \,,$$

where the scaling factor is taken into account inside the 𝛼 ′ B 𝑙 𝑥 𝛼 𝑎 . However, utilizing this equation directly as the loss function presents significant challenges in solving the problem. To understand why this is the case, let us consider two extremes for different regions of the atmosphere: 𝛼 ′ ≪ 1 and 𝛼 ′ ≫ 1.

For an extremely high value of 𝛼 ′ , the RT equation effectively becomes 𝛼 ′ · 𝑢 = 0 with the simple solution 𝑢 = 0. Now consider small deviations of 𝑢 from 0; Because of the high 𝛼 ′ they can result in high values of the loss function and proportionally steep gradients. This means that the PINN effectively overfits these regions to be as close to zero as possible, even though other regions might be of more interest. We can solve this by dividing the whole equation by 𝛼 ′ . This gives us approximately 𝑢 = 0 as our loss function in such regions, resulting in the PINN learning the right solution without overfitting.

In the second case where 𝛼 ′ ≪ 1, the form of Eq. (10) would be preferable, since it would effectively be d 𝑢 / d 𝑥 = 0. Again, the neural network would learn the right solution of 𝑢 = const , but would not . overfit the region as would be the case if we would divide the gradient by 𝛼 ′ .

Wecan combine these two forms of the same equation by calculating both and taking the (element-wise) minimum. The RT equation, which only considers Rayleigh scattering, is normalized in a similar fashion, details of which can be found in Appendix C.

Meanwhile, the boundary condition of the scattering PINN has more complexity to it. On the left side, the light from the central star enters the atmosphere, this time with an angular extent. Because of this angular extent, also light from the upper boundary could enter the atmosphere, while from the right boundary, no light should enter but only leave the considered slab:

$$u ( x = - 1, y, | \phi | \leq \frac { \pi } { 2 } ) & = 1 \text{ if } | \phi | \leq \Delta _ { * }, \text{else} = 0 \\ \text{i nitially } u ( x, y = 1, \phi \leq 0 ) & = 1 \text{ if } | \phi | \leq \Delta _ { * }, \text{else} = 0 \quad \quad ( 1 1 ) \\ s. \text{ The } u ( x = 1, y, | \phi | \geq \frac { \pi } { 2 } ) & = 0 \,. \\ \cdot \quad \, \quad \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real ^ { \, \real } \, \real.$$

For the lower boundary, defining an exact boundary condition poses a challenge due to the uncertainty in the amount of light that could be scattered into the observed atmospheric slab from below. However, the absence of a boundary condition can potentially lead

Figure 3. The architecture of the scattering PINN, where the numbers in the (fully-connected) layers indicate their width. The 4th layer connects both to an output 𝑢 𝑎 and to two subsequent layers that have an output 𝑢 𝑠 . Because 𝑢 𝑎 models the part of the solution that fulfills the RTE without scattering, a filter is applied, dependent on 𝜙 , such that 𝑢 𝑎 = 0 where the light cannot get without changing directions. The combined output 𝑢 𝑎 + 𝑢 𝑠 should fulfill the complete RTE, such that 𝑢 𝑠 only has to model the scattered light component of the solution. Since the diffuse light only becomes relevant for large enough Δ ∗ , 𝑢 𝑠 is multiplied with the solid angle subtended by the star.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00084v1-figure-d0ae8d23-425f-4cb5-a318-9def71fd2ec1.jpg)

to erroneous solutions. To mitigate this issue, we approximate the lower boundary condition by setting the scattered light component to zero. This implies that no additional scattered light enters the system from this boundary. This approximation is expected to yield anaccurate solution, provided that the optical depth of the atmosphere is sufficiently high.

For the absorption component, an accurate boundary condition can be found by recognizing the spherical symmetry of the system and rotating the cord of the light ray along a path of constant radius until the point of interest reaches the x-coordinate of 𝑥 = 0:

$$& u _ { a } ( x _ { 0 } \leq 0, y = - 1, \phi \geq 0 ) = u _ { a } ( x = 0, y ^ { \prime }, \phi + \Delta \phi ), \quad \text{$\mod$} \\ & \text{where $y^{\prime} = y \left ( x _ { 0 }, r ( x _ { 0 }, y = - 1 ) \right ), \, \tan ( \Delta \phi ) = \frac { l _ { x } x _ { 0 } } { r _ { 9 9 } }. \quad \text{$(12)} \quad \text{$\text{$thephi$}$} \\ & \text{to opti}$$

Note that this only works for 𝑥 ≤ 0. Otherwise, this boundary condition is set to zero, which again should hold, if the optical depth of the atmosphere is high enough at the highest pressures. The other edge case that we have to consider is if abs ( 𝜙 + Δ 𝜙 ) &gt; Δ ∗ , max ( = 0 3 . 𝜋 ) , for which 𝑢 𝑎 is zero because of the inherent architecture of the PINN. As this only happens for particularly low values of 𝑟 99 and large values of | 𝑙 𝑥 𝑥 0 | (i.e., for the outer layers of the atmosphere), the boundary condition is set to 𝑢 𝑎 = 1 for these rare instances.

The final boundary condition is to ensure the periodicity of 𝜙 : 𝑢 𝑥, 𝑦, 𝜙 ( = 𝜋 ) = 𝑢 𝑥, 𝑦, 𝜙 ( = - ) 𝜋 .

Because of its boundary conditions, the scattering PINN has a special architecture, as shown in Fig. 3. This architecture has two outputs: One ( 𝑢 𝑎 ) should fulfill the RTE without scattering and utilizes a hard-coded filter, such that the neural network does not have to learn the non-differentiable boundary at | 𝜙 | = Δ ∗ , where the value of 𝑢 can jump from 1 to 0. The other output ( 𝑢 𝑠 ) should then model the diffuse component of the light such that when adding the two outputs 𝑢 = ( 𝑢 𝑎 + 𝑢 𝑠 ) , they together fulfill the complete RTE with scattering.

If we would now sample 𝜙 uniformly in [-𝜋, 𝜋 ] , most of the "absorption output" 𝑢 𝑎 would equate to zero only because of the mentioned filter. This not only wastes a lot of training points, but also effectively lowers the weight (or importance) of the corresponding residual in the loss term. We mitigate this problem by sampling half of 𝜙 within [-Δ ∗ , Δ ∗ ] .

Both the boundary conditions and the residuals enter the loss function with their mean square. To make sure that the boundary conditions are fulfilled first, a weight of 0.5 was given to the residual term ( 𝜔 F = 0 5 and . 𝜔 B = 1 in Eq. (3)). The lower boundary condition is given a weight of 0.1, as it is reasonable to fulfill only if 𝑢 𝑎 at 𝑥 = 0 is already sufficiently accurate. The final loss function is then the logarithm of their weighted sum.

The used optimizer is the higher-order quasi-Newton L-BFGS. We found that initially using Adam and then switching to L-BFGS yielded no benefit, so L-BFGS was used from the start. The parameters of the optimizer differ from the default PyTorch parameters in a learning rate of 0.8, the maximal number of both function evaluations and iterations per optimization steps of 50 (or 20), the tolerance change is set to the machine epsilon and the Strong Wolfe condition is chosen for the line search. As described in many PINN papers before (e.g. Sankaran et al. 2022), we confirmed that having a larger batch size generally improves the final model. Accordingly, we always used the maximum batch size possible on an NVIDIA A100-SXM4-40GB GPU,withdoublethe number of points for the residual loss compared to the boundary loss. To mitigate overfitting, the batch was resampled every 50 training steps (every 20 steps for the first 8000 steps in the case of the scattering PINN). The parameters that are computationally more expensive to sample (mainly the 𝛼 𝑎 -profiles, but also the two radius profile parameters) were sampled pre-training with a size of 10 7 . During training, these parameters were then sampled from this saved array and loaded onto the GPU, thereby accelerating the process.

Throughout the training of the absorption PINNs, the inputs were resampled 100 times. Considering the number of points for which the residual loss is computed (2 21 or 2 22 ), the total number of input combinations amounts to roughly 2 1 . × 10 8 or 4 2 . × 10 8 , depending on the network size. For the scattering PINN, which was resampled 400 times and repeated thereafter, with 2 19 residual points, the model was similarly trained on approximately 2 1 . × 10 8 unique input combinations.

For the absorption PINN, a hyperparameter grid search is used to optimize the number of layers (4 or 6) as well as the width of each layer (32, 64, or 128). A PINN with more free parameters may improve our result even more, but the 6x128 network already needed more time to calculate a spectrum than the classical numerical method, so we decided against trying out a deeper or wider network. For each configuration, 7 PINNs were trained with different random Xavier initialization. They were evaluated on 95 test spectra with the help of petitRADTRANS (Mollière et al. 2019).

The number of subsequent layers of the scattering PINN were not optimized, although one layer was insufficient to accurately model the correct solution.

## 3 RESULTS

## 3.1 Absorption PINN

After trying out what maximum batch size is possible for each hyperparameter configuration, we found that the three models with the most adjustable parameters (6x64, 4x128, 6x128) are able to utilize a batch size of 3 × 2 20 while the three smaller models (4x32, 6x32, 4x64) could use double the batch size (3 × 2 21 ).

The distribution of the root mean squared errors with respect to each example spectrum can be seen in Fig. 4. It seems that the 4x64 model consistently gives us the best result, although no architecture performs significantly worse.

Given that the errors of the smallest networks with a width of 32 were not significantly higher than those of larger networks, we experimented with an even smaller architecture of 4x16 and 4x8. For

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00084v1-figure-b72e8223-ef70-4b7e-bd0c-86ef38b14f0e.jpg)

Figure 4. Violin plot of the 95 root mean squared errors of each test spectrum for each trained absorption PINN. The maximum RMSE of the 4x8 PINN lies outside the shown range at 7.7%.

Figure 5. Example for the evolution of the different loss components, for the 3rd 4x64 PINN. The total loss L is calculated from the boundary ( LB ) and residual loss ( LF ) as log 10 ( L) = log 10 ( L B + 0 5 . LF ) , as defined in Eq. 3 5.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00084v1-figure-8bedb501-ca60-4be3-bb88-c90c6ac2484d.jpg)

Figure 6. Example spectrum for one of the best-performing absorption PINNs, the 3rd 4x64 PINN. A hydrogen- and helium-dominated atmosphere, containing H 2 O, CO, CO 2 , CH 4 , Na, and K.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00084v1-figure-68da27aa-23c6-4a91-b478-4ea023ff2694.jpg)

the 4x16 PINN, this adjustment resulted in a slight increase in error, similar to the worst-performing 6x32 network, and a considerably larger error in the case of the 4x8 PINN, but it produced diminishing returns in terms of speed.

For one of the 4x64 models, we also show the evolution of the different loss components during training in Fig 5. Here we see how the weight 𝜔 F of 0.5 allows the network to quickly fulfill the boundary condition, while still giving the residual loss component sufficient weight to decrease over the epochs. The final value of this residual loss log 10 (L F ) = -4 0 already gives us a first order of estimate for . the errors of the PINN of 10 -2 .

We also checked for possible correlations of the error with any atmospheric parameters. We only found that for really low radii, our formula for the radius profile 𝑟 ( 𝑃 ) breaks down for the defined pressure range, yielding inaccurate results.

Fig. 6 presents an example spectrum, notably incorporating sodium (Na) and potassium (K). These elements serve as robust indicators for assessing the PINN's capacity to generalize beyond its initial training data, due to their potential for exhibiting abrupt increases in absorption coefficients at specific atmospheric layers. During the training phase, the 𝛼 𝑎 ( 𝑃 ) profiles we generated typically did not showcase logarithmic gradients d ln ( 𝛼 𝑎 ) d ln ( 𝑃 ) &gt; 2. However, for these atomic species, the gradient between two layers can surge to values well above 2, occasionally surpassing 10. Such dramatic shifts can result in the optical depths of subsequent layers abruptly increasing from negligible to substantial values. Although 𝛼 𝑎 ( 𝑃 ) profiles with such non-differentiable characteristics were not included in the training due to the complexity they introduced in minimizing the residual loss, the PINN demonstrates a commendable ability to approximate accurate solutions, even though this results in slightly larger errors than for known 𝛼 𝑎 ( 𝑃 ) profiles. This success underscores the network's adaptability and effectiveness, showcasing its capability to handle complex scenarios that extend beyond the initial training data.

While the observed accuracy gives us confidence in our method, it does not offer any improvement in performance. Using the larger models of width 128 can even significantly slow down the simulation of a transit spectrum compared to the numerical method.

This contrasts with the results of Himes et al. (2022), who could significantly speed up their forward model using a neural network. But their method directly calculates the spectrum (in a fixed spectral range and spectral resolution) from the atmospheric parameters (also fixed in their number and range). Our PINN, on the other hand, calculates the transmission function 𝑢 𝑦 ( ) for every individual 𝛼 𝑎 ( 𝑃 ) profile, so for every wavelength and every source of absorption. Depending on the number of wavelength points and absorbing species, this can mean that the PINN has to calculate the solution in multiple batches, significantly reducing its speed. The advantage of our model lies within its versatility, as it is not trained on, e.g., specific molecules.

## 3.2 Scattering PINN

Since we cannot easily compare the spectra obtained from the scattering PINN to the ones from a traditional numerical method, which do not take scattering into account, evaluating its performance is harder to do.

In Fig. 7, we can see an example solution with the two components, the absorbed 𝑢 𝑎 and scattered light 𝑢 𝑠 separated, as well as the corresponding residuals. These residuals show how well the RT equation is fulfilled and give us an estimate of the error. In the example, we can see residuals on the order of 10 -2 spanning a significant fraction of the considered slab.

Additionally, we compared our method to a fixed PINN for the same example, the result of which can be seen in Fig. 8. By fixed PINN, we mean that all the parameters of this PINN have been fixed to certain values during training, apart from the spacial coordinates

Figure 7. Example solution for the 𝛼 𝑥, 𝑦 ( ) shown in Fig. 1 for the direction of the light in positive x-direction ( 𝜙 = 0).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00084v1-figure-fcb42574-99c5-4a48-a4a4-b80103872f53.jpg)

Top: The complete solution (left), the solution if scattering would be treated as absorption (middle) and the difference between the two: the scattered light component (right), all in units of 𝑢 0 : = 𝑢 𝑥 ( = -1 , 𝜙 = 0 ) = 1.

Bottom: The logarithm (log 10 ) of the respective residuals of the differential equations. The middle plot shows the residual of 𝑢 𝑎 with respect to the RTE without scattering and the right one the residual of ( 𝑢 𝑎 + 𝑢 𝑠 ) with respect to the complete RTE. The left plot is calculated from the two residuals by taking their root mean square.

𝑥, 𝑦, 𝜙 (and Δ ∗ , which cannot be fixed because of the boundary condition at 𝑦 = -1). Leaving the batch size the same should yield a solution of higher accuracy, as the fixed PINN has to explore a significantly smaller parameter space. Similarly to the absorption PINN, we can see errors of around 1%.

Finally, we want to explore the usefulness of this method by plotting an example spectrum (see Fig. 9). Even for the extreme case of Δ ∗ = 0 15 . 𝜋 , which would correspond to an exoplanet orbiting a sun-like star ( 𝑅 ∗ = 𝑅 ⊙ ) at a distance of 𝑎 = 0 01 AU, the Rayleigh . scattering barely influences the spectrum. This is because our assumptions for calculating the scattering integral only hold for small angles. Thus, while the scattering PINN is able to find an accurate solution, the approximation of parameterizing our problem with only one angle is not sufficient to accurately model Rayleigh scattering for large angular extents of its host star, which is when it becomes relevant.

## 4 LIMITATIONS

The assumptions under which both the absorption and scattering PINN operate include isothermality of the atmosphere and a radius range of 𝑅 0 ∈ [ 0 2 . , 2 0 . ] 𝑅 J .

In the case of the scattering PINN we further assume an orbital inclination of 𝑖 = 90 . It also focuses narrowly on pressure-independent ◦ scattering cross-sections, characteristic of Rayleigh scattering from hydrogen and helium, and does not consider any absorbing species simultaneously being present in the atmosphere. Additionally, it lacks an exact lower boundary condition. The most significant shortcoming is the small angle approximation for the calculation of the scattering integral, which is ultimately insufficient to accurately model Rayleigh scattering.

Lastly, the fully-connected layers of both PINNs offer no interpretability, although the special architecture of the scattering PINN at least simplifies differentiating the diffuse light component from the total solution.

Figure 8. Comparison of the solution from the scattering PINN shown in Fig. 7 with a higher-accuracy PINN with fixed parameters (apart from 𝑥, 𝑦, 𝜃 and Δ ∗ ). Grey indicates differences of less than 0.1%.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00084v1-figure-0d995e71-5f6f-42c4-baa3-bf5c6c6b24f6.jpg)

Figure 9. Example spectrum similar to Fig. 6 (without Na and K) with the solutions from petitRADTRANS (Mollière et al. 2019) and the scattering PINN, each for different values of Δ ∗ . We can see that even in the extreme case of Δ ∗ = 0 15 . 𝜋 Rayleigh scattering barely influences the spectrum, when using our approximations.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00084v1-figure-a4631081-2b58-4fe8-8978-ad9e1bab9739.jpg)

## 5 CONCLUSIONS

In this study, we have assessed the effectiveness of PINNs in simulating complex atmospheric phenomena like scattering, which are traditionally challenging and resource-intensive to model.

The absorption PINN demonstrated in this research exhibits a high level of precision in forecasting atmospheric transmission for transiting exoplanets, achieving an error margin within a few percent, even when confronted with unfamiliar non-differentiable 𝛼 𝑎 ( 𝑃 ) profiles. This accuracy underscores the model's robustness and reliability in complex atmospheric simulations. However, it is worth noting that, in terms of computational speed, this method does not yet surpass the efficiency of existing numerical techniques. This opens up opportunities for further optimization and development in the realm of PINN methodologies to enhance their speed performance while maintaining their high accuracy. The use of other architectures than the simple MLP chosen in this work, could aid in achieving this. A 1D-CNN for the 𝛼 𝑎 inputs could, for example, significantly reduce the number of weights in the network.

The Rayleigh scattering PINN in this study, while showcasing po-

tential, operates within certain constraints, as summarized in Section 4, that currently limit its broader applicability. Having these caveats in mind, nevertheless, a parameterized PINN could prove effective in modeling scattering processes in the atmosphere with a minimal additional computational cost beyond the training phase.

In this work, we have laid out one possible path of how to present this complex problem to a neural network. Here we want to emphasize the key ingredients for successfully training a PINN:

- · especially when dealing with a complicated differential equation, one has to carefully consider all possible configurations of parameters, how they affect the absolute values of the respective residual and normalize all cases accordingly
- · similar to the insight from Wang et al. (2022), it is important to respect causality when choosing weights for different parts of the loss function
- · having a larger batch size generally improves the accuracy of a PINN
- · incorporating physical knowledge into, e.g., the architecture of the PINN can significantly simplify finding the correct solution, especially if any non-differentiable boundaries are involved

However, further enhancements to this approach are necessary to calculate the correct solution to the RTE with scattering, where it becomes relevant for transit spectra. E.g., auxiliary PINNs (Yuan et al. 2022; Riganti &amp; Negro 2023) could aid in expanding the problem to a second angular dimension ( 𝜃 ), while making it possible to evaluate the scattering integral without any numerical integration. Experimenting with alternative architectures like FNOs (Li et al. 2021) could also make addressing these challenges more feasible.

In the future, the inclusion of specific pressure-temperature profiles can be achieved by simply extending the input parameters, though this might slightly impact the model's accuracy due to the increased dimensionality of the problem. Integrating other scattering sources, such as clouds with their variable 𝛼 𝑃 ( ) profiles and complex scattering phase functions Φ ( 𝜃, 𝜃 ′ ) , presents a more significant challenge. One major concern is the possible non-differentiability of cloud decks. But, as the absorption model's predictions for Na and K show, a PINN might be able to overcome such difficulties.

## ACKNOWLEDGEMENTS

We thank the referee for the helpful comments. We also thank Lukas Heinrich, Nicole Hartman and the ORIGINS Data Science Lab (ODSL) for valuable discussions and hardware support during the development of the project, which was made possible through the funding of the ORIGINS Excellence Cluster by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC 2094 - 390783311; in particular through ORIGINS Excellence Cluster's Seed Money Program.

This research has made use of the NASA Astrophysics Data System, and the Python packages NumPy, SciPy, Matplotlib, and PyTorch.

## DATA AVAILABILITY

The code, based on Mishra &amp; Molinaro (2021), is available at https://github.com/DavidDahlbudding/ AtmosphericScatteringPinn (commit 3c805f5). The respective models and other data can be downloaded from https://zenodo.org/doi/10.5281/zenodo.10727497 .

## REFERENCES

Ardévol Martínez F., Min M., Kamp I., Palmer P. I., 2022, A&amp;A, 662, A108 Basri R., Galun M., Geifman A., Jacobs D., Kasten Y., Kritchman S., 2020, in III H. D., Singh A., eds, Proceedings of Machine Learning Research Vol. 119, Proceedings of the 37th International Conference on Machine Learning. PMLR, pp 685-694, https://proceedings.mlr.press/ v119/basri20a.html

| Borysow A., 2002, A&A, 390, 779                                                                                                                                                                   |
|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Borysow A., Frommhold L., 1989, ApJ, 341, 549                                                                                                                                                     |
| Borysow J., Frommhold L., Birnbaum G., 1988, ApJ, 326, 509                                                                                                                                        |
| Borysow A., Frommhold L., Moraldi M., 1989, ApJ, 336, 495                                                                                                                                         |
| Borysow A., Jorgensen U. G., Fu Y., 2001, J. Quant. Spectrosc. Radiative Transfer, 68, 235                                                                                                        |
| Chan Y. M., Dalgarno A., 1965, Proceedings of the Physical Society, 85, 227 Chubb K. L., et al., 2021, A&A, 646, A21                                                                              |
| Cobb A. D., et al., 2019, AJ, 158, 33                                                                                                                                                             |
| Cuomo S., Schiano di Cola V., Giampaolo F., Rozza G., Raissi M.,                                                                                                                                  |
| Piccialli F., 2022, Journal of Scientific Computing,                                                                                                                                              |
| Dalgarno A., Williams D. A., 1962, ApJ, 136, 690 Fisher C., Hoeĳmakers H. J., Kitzmann D., Márquez-Neila P., Grimm S. L., Sznitman R., Heng K., 2020, AJ, 159, 192                                |
| Gebhard T. D., Wildberger J., Dax M., Angerhausen D., Quanz S. P., Schölkopf B., 2023, arXiv e-prints, p. arXiv:2312.08295 Goodfellow I., Bengio Y., Courville A., 2016, Deep Learning. MIT Press |
| Hendrix J. L. A. M., Louca A. J., Miguel Y., 2023, MNRAS, 524, 643 Himes M. D., et al., 2022, PSJ, 3, 91                                                                                          |
| Hou Yip K., Changeat Q., Al-Refaie A., Waldmann I., 2022, arXiv e-prints, p. arXiv:2205.07037                                                                                                     |
| Li Z., Kovachki N. B., Azizzadenesheli K., liu B., Bhattacharya K., Stuart A., Anandkumar A., 2021, in International Conference on Learning Repre-                                                |
| sentations. https://openreview.net/forum?id=c8P9NQVtmnO                                                                                                                                           |
| Márquez-Neila P., Fisher C., Sznitman R., Heng K., 2018, Nature Astronomy, 2, 719                                                                                                                 |
| Mishra S., Molinaro R., 2021, J. Quant. Spectrosc. Radiative Transfer, 270, 107705                                                                                                                |
| Mollière P., Wardenier J. P., van Boekel R., Henning T., Molaverdikhani K., Snellen I. A. G., 2019, A&A, 627, A67                                                                                 |
| Nixon M. C., Madhusudhan N., 2020, MNRAS, 496, 269 Piskunov N. E., Kupka F., Ryabchikova T. A., Weiss W. W., Jeffery C. 1995, A&AS, 112, 525                                                      |
| S., Raissi M., Perdikaris P., Karniadakis G. E., 2019, Journal of Computational                                                                                                                   |
| Physics, 378, 686 Riganti R., Negro L. D., 2023, Applied Physics Letters, 123, 171104 Rocchetto M., Waldmann I. P., Venot O., Lagage P. O., Tinetti G., 2016,                                     |
| 833, 120                                                                                                                                                                                          |
| Rothman L. S., et al., 2010, J. Quant. Spectrosc. Radiative Transfer,                                                                                                                             |
| Sankaran S., Wang H., Guilhoto L. F., Perdikaris P., 2022, in The of Deep Learning and Differential Equations II.                                                                                 |
| net/forum?id=THCvohg1RV                                                                                                                                                                           |
| S., Chakrabarty A., Tinetti G., 2020, ApJ, 889, 181 J., 2004, in Fischer R., Preuss R., Toussaint U. V.,                                                                                          |
| Sengupta Skilling eds, Institute of Physics Conference Series Vol. 735, Bayesian Inference                                                                                                        |
| American and Maximum Entropy Methods in Science and Engineering: 24th Interna- tional Workshop on Bayesian Inference and Maximum Entropy Methods                                                  |
| in Science and Engineering. pp 395-405, doi:10.1063/1.1835238 F., et al., 2018, arXiv e-prints, p. arXiv:1811.03390                                                                               |
| M., Rozet F., Absil O., Mollière P., Nasedkin E., Louppe G., A&A, 672, A147 I. P., 2016, ApJ, 820, 107                                                                                            |
| Waldmann S., Sankaran S., Perdikaris P., 2022, arXiv e-prints,                                                                                                                                    |
| arXiv:2203.07404                                                                                                                                                                                  |
| L., Ni Y.-Q., Deng X.-Y., Hao S., 2022, Journal of                                                                                                                                                |
| Yuan Physics, 462, 111260                                                                                                                                                                         |
| Zingales T., Waldmann I. P., 2018, AJ, 156, 268                                                                                                                                                   |
| ApJ,                                                                                                                                                                                              |
| 111,                                                                                                                                                                                              |
| 2139                                                                                                                                                                                              |
| https://openreview.                                                                                                                                                                               |
| Symbiosis                                                                                                                                                                                         |
| Soboczenski Vasist                                                                                                                                                                                |
| Wang                                                                                                                                                                                              |
| 2023,                                                                                                                                                                                             |
| p.                                                                                                                                                                                                |
| Computational                                                                                                                                                                                     |

## APPENDIX A: THE MULTI-LAYER PERCEPTRON

A feedforward neural network or multilayer perceptron (MLP) is a versatile approximator that maps 𝒙 to 𝒚 = 𝑓 ( 𝒙 𝜽 , ) , where 𝜽 defines the free network parameters that are adjusted during training.

What makes these standard MLPs so versatile is just simple linear algebra: a linear mapping 𝑾 (the weights), a displacement 𝒃 (the biases) and, most importantly, a non-linear activation function 𝜎 . The output 𝒚 of one such a layer 𝑓 1 can then be computed with:

$$y & = f _ { 1 } ( x ) = \sigma ( W x + b ), & ( A 1 ) & \quad \text{For}$$

where the activation function 𝜎 is applied element-wise: 𝜎 (( 𝑥 1 , 𝑥 2 ) 𝑇 ) = ( 𝜎 𝑥 ( 1 ) , 𝜎 ( 𝑥 2 )) 𝑇 . Multiple of these layers can be stacked on top of each other to be able to represent highly complex functions between the input 𝒙 and the output 𝒚 = 𝑓 𝑛 ◦· · · ◦ 𝑓 2 ◦ 𝑓 1 ( 𝒙 ) (where ◦ denotes a function composition: 𝑓 2 ◦ 𝑓 1 ( 𝒙 ) = 𝑓 2 ( 𝑓 1 ( 𝒙 )) ) (e.g. Goodfellow et al. 2016). The layers between the input and the output layer are called the hidden layers and the number of hidden layers is called the depth of the network, while the number of nodes per layer is called the width. This simple architecture is also mostly used in this paper.

To train a neural network, we need another important function, the so-called objective or loss function. This function should be minimized to achieve a desired goal (e.g., the mean squared error when fitting a function to given data points). Then, the gradient of the loss function with respect to the weights and biases can be cumulatively calculated for a batch of data using automatic differentiation, and the neural network can "move" towards a minimum by adjusting its parameters in small steps in the negative direction of this gradient.

## APPENDIX B: THE SCATTERING INTEGRAL

As mentioned in Section 2.2 we need to make a few approximations valid for small angular extents of the star for the problem to be parameterizable by only one angle ( 𝜙 ). Here, we choose different symmetries for the two components of the radiant flux: 𝑢 𝑎 and 𝑢 𝑠 . Because they are added together to get the total solution to the RTE, we can calculate the total integral by adding the two integrals for each component:

$$\int \Phi \, u \, d \Omega = \int \Phi \, u _ { a } \, d \Omega + \int \Phi \, u _ { s } \, d \Omega \quad & & ( B 1 )$$

For the absorption component 𝑢 𝑎 we choose the conventional spherical (local) coordinate system, such that the angle 𝜙 ∈ [-𝜋, 𝜋 ] lies in the x-y-plane and 𝜃 ∈ [ 0 , 𝜋 ] is the angle relative to the z-axis. If we now take a fixed 𝜙 , let 𝜃 vary and expand this into a plain, the intersection of this plain with the sphere of an atmospheric layer is a circle. For small angular extents of the star each light ray for varying 𝜃 travels roughly the same path through this circle, which is why we assume the symmetry of 𝑢 𝑎 ( 𝜙, 𝜃 ) ≈ 𝑢 𝑎 ( 𝜙, 𝜃 = 𝜋 2 ) .

We can parameterize 𝜃 ′ in terms of 𝜙 ′ such that we can integrate over the complete solid angle of the star with angular extent 2 Δ ∗ , for which 𝑢 𝑎 ≠ 0, by using the law of cosines with 𝜃 = 𝜋 2 :

$$\cos \Delta & = \sin \theta ^ { \prime } \cos ( \phi - \phi ^ { \prime } ). & ( B 2 ) \quad \inf _ { \phi \, 1 0 } \,,$$

Inserting this into the scattering phase function Φ Δ ( ) (7) and the integration limits, yields:

$$P \right ) & \text{ is a} \\ \text{defines} \quad & \int _ { - \Delta _ { * } } ^ { \Delta _ { * } } d \phi ^ { \prime } \, u _ { a } ( \phi ^ { \prime } ) \int _ { - \theta ^ { \prime } ( \phi ^ { \prime } ) } ^ { \theta ^ { \prime } ( \phi ^ { \prime } ) } d \theta ^ { \prime } \, | \sin \theta ^ { \prime } | \, \Phi ( \Delta ) = \\ \colon \text{linear} \quad & \int _ { - \Delta _ { * } } ^ { \Delta _ { * } } u _ { a } ( \phi ^ { \prime } ) \, \frac { 3 } { 2 } \sqrt { 1 - a ^ { 2 } } \left ( \cos ^ { 2 } ( \phi - \phi ^ { \prime } ) \left ( \frac { 1 - a ^ { 2 } } { 3 } - 1 \right ) - 1 \right ) d \phi ^ { \prime }, \\ \text{ion} \, \sigma. \quad & \text{ where } a \coloneqq \frac { \cos \Delta _ { * } } { \cos \phi ^ { \prime } } \\ \text{h:} \quad & \int _ { \Delta 1 \right \rangle } ^ { \quad } \quad \text{ For exnolitino} \, a \, \text{vmmetrix} v \, \text{ of the scatterino} \, \text{ commonent we have to}$$

For exploiting a symmetry of the scattering component we have to rotate our local coordinate system. The 𝜙 𝑠 , coordinate now lies in the y-z-plane and 𝜃 𝑠 is the angle relative to the x-axis. The absorption component 𝑢 𝑎 therefore is a concentrated beam of light around 𝜃 𝑠 = 0 for small Δ ∗ . Since the scattering phase function is only dependent on the angular distance Δ , relative to 𝜃 = 0, we get Δ = 𝜃 . Hence we can approximate, while allowing for some slight asymmetry: 𝑢 𝑠 ( 𝜃 𝑠 , 𝜙 𝑠 ≶ 0 ) ≈ 𝑢 𝑠 ( 𝜃 𝑠 , 𝜙 = ± 𝜋 2 ) . The integral can therefore be calculated using:

$$\colon n o d e s \, \text{mostly} \quad \iint \Phi ( \Delta ) \, u _ { s } ( \theta _ { s } ^ { \prime }, \phi _ { s } ^ { \prime } ) \, | \sin \theta _ { s } ^ { \prime } | \, d \theta _ { s } ^ { \prime } \, d \phi _ { s } ^ { \prime } = \\ \colon \min _ { 0 } \quad \int _ { 0 } ^ { \pi } | \sin \theta _ { s } ^ { \prime } | \, d \theta _ { s } ^ { \prime } \left [ u _ { s } ( \theta _ { s } ^ { \prime }, \phi _ { s } ^ { \prime } = - \frac { \pi } { 2 } ) \int _ { - \pi } ^ { 0 } \Phi ( \Delta ) d \phi _ { s } ^ { \prime } + \quad \quad ( B 4 ) \\ \text{he loss} \quad \text{$u_{s}(\theta _ { s } ^ { \prime }, \phi _ { s } ^ { \prime } = \frac { \pi } { 2 } ) \int _ { 0 } ^ { \pi } \Phi ( \Delta ) d \phi _ { s } ^ { \prime } \right ] \,. \\ \text{n, and} \quad \text{$The interval over $\phi^{\prime}$ can again be incorporated analyticallv hv insert-}$$

The integral over 𝜙 ′ 𝑠 can again be integrated analytically by inserting the law of cosines (assuming 𝜙 𝑠 = ± 𝜋 2 ) into the scattering phase function Φ Δ ( ) (7):

$$& \int _ { 0 } ^ { \pm \pi } \Phi ( \Delta ) d \phi _ { s } ^ { \prime } = \\ \intertext { a t o b e } \text{ferent} & \quad \pm \frac { 3 } { 4 } \left [ \pi \left ( 1 + \cos \theta _ { s } \cos \theta _ { s } ^ { \prime } \pm \frac { 1 } { 2 } \sin \theta _ { s } \sin \theta _ { s } ^ { \prime } \right ) & & ( B 5 ) \\ \text{ferent} & \quad \pm \text{sgn} ( \phi _ { s } ) \sin ( 2 \theta _ { s } ) \sin ( 2 \theta _ { s } ^ { \prime } ) \right ]. \\ \text{nd } u. & \quad \text{the numerical information above the Transactual rule with } a$$

The numerical integration employs the Trapezoidal rule with a total of 19 sample points in the interval [-𝜋, 𝜋 ] , where 5 of these points are equally spaced in the interval [-Δ ∗ , Δ ∗ ] and the remaining intervals [-𝜋, -Δ ∗ ] and [ Δ ∗ , 𝜋 ] each have another 7 equidistant sample points.

## APPENDIX C: RESIDUAL LOSS FOR SCATTERING PINN

As mentioned in Section 2.3, we have to divide the gradients in the RTE (Eq. 6) with the respective scale lengths 𝑙 𝑥 and 𝑙 𝑦 . Given in cm, these parameters are quite large and can therefore diminish the value of the residual loss, especially for low values of 𝛼 . Hence, we multiply the whole equation with 𝑙 𝑥 :

$$\begin{smallmatrix} \ell ^ { 2 } \\ \text{tegrate} \\ \Delta _ { s }, \text{ for} \\ \end{smallmatrix} \quad \cos ( \phi ) \frac { d u } { d x } + \sin ( \phi ) \frac { l _ { x } } { l _ { y } } \frac { d u } { d y } + l _ { x } \alpha _ { s } \cdot u \\ - \frac { l _ { x } \alpha _ { s } } { 4 \pi } \cdot \int \Phi ( \theta, \phi, \theta ^ { \prime }, \phi ^ { \prime } ) u ( \theta ^ { \prime }, \phi ^ { \prime } ) d \Omega ^ { \prime } = \text{res} _ { 1 }.$$

The remaining problem is that the value of 𝑙 𝑥 𝑙 𝑦 can be on the order of 10, possibly inflating the value of the residual loss for a certain set of parameters. Because of this, we further normalize this formulation of the residual with norm = | cos ( 𝜙 )| + 𝑙 𝑥 𝑙 𝑦 | sin ( 𝜙 )| .

## 10 D. Dahlbüdding et al.

The second formulation is similar to the one for the absorption PINN, where we divide the RTE by 𝛼 𝑠 :

$$\frac { \cos ( \phi ) } { l _ { x } \alpha _ { s } } \frac { d u } { d x } + \frac { \sin ( \phi ) } { l _ { y } \alpha _ { s } } \frac { d u } { d y } + u \\ - \frac { 1 } { 4 \pi } \cdot \int \Phi ( \theta, \phi, \theta ^ { \prime }, \phi ^ { \prime } ) u ( \theta ^ { \prime }, \phi ^ { \prime } ) d \Omega ^ { \prime } = \text{res} _ { 2 }.$$

The complete residual is then again the minimum of the two formulations:

$$\text{res} = \min \left ( \frac { \text{res} _ { 1 } } { \text{norm} }, \text{res} _ { 2 } \right ).$$

This paper has been typeset from a T E X/L A T E X file prepared by the author.
David Dahlbüdding Karan Molaverdikhani Barbara Ercolano Tommaso Grassi
Approximating Rayleigh Scattering in Exoplanetary Atmospheres using Physics-informed Neural Networks (PINNs)
This research introduces an innovative application of physics-informed neural
networks (PINNs) to tackle the intricate challenges of radiative transfer (RT)
modeling in exoplanetary atmospheres, with a special focus on efficiently
handling scattering phenomena. Traditional RT models often simplify scattering
as absorption, leading to inaccuracies. Our approach utilizes PINNs, noted for
their ability to incorporate the governing differential equations of RT
directly into their loss function, thus offering a more precise yet potentially
fast modeling technique. The core of our method involves the development of a
parameterized PINN tailored for a modified RT equation, enhancing its
adaptability to various atmospheric scenarios. We focus on RT in transiting
exoplanet atmospheres using a simplified 1D isothermal model with
pressure-dependent coefficients for absorption and Rayleigh scattering. In
scenarios of pure absorption, the PINN demonstrates its effectiveness in
predicting transmission spectra for diverse absorption profiles. For Rayleigh
scattering, the network successfully computes the RT equation, addressing both
direct and diffuse stellar light components. While our preliminary results with
simplified models are promising, indicating the potential of PINNs in improving
RT calculations, we acknowledge the errors stemming from our approximations as
well as the challenges in applying this technique to more complex atmospheric
conditions. Specifically, extending our approach to atmospheres with intricate
temperature-pressure profiles and varying scattering properties, such as those
introduced by clouds and hazes, remains a significant area for future
development.