# Exact Solution of the Einstein Field Equations
for a
Spherical Shell of Fluid Matter

Jorge L. deLyra111Email<EMAIL_ADDRESS>
Rodrigo de A. Orselli222Email<EMAIL_ADDRESS>
C. E. I. Carneiro333Email<EMAIL_ADDRESS>
Universidade de São Paulo
Instituto de Física
Rua do Matão, 1371,
05508-090 São Paulo, SP, Brazil

(April 5, 2021)

###### Abstract

We determine the exact solution of the Einstein field equations for the case
of a spherically symmetric shell of liquid matter, characterized by an energy
density which is constant with the Schwarzschild radial coordinate $r$ between
two values $r_{1}$ and $r_{2}$. The solution is given in three regions, one
being the well-known analytical Schwarzschild solution in the outer vacuum
region, one being determined analytically in the inner vacuum region, and one
being determined mostly analytically but partially numerically, within the
matter region. The solutions for the temporal coefficient of the metric and
for the pressure within this region are given in terms of a non-elementary but
fairly straightforward real integral. For some values of the parameters this
integral can be written in terms of elementary functions.

We show that in this solution there is a singularity at the origin, and give
the parameters of that singularity in terms of the geometrical and physical
parameters of the shell. This does not correspond to an infinite concentration
of matter, but in fact to zero energy density at the center. It does, however,
imply that the spacetime within the spherical cavity is not flat, so that
there is a non-trivial gravitational field there, in contrast with Newtonian
gravitation. This gravitational field is repulsive with respect to the origin,
and thus has the effect of stabilizing the geometrical configuration of the
matter, since any particle of the matter that wanders out into either one of
the vacuum regions tends to be brought back to the bulk of the matter by the
gravitational field.

##  1 Introduction

The Schwarzschild external solution [1, 2] of the Einstein field equations has
played a major role in General Relativity. It describes the effects of
gravitation in the vacuum outside a time-independent spherically symmetric
distribution of matter. One of the reasons for its importance is its
generality — it only depends on the spherical symmetry and on the total energy
of the matter distribution. Jebsen and Birkhoff [3, 4] have shown that this
solution is still valid even in time-dependent situations, provided that the
spherical symmetry is preserved. Another reason for its popularity is the
association of the coordinate singularity of this solution, which occurs for a
certain value of the radial coordinate, with the presence of an event horizon,
thus leading to the concept of black holes.

Less known — even absent in many standard textbooks on General Relativity — is
the interior Schwarzschild solution [5, 2]. It gives the metric of the space
inside a spherically symmetric matter distribution with an energy density
which is constant with the radial coordinate. This other solution can be
continuously joined with the Schwarzschild vacuum solution that is valid
outside the matter distribution. It is less general in that it only describes
matter distributions with energy densities that do not depend on the radial
coordinate $r$. In addition, it does not contain any singularities. This point
is emphasized in many texts, for example in [2, 6]. Basically, in order to
avoid singularities at the center of the matter distribution a certain
integration constant is set equal to zero.

For a spherical matter shell characterized by an inner radius $r_{1}$, an
outer radius $r_{2}$ and an energy density constant with $r$ the situation is
more involved. In the inner vacuum region, where $r<r_{1}$, the solution of
the Einstein equations leads to an integration constant, heretofore denoted by
$r_{\mu}$, which determines the singularities in the entire inner vacuum
region. There are no singularities only if $r_{\mu}=0$. In analogy with what
is done for the interior Schwarzschild solution one may feel tempted to set
$r_{\mu}=0$ by hand and eliminate all singularities. However, as we are going
to show in this paper, the correct approach is to start in the outer vacuum
region ($r>r_{2}$), where the Schwarzschild external solution holds, and use
the continuity of the solution in the two boundaries of the three regions to
determine the constant $r_{\mu}$. The rather surprising result is that the
imposition of the surface boundary conditions implies that $r_{\mu}>0$, so
that the solutions do contain a singularity at the origin. In addition, one
can prove that this condition has to be satisfied in order to produce
solutions with non-negative pressure inside the matter shell.

It is remarkable that the boundary conditions on matter interfaces for the
Einstein field equations seem to play a smaller than expected role in the
literature. A rare example in which the role of these boundary conditions is
emphasized can be found in [7], although the author of that paper only
obtained solutions containing a negative pressure region inside the matter
shell. By analyzing these negative pressure solutions the author concluded
that matter cannot collapse towards the center of black holes in general
relativity. We are going to show in this paper that it is possible to obtain
physically reasonable matter shell solutions of the Einstein equations with
non-negative and finite pressure inside the shell. It is important to
emphasize that the singularity at the origin in the inner vacuum region does
not lead to any divergence of the matter quantities, and in fact stabilizes
the matter shell structure. This is so because the gravitational field within
the inner vacuum region turns out to be repulsive with respect to the origin.
Our solutions for matter shells are expressed in terms of a single integral
which for some values of the physical parameters can be written in terms of
elementary functions and constitute a new class of exact solutions of the
Einstein field equations.

Results similar to the ones we present here were obtained for the case of
neutron stars, with the Chandrasekhar equation of state [8], by Ni [9],
including the presence of inner and outer matter-vacuum interfaces. However,
the crucial consideration of the interface boundary conditions was missing
from that analysis, thus leading to incomplete results. The discussion of the
interface boundary conditions was subsequently introduced by Neslušan [10],
thus completing the analysis of the case of the neutron stars. Just as in the
present work, the discussion of the interface boundary conditions led, also in
that case, to an inner vacuum region containing a singularity at the origin
and a gravitational field pointing away from the origin, that is, repulsive
with respect to the origin. The present work can be considered as an exactly
solvable laboratory model that illustrates some of the properties of that
solution. It also shows that the properties of the inner vacuum region are not
artifacts of that particular problem or of that particular type of equation of
state.

This paper is organized as follows. In Section 2 we state and solve the
problem; in Section 3 we derive the main physical properties of the solution;
in Section 4 we present a two-parameter family of explicit solutions and a few
numerical examples; and in Section 5 we present our conclusions.

##  2 The Problem and its Solution

We will present, in the case of a spherically symmetric shell of liquid fluid
with constant energy density, the exact solution of the Einstein field
equations of General Relativity [11],

$R_{\mu}^{\;\nu}-\frac{1}{2}\,R\,g_{\mu}^{\;\nu}=-\kappa\,T_{\mu}^{\;\nu},$
(1)

where $\kappa=8\pi G/c^{4}$, $G$ is the universal gravitational constant and
$c$ is the speed of light. Under the conditions of time independence and of
spherical symmetry around the origin of a spherical system of coordinates
$(t,r,\theta,\phi)$, the Schwarzschild system of coordinates, the most general
possible metric is given by the invariant interval, written in terms of this
spherical system of coordinates,

$ds^{2}=\,{\rm e}^{2\nu(r)}c^{2}dt^{2}-\,{\rm
e}^{2\lambda(r)}dr^{2}-r^{2}\left[d\theta^{2}+\sin^{2}(\theta)d\phi^{2}\right],$
(2)

where $\exp[\nu(r)]$ and $\exp[\lambda(r)]$ are two positive functions of only
$r$. As one can see, in this work we will use the time-like signature
$(+,-,-,-)$, following [11]. Under these conditions the matter energy-momentum
tensor density $T_{\mu}^{\;\nu}$ on the right-hand side of the equation is
diagonal, and given by the four diagonal components $T_{0}^{\;0}(r)=\rho(r)$,
where $\rho(r)$ is the energy density of the matter, and
$T_{1}^{\;1}(r)=T_{2}^{\;2}(r)=T_{3}^{\;3}(r)=-P(r)$, where $P(r)$ is the
pressure, which is isotropic, thus characterizing a fluid.

Since under these conditions $R_{\mu}^{\;\nu}$ and $T_{\mu}^{\;\nu}$ are both
diagonal, there are just four non-trivial field equations contained in
Equation (1). In addition to these four field equations we have the
consistency condition

$D_{\nu}T_{\mu}^{\;\nu}=0,$ (3)

which is due to the fact that the combination of tensors that constitutes the
left-hand side of the Einstein field equation satisfies the Bianchi identity
of the Ricci curvature tensor. Writing these equations explicitly in the
chosen coordinate system, one finds that the component equations involving
$T_{2}^{\;2}(r)$ and $T_{3}^{\;3}(r)$ turn out to be identical, so that we are
left with the set of four equations, including the consistency condition,

$\displaystyle\left\\{\rule{0.0pt}{12.91663pt}1-2\left[r\lambda^{\prime}(r)\right]\right\\}\,{\rm
e}^{-2\lambda(r)}$ $\displaystyle=$ $\displaystyle 1-\kappa r^{2}\rho(r),$ (4)
$\displaystyle\left\\{\rule{0.0pt}{12.91663pt}1+2\left[r\nu^{\prime}(r)\right]\right\\}\,{\rm
e}^{-2\lambda(r)}$ $\displaystyle=$ $\displaystyle 1+\kappa r^{2}P(r),$ (5)
$\displaystyle\left\\{\rule{0.0pt}{12.91663pt}r^{2}\nu^{\prime\prime}(r)-\left[r\lambda^{\prime}(r)\right]\left[r\nu^{\prime}(r)\right]\right.\hskip
70.0001pt$
$\displaystyle\left.\rule{0.0pt}{12.91663pt}+\left[r\nu^{\prime}(r)\right]^{2}+\left[r\nu^{\prime}(r)\right]-\left[r\lambda^{\prime}(r)\right]\right\\}\,{\rm
e}^{-2\lambda(r)}$ $\displaystyle=$ $\displaystyle\kappa r^{2}P(r),$ (6)
$\displaystyle\left[\rho(r)+P(r)\right]\nu^{\prime}(r)$ $\displaystyle=$
$\displaystyle-P^{\prime}(r),$ (7)

where the primes indicate differentiation with respect to $r$. Next, it can be
shown that Equation (2) can be obtained from the others, being in fact a
linear combination of the derivative of Equation (5) and of Equations (4), (5)
and (7). If we denote Equations (4) through (7) respectively by $E_{t}$,
$E_{r}$, $E_{\theta}$ and $E_{c}$, we have that

$E_{\theta}=\frac{1}{2}\left[-r\nu^{\prime}(r)\left(E_{t}-E_{r}\right)+rE^{\prime}_{r}+\kappa
r^{2}E_{c}\right].$ (8)

This leaves us with a set of just three differential equations to solve. In
addition to this, we will assume that we have an energy density
$\rho(r)=\rho_{0}$ which is constant as a function of $r$ within the shell of
fluid matter, thus characterizing a liquid fluid. The equations that we
propose to solve are therefore those given in Equations (4), (5) and (7). It
is important to note that, in this way, we are left with a system of just
three first-order differential equations. Therefore, the discussion of
boundary conditions can be limited to the discussion of the behavior of the
functions involved, thus eliminating the need for any discussion of the
behavior of their derivatives.

We will assume that the matter consists of a spherical shell of liquid,
located between the radial positions $r_{1}$ and $r_{2}$, meaning that we will
have an inner vacuum region within $(0,r_{1})$, a matter region within
$(r_{1},r_{2})$, and an outer vacuum region within $(r_{2},\infty)$. This
means that we will have for $\rho(r)$ and $P(r)$

$\displaystyle\rho(r)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{lcl}0&\mbox{for}&0\;\leq
r<r_{1},\\\\[12.91663pt] \rho_{0}&\mbox{for}&r_{1}<r<r_{2},\\\\[12.91663pt]
0&\mbox{for}&r_{2}<r<\infty,\end{array}\right.$ (12) $\displaystyle P(r)$
$\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{lcl}0&\mbox{for}&0\;\leq r\leq
r_{1},\\\\[12.91663pt] 0&\mbox{for}&r_{2}\leq r<\infty.\end{array}\right.$
(15)

The function $P(r)$ within the matter region is, of course, one of the
unknowns of our problem. In addition to this, we have the boundary conditions
for $P(r)$ at the two interfaces, in the limits coming from within the liquid,

$\displaystyle P(r_{1})$ $\displaystyle=$ $\displaystyle 0,$ $\displaystyle
P(r_{2})$ $\displaystyle=$ $\displaystyle 0,$ (16)

since these constitute a requirement in any interface between fluid matter and
a vacuum. The remaining boundary conditions are those requiring the continuity
of $\lambda(r)$ and $\nu(r)$ across the interfaces, and the asymptotic
conditions leading to the Newtonian limit at radial infinity.

###  2.1 Solutions in the Vacuum Regions

Within either vacuum region the consistency condition in Equation (7) becomes
a mere identity, so that we are left with only two equations, in which we
replace both $\rho(r)$ and $P(r)$ by zero,

$\displaystyle 1-2\left[r\lambda^{\prime}(r)\right]$ $\displaystyle=$
$\displaystyle\,{\rm e}^{2\lambda(r)},$ $\displaystyle
1+2\left[r\nu^{\prime}(r)\right]$ $\displaystyle=$ $\displaystyle\,{\rm
e}^{2\lambda(r)}.$ (17)

This immediately implies that $\lambda^{\prime}(r)+\nu^{\prime}(r)=0$, and
hence that $\lambda(r)+\nu(r)=A$, where $A$ is a dimensionless integration
constant. The first of these two equations involves only $\lambda(r)$, and can
also be written as

$\left[r\,{\rm e}^{-2\lambda(r)}\right]^{\prime}=1,$ (18)

which can be immediately integrated to

$\,{\rm e}^{-2\lambda(r)}=1-\frac{R}{r},$ (19)

where $R$ is an integration constant with dimensions of length.

We must now discriminate between the inner and outer vacuum regions. In the
outer vacuum region we must get flat space at radial infinity, which requires
that both $\lambda(r)$ and $\nu(r)$ go to zero for $r\to\infty$. This in turn
implies that $A=0$ in the outer vacuum region, thus leading to
$\nu(r)=-\lambda(r)$. As is well known, the condition that the Newtonian limit
be realized at radial infinity requires that $R=r_{M}$, the Schwarzschild
radius $r_{M}=2MG/c^{2}$ associated to the asymptotic gravitational mass $M$
of the system. Thus we arrive at the time-honored Schwarzschild solution [1,
2] in the outer vacuum region,

$\displaystyle\lambda_{s}(r)$ $\displaystyle=$
$\displaystyle-\,\frac{1}{2}\,\ln\\!\left(\frac{r-r_{M}}{r}\right),$
$\displaystyle\nu_{s}(r)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\ln\\!\left(\frac{r-r_{M}}{r}\right),$ (20)

where the subscript $s$ denotes the outer vacuum region. Note that there is a
limitation on the values of the parameters $r_{2}$ and $r_{M}$ describing the
distribution of matter, because these expressions have a singular behavior at
$r=r_{M}$. We must have $r_{M}<r_{2}$ to ensure that there is no event horizon
formed outside the distribution of matter.

In the inner vacuum region there are no asymptotic conditions to be applied,
and thus the integration constants $A$ and $R$ will have to be left
undetermined, to be determined later on via the boundary conditions at the
interfaces between the vacuum and the matter, as we come in from radial
infinity towards the origin. For convenience we will put $R=-r_{\mu}$, and
write the solution in the inner vacuum region as

$\displaystyle\lambda_{i}(r)$ $\displaystyle=$
$\displaystyle-\,\frac{1}{2}\,\ln\\!\left(\frac{r+r_{\mu}}{r}\right),$
$\displaystyle\nu_{i}(r)$ $\displaystyle=$ $\displaystyle
A+\frac{1}{2}\,\ln\\!\left(\frac{r+r_{\mu}}{r}\right),$ (21)

where the subscript $i$ denotes the inner vacuum region. Note that the value
of $r_{\mu}$ determines the singularity structure of this solution within the
inner vacuum region. If $r_{\mu}<0$ then there is a singularity at the
strictly positive radial position $r=-r_{\mu}$, corresponding to the formation
of an event horizon at that position. If $r_{\mu}=0$ then there are no
singularities at all within this region. If $r_{\mu}>0$ then there is only one
singularity, located at the origin $r=0$. We will show later on that we do
indeed have that $r_{\mu}>0$.

We therefore have the complete analytical solutions in the inner and outer
vacuum regions, which contain one input parameter of the problem, the mass $M$
associated to the Schwarzschild radius $r_{M}$, and two integration constants
still to be determined, $A$ and $r_{\mu}$.

###  2.2 Solution in the Matter Region

In the matter region Equation (4) for $\lambda(r)$ can be written as

$\left[r\,{\rm e}^{-2\lambda(r)}\right]^{\prime}=1-\kappa\rho_{0}r^{2},$ (22)

which can be immediately integrated to

$\,{\rm e}^{-2\lambda(r)}=1+\frac{B}{r}-\frac{\kappa\rho_{0}}{3}\,r^{2},$ (23)

where $B$ is an integration constant with dimensions of length, thus leading
to the general solution for $\lambda(r)$ in the matter region,

$\lambda_{m}(r)=-\,\frac{1}{2}\,\ln\\!\left(1+\frac{B}{r}-\frac{\kappa\rho_{0}}{3}\,r^{2}\right),$
(24)

where the subscript $m$ denotes the matter region. This solution contains one
integration constant, the constant $B$, and one parameter characterizing the
system, namely $\rho_{0}$, which is not, however, a free input parameter of
the problem, since it will depend on $M$ and thus on $r_{M}$.

In order to deal with $\nu(r)$ in the matter region, we consider the
consistency condition given in Equation (7), which can be written in this
region as

$\nu^{\prime}(r)=-\,\frac{P^{\prime}(r)}{\rho_{0}+P(r)},$ (25)

thus allowing us to separate variables and hence to write $\nu(r)$ in terms of
$P(r)$,

$\displaystyle d\nu$ $\displaystyle=$ $\displaystyle-\,\frac{dP}{\rho_{0}+P}$
(26) $\displaystyle=$ $\displaystyle-d\ln\\!\left(\rho_{0}+P\right).$

If we integrate from the left end $r_{1}$ of the matter interval to a generic
point $r$ within that interval, we get

$\nu(r)-\nu(r_{1})=-\ln\\!\left[\frac{\rho_{0}+P(r)}{\rho_{0}+P(r_{1})}\right].$
(27)

However, the boundary conditions for $P(r)$ at the interfaces tell us that we
must have $P(r_{1})=0$, and hence we get the general solution for $\nu(r)$
within the matter region, written in terms of $P(r)$,

$\nu_{m}(r)=\nu_{1}-\ln\\!\left[\frac{\rho_{0}+P(r)}{\rho_{0}}\right],$ (28)

where $\nu_{1}=\nu(r_{1})$. The solutions for $\lambda(r)$ and $\nu(r)$ within
the matter region involve therefore two integration constants, $B$ and
$\nu_{1}$. The solution for $\nu(r)$ is not yet completely determined, since
it is given in terms of $P(r)$, which is also as yet undetermined. However,
the information obtained so far already allows us to impose the boundary
conditions at the interfaces, in order to determine the integration constants,
which is what we turn to now.

###  2.3 Interface Boundary Conditions

The condition of the continuity of $\lambda(r)$ at the interface $r_{1}$
implies that we must have that $\lambda_{i}(r_{1})=\lambda_{m}(r_{1})$, which
from Equations (2.1) and (24) gives us the following relation between the
parameters

$B-r_{\mu}=\frac{\kappa\rho_{0}}{3}\,r_{1}^{3}.$ (29)

In addition to this, the condition of the continuity of $\lambda(r)$ at the
interface $r_{2}$ implies that we must have
$\lambda_{m}(r_{2})=\lambda_{s}(r_{2})$, which from Equations (2.1) and (24)
gives us the following relation between the parameters

$B+r_{M}=\frac{\kappa\rho_{0}}{3}\,r_{2}^{3}.$ (30)

This last condition already determines the integration constant $B$ in terms
of the parameters of the problem,

$B=-r_{M}+\frac{\kappa\rho_{0}}{3}\,r_{2}^{3},$ (31)

and the difference of the two conditions just obtained determines the
integration parameter $r_{\mu}$ in terms of the parameters of the problem,

$r_{\mu}=-r_{M}+\frac{\kappa\rho_{0}}{3}\left(r_{2}^{3}-r_{1}^{3}\right).$
(32)

We have therefore the solution for $\lambda(r)$ in the matter region, in terms
of the parameters of the problem,

$\lambda_{m}(r)=-\,\frac{1}{2}\,\ln\\!\left[\frac{\kappa\rho_{0}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)}{3r}\right].$
(33)

Let us point out that there is a consistency condition to be applied to this
result, since we must have that the cubic polynomial appearing in the argument
of the logarithm be strictly positive for all values of $r$ within the matter
region, that is

$\kappa\rho_{0}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)>0,$ (34)

for all $r\in[r_{1},r_{2}]$. Note that the term with the cubes is necessarily
non-negative, but that the other term may be negative, if $r_{M}$ is not
smaller than $r_{1}$. Therefore, so long as $r_{M}<r_{1}$, this strict
positivity condition is automatically satisfied. If, however, we have that
$r_{1}<r_{M}<r_{2}$, then the condition must be actively verified for all
$r\in[r_{M},r_{2}]$. If it fails, then there is no solution for that
particular set of input parameters.

Since we have $\nu_{m}(r)$ written in terms of $P(r)$, and since we know the
interface boundary conditions for $P(r)$ in limits from within the matter
region, we are in a position to impose the boundary conditions on $\nu(r)$
across the interfaces, even without having available the complete solution for
$\nu_{m}(r)$. To this end, let us note that from Equation (28) we have that
$\nu_{m}(r_{1})=\nu_{m}(r_{2})=\nu_{1}$. At the interface $r_{1}$ the
condition of the continuity of $\nu(r)$ implies that we must have
$\nu_{i}(r_{1})=\nu_{m}(r_{1})$, which from Equations (2.1) and (28) gives us
the following relation between the parameters,

$\nu_{1}=A+\frac{1}{2}\,\ln\\!\left(\frac{r_{1}+r_{\mu}}{r_{1}}\right).$ (35)

In addition to this, the condition of the continuity of $\nu(r)$ at the
interface $r_{2}$ implies that we must have $\nu_{m}(r_{2})=\nu_{s}(r_{2})$,
which from Equations (2.1) and (28) gives us the following relation between
the parameters,

$\nu_{1}=\frac{1}{2}\,\ln\\!\left(\frac{r_{2}-r_{M}}{r_{2}}\right).$ (36)

This last condition gives us the integration constant $\nu_{1}$ in terms of
the parameters of the problem, and its difference with the previous one
determines the integration constant $A$,

$A=\frac{1}{2}\,\ln\\!\left(\frac{r_{1}}{r_{2}}\,\frac{r_{2}-r_{M}}{r_{1}+r_{\mu}}\right).$
(37)

This completes the determination of the solution for both $\nu(r)$ and
$\lambda(r)$ in the inner vacuum region, for which we now have

$\displaystyle\lambda_{i}(r)$ $\displaystyle=$
$\displaystyle-\,\frac{1}{2}\,\ln\\!\left(\frac{r+r_{\mu}}{r}\right),$
$\displaystyle\nu_{i}(r)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\ln\\!\left(\frac{r_{1}}{r_{2}}\,\frac{r_{2}-r_{M}}{r_{1}+r_{\mu}}\right)+\frac{1}{2}\,\ln\\!\left(\frac{r+r_{\mu}}{r}\right),$
(38)

with $r_{\mu}$ given by Equation (32). We also have the following form for the
solution for $\nu(r)$ within the matter region, still in terms of $P(r)$,

$\nu_{m}(r)=\frac{1}{2}\,\ln\\!\left(\frac{r_{2}-r_{M}}{r_{2}}\right)-\ln\\!\left[\frac{\rho_{0}+P(r)}{\rho_{0}}\right].$
(39)

At this point the situation is as follows, in regard to the complete solution
of the problem. Given values of $r_{1}$, $r_{2}$ and $r_{M}$, which completely
characterize the geometrical and physical nature of the object under study, we
have the complete solution for both $\lambda(r)$ and $\nu(r)$ in the outer
vacuum region. We also have the complete solution for both $\lambda(r)$ and
$\nu(r)$ in the inner vacuum region, except for the determination of the
parameter $\rho_{0}$. We have as well the complete solution for $\lambda(r)$
in the matter region, again up to the determination of the parameter
$\rho_{0}$. The one element of the solution still missing is the complete
solution for $\nu(r)$ in the matter region. However, since we have $\nu(r)$
determined in terms of $P(r)$ in this region, this can also be accomplished by
the complete determination of $P(r)$ in this region, which is the task we
tackle next. Let us emphasize that the parameter $\rho_{0}$ is not a free
input parameter of the problem, since it must be chosen so that the given
value of $r_{M}$ results, that is, the local value of the energy density must
be chosen so that the given value of the asymptotic gravitational mass $M$
results at radial infinity.

###  2.4 The Equation for the Pressure

The equation determining the pressure $P(r)$ in the matter region can be
obtained by eliminating $\nu^{\prime}(r)$ from Equations (5) and (7), which
gives us

$\rho_{0}+P(r)-2\left[rP^{\prime}(r)\right]=\,{\rm
e}^{2\lambda_{m}(r)}\left[1+\kappa
r^{2}P(r)\right]\left[\rho_{0}+P(r)\right].$ (40)

In this equation the quantity $\exp[2\lambda_{m}(r)]$ is a known function,
since we have already determined $\lambda(r)$ in the matter region. This is a
first-order non-linear differential equation determining $P(r)$, with the
boundary conditions $P(r_{1})=0$ and $P(r_{2})=0$. Since the equation is
first-order and there are two boundary conditions to be satisfied, it is clear
that the parameter $\rho_{0}$ will have to be adjusted so that the second
condition can be satisfied. This will therefore determine the parameter
$\rho_{0}$ in terms of the other parameters of the problem. This equation can
be simplified by a series of transformations on the variables and parameters.
First we define the parameter $\Upsilon_{0}$, which has dimensions of inverse
length and is such that

$\Upsilon_{0}^{2}=\kappa\rho_{0},$ (41)

and the dimensionless pressure

$p(r)=\frac{P(r)}{\rho_{0}},$ (42)

in terms of which Equation (40) becomes

$\left[rp^{\prime}(r)\right]=\frac{1}{2}\left[1+p(r)\right]\left\\{1-\,{\rm
e}^{2\lambda_{m}(r)}\left[1+\Upsilon_{0}^{2}r^{2}p(r)\right]\right\\}.$ (43)

Substituting the known value of $\lambda_{m}(r)$ from Equation (33) we get

$p^{\prime}(r)=\frac{1}{2r}\left[1+p(r)\right]\frac{\Upsilon_{0}^{2}\left(r_{2}^{3}-r^{3}\right)-3r_{M}-3\Upsilon_{0}^{2}r^{3}p(r)}{\Upsilon_{0}^{2}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)}.$
(44)

This has the form of a Riccati equation, and can be linearized by the
transformation of variables

$p(r)=\frac{1}{z(r)}-1,$ (45)

thus resulting in the equation for $z(r)$,

$z^{\prime}(r)+\frac{\Upsilon_{0}^{2}\left(r_{2}^{3}+2r^{3}\right)-3r_{M}}{2r\left[\Upsilon_{0}^{2}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)\right]}\,z(r)=\frac{3\Upsilon_{0}^{2}r^{3}}{2r\left[\Upsilon_{0}^{2}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)\right]}.$
(46)

This equation has an integrating factor given by $\exp[F(r)]$, where $F(r)$ is
defined as an integral of the coefficient of the second term from $r_{2}$ to
some arbitrary $r$ within $[r_{1},r_{2}]$,

$\displaystyle F(r)$ $\displaystyle=$
$\displaystyle\int_{r_{2}}^{r}ds\,\frac{\Upsilon_{0}^{2}\left(r_{2}^{3}+2s^{3}\right)-3r_{M}}{2s\left[\Upsilon_{0}^{2}\left(r_{2}^{3}-s^{3}\right)+3\left(s-r_{M}\right)\right]}$
(47) $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{r_{2}}^{r}ds\,\frac{1}{s}-\frac{1}{2}\int_{r_{2}}^{r}ds\,\frac{-3\Upsilon_{0}^{2}s^{2}+3}{\Upsilon_{0}^{2}\left(r_{2}^{3}-s^{3}\right)+3\left(s-r_{M}\right)}.$

One can see now that both integrals can be done, and thus we obtain

$e^{F(r)}=\sqrt{\frac{r}{r_{2}}}\;\sqrt{\frac{3\left(r_{2}-r_{M}\right)}{\Upsilon_{0}^{2}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)}},$
(48)

in terms of which the equation for $z(r)$ can be written as

$\left[\,{\rm
e}^{F(r)}z(r)\right]^{\prime}=\frac{3}{2}\,\frac{\Upsilon_{0}^{2}r^{2}\,{\rm
e}^{F(r)}}{\Upsilon_{0}^{2}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)},$
(49)

which can then be integrated over the interval $[r_{2},r]$ giving

$z(r)=e^{-F(r)}+\frac{3}{2}\,e^{-F(r)}\\!\int_{r_{2}}^{r}ds\,\frac{\Upsilon_{0}^{2}s^{2}e^{F(s)}}{\Upsilon_{0}^{2}\left(r_{2}^{3}-s^{3}\right)+3\left(s-r_{M}\right)},$
(50)

where we used the fact that by definition $F(r_{2})=0$, and the fact that
$P(r_{2})=0$ implies $z(r_{2})=1$.

Note that once more the existence of the solutions for $F(r)$ and for $z(r)$
is conditioned by the strict positivity of the same cubic polynomial that we
discussed before in Equation (34), which we can now write as

$\Upsilon_{0}^{2}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)>0,$ (51)

for all $r\in[r_{1},r_{2}]$. Substituting the value of $\exp[F(r)]$ we have
the solution for $z(r)$ written in terms of a real integral,

$\displaystyle z(r)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\Upsilon_{0}^{2}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)}{r}}$
(52)
$\displaystyle\times\left\\{\sqrt{\frac{r_{2}}{3\left(r_{2}-r_{M}\right)}}+\frac{3}{2}\,\int_{r_{2}}^{r}ds\,\frac{\Upsilon_{0}^{2}s^{5/2}}{\left[\Upsilon_{0}^{2}\left(r_{2}^{3}-s^{3}\right)+3\left(s-r_{M}\right)\right]^{3/2}}\right\\}.$

In most cases this remaining integral is elliptic and therefore cannot be
written in terms of elementary functions, so that in general this remaining
last step of the resolution procedure has to be performed by numerical means.
However, as we are going to show in Section 4, for some values of the
parameters it is possible to express this integral in terms of elementary
functions.

After determining $z(r)$ in the matter region, Equations (45) allows us to
calculate the dimensionless pressure $p(r)$ which, according to Equation (42),
is equal to the pressure divided by the energy density $\rho_{0}$,

$\displaystyle p(r)$ $\displaystyle=$
$\displaystyle\frac{1}{z(r)}-1\;\;\;\Longrightarrow$ $\displaystyle P(r)$
$\displaystyle=$ $\displaystyle\frac{\rho_{0}}{z(r)}-\rho_{0}.$ (53)

Note that $z(r)$ also determines $\nu(r)$ in the matter region, since in
Equation (39) we have $\nu_{m}(r)$ in terms of $P(r)$, and therefore we have
for the exponential of $\nu_{m}(r)$,

$\,{\rm
e}^{\nu_{m}(r)}=\sqrt{\frac{r_{2}-r_{M}}{r_{2}}}\,\frac{\rho_{0}}{\rho_{0}+P(r)},$
(54)

which, using Equation (2.4), implies that

$\,{\rm e}^{\nu_{m}(r)}=\sqrt{\frac{r_{2}-r_{M}}{r_{2}}}\,z(r),$ (55)

so that, up to a constant factor, $z(r)$ turns out to be the square root of
the temporal coefficient of the metric. This completes the determination of
the solution in all three regions, in terms of the parameters of the problem.
Given certain values of $r_{1}$, $r_{2}$ and $r_{M}$, one must still find a
value of the parameter $\rho_{0}$, and hence of $\Upsilon_{0}$, such that the
boundary conditions for $P(r)$ at the two interfaces are satisfied. One can
obtain an equation determining this value of $\Upsilon_{0}$ by considering the
value of $z(r_{1})$. Since $P(r_{1})=0$, we have that $z(r_{1})=1$, so that
from Equation (52) we get

$\displaystyle\sqrt{\frac{r_{2}}{3\left(r_{2}-r_{M}\right)}}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{r_{1}}{\Upsilon_{0}^{2}\left(r_{2}^{3}-r_{1}^{3}\right)+3\left(r_{1}-r_{M}\right)}}$
(56)
$\displaystyle+\frac{3}{2}\int_{r_{1}}^{r_{2}}dr\,\frac{\Upsilon_{0}^{2}r^{5/2}}{\left[\Upsilon_{0}^{2}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)\right]^{3/2}}.$

The solution of this algebraic equation gives the value of $\Upsilon_{0}$, and
hence the value of $\rho_{0}$, for which the two interface boundary conditions
for $P(r)$ will be satisfied. The solution of this equation necessarily
includes the consistency check of the solution obtained, since the calculation
of the integral is dependent on the strict positivity of the polynomial in
Equation (51), for all $r$ within $[r_{1},r_{2}]$. This is the same condition
that guarantees the consistency of the results for $F(r)$ and $z(r)$, and
hence the consistency of the results for $P(r)$ and $\nu(r)$ within the matter
region.

##  3 Main Properties of the Solution

In this section we will state and prove a few important properties of the
solution. We will assume that, given certain values of $r_{1}$, $r_{2}$ and
$r_{M}$, the corresponding solution exists. In other words, we are assuming
that a solution of Equation (56) for $\Upsilon_{0}$ can be found, thus
determining $\rho_{0}$, which includes establishing the strict positivity of
the cubic polynomial within the square roots in the denominators, and that a
corresponding function $z(r)$ is therefore determined via Equation (52). This
then implies that the solutions for both $\lambda(r)$ and $\nu(r)$, as well as
for $P(r)$, are all determined, with all the boundary conditions duly
satisfied. A simpler way to put this is to say that we are establishing the
most important properties of all existing solutions of the problem. For easy
reference, we state the complete solution explicitly in Table 1, where we have
that $\rho_{0}$ is determined algebraically via Equation (56), $z(r)$ is
determined by Equation (52), and $r_{\mu}$ is given by Equation (32). We will
start by the discussion of the presence of the singularity at the origin.

Table 1: Summary of the results.

$\displaystyle\lambda(r)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{lcl}-\,{\displaystyle\frac{\displaystyle
1}{\displaystyle 2}}\,\ln\\!\left({\displaystyle\frac{\displaystyle
r+r_{\mu}}{\displaystyle r}}\right)&\mbox{for}&0\;\leq r\leq
r_{1},\\\\[12.91663pt] -\,{\displaystyle\frac{\displaystyle 1}{\displaystyle
2}}\,\ln\\!\left[{\displaystyle\frac{\displaystyle\kappa\rho_{0}\left(r_{2}^{3}-r^{3}\right)+3\left(r-r_{M}\right)}{\displaystyle
3r}}\right]&\mbox{for}&r_{1}\leq r\leq r_{2},\\\\[12.91663pt]
-\,{\displaystyle\frac{\displaystyle 1}{\displaystyle
2}}\,\ln\\!\left({\displaystyle\frac{\displaystyle r-r_{M}}{\displaystyle
r}}\right)&\mbox{for}&r_{2}\leq r<\infty,\end{array}\right.$
$\displaystyle\nu(r)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{lcl}{\displaystyle\frac{\displaystyle
1}{\displaystyle 2}}\,\ln\\!\left({\displaystyle\frac{\displaystyle
r_{1}}{\displaystyle r_{2}}}\,{\displaystyle\frac{\displaystyle
r_{2}-r_{M}}{\displaystyle
r_{1}+r_{\mu}}}\right)+{\displaystyle\frac{\displaystyle 1}{\displaystyle
2}}\,\ln\\!\left({\displaystyle\frac{\displaystyle r+r_{\mu}}{\displaystyle
r}}\right)&\mbox{for}&0\;\leq r\leq r_{1},\\\\[12.91663pt]
{\displaystyle\frac{\displaystyle 1}{\displaystyle
2}}\,\ln\\!\left({\displaystyle\frac{\displaystyle r_{2}-r_{M}}{\displaystyle
r_{2}}}\right)+\ln\\!\left[z(r)\right]&\mbox{for}&r_{1}\leq r\leq
r_{2},\\\\[12.91663pt] {\displaystyle\frac{\displaystyle 1}{\displaystyle
2}}\,\ln\\!\left({\displaystyle\frac{\displaystyle r-r_{M}}{\displaystyle
r}}\right)&\mbox{for}&r_{2}\leq r<\infty.\end{array}\right.$

###  3.1 Existence of the Singularity at the Origin

The existence of the singularity at the origin is equivalent to the statement
that $r_{\mu}>0$, because the only way to avoid that singularity would be to
have $r_{\mu}=0$. If we put $r_{\mu}=0$ and take the limit $r_{1}\rightarrow
0$ we no longer have a matter shell, and we obtain instead the Schwarzschild
interior solution.

We start with a preliminary lemma, in which we will prove that the following
combination of parameters

$\frac{1}{3}\,\Upsilon_{0}^{2}\left(r_{2}^{3}-r_{e}^{3}\right)-r_{M}>0,$ (59)

is strictly positive, where $r_{e}$ is the position of the maximum of the
dimensionless pressure $p(r)$ within the interval $[r_{1},r_{2}]$. In order to
do this, we consider the equation for $p(r)$ given in Equation (44). Applying
that equation at $r_{2}$, since we have that $p(r_{2})=0$, we get for the
derivative at the right end of the matter interval,

$p^{\prime}(r_{2})=-\,\frac{r_{M}}{2r_{2}\left(r_{2}-r_{M}\right)}.$ (60)

Since by hypothesis we have that $r_{2}>r_{M}$ and that $r_{M}>0$, we conclude
that the derivative $p^{\prime}(r_{2})$ is strictly negative. In addition to
this, since $p(r)$ is a positive function that is the solution of a first-
order differential equation within $(r_{1},r_{2})$, it must be a continuous
and differentiable function. Therefore, given that it is zero at both ends and
always increases as we go to the interior of the interval, it must have a
point of maximum $r_{e}$ somewhere in the interior of the interval, where we
will have that $p^{\prime}(r_{e})=0$. Using the differential equation for
$p(r)$ given by Equation (44) at this point we thus obtain

$\frac{1}{2r_{e}}\left[1+p(r_{e})\right]\frac{\Upsilon_{0}^{2}\left(r_{2}^{3}-r_{e}^{3}\right)-3r_{M}-3\Upsilon_{0}^{2}r_{e}^{3}p(r_{e})}{\Upsilon_{0}^{2}\left(r_{2}^{3}-r_{e}^{3}\right)+3\left(r_{e}-r_{M}\right)}=0.$
(61)

This can only be zero if the numerator is zero, so we have that

$\Upsilon_{0}^{2}r_{e}^{3}p(r_{e})=\frac{1}{3}\,\Upsilon_{0}^{2}\left(r_{2}^{3}-r_{e}^{3}\right)-r_{M}.$
(62)

Since $\Upsilon_{0}^{2}>0$ and at its maximum we must have $p(r_{e})>0$ for
the dimensionless pressure, we conclude that our lemma holds,

$\frac{1}{3}\,\Upsilon_{0}^{2}\left(r_{2}^{3}-r_{e}^{3}\right)-r_{M}>0.$ (63)

Let us now consider the result for $r_{\mu}$ in terms of the given parameters
of the problem, as shown in Equation (32), which we can write as

$r_{\mu}=\frac{1}{3}\,\Upsilon_{0}^{2}\left(r_{2}^{3}-r_{1}^{3}\right)-r_{M}.$
(64)

By adding and subtracting terms to this equation, we can write it as

$r_{\mu}=\left[\frac{1}{3}\,\Upsilon_{0}^{2}\left(r_{2}^{3}-r_{e}^{3}\right)-r_{M}\right]+\frac{1}{3}\,\Upsilon_{0}^{2}\left(r_{e}^{3}-r_{1}^{3}\right).$
(65)

The quantity within square brackets is the one we just proved to be strictly
positive in our lemma. The other term is also strictly positive because we
certainly have that $r_{e}>r_{1}$. Therefore, we have our theorem,

$r_{\mu}>0.$ (66)

Therefore, every solution of the problem that exists at all is bound to have a
singularity at the origin, which is characterized by the factor

$\ln\\!\left(\frac{r+r_{\mu}}{r}\right),$ (67)

that appears with a negative sign in $\lambda_{i}(r)$ and with a positive sign
in $\nu_{i}(r)$. This implies that at this singular point we have that

$\displaystyle\lim_{r\to 0}\lambda_{i}(r)$ $\displaystyle=$
$\displaystyle-\infty,$ $\displaystyle\lim_{r\to 0}\,{\rm e}^{\lambda_{i}(r)}$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\lim_{r\to 0}\nu_{i}(r)$
$\displaystyle=$ $\displaystyle\infty,$ $\displaystyle\lim_{r\to 0}\,{\rm
e}^{\nu_{i}(r)}$ $\displaystyle=$ $\displaystyle\infty.$ (68)

Note that this singularity does not have any disastrous consequences, since it
does not imply infinite concentrations of matter. In fact, we have $\rho(r)=0$
in the whole inner vacuum region, including at the origin. For the proper
lengths in the radial direction, it just implies that they get progressively
more contracted as we approach the origin, rather than being expanded with
respect to the corresponding variations of the radial coordinate $r$, as is
the case in the outer vacuum region. For the proper times it just means that
we get progressively more severe red shifts as we approach the origin, rather
than the blue shifts that we get as we approach the event horizon from the
outer vacuum region.

As a corollary to the proof that $r_{\mu}>0$, note that this fact guarantees
the positivity of the cubic polynomial in Equation (34). This is so because
the second derivative of that polynomial is given by $-6\kappa\rho_{0}r$,
being therefore negative for all $r\in[r_{1},r_{2}]$. This means that the
graph of the cubic polynomial has a concavity turned downward throughout this
interval. In addition to this, it is easy to see that at $r=r_{2}$ the
polynomial is given by $3\left(r_{2}-r_{M}\right)$, which is strictly positive
so long as $r_{2}>r_{M}$. Finally, at $r=r_{1}$ the polynomial is given by

$\kappa\rho_{0}\left(r_{2}^{3}-r_{1}^{3}\right)+3\left(r_{1}-r_{M}\right)=3\left(r_{1}+r_{\mu}\right),$
(69)

where we used Equation (32), which is also strictly positive since
$r_{\mu}>0$. As a consequence of this, we may conclude that, so long as the
conditions $r_{2}>r_{M}$ and $r_{\mu}>0$ hold, as they must for physically
sensible solutions, the polynomial is strictly positive for all
$r\in[r_{1},r_{2}]$.

###  3.2 Nature of the Inner Gravitational Field

The physical interpretation of the function $\nu(r)$ is that the proper time
interval at the radial position $r$, between two events occurring at the same
spatial point, is given by $d\tau=\exp[\nu(r)]dt$, where $dt$ is the time
interval between the two events as seen at spatial infinity, where spacetime
is flat. If we consider a photon traveling in the radial direction, either
inwards or outwards, this means that the proper frequency $f(r)$ of the photon
changes with position, between a first point $r_{a}$ and a second point
$r_{b}$, according to

$\displaystyle f(r_{a})$ $\displaystyle=$ $\displaystyle\,{\rm
e}^{-\nu(r_{a})}f_{\infty},$ $\displaystyle f(r_{b})$ $\displaystyle=$
$\displaystyle\,{\rm e}^{-\nu(r_{b})}f_{\infty},$ (70)

where $f_{\infty}$ is the frequency of the photon at radial infinity. Dividing
these two equations and making the two points very close together, so that
$r_{a}=r$ and $r_{b}=r_{a}+\delta r$, we have

$\frac{f(r+\delta r)}{f(r)}=\,{\rm e}^{-[\nu(r+\delta r)-\nu(r)]}.$ (71)

For sufficiently small $\delta r$ we may write the variation of the function
$\nu(r)$ in terms of its derivative $\nu^{\prime}(r)$, so that we get

$\frac{f(r+\delta r)}{f(r)}\simeq\,{\rm e}^{-\delta r\,\nu^{\prime}(r)}.$ (72)

Since the energy $hf(r)$ of a photon, $h$ being the Planck constant, is
proportional to its frequency, we have an interpretation of the red and blue
shifts of the frequency of the photons as decreases or increases in their
energies, respectively. We thus observe that, if a photon is going outward, so
that $\delta r>0$, and if the derivative $\nu^{\prime}(r)$ is positive, then
we will have that $f(r+\delta r)<f(r)$, and therefore a red shift in the
frequency. If it is going outward but the derivative is negative, then we will
have that $f(r+\delta r)>f(r)$ and hence a blue shift. On the other hand, if
the photon is going inward, so that $\delta r<0$, and the derivative is
positive, then we will have a blue shift, and finally, if it is going inward
and the derivative is negative, then we will have a red shift. Let us write
down the derivative of $\nu(r)$ in the inner and outer vacuum regions,

$\displaystyle\nu^{\prime}(r)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{rcl}-\,{\displaystyle\frac{\displaystyle
1}{\displaystyle 2}}\,{\displaystyle\frac{\displaystyle r_{\mu}}{\displaystyle
r(r+r_{\mu})}}&\mbox{for}&0\;\leq r\leq r_{1},\\\\[12.91663pt]
{\displaystyle\frac{\displaystyle 1}{\displaystyle
2}}\,{\displaystyle\frac{\displaystyle r_{M}}{\displaystyle
r(r-r_{M})}}&\mbox{for}&r_{2}\leq r<\infty.\end{array}\right.$ (75)

Let us now consider the consequences of Equation (72) in more detail in each
one of these two regions, starting with the outer vacuum region. As one can
see above, in the outer vacuum region, since we have that $r>r_{2}>r_{M}>0$,
the derivative $\nu^{\prime}(r)$ is always positive. Therefore, photons
traveling outward undergo red shifts, while those traveling inward undergo
blue shifts. This can be interpreted in energetic terms as the statement that
when traveling inward the photons gain energy from the gravitational field,
and when traveling outward they lose energy to it. This is characteristic of a
gravitational field that is attractive towards the origin.

However, in the inner vacuum region the situation is reversed. Since we have
that $r_{\mu}>0$, the derivative is everywhere negative in that region. This
means that photons traveling outward within this region are blue shifted, and
therefore gain energy from the gravitational field, while photons traveling
inward within this region are red shifted, and therefore lose energy to the
gravitational field. This is characteristic of a gravitational field that is
repulsive, driving matter and energy away from the origin. This is the exact
opposite of what happens in the outer vacuum region. It is important to note
that this repulsion is not from the matter in itself, but from the origin,
consisting therefore of an outward attraction towards the shell of matter.

##  4 Examples of Specific Solutions

In order to calculate $z(r)$ either analytically or numerically it is
convenient to define a dimensionless variable $x$ such that

$\displaystyle x$ $\displaystyle\equiv$
$\displaystyle\Upsilon_{0}\,r\;\;\;\Longrightarrow$
$\displaystyle\frac{d}{dr}$ $\displaystyle=$
$\displaystyle\Upsilon_{0}\frac{d}{dx}.$ (76)

In terms of $x$, Equation (46), that determines $z(r)$, becomes

$z^{\prime}(x)+\frac{\eta+2x^{3}}{2x(\eta+3x-x^{3})}z(x)=\frac{3x^{3}}{2x(\eta+3x-x^{3})},$
(77)

where the primes indicate now derivatives with respect to $x$, and where we
define

$\eta\equiv
x_{2}^{3}-3x_{M},\;\;x_{1}\equiv\Upsilon_{0}\,r_{1},\;\;x_{2}\equiv\Upsilon_{0}\,r_{2},\;\;x_{M}\equiv\Upsilon_{0}\,r_{M}.$
(78)

Thus $x_{1}$, $x_{2}$ and $x_{M}$ correspond respectively to the internal
radius $r_{1}$, the external radius $r_{2}$ and the Schwarzschild radius
$r_{M}$, expressed in terms of the new variable $x$. The solution of Equation
(77) is obtained by writing Equation (52) in terms of $x$,

$z(x)=\sqrt{\frac{\eta+3x-x^{3}}{x}}\left[\sqrt{\frac{x_{2}}{3(x_{2}-x_{M})}}+\frac{3}{2}\,\int_{x_{2}}^{x}dy\,\frac{y^{5/2}}{(\eta+3y-y^{3})^{3/2}}\right],$
(79)

where, in order to remain within the matter region, we must have $x_{1}\leq
x\leq x_{2}$. If we multiply both the numerator and the denominator of the
integral in Equation (79) by $y^{3/2}$, define the polynomial
$Q(y)=y\left(\eta+3y-y^{3}\right)$ and the rational function $S(y,Q)\equiv
y^{4}/Q^{3}$, then the integral in Equation (79) can be rewritten as

$\int_{x_{2}}^{x}dy\,\frac{y^{5/2}}{(\eta+3y-y^{3})^{3/2}}=\int_{x_{2}}^{x}S\\!\left[y,\sqrt{Q(y)}\,\right]dy.$
(80)

The expression on the right-hand side of Equation (80) is by definition an
elliptic integral [12] and cannot be expressed in terms of elementary
functions except in two cases: 1) $S\\!\left(y,Q^{1/2}\right)$ contains no odd
powers of $y$; in our case this happens when $\eta=0$ and leads to the
Schwarzschild interior solution; 2) the polynomial $Q(y)$ has two equal roots;
this leads to the explicit solutions that we discuss next.

###  4.1 A Family of Explicit Solutions

The integral in Equation (79) contains a cubic polynomial. The nature of its
three roots depends on the value of its discriminant $\Delta$ [13]. For cubic
polynomials of the form $ax^{3}+cx+d$ we have $\Delta=-4ac^{3}-27a^{2}d^{2}$.
If $\Delta>0$ the polynomial has three distinct real roots, if $\Delta=0$ it
has three real roots but two of them are equal, and if $\Delta<0$ it has one
real and two complex roots which are conjugate to each other. In our case we
have $a=-1$, $c=3$, $d=\eta$ and therefore $\Delta=27(4-\eta^{2})$.

Figure 1: Comparison between the dimensionless pressure $p(x)$ calculated
analytically and numerically using the Runge-Kutta fourth-order algorithm for
$\eta=2.0$, $x_{2}=5^{1/3}$, $x_{M}=1.0$ and $x_{1}=0.594881$.

Figure 2: The functions $\nu(x)$ and $\lambda(x)$ for $\eta=2.0$,
$x_{2}=5^{1/3}$ and $x_{M}=1.0$. The shaded area indicates the matter region,
to its right is the outer vacuum and to its left is the inner vacuum.

Figure 3: The dimensionless pressure $p$ calculated numerically for
$\eta=5.0$, $x_{2}=2.0$ and $x_{M}=1.0$.

Figure 4: The functions $\nu(x)$ and $\lambda(x)$ for $\eta=5.0$, $x_{2}=2.0$
and $x_{M}=1.0$. The shaded area indicates the matter region, to its right is
the outer vacuum and to its left is the inner vacuum.

The value $\Delta=0$ corresponds to the case where the solution for $z(x)$ can
be expressed in terms of elementary functions. Note that we have $\Delta=0$
when $\eta=\pm 2$, which corresponds to $x_{2}^{3}=\pm 2+3x_{M}$. For
$\eta=-2$ the polynomial in the integral of Equation (79) is non-positive for
$x\geq 0$. Therefore, we must choose $\eta=2$. For this value of $\eta$ the
polynomial is strictly positive in the interval $[0,2)$ and can be factored as

$2+3y-y^{3}=(2-y)(y+1)^{2}.$ (81)

In this case we can express the integral in Equation (79) in terms of
elementary functions. The calculation can be considerably simplified using a
new integration variable $u$ defined by $u=\sqrt{y/(2-y)}$. The final result,
up to an integration constant, is

$\displaystyle\mathcal{I}(y)$ $\displaystyle\equiv$
$\displaystyle\int\,dy\,\frac{y^{5/2}}{(2-y)^{3/2}(y+1)^{3}}$ (82)
$\displaystyle=$
$\displaystyle\frac{2y^{2}+15y+10}{18\,(y+1)^{2}}\sqrt{\frac{y}{2-y}}-\,\frac{5\sqrt{3}}{27}\arctan\\!\left(\sqrt{\frac{3y}{2-y}}\,\right).$

Thus, in terms of $\mathcal{I}(y)$ Equation (79) reads

$z(x)=\sqrt{\frac{2+3x-x^{3}}{x}}\left\\{\sqrt{\frac{x_{2}}{3(x_{2}-x_{M})}}+\frac{3}{2}\left[\rule{0.0pt}{10.76385pt}\mathcal{I}(x)-\mathcal{I}(x_{2})\right]\right\\}.$
(83)

Note that, in order to guarantee that the cubic polynomial for $\eta=2$ shown
in Equation (81) is always positive, we need to have $y<2$. Therefore, since
we already know that the polynomial is positive, the arguments of the square
roots in Equation (82) are always positive.

###  4.2 Examples of Numerical Solutions

In our numerical approach here, we assume that the external radius
$x_{2}=\Upsilon_{0}r_{2}$ is given. In order to complete the calculation we
have to determine the interior radius $x_{1}$. This can be done recalling that
the dimensionless pressure $p(x)$ is zero for $x=x_{2}$ and $x=x_{1}$. Since
according to Equation (2.4) $p(x)=1/z(x)-1$, this is equivalent to the
determination of the values of $x$ for which $z(x)=1$. By the determination of
$x_{1}$ we would have solved the problem in the entire matter region. Note
that since $x=\Upsilon_{0}r=\sqrt{\kappa\rho_{0}}\;r$ we have obtained a
family of solutions parametrized by two parameters, the external radius
$r_{2}$ and the parameter $\eta$.

If the discriminant $\Delta\neq 0$ the integral in Equation (79) is expressed
in terms of elliptic integrals and the result is not very transparent. It is
more convenient to integrate the differential Equation (77) using the fourth-
order Runge-Kutta algorithm (RK4) [14]. We start by choosing a value of
$x=x_{2}$ for which the cubic polynomial is positive and we put $z(x_{2})=1$.
This determines the outer radius of the matter shell. We then iterate the
differential equation given in Equation (77) in the decreasing $x$ direction
until we reach the first point for which the value of $z$ returns to $1$. This
point is chosen as $x_{1}$. If a value for $x_{1}$ cannot be found, we
conclude that there is no solution to the problem with the given values of
$x_{2}$ and $x_{M}$. A good test for the efficiency of the algorithm is to
compare the exact analytic result given in Equation (83) with the result from
the numerical integration in that same case. These results are shown in Figure
1. On any current $64$-bit desktop computer one can easily reach a high degree
of precision with little numerical effort. After iterating the RK4 algorithm
from $x_{2}$ to $x_{1}$ the difference between the exact and the numerical
results for $z(x)$ stays below $1.03536\times 10^{-29}$ for an iteration step
of $\delta x\approx 10^{-7}$.

In the comments that follow $x_{\mu}\equiv\Upsilon_{0}r_{\mu}$, where
$r_{\mu}$ is the integration constant that results from the solution of the
Einstein equations in the inner vacuum region, given in Equation (32). In the
matter region the input parameters are $\eta$ and $x_{2}$. The parameter
$x_{1}$ is obtained from the iteration of Equation (77). The value of $x_{M}$
that is necessary for plotting the curves is given in Equation (78). The
expressions for $\lambda(x)$ and $\nu(x)$ are given in Table 1. Figure 2 shows
the plots of the functions $\nu(x)$ and $\lambda(x)$ for $\eta=2.0$ and
$x_{2}=5^{1/3}$. The curves were obtained analytically using Equation (83) and
the expressions in Table 1, but using the numerically calculated parameters
$x_{1}=0.594881$ and $x_{\mu}=0.596494$.

In Figure 3 we plot the dimensionless pressure $p(x)$ as a function of $x$, in
a case in which there is no analytic expression and the calculation is
performed numerically. The parameters are $x_{1}=1.24050$ and
$x_{\mu}=1.03035$. Comparing Figures 1 and 3, that depict the dimensionless
pressure $p(x)$ as a function of $x$ for $\eta=2.0$ and $\eta=5.0$, one notes
that the two graphs are similar but for larger values of $\eta$ the graph
becomes less symmetric.

Figure 4 shows the plots of the functions $\nu(x)$ and $\lambda(x)$, for
$\eta=5.0$ and $x_{2}=2.0$. In this case there are no analytical solutions
available in the matter region and the values of $\nu(x)$ and $\lambda(x)$
were obtained numerically. In the vacuum regions we used the analytical
expressions given in Table 1 with the parameters $x_{1}=1.24050$ and
$x_{\mu}=1.03035$.

##  5 Conclusions

In this paper we have given the complete and exact solution of the Einstein
field equations for the case of a shell of liquid matter. Although this
particular problem can be seen as having a somewhat academic nature, it does
lead us to two important and unexpected conclusions. One of them is that all
solutions for shells of liquid matter have a singularity at the origin, within
the inner vacuum region, that does not, however, lead to any kind of
pathological behavior involving the matter. The other is that, contrary to
what is usually thought, a non-trivial gravitational field does exist within a
spherically symmetric central cavity, namely the inner vacuum region.

The geometry within the cavity is associated with a spacetime that is
contracted in the radial direction, rather than expanded. It is easy to verify
that, unlike what happens in the outer vacuum region, the proper radial
length, $\ell_{1}$, say from $r=0$ to $r=r_{1}$, is in fact smaller than the
corresponding radial coordinate $r_{1}$. We have that
$d\ell_{1}=\sqrt{g_{11}}\,dr$, and therefore

$\displaystyle\ell_{1}$ $\displaystyle=$
$\displaystyle\int_{0}^{r_{1}}dr\,\sqrt{\frac{r}{r+r_{\mu}}}$ (84)
$\displaystyle<$ $\displaystyle\int_{0}^{r_{1}}dr$ $\displaystyle=$
$\displaystyle r_{1},$

given that $r_{\mu}>0$. This illustrates the fact that the radial lengths
within the inner vacuum region are contracted rather than expanded. The true
physical volume of the inner vacuum region is therefore correspondingly
smaller than the apparent coordinate volume. This renders this inner geometry
not embeddable in the illustrative way that is usually employed in the case of
the outer vacuum region.

The gravitational field associated to this geometry, inside the inner vacuum
region, can be interpreted as a repulsive field with respect to the origin.
This can be ascertained from an examination of the sign of the derivative of
$\nu(r)$ in the inner and outer vacuum regions, and its interpretation in
terms of the energy of a photon traveling in the radial direction. This sign
is positive in the outer vacuum region, corresponding to an attractive field
towards the origin, and negative in the inner vacuum region, corresponding to
an repulsive field away from the origin. Of course, since $\nu^{\prime}(r)$ is
a continuous function, and since we enter the matter region from the outer
vacuum region with a positive derivative, and exit it into the inner vacuum
region with a negative derivative, there must be a point within the matter
region where $\nu^{\prime}(r)=0$, and where the derivative flips sign. This is
clearly the point $r_{e}$ of minimum of $\nu(r)$, which is also the point of
minimum of $z(r)$, and hence the point of maximum of the pressure $P(r)$, a
point which already had a role to play in our arguments.

One can acquire an intuitive understanding of the unexpected situation in the
inner vacuum region by observing that such a situation can arise even within
the Newtonian framework, if we use a slightly modified potential. We can do
this if we consider the Newtonian argument for the gravitational force within
a hollow thin spherical shell of matter, but with a potential that behaves as
$1/r^{1+\epsilon}$ for some $|\epsilon|\ll 1$, thus leading to a force that
behaves as $1/r^{2+\epsilon}$. If one considers a test mass at a point in the
interior of the hollow shell, at the position $\vec{r}$ with respect to the
center, it is not difficult to use the usual Newtonian argument to show that,
if $\epsilon>0$, then the resulting gravitational force at that point is
oriented outward, in the direction of $\vec{r}$, towards the shell of matter.
In other words, the attraction by the part of the shell that is closer to the
point $\vec{r}$ outweighs the attraction from the opposite side, thus leading
to a resulting force that repels particles away from the origin. Note that
this argument involving a potential behaving in a way other than $1/r$ is the
same that can be used to model the precession of the perihelion of orbits in
General Relativity using this semiclassical Newtonian framework. That
precession is prograde precisely if $\epsilon>0$.

It is interesting to note that this configuration of the gravitational field
tends to stabilize the shell of liquid matter, since any particle of matter
that detaches from the liquid and wanders into one of the vacuum regions will
be driven back to the bulk of the liquid by the gravitational field. This can
be interpreted as a successful stability test satisfied by all the solutions.
The general tendency of the gravitational field is therefore that of
compressing the shell of fluid matter, from both sides. This suggests that the
same interpretation should be valid in the case of a gaseous fluid.

The singularity at the origin is usually thought to be associated with an
infinite concentration of matter there, and thus considered to be an evil that
must be avoided at any cost. However, this only makes any sense at all if one
thinks of that singularity as a point of gravitational attraction, rather than
as a point of repulsion of matter. Here we do have the singularity, but not
the infinite concentration of matter at the origin, due to the repulsive
character of the gravitational field around the origin. In any case, the
existence of the singularity is not a question of choice, of course, since it
is required by the field equations and by the interface boundary conditions
that follow from them. One is not at liberty to impose that $r_{\mu}=0$ in
order to avoid this singularity.

## References

  * [1] K. Schwarzschild, “Über das gravitationsfeld eines massenpunktes nach der einsteinschen theorie (on the gravitational field of a mass point according to einstein’s theory),” Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften, vol. 7, pp. 189–196, 1916.
  * [2] R. Wald, General Relativity.  University of Chicago Press, 2010.
  * [3] J. T. Jebsen, “Über die allgemeinen kugelsymmetrischen lösungen der einsteinschen gravitationsgleichungen im vakuum (on the general spherically symmetric solutions of einstein’s gravitational equations in vacuo),” Arkiv för Matematik, Astronomi och Fysik, vol. 15, pp. 1–9, 1921.
  * [4] G. D. Birkhoff, Relativity and Modern Physics.  Harvard University Press, 1923.  Cambridge, Massachusetts, LCCN 23008297.
  * [5] K. Schwarzschild, “Über das gravitationsfeld einer kugel aus inkompressibler flüssigkeit nach der einsteinschen theorie (on the gravitational field of a ball of incompressible fluid following einstein’s theory),” Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften, vol. 7, pp. 424–434, 1916.
  * [6] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation.  San Francisco: W.H. Freeman and Co., 1973.
  * [7] X. Mei, “The precise inner solutions of gravity field equations of hollow and solid spheres and the theorem of singularity,” International Journal of Astronomy and Astrophysics, vol. 1, pp. 109–116, 2011.
  * [8] S. Weinberg, Gravitation and Cosmology.  New York: John Wiley and Sons, 1972.
  * [9] J. Ni, “Solutions without a maximum mass limit of the general relativistic field equations for neutron stars,” Science China, vol. 54, no. 7, pp. 1304–1308, 2011.
  * [10] L. Neslušan, “Solutions without a maximum mass limit of the general relativistic field equations for neutron stars,” Journal of Modern Physics, vol. 6, pp. 2164–2183, 2015.
  * [11] P. A. M. Dirac, General Theory of Relativity.  John Wiley & Sons, Inc., 1975.  ISBN 0-471-21575-9.
  * [12] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.  Applied mathematics series, U.S. Government Printing Office, 1965.
  * [13] “Cubic equation.” Wikipedia.  https://en.wikipedia.org/wiki/Cubic_equation.
  * [14] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling, Numerical Recipes in FORTRAN 77: Volume 1, Volume 1 of Fortran Numerical Recipes: The Art of Scientific Computing.  Cambridge University Press, 1992.


# Physics-Informed Learning of Aerosol Microphysics

Paula Harder  Atmospheric, Oceanic and Planetary Physics, Department of
Physics, Universty of Oxford  Fraunhofer Center High-Performance Computing,
Fraunhofer ITWM  Fraunhofer Center Machine Learning, Fraunhofer Society  Chair
for Scientific Computing, TU Kaiserslautern  Duncan Watson-Parris
Atmospheric, Oceanic and Planetary Physics, Department of Physics, Universty
of Oxford  Philip Stier  Atmospheric, Oceanic and Planetary Physics,
Department of Physics, Universty of Oxford  Dominik Strassel  Fraunhofer
Center High-Performance Computing, Fraunhofer ITWM  Nicolas R. Gauger  Chair
for Scientific Computing, TU Kaiserslautern  Janis Keuper  Fraunhofer Center
High-Performance Computing, Fraunhofer ITWM  Fraunhofer Center Machine
Learning, Fraunhofer Society  Institute for Machine Learning and Analytics,
Institute for Machine Learning and Analytics

Aerosol particles play an important role in the climate system by absorbing
and scattering radiation and influencing cloud properties. They are also one
of the biggest sources of uncertainty for climate modeling. Many climate
models do not include aerosols in sufficient detail due to computational
constraints. In order to represent key processes, aerosol microphysical
properties and processes have to be accounted for. This is done in the ECHAM-
HAM global climate aerosol model using the M7 microphysics, but high
computational costs make it very expensive to run with finer resolution or for
a longer time. We aim to use machine learning to emulate the microphysics
model at sufficient accuracy and reduce the computational cost by being fast
at inference time. The original M7 model is used to generate data of input-
output pairs to train a neural network on it. We are able to learn the
variables’ tendencies achieving an average $R^{2}$ score of $77.1\%$. We
further explore methods to inform and constrain the neural network with
physical knowledge to reduce mass violation and enforce mass positivity. On a
GPU we achieve a speed-up of up to over 64x compared to the original model.

##  1 Introduction

Aerosol forcing remains the largest source of uncertainty in the anthropogenic
effect on the current climate [1]. The aerosol cooling effect hides some of
the positive radiative forcing caused by greenhouse gas emissions and future
restrictions to lower air pollution might result in stronger observed warming.
Aerosols impact climate change through aerosol-radiation interactions and
aerosol-cloud interactions [2]. They can either scatter or absorb radiation,
which depends on the particles’ compounds. Black carbon aerosols from fossil
fuel burning for example have a warming effect by strongly absorbing
radiation, whereas sulphate from volcanic eruptions has a cooling effect by
being less absorbing and primarily scattering radiation. Clouds influence the
earth’s radiation budget by reflecting sunlight and aerosols can change cloud
properties significantly by acting as cloud condensation nuclei (CCN). A
higher concentration of aerosols leads to more CCN which, for a fixed amount
of water, results in more but smaller cloud droplets. Smaller droplets
increase the clouds’ albedo [3] and can enhance the clouds’ lifetime [4].

Many climate models consider aerosols only as external parameters, they are
read once by the model, but then kept constant throughout the whole model run.
In the case where some aerosol properties are modeled, there might be no
distinction between aerosol types, and just an overall mass is considered. To
incorporate a more accurate description of aerosols, aerosol-climate modeling
systems like ECHAM-HAM [5] have been introduced. It couples the ECHAM General
Circulation Model (GCM) with a complex aerosol model called HAM. The
microphysical core of HAM is either SALSA [6] or the M7 model [7]. We consider
the latter here. M7 uses seven log-normal modes to describe aerosol
properties, the particle sizes are represented by four size modes, nucleation,
aitken, accumulation, and coarse, of which the aitken, accumulation, and
coarse can be either soluble or insoluble 11footnotetext: In the paper we will
use following abbreviations: NS=nucleation soluble, KS/KI=aitken
soluble/insoluble, AS/AI=accumulation soluble/insoluble, CS/CI=coarse
soluble/insoluble. It includes processes like nucleation, coagulation,
condensation, and water uptake, which lead to the redistribution of particle
numbers and mass among the different modes. In addition, M7 considers five
different components - Sea salt (SS), sulfate (SO4), black carbon (BC),
primary organic carbon (OC), and dust (DU). M7 is applied to each grid box
independently, it does not model any spatial relations.

Change in H2SO4 concentration

Figure 1: The change in concentration modeled by the M7 module for the first
time step of the test data is plotted on the left. The predicted change is
plotted on the right. Both plots show the change in concentration on a
logarithmic scale

More detailed models come with the cost of increased computational time:
ECHAM-HAM can be run at 150km resolution for multiple decades. But to run
storm-resolving models, for example, the goal is ideally a 1km horizontal grid
resolution and still be able to produce forecasts up to a few decades. If we
want to keep detailed aerosol descriptions for this a significant speedup of
the aerosol model is needed.

Replacing climate model components with machine learning approaches and
therefore decreasing a model’s computing time has shown promising results in
the past. There are several works on emulating convection, both random forest
approaches, [8] and deep learning models [9], [10], [11] have been explored.
Recently, multiple consecutive neural networks have been used to emulate a bin
microphysical model for warm rain processes [12]. Silva et al. [13] compare
several methods, including deep neural networks, XGBoost, and ridge regression
as physically regularized emulators for aerosol activation. In addition to the
aforementioned approaches random forest approaches have been used to
successfully derive the CCN from atmospheric measurements [14]. Next to many
application cases there now exist two benchmark datasets for climate model
emulation [15], [16].

Physics-Informed machine learning recently gained attention in machine
learning research, including machine learning for weather and climate modeling
[17]. Whereas a large amount of work focuses on soft-constraining neural
networks by adding equations to the loss term, there was also the framework
DC3 [18] developed, that incorporates hard constraints. Using emulators in
climate model runs can require certain physical constraints to be met, that
are usually not automatically learned by a neural network. For convection
emulation there has been work on both hard and soft constraints to enforce
conversation laws, adding a loss term or an additional network layer [19].

Building on our previous work [20], we demonstrate a machine learning approach
to emulate the M7 microphysics module. We investigated different approaches,
neural networks as well as ensemble models like random forest and gradient
boosting, with the aim of achieving the desired accuracy and computational
efficiency, finding the neural network (NN) appeared to be the most
successful. We use data generated from a realistic ECHAM-HAM simulation and
train a model offline to predict one time step of the aerosol microphysics.
The underlying data distribution is challenging, as the changes in the
variables are often zero or very close to zero. We do not predict the full
values, but the tendencies. An example of our global prediction is shown in
Figure 1. To incorporate physics into our network we explore both soft
constraints for our emulator by adding regularization terms to the loss
function and hard constraints by adding completion and correction layers.

##  2 Methodology

###  2.1 Data

####  2.1.1 Data Generation

To easily generate data we extract the aerosol microphysics model from the
global climate model ECHAM-HAM and develop a stand-alone version. To obtain
input data ECHAM-HAM is run for four days within a year. We use a horizontal
resolution of 150km at the equator, overall 31 vertical levels, 96 latitudes,
and 192 longitudes. A time step length of 450s is chosen. This yields a data
set of over 100M points for each day. We only use a subset of 5 times a day,
resulting in about 2.85M points per day. To test if our emulator generalizes
well in an unknown setting, we use the days in January and April for training,
a day in July for validation, and the October data for testing. Because the
two sea salt variables only change in $0.2\%$ of the time as they undergo
little microphysical changes we do not model them here. The masses and
particle concentrations for different aerosol types and different modes are
both inputs and outputs, the output being the value one time step later.
Atmospheric state variables like pressure, temperature, and relative humidity
are used as inputs only. The dataset is available at [21]. A full list of
input and output values can be found in the supplementary material.

####  2.1.2 Data distribution and transformation

Compared to the full values the tendencies (changes) are very small, therefore
we aim to predict the tendencies, not the full values, where it is possible.
Depending on the variable we have very different size scales, but also a
tendency for a specific variable may span several orders of magnitude. In some
modes, variables can only grow, in some only decrease, and in others do both.
Often the majority of tendencies are either zero or very close to zero, but a
few values might be very high. The absolute tendencies are roughly log-normal
distributed. A logarithmic transformation has been used in [12] to give more
weight to the prominent case, the smaller changes. This then transforms our
inputs close to a normal distribution, which is favorable for the neural
network to be trained on. On the other hand, a linear scale might be more
representative, as it gives more weight to values with the largest actual
physical impact. Results and a discussion for training with logarithmically
transformed variables can be found in the appendix.

###  2.2 General Network Architecture

For emulating the M7 microphysics model we explored different machine learning
approaches, including random forest regression, gradient boosting, and a
simple neural network (see suppl.). Providing more expressivity, the neural
network approach appears to be the most successful for this application and
will be shown here. We do not only present one neural model here, but
different versions building on the same base architecture and discuss the
advantages and disadvantages of these different designs.

We employ a fully connected neural network, where the values are propagated
from the input layer through multiple hidden to the output layer, using a
combination of linear operations and non-linear activations. We use a ReLU
activation function (a comparison of scores for different activation functions
can be found in the supplementary material) and three hidden layers, each
hidden layer containing 128 nodes. Using zero hidden layers results in linear
regression and does not have the required expressivity for our task, one
hidden layer already achieves a good performance, after two layers the model
we could not see any significant further improvements.

####  2.2.1 Training

We train all our NNs using the Adam optimizer with a learning rate of
$10^{-3}$, a weight decay of $10^{-9}$, and a batch size of 256. Our objective
to optimize during the training of the neural networks is specified in
Equation 2, using an MSE loss in every version and activating the additional
loss terms depending on the version. WE train our model for 100 epochs. The
training on a single NVIDIA Titan V GPU takes about 2 hours.

For the classification network, we use a similar setup as described above. As
a loss function be use the Binary-Cross-Entropy loss, the training takes about
30 minutes.

###  2.3 Physics-Constrained Learning

In this work, we explore two different ways of physics-informing similar to
[19]: Adding regularizer terms to the loss function and adding hard
constraints as an additional network layer. The physical constraints that
naturally come up for our setting are mass conservation, the mass of one
aerosol species has to stay constant, and mass or concentration positivity, a
mass or concentration resulting from our predicted tendencies has to be
positive.

####  2.3.1 Soft Constraining

To encourage the neural network to decrease mass violation or negative
mass/concentration values we add additional loss terms to the objective:

$\min_{\theta}\mathcal{L}(\tilde{y},y_{\theta})+\lambda\mathcal{L}^{\text{(mass)}}(y_{\theta})+\mu\mathcal{L}^{\text{(pos)}}(y_{\theta}).$
(1)

Where $y_{\theta}=f_{\theta}(x)$ is the network output for an input
$x\in\mathbb{R}^{32}$, parameterized by the network’s weights $\theta$. In our
case $\mathcal{L}$ is a MSE error, $\mu,\ \lambda\in\\{0,1\\}$ depending if a
regularizer is activated or not. The term $\mathcal{L}^{\text{(mass)}}$
penalizes mass violation:

$\mathcal{L}^{\text{(mass)}}(y_{\theta}):=\sum_{s\in S}\alpha_{s}|\sum_{i\in
I_{s}}y_{\theta}^{(i)}|,$ (2)

Here $S$ is the set of different species and $I_{s}$ are the indices for a
specific species $s$. The parameters $\alpha_{i}$ need to be tuned, chosen too
small the mass violation will not be decreased, if a factor is too large, the
network can find the local optimum that is constantly zero, the values can be
found in the supplementary material. The other term penalizes negative mass
values and is given by 22footnotetext: Skipping the back transformation in to
original scale for a clearly notation, details on that are provided in the
supplementary material.:

$\mathcal{L}^{\text{(pos)}}(y_{\theta}):=\sum_{i=0}^{n-1}\beta_{i}\text{ReLU}(-(y_{\theta}^{(i)}+x_{i}))^{2}$
(3)

Using a ReLU function, all negative values are getting penalized. In our case
$n=28$, $x_{i}$ is the corresponding full input variable to predicted tendency
$y_{\theta}^{(i)}$. The factors $\beta_{i}$ need to be tuned, similar to
$\alpha_{i}$.

####  2.3.2 Hard Constraining

In order to have a guarantee for mass conservation or positive masses, we
implement a correction and a completion method. This is done at test time as
an additional network layer, a version where this is implemented within the
training of the neural network did not improve the performance.

##### Inequality Correction

For the correction method, predicted values that result in an overall negative
mass or number are set to a value such that the full value is zero. For the
neural networks intermediate output $\tilde{y}_{\theta}$, the final output is
given by:

$y_{\theta}^{(i)}=\text{ReLU}((\tilde{y}_{\theta}^{(i)}+x_{i}))-x_{i}$ (4)

##### Equality Completion

The completion method addresses mass violation, for each species one
variable’s prediction is replaced by the negative sum of the other same
species’ tendencies. We obtained the best results by replacing the worst-
performing variable. This completion results in exact mass conservation. For a
species $S$ we choose index $j\in I_{S}$, the completion layer is defined as
follows:

$y_{\theta}^{(j)}=-\sum_{i\in I_{s}\setminus\\{j\\}}y_{\theta}^{(i)}$ (5)

##  3 Results

###  3.1 Predictive Performance

####  3.1.1 Metrics

We consider multiple metrics to get an insight into the performance of the
different emulator designs, covering overall predictive accuracy, mass
violation, and predicting nonphysical values. We look at the $R^{2}$ score and
the MSE. To understand mass violation in we look at the mass biases for the
different species and the overall mass violation, where all scores are
normalized by the mean over the respective species. The metrics are completed
with two scores about negative value predictions: An overall fraction of
negative and therefore nonphysical predictions and the average negative extend
per predicted value. For all the different scores and architectures we take
the mean over five different random initializations of the underlying neural
network.

####  3.1.2 Evaluation and Comparison

Table 1: Test metrics for different architectures and transformations. Base means the usage of only the base NN, correct adds the correction procedure and complete the completion procedure. Mass reg. and positivity reg. include the regularization terms. Best scores are bold and second best scores are highlighted in blue . Architecture | Base | +Correct | +Complete | +Mass Loss | +Positivity Loss
---|---|---|---|---|---
R2 | $0.763$ | 0.771 | $0.738$ | 0.730 | 0.709
MSE | 0.162 | 0.161 | 0.162 | 0.187 | 0.211
Mass Bias SO4 | $1.1\cdot 10^{-5}$ | $8.5\cdot 10^{-5}$ | 0.00 | $8.6\cdot 10^{-6}$ | $1.0\cdot 10^{-3}$
Mass Bias BC | $3.8\cdot 10^{-5}$ | $1.4\cdot 10^{-4}$ | 0.00 | $3.4\cdot 10^{-5}$ | $3.6\cdot 10^{-4}$
Mass Bias OC | $3.3\cdot 10^{-5}$ | $6.0\cdot 10^{-5}$ | 0.00 | $1.1\cdot 10^{-5}$ | $6.4\cdot 10^{-4}$
Mass Bias DU | $1.0\cdot 10^{-6}$ | $3.9\cdot 10^{-5}$ | 0.00 | $2.8\cdot 10^{-7}$ | $1.5\cdot 10^{-5}$
Mass Violation | $3.7\cdot 10^{-4}$ | $1.1\cdot 10^{-3}$ | 0.00 | $1.4\cdot 10^{-4}$ | $2.4\cdot 10^{-3}$
Neg. Fraction | 0.134 | 0.00 | 0.146 | 0.144 | 0.0894
Neg. Mean | 0.00151 | 0.00 | 0.00170 | 0.00169 | 0.000081

In Table LABEL:many_scores we display the scores for our emulator variants. We
achieve good performance with a $77.1\%$ $R^{2}$ score and an MSE of $0.16$.
Using the correction method results per construction in no negative
predictions, but increases the mass violation. The accuracy scores are not
negatively affected by the correction operation. With the completion method,
we achieve perfect mass conservation and a very slight worsening of the other
metrics. The mass regularization decreases the overall mass violation and
decreases the mass biases for most cases. The positivity loss term decreases
the negative fraction and negative mean by a large amount. A few examples are
plotted in Figure 2. Overall a the architecture with an additional correction
procedure could be a good choice to use in a GCM run, having the guarantee of
no negative masses and good accuracy in the original units. The mass violation
is relatively low for all cases.

Although the emulator performs well in an offline evaluation, it still remains
to be shown how it performs when implemented in the GCM and run for multiple
time steps. Good offline performance is not necessarily a guarantee for good
online performance [22] in the case of convection parameterization, where a
model crash could occur. In our case, it is likely that offline performance is
a good indicator of online performance, as a global climate model crash is not
expected because aerosols do not strongly affect large-scale dynamics.

Figure 2: This figure shows the test prediction of our emulators against the
true M7 values.For each type (species, number particles, water) we plot the
performance of one variable (using the median or worse performing, see all in
suppl.)

###  3.2 Runtime

We conduct a preliminary runtime analysis by comparing the Python runtime for
the emulator with the Fortran runtime of the original aerosol model. We take
the time for one global time step, which is about 570,392 data points to
predict. For the M7 model, the 31 vertical levels are calculated
simultaneously and for the emulator, we predict the one time step at once,
using a batch size of 571,392 and taking into account the time for
transforming variables and moving them on and off the GPU. We use a single
NVIDIA Titan V GPU and a single Intel Xeon 4108 CPU. As shown in Table 2 we
can achieve a large speed-up of over 11,000x in a pure GPU setting. Including
the time it takes to move the data from the CPU onto the GPU and back the
acceleration is 64x compared to the original model. In case of no available
GPU, the NN emulator is still 2.8 times faster. Here, further speed-ups will
be achieved by using multiple CPUs, a smaller network architecture, and
efficient implementation in Fortran.

Table 2: Runtime comparison for the original M7 model and the NN emulator. NN
pure GPU includes only the transformation of the variables and the prediction,
whereas NN CPU-GPU-CPU also includes the time to transfer the data from the
CPU to the GPU and back to the CPU.
Model | M7 | NN pure GPU | NN CPU-GPU-CPU | NN CPU
---|---|---|---|---
time (s) | 5.781 | 0.000517 | 0.0897 | 2.042
speed-up | - | 11181.8 | 64.4 | 2.80

##  4 Conclusion and Future Work

To enable accurate climate forecasts aerosols need to be modeled with lower
computational costs. This work shows how neural networks can be used to learn
the mapping of an aerosol microphysics model and how simple physical
constraints can be included in the model. Our neural models approximate the
log-transformed tendencies excellently and the original units’ tendencies
well. On the test set, an overall regression coefficient of $77.1\%$ as well
as an MSE of $16.1\%$ are achieved. Using a GPU we accomplish a large speed-up
of 11,1181x compared to the original M7 model in a pure GPU setting, with the
time to move from and to the CPU we are still significantly faster having a
speed-up factor of 64. On a single CPUm the speed-up is 2.8x. By adding
completion and correction mechanisms we can remove mass violation or negative
predictions completely and make our emulator feasible for stable a GCM run.

How much of a speed-up can be achieved in the end remains to be shown, when
the machine learning model is used in a GCM run. Different versions of our
emulator need to be run within the global climate model to show which setup
performs the best online. A further step would be the combination of the
different methods for physics constraining, to achieve both mass conservation
and mass positivity at the same time.

## References

  * [1] N. Bellouin, J. Quaas, E. Gryspeerdt, S. Kinne, P. Stier, D. Watson-Parris, O. Boucher, K. S. Carslaw, M. Christensen, A.-L. Daniau, and el al. Dufresne.  Bounding global aerosol radiative forcing of climate change.  Reviews of Geophysics, 58(1):e2019RG000660, 2020.  e2019RG000660 10.1029/2019RG000660.
  * [2] IPCC.  Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change.  Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 2013.
  * [3] S. Twomey.  Pollution and the planetary albedo.  Atmospheric Environment (1967), 8(12):1251–1256, 1974.
  * [4] B. A. Albrecht.  Aerosols, cloud microphysics, and fractional cloudiness.  Science, 245(4923):1227–1230, 1989.
  * [5] P. Stier, J. Feichter, S. Kinne, S. Kloster, E. Vignati, J. Wilson, and et al.  The aerosol-climate model echam5-ham.  Atmospheric Chemistry and Physics, 5(4):1125–1156, 2005.
  * [6] H. Kokkola, H. Korhonen, K. E. J. Lehtinen, R. Makkonen, A. Asmi, S. Järvenoja, T. Anttila, A.-I. Partanen, M. Kulmala, H. Järvinen, A. Laaksonen, and V.-M. Kerminen.  Salsa: a sectional aerosol module for large scale applications.  Atmospheric Chemistry and Physics, 8(9):2469–2483, 2008.
  * [7] Elisabetta Vignati, Julian Wilson, and Philip Stier.  M7: An efficient size-resolved aerosol microphysics module for large-scale aerosol transport models.  Journal of Geophysical Research: Atmospheres, 109(D22), 2004.
  * [8] Paul A. O’Gorman and John G. Dwyer.  Using machine learning to parameterize moist convection: Potential for modeling of climate, climate change, and extreme events.  Journal of Advances in Modeling Earth Systems, 10(10):2548–2563, 2018.
  * [9] Stephan Rasp, Michael S. Pritchard, and Pierre Gentine.  Deep learning to represent subgrid processes in climate models.  Proceedings of the National Academy of Sciences, 115(39):9684–9689, 2018.
  * [10] P. Gentine, M. Pritchard, S. Rasp, G. Reinaudi, and G. Yacalis.  Could machine learning break the convection parameterization deadlock?  Geophysical Research Letters, 45(11):5742–5751, 2018.
  * [11] Tom Beucler, Michael Pritchard, Pierre Gentine, and Stephan Rasp.  Towards physically-consistent, data-driven models of convection, 2020\.
  * [12] A. Gettelman, D. J. Gagne, C.-C. Chen, M. W. Christensen, Z. J. Lebo, H. Morrison, and G. Gantos.  Machine learning the warm rain process.  Journal of Advances in Modeling Earth Systems, 13(2):e2020MS002268, 2021.  e2020MS002268 2020MS002268.
  * [13] S. J. Silva, P.-L. Ma, J. C. Hardin, and D. Rothenberg.  Physically regularized machine learning emulators of aerosol activation.  Geoscientific Model Development Discussions, 2020:1–19, 2020.
  * [14] A. A. Nair and F. Yu.  Using machine learning to derive cloud condensation nuclei number concentrations from commonly available measurements.  Atmospheric Chemistry and Physics, 20(21):12853–12869, 2020.
  * [15] Duncan Watson-Parris, Yuhan Rao, Dirk Olivié, Øyvind Seland, Peer J Nowack, Gustau Camps-Valls, Philip Stier, Shahine Bouabid, Maura Dewey, Emilie Fons, and et al.  Climatebench: A benchmark dataset for data-driven climate projections.  Earth and Space Science Open Archive, page 36, 2021.
  * [16] Salva Rühling Cachay, Venkatesh Ramesh, Jason N. S. Cole, Howard Barker, and David Rolnick.  ClimART: A benchmark dataset for emulating atmospheric radiative transfer in weather and climate models.  In Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track, 2021.
  * [17] Karthik Kashinath, M Mustafa, Adrian Albert, J.-L Wu, C Jiang, Soheil Esmaeilzadeh, Kamyar Azizzadenesheli, R Wang, and et al.  Physics-informed machine learning: Case studies for weather and climate modelling.  Philosophical transactions. Series A, Mathematical, physical, and engineering sciences, 379:20200093, 04 2021.
  * [18] Priya Donti, David Rolnick, and J Zico Kolter.  Dc3: A learning method for optimization with hard constraints.  In International Conference on Learning Representations, 2021.
  * [19] Tom Beucler, Michael Pritchard, Stephan Rasp, Jordan Ott, Pierre Baldi, and Pierre Gentine.  Enforcing analytic constraints in neural networks emulating physical systems.  Physical Review Letters, 126(9), Mar 2021.
  * [20] Paula Harder, Duncan Watson-Parris, Dominik Strassel, Nicolas Gauger, Philip Stier, and Janis Keuper.  Emulating aerosol microphysics with machine learning.  arXiv preprint arXiv:2109.10593, 2021.
  * [21] Paula Harder and Duncan Watson-Parris.  Aerosol microphysics emulation dataset, January 2022.
  * [22] S. Rasp.  Coupled online learning as a way to tackle instabilities and biases in neural network parameterizations: general algorithms and lorenz 96 case study (v1.0).  Geoscientific Model Development, 13(5):2185–2196, 2020.

##  5 Supplementary Material

###  5.1 Modelled Variables

Table 3 shows all the variables used for our emulator. Overall we consider 36
quantities, the first eight shown in the Table are only input variables and
are not changed by M7/our emulator. Masses and concentrations of different
aerosol species are both inputs and outputs for the model. Additionally, we
output the water content of different aerosol size bins.

Table 3: Modelled variables. Variable | Unit | Input | Output
---|---|---|---
Pressure | Pa | $\surd$ |
Temperature | K | $\surd$ |
Rel. Humidity | - | $\surd$ |
Ionization Rate | - | $\surd$ |
cloud cover | - | $\surd$ |
Boundary layer | - | $\surd$ |
Forest fraction | - | $\surd$ |
H2SO4 prod. rate | $cm^{-3}s^{-1}$ | $\surd$ |
H2SO4 mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
SO4 ns mass | $molec.\ m^{-3}$ | $\surd$ | $\surd$
SO4 ks mass | $molec.\ m^{-3}$ | $\surd$ | $\surd$
SO4 as mass | $molec.\ m^{-3}$ | $\surd$ | $\surd$
SO4 cs mass | $molec.\ m^{-3}$ | $\surd$ | $\surd$
bc ks mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
bc as mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
bc cs mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
bc ki mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
oc ks mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
oc as mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
oc cs mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
oc ki mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
du as mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
du cs mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
du ai mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
du ci mass | $\mu g\ m^{-3}$ | $\surd$ | $\surd$
ns concentration | $cm^{-3}$ | $\surd$ | $\surd$
ks concentration | $cm^{-3}$ | $\surd$ | $\surd$
as concentration | $cm^{-3}$ | $\surd$ | $\surd$
cs concentration | $cm^{-3}$ | $\surd$ | $\surd$
ki concentration | $cm^{-3}$ | $\surd$ | $\surd$
ai concentration | $cm^{-3}$ | $\surd$ | $\surd$
ci concentration | $cm^{-3}$ | $\surd$ | $\surd$
ns water | $kg\ m^{-3}$ |  | $\surd$
ks water | $kg\ m^{-3}$ |  | $\surd$
as water | $kg\ m^{-3}$ |  | $\surd$
cs water | $kg\ m^{-3}$ |  | $\surd$

###  5.2 Logarithmic Transformation

Training and evaluating on a logarithmic scale we can achieve much better
scores than in the original units. On the test set, we achieve an
$R^{2}$-score of $97.4\%$. The predictive quality on the logarithmic scale
though does not transfer to the original scale, as we can see in the scores in
Table LABEL:many_scores2.

Table 4: Test metrics for different architectures and transformations. Standard refers to learning with the standard transformation, log-scaled to the learning on logarithmitically transformed values. Base means the usage of only the base NN, correct adds the correction procedure and complete the completion procedure. Mass reg. and positivity reg. include the regularization terms. Best scores are bold and second best scores are highlighted in blue . | Standard | Log-Scaled
---|---|---
Architecture | Base | +Correct | +Complete | +Mass Loss | +Positivity Loss | Base | +Correct | +Complete | +Mass Loss
R2 | $0.763$ | 0.771 | $0.738$ | 0.730 | 0.709 | $-$ | $-$ | $-$ | $-$
R2 log | $-$ | $-$ | $-$ | $-$ | $-$ | $0.974$ | $0.974$ | 0.959 | 0.984
RMSE | 0.402 | 0.401 | 0.403 | 0.433 | 0.459 | 115431 | 115431 | 260706 | 21057757
Mass Bias SO4 | $1.1\cdot 10^{-5}$ | $8.5\cdot 10^{-5}$ | 0.00 | $8.6\cdot 10^{-6}$ | $1.0\cdot 10^{-3}$ | $1.2\cdot 10^{0}$ | $1.2\cdot 10^{0}$ | 0.00 | $1.1\cdot 10^{-3}$
Mass Bias BC | $3.8\cdot 10^{-5}$ | $1.4\cdot 10^{-4}$ | 0.00 | $3.4\cdot 10^{-5}$ | $3.6\cdot 10^{-4}$ | $2.6\cdot 10^{-4}$ | $2.6\cdot 10^{-4}$ | 0.00 | $3.2\cdot 10^{-4}$
Mass Bias OC | $3.3\cdot 10^{-5}$ | $6.0\cdot 10^{-5}$ | 0.00 | $1.1\cdot 10^{-5}$ | $6.4\cdot 10^{-4}$ | $3.8\cdot 10^{-4}$ | $1.4\cdot 10^{-4}$ | 0.00 | $8.0\cdot 10^{-5}$
Mass Bias DU | $1.0\cdot 10^{-6}$ | $3.9\cdot 10^{-5}$ | 0.00 | $2.8\cdot 10^{-7}$ | $1.5\cdot 10^{-5}$ | $1.0\cdot 10^{-3}$ | $2.2\cdot 10^{-4}$ | 0.00 | $8.5\cdot 10^{-5}$
Mass Violation | $3.7\cdot 10^{-4}$ | $1.1\cdot 10^{-3}$ | 0.00 | $1.4\cdot 10^{-4}$ | $2.4\cdot 10^{-3}$ | $4.0\cdot 10^{2}$ | $4.0\cdot 10^{2}$ | 0.00 | $3.2\cdot 10^{-1}$
Neg. Fraction | 0.134 | 0.00 | 0.146 | 0.144 | 0.0894 | 0.0807 | 0.00 | 0.0812 | 0.061
Neg. Mean | 0.00151 | 0.00 | 0.00170 | 0.00169 | 0.000081 | 0.00921 | 0.00 | 0.108 | 0.00115

###  5.3 Other ML Approaches

Apart from a neural network approach, we investigated other ML models to
emulate the M7 module. We looked at linear regression (LR) and different
ensemble methods such as a random forest regressor t (RF) and a gradient
boosting model (GB). In Table LABEL:other_ML_approaches we report the test
scores for linear regression, a random forest, and a gradient boosting
approach. Being linear combinations of training points, linear regression and
random forest both automatically conserve mass. Overall the random forest
shows the best accuracy among these models. We note that the random forest
performs worse than the neural network in terms of predictive accuracy, but
could still be considered for future work because of the mass conserving
property.

Table 5: Test metrics for different other ML models. | Accuracy | Mass Conservation | Positivity
---|---|---|---
Architecture | R2 | R2 log | RMSE | Bias SO4 | Bias BC | Bias OC | Bias DU | RMSE | Neg. Fraction | Neg. Mean
LR | +2.10e-1 | -4.19e+0 | +6.95e-1 | -4.21e-18 | -1.03e-16 | -1.08e-16 | +1.03e-17 | +4.16e-16 | +1.92e-1 | -1.32e-2
RF | +7.07e-1 | -7.30e-1 | +4.66e-1 | +9.90e-17 | -4.78e-16 | -6.19e-17 | -9.27e-19 | +3.85e-16 | +1.31e-1 | -7.28e-4
GB | +6.66e-1 | -2.95e-1 | +4.76e-1 | +1.90e-5 | +1.75e-5 | +1.45e-5 | -8.56e-5 | +1.14e-2 | +1.27e-1 | -1.25e-3

###  5.4 Network Architectures and Losses

The network architectures used all build on the same base neural network
containing two hidden layers and ReLU activations. For enforcing hard
constraints we add an additional layer at the end, either the correction layer
for the inequality constraint of positive masses or the completion layer of
the equality constraint for mass conservation. Schematics and the combination
of architectures and loss functions are shown in Figure 3.

Figure 3: Different architectures and loss functions used for comparison.

###  5.5 Activation Functions

We explored different activation functions for our neural architecture,
Sigmoid, Tanh, ReLU, and Leaky ReLU. We evaluated their performance (see Table
6) on the validation set and find that the ReLU-based activations result in
better MSE and $R^{2}$-scores.

Table 6: Validation scores for different activation functions.
Activation Function | MSE | $R^{2}$
---|---|---
sigmoid | 0.442 | 74.2
tanh | 0.446 | 73.4
relu | 0.395 | 78.9
leaky relu | 0.395 | 78.7

###  5.6 Parameters for Regularization

The neural network’s success is very sensitive to changes in the parameters of
both mass and positivity regularizers. The parameters we found to work the
best in this scenario are $\alpha=[10^{-7},2\cdot 10^{4},2\cdot
10^{3},10^{-1}]$ and $\beta=[10^{-11},10^{7},10^{7},10^{3},10^{-8},10^{1}]$
for the training with linear transformed values and
$\alpha=[10^{-8},10^{3},10^{4},10^{5}]$ for log-transformation. $\beta$ was
tuned for the different variable groups (SO4, BC, OC, DU, NUM, WAT) and not
each variable independently, which could lead to further improvement.

###  5.7 Transformation for Constraining

For clarity and simplicity, we omit the fact that we need to back-transform
the variables before calculating masses and negative fraction in the main
paper. The more technically precise formulation is introduced here. Let
$\mu_{x},\sigma_{x}\in\mathbb{R}^{32}$ be the mean and standard deviation
vectors over all training input data and
$\mu_{y},\sigma_{y}\in\mathbb{R}^{28}$ the mean over all training output data.
Then we set the function $g:\mathbb{R}^{28}\mapsto\mathbb{R}^{28}$ as the back
transformation for the output data, $g(y)=y*\sigma_{y}+\mu_{y}$ and
$h:\mathbb{R}^{32}\mapsto\mathbb{R}^{32}$ as the back transformation for the
input data, $h(x)=x*\sigma_{x}+\mu_{x}$. With that, the loss terms are defined
as:

$\mathcal{L}^{\text{(mass)}}(y_{\theta}):=\sum_{s\in S}\alpha_{s}|\sum_{i\in
I_{s}}g(y_{\theta}^{(i)})|,$ (6)

for the mass conservation and

$\mathcal{L}^{\text{(pos)}}(y_{\theta}):=\sum_{i=0}^{n}\beta_{i}ReLU(-(g(y_{\theta}^{(i)})+h(x_{i})))^{2}$
(7)

for the positivity term.

The same applies for hard-constraining. The equations for that are:

$y_{\theta}^{(i)}=\text{ReLU}((g(\tilde{y}_{\theta}^{(i)})+h(x_{i})))-h(x_{i})$
(8)

and

$y_{\theta}^{(j)}=(-\sum_{i\in
I_{s}\setminus\\{j\\}}g(y_{\theta}^{(i)})-\mu_{y_{j}})/\sigma_{y_{j}}.$ (9)

###  5.8 Metrics Development

In Figure 4 we show the development of MSE, $R^{2}$-score, average mass
conservation violation, and overall negative fraction for different losses
throughout training. We can observe that the accuracies expressed in MSE and
$R^{2}$-score develop very similar for all loss types. Adding the mass
violation term we can see an improvement in mass conservation, especially for
later layers, adding the positivity term we can see a decrease in negative
fraction for all the epochs. It is also noticeable, that the negative fraction
shows high variability between different epochs.

Figure 4: Development of MSE, $R^{2}$-score, mass conservation violation and
negative fraction over the course of training epochs, measured on the
validation set.

###  5.9 Results for Classification

To enable log-transformed learning we need to combine our regression neural
network with a classification network. The classification network predicts
whether a tendency is positive or negative. The scores for this binary
classification task can be found in Table LABEL:class_scores. Note that this
classification problem is highly imbalanced and could be improved by using
techniques for this case.

Table 7: Scores for classification neural network, calculated on test set. Variable | SO4 NS | SO4 KS | SO4 AS | BC KS | BC AS | OC KS | OC AS | DU AS | NUM NS | NUM KS | NUM AS
---|---|---|---|---|---|---|---|---|---|---|---
Accuracy | 0.986 | 0.962 | 0.999 | 0.960 | 0.999 | 0.962 | 0.999 | 0.999 | 0.989 | 0.984 | 0.959
Precision | 0.954 | 0.931 | 0.999 | 0.914 | 0.999 | 0.908 | 0.999 | 1.000 | 0.972 | 0.755 | 0.919
Recall | 0.990 | 0.988 | 1.000 | 0.972 | 1.000 | 0.968 | 1.000 | 1.000 | 0.989 | 0.844 | 0.768
Table 8: $R^{2}$ scores on the test set for each variable. For our neural network that learns the log-transformed tendencies we consider the $R^{2}$ score on the log-transformed values, for the neural network trained on standardized values we consider the $R^{2}$ in original units. | SO4 | Black Carbon | Organic Carbon
---|---|---|---
Variable | H2SO4 | NS | KS | AS | CS | KS | AS | CS | CI | KS | AS | CS | CI
Log. ($R^{2}$) | 0.999 | 0.993 | 0.981 | 0.993 | 0.992 | 0.979 | 0.995 | 0.990 | 0.997 | 0.984 | 0.997 | 0.994 | 0.998
Stand. ($R^{2}$) | 0.983 | 0.877 | 0.652 | 0.875 | 0.292 | 0.837 | 0.844 | 0.000 | 0.867 | 0.616 | 0.669 | 0.200 | 0.710
Table 9: Same as table above for additional variables. | Dust | Number particles | Water content
---|---|---|---
Variable | AS | CS | AI | CI | NS | KS | AS | CS | KI | AI | CI | NS | KS | AS | CS
Log. ($R^{2}$) | 0.995 | 0.996 | 0.995 | 0.996 | 0.993 | 0.997 | 0.998 | 0.997 | 0.997 | 0.977 | 0.983 | 0.979 | 0.982 | 0.979 | 0.982
Stand. ($R^{2}$) | 0.872 | 0.966 | 0.945 | 0.966 | 0.900 | 0.663 | 0.874 | 0.932 | 0.949 | 0.936 | 0.940 | 0.979 | 0.978 | 0.928 | 0.954

###  5.10 Individual Scores and Plots

In Table LABEL:scores_var1 and LABEL:scores_var2 we present the $R^{2}$ scores
for all predicted variables individually, using the base neural network for
logarithmic transformed or standardized values. Figure 5 and 6 show plots for
all predicted variables.

Figure 5: Predicted test values by emulator versus true values from M7 model.
Trained on logarithmically transformed quantities and plotted in the
logarithmic scale

Figure 6: Predicted test values by emulator versus true values from M7 model.
Trained on standardized quantities and plotted standardized