# Spectral energy density in an axisymmetric galaxy as predicted by an
analytical model for the Maxwell field

Mayeul Arminjon
Univ. Grenoble Alpes, CNRS, Grenoble INP, 3SR, F-38000 Grenoble, France
E-mail<EMAIL_ADDRESS>Orcid: 0000-0002-7035-351X

###### Abstract

An analytical model for the Maxwell radiation field in an axisymmetric galaxy,
proposed previously, is first checked for its predictions of the spatial
variation of the spectral energy distributions (SEDs) in our Galaxy. First,
the model is summarized. It is now shown how to compute the SED with this
model. Then the model is adjusted by asking that the SED predicted at our
local position in the Galaxy coincide with the available observations. Finally
the first predictions of the model for the spatial variation of the SED in the
Galaxy are compared with those of a radiation transfer model. We find that the
two predictions do not differ too much. This indicates that, in a future work,
it should be possible with the present model to check if the “interaction
energy” predicted by an alternative, scalar theory of gravitation, contributes
to the dark matter.

Keywords: Disk galaxy; interstellar radiation field; axial symmetry;
electromagnetic field; Milky Way.

Short title: Spectral energy density in a galaxy from an exact Maxwell field

##  1 Introduction

The interstellar radiation field (ISRF) is a very important physical
characteristic of a galaxy, and there is indeed a rich astrophysical
literature that studies the modeling of the ISRF; see e.g. [1, 2, 3, 4, 5, 6]
and references therein. This subject is of importance not only in itself, but
also because the radiation field produced by the stars and other bright
objects has a strong interaction with the cosmic rays, which are another
important component of the interstellar medium. The necessity of a modeling is
obvious, since we can directly measure the radiation field only in the solar
system, in fact in the neighborhood of the Earth [7, 8, 9]. The ISRF is
significantly affected by its interaction with interstellar dust (e.g. [10,
11]). In order to account for this interaction, equations for radiative
transfer are solved: either by using Monte-Carlo methods, that simulate in a
probabilistic framework the propagation of photons, and their absorption and
scattering by dust, in a discretized spatial domain (e.g. [1, 4]); or, by
using “ray-tracing” algorithms that calculate the variation, due to absorption
and scattering by dust, of the radiation’s specific intensity along a finite
set of directions (e.g. [3, 4, 6]). Thus, until now, the electromagnetic (EM)
field, made of coupled electric and magnetic fields, has not been considered
in the studies of the ISRF in galaxies, even less an EM field that would be a
solution of the Maxwell equations — as should actually be the case for the
real ISRF, however.

Obtaining for the ISRF such an exact solution of the Maxwell equations is
somewhat demanding from both the mathematical and the numerical points of
view, and has been the main purpose of a recent work [12]. A decisive
motivation for that recent work was to prepare for checking an exotic
prediction of an alternative, scalar theory of gravitation: namely, that the
total energy(-momentum-stress) tensor, which provides the source of the
gravitational field, is not in general the sum of the energy tensor of the EM
field and the energy tensor of the material medium producing it. Instead, each
time that an EM field is present together with a variable gravitational field
(which is the general case), there should appear an additional energy tensor:
the “interaction” energy tensor, ${{\boldsymbol{T}}}_{\mathrm{inter}}$ [13].
This speculative energy tensor ${{\boldsymbol{T}}}_{\mathrm{inter}}$, being
thus present virtually everywhere according to that theory, and having clearly
an exotic character, could contribute to the dark matter. In order to
calculate ${{\boldsymbol{T}}}_{\mathrm{inter}}$, and to see if its
distribution could be similar to dark matter halos, it is necessary to be able
to calculate the EM field and its first order derivatives [13], whence the
work [12]. In the latter work, a model has been formulated. It is based on an
analytical representation of the general solution of the source-free Maxwell
equations in the axisymmetric case [14]. Moreover, the model has been
numerically implemented and has passed a numerical validation test based on an
exact solution [12].

The aim of the present paper is to make a first check of the predictions of
that “Maxwell model of the interstellar radiation field” for the spectral
energy distribution (SED) and its spatial variation in a galaxy. Such SEDs are
a very important output of the mainstream models for the analysis of the ISRF
and, as was recalled at the beginning of this paper, those models depend
essentially of the radiation transfer (mainly through dust). Note that, in
contrast, the present model [12] does not directly take into account any
radiation transfer. Section 2 summarizes that model. Section 3 details the
calculation of the SED in that model. Section 4 presents the results obtained
with that model for a simulation of the Milky Way. To obtain these results,
the model is adjusted on the locally observed SED. In this way, it does
indirectly take into account the radiation transfer, because the observed SED
is indeed affected by the latter.

##  2 The Maxwell model of the interstellar radiation field

In this section, we give a new presentation of the model which was built in
Ref. [12], and that uses results of Ref. [14] and references therein. This
presentation uses an easy extension of the result [14]. Such an extension was
used implicitly in Ref. [12], and is explicitly formulated in Appendix A.

###  2.1 General assumptions for the Maxwell field

We make two assumptions regarding the form of the Maxwell radiation field
$({\bf E},{\bf B})$ (with ${\bf E}$ the electric field and ${\bf B}$ the
magnetic field):

i) We assume that it is an axially symmetric field. This is a relevant first
approximation for many galaxies, even though it neglects some aspects like the
arm structure of a spiral galaxy. Axial symmetry is often assumed in the
literature on modeling the ISRF in galaxies and/or its interaction with dust
(e.g. [2, 3, 6, 10]).

ii) We search for a source-free solution of the Maxwell equations. This is
because we seek to describe the EM field at the galactic scale, not the local
EM field in the neighborhood of the stars that are the primary sources of that
field.

###  2.2 The explicit representation for axisymmetric free Maxwell fields

This is based on the following recent Theorem [14]: any time-harmonic
axisymmetric solution $({\bf E},{\bf B})$ of the source-free Maxwell equations
has a unique decomposition

${\bf E}={\bf E}_{1}+{\bf E}^{\prime}_{2},\qquad{\bf B}={\bf B}_{1}+{\bf
B}^{\prime}_{2},$ (1)

where $({\bf E}_{1},{\bf B}_{1})$ has the form (in cylindrical coordinates
$\rho,\phi,z$ corresponding with the symmetry axis $Oz$):

$\displaystyle{\color[rgb]{0.,0.,0.}B_{\phi}}$ $\displaystyle=$
$\displaystyle{\color[rgb]{0.,0.,0.}-\frac{\partial\Psi}{\partial\rho}},\qquad{\color[rgb]{0.,0.,0.}E_{\phi}=0},$
(2) $\displaystyle{\color[rgb]{0.,0.,0.}E_{\rho}}$ $\displaystyle=$
$\displaystyle{\color[rgb]{0.,0.,0.}\mathrm{i}\frac{c^{2}}{\omega}\frac{\partial^{2}\Psi}{\partial\rho\,\partial
z}},\qquad{\color[rgb]{0.,0.,0.}B_{\rho}=0},$ (3)
$\displaystyle{\color[rgb]{0.,0.,0.}E_{z}}$ $\displaystyle=$
$\displaystyle{\color[rgb]{0.,0.,0.}\mathrm{i}\frac{c^{2}}{\omega}\frac{\partial^{2}\Psi}{\partial
z^{2}}+\mathrm{i}\omega\Psi},\qquad{\color[rgb]{0.,0.,0.}B_{z}=0};$ (4)

and where $({\bf E}^{\prime}_{2},{\bf B}^{\prime}_{2})$ is deduced from a
solution $({\bf E}_{2},{\bf B}_{2})$ having just that same form (2)–(4), by
the EM duality — i.e.,

${\bf E}^{\prime}=c{\bf B},\quad{\bf B}^{\prime}=-{\bf E}/c.$ (5)

(Here $c$ is the velocity of light.) Note that the dual solution (5) has the
complementary zero and non-zero components as compared with the direct
solution (2)–(4), i.e.,
$B^{\prime}_{\phi}=E^{\prime}_{\rho}=E^{\prime}_{z}=0$, while
$E^{\prime}_{\phi},B^{\prime}_{\rho},B^{\prime}_{z}\neq 0$ in general. Both
solutions $({\bf E}_{1},{\bf B}_{1})$ and $({\bf E}^{\prime}_{2},{\bf
B}^{\prime}_{2})$ are time-harmonic with the same angular frequency $\omega$
as has $({\bf E},{\bf B})$. In Eqs. (2)–(4), the potential
$\Psi(t,\rho,z)=e^{-\mathrm{i}\omega t}\hat{\Psi}(\rho,z)$ is a time-harmonic
axisymmetric solution of the standard scalar wave equation, which also has the
same frequency $\omega$ as has $({\bf E},{\bf B})$; in general, the two
solutions $({\bf E}_{1},{\bf B}_{1})$ and $({\bf E}_{2},{\bf B}_{2})$
correspond with two different scalar potentials $\Psi$ and $\Psi^{\prime}$
[14]. The solution (2)–(4) actually derives from a vector potential having
just its axial component $A_{z}=\Psi$ non-zero, i.e., ${\bf A}=\Psi{\bf
e}_{z}$ [14, 15]. The Theorem [14] states the existence of the two potentials
$\Psi$ and $\Psi^{\prime}$, not their uniqueness. In Appendix A, we extend
this result to EM fields having a finite spectrum, which is easy, and allows
us to give a more rigorous presentation of the model in Subsect. 2.3.

If a time-harmonic, axisymmetric solution of the scalar wave equation is
totally propagating, it has the following explicit form [15, 16]:

${\color[rgb]{0.,0.,0.}\psi_{\omega\ S}\,(t,\rho,z)=e^{-\mathrm{i}\omega
t}\int_{-K}^{+K}\ J_{0}\left(\rho\sqrt{K^{2}-k^{2}}\right)\
e^{\mathrm{i}k\,z}\,S(k)\,\mathrm{d}k},$ (6)

with $\omega$ the angular frequency, $K:=\omega/c$, $J_{0}$ the first-kind
Bessel function of order $0$, and where $S$ is some (generally complex)
function of $k\in[-K,+K]$. “Totally propagating” means, in short, that $\Psi$
has no “evanescent” component [15]; a detailed explanation is given in Ref.
[12]. In turn, we shall define simply that a time-harmonic, axisymmetric EM
field ${{\boldsymbol{F}}}=({\bf E},{\bf B})$ is totally propagating if the two
potentials $\Psi$ and $\Psi^{\prime}$, whose existence is guaranteed by the
Theorem recalled above, can themselves be chosen totally propagating, i.e.,
having the form (6).

###  2.3 Statement of the model

We model a disk galaxy as a disk-like set of point-like “stars”. 111 The
corresponding points are got by pseudo-random generation of their cylindrical
coordinates. Here, the same set of $16\times 16\times 36$ triplets
$(\rho,z,\phi)$ was used as in Ref. [12], thus $9216$ points. The distribution
of $\rho$ and $z$ is approximately that valid for the star distribution in the
Milky Way, and the distribution of $\phi$ is uniform, thus ensuring
approximate axial symmetry.  Each of them is assumed to produce an EM field
with a finite spectrum: that field derives, in the way defined in the previous
section, from potentials that obey the scalar wave equation, more precisely
ones with spherical symmetry:

$\psi_{\omega_{j}}\ (t,{\bf
x})=\frac{e^{\mathrm{i}(K_{j}r-\omega_{j}t)}}{K_{j}r},\qquad
K_{j}:=\frac{\omega_{j}}{c},\qquad r:=\left|{\bf x}\right|.$ (7)

Thus, the sum of the potentials emitted by the different stars at the
different frequencies is:

$\Psi_{({\bf x}_{i})\ (\omega_{j})\ (w_{j})}\,(t,{\bf
x})=\sum_{i=1}^{i_{\mathrm{max}}}\sum_{j=1}^{N_{\omega}}\,w_{j}\,\frac{e^{\mathrm{i}(K_{j}\,r_{i}-\omega_{j}t)}}{K_{j}\
r_{i}}.$ (8)

Here $\,r_{i}:=\left|{\bf x}-{\bf x}_{i}\right|$, with $\,{\bf x}\,$ the
spatial position at which the function is evaluated and $\,{\bf x}_{i}$ the
spatial position of the $i$-th star; the numbers $\,w_{j}>0$, with
$\,\sum_{j}w_{j}=1$, are the weights affected to the different frequencies.
222 In Ref. [12], these weights were noted $S^{\prime}_{j}$, which now would
clash with the functions $S^{\prime}_{j}$ mentioned at the end of the present
paper.  Since the disk-like distribution of the “stars” is built
(approximately) axisymmetric [12], the function (8) is also (approximately)
axisymmetric. Moreover, this is a solution of the scalar wave equation, and it
is clearly a totally propagating one, as is each of the individual spherical
scalar waves in (8). It should thus be possible to put the function (8) — at
least approximately — in the form of a sum, over the frequencies, of functions
having the form (6). This is done by the least-squares method applied on a
finite set of events $\,(t,{\bf x})$ (with ${\bf x}={\bf x}(\rho,\phi=0,z)$
restricted to the plane $\phi=0$, due to the axial symmetry), making a regular
spatio-temporal grid $\\{(t_{l},\rho_{m},z_{p})\\}$. We thus want to determine
complex functions $S_{j}=S_{j}(k)\ (-K_{j}\leq k\leq+K_{j})$, such that

$\Psi_{({\bf x}_{i})\ (\omega_{j})\
(w_{j})}\,(t,\rho,z)\cong\sum_{j}\,\psi_{\omega_{j}\
S_{j}}\,(t,\rho,z):=\Psi_{(\omega_{j})\ (S_{j})}\,(t,\rho,z),$ (9)

where the sign $\cong$ indicates that the equality is in the sense of the
least squares (on the relevant spatio-temporal grid, just mentioned). To the
unique function on the r.h.s. of (9): $\Psi=\sum_{j}\,\psi_{\omega_{j}\
S_{j}}$, we may then associate:

  * •

either the EM field obtained by summing the time-harmonic contributions given
by Eqs. (2)–(4), with $A_{z}:=\psi_{\omega_{j}\ S_{j}}$ — that is,

${{\boldsymbol{F}}}_{1}:=Z_{1}(\Psi)$ (10)

(cf. Eq. (35) in the Appendix);

  * •

or, the dual (5) of that field — that is,

${{\boldsymbol{F}}}_{2}:=Z_{2}(\Psi);$ (11)

  * •

or still, the sum of these two fields, that is

${{\boldsymbol{F}}}={{\boldsymbol{F}}}_{1}+{{\boldsymbol{F}}}_{2}=Z_{1}(\Psi)+Z_{2}(\Psi):=Z(\Psi,\Psi)$
(12)

(cf. Eq. (36) in the Appendix).

###  2.4 Discretization

To calculate the integrals in Eq. (6) and in Eqs. (2)–(4) applied with
functions $\psi_{\omega_{j}\ S_{j}}$ of the form (6), we discretize the
integration intervals $[-K_{j},+K_{j}]$ and the functions $S_{j}$, and we use
the “Simpson $\frac{3}{8}$ composite rule” [12]. We set $K_{0}:=\omega_{0}/c$,
with $\omega_{0}$ some (arbitrary) reference frequence, and

$k_{n}:=-K_{0}+n\delta_{0}\quad(n=0,...,N),\qquad\delta_{0}:=2K_{0}/N,$ (13)
$S_{nj}:=S_{j}\left(\frac{\omega_{j}}{\omega_{0}}k_{n}\right)\quad(n=0,...,N,\quad
j=1,...,N_{\omega}),$ (14) $\displaystyle a_{n}$ $\displaystyle:=$
$\displaystyle(3/8)\,\delta_{0}\quad(n=0\ \mathrm{or}\ n=N),$ (15)
$\displaystyle a_{n}$ $\displaystyle:=$ $\displaystyle
2\times(3/8)\,\delta_{0}\quad(\mathrm{mod}(n,3)=0\ \mathrm{and}\ n\neq 0\
\mathrm{and}\ n\neq N),$ (16) $\displaystyle a_{n}$ $\displaystyle:=$
$\displaystyle 3\times(3/8)\,\delta_{0}\quad\mathrm{otherwise}.$ (17)

With this discretization, the unknown functions $S_{j}\ (j=1,...,N_{\omega})$
in Eq. (9) are characterized by the complex numbers $S_{nj}\
(n=0,...,N;\,j=1,...,N_{\omega})$. Thus, the equality in the sense of the
least squares, schematized by Eq. (9), actually becomes [12]

$\Psi_{({\bf x}_{i})\ (\omega_{j})\
(w_{j})}\,(t,\rho,z)\cong\sum_{n=0}^{N}\sum_{j=1}^{N_{\omega}}f_{nj}(t,\rho,z)\,S_{nj},$
(18)

with

$f_{nj}(t,\rho,z)=\frac{\omega_{j}}{\omega_{0}}a_{n}\,J_{0}\left(\rho\frac{\omega_{j}}{\omega_{0}}\sqrt{K_{0}^{2}-k_{n}^{2}}\right)\exp\left[\mathrm{i}\left(\frac{\omega_{j}}{\omega_{0}}k_{n}z-\omega_{j}\,t\right)\right].$
(19)

The $S_{nj}$ ’s are the solved-for parameters in the least-squares problem
(18).

This gives us for the EM field (10) [12], as also for the field (12):

$B_{\phi}(t,\rho,z)=\sum_{n=0}^{N}\sum_{j=1}^{N_{\omega}}R_{n}\,J_{1}\left(\rho\frac{\omega_{j}}{\omega_{0}}R_{n}\right)\,{\mathcal{R}e}\left[F_{nj}(t,z)\right]+O\left(\frac{1}{N^{4}}\right),$
(20)
$E_{\rho}(t,\rho,z)=\sum_{n=0}^{N}\sum_{j=1}^{N_{\omega}}\frac{c^{2}}{\omega_{0}}k_{n}R_{n}\,J_{1}\left(\rho\frac{\omega_{j}}{\omega_{0}}R_{n}\right)\,{\mathcal{R}e}\left[F_{nj}(t,z)\right]+O\left(\frac{1}{N^{4}}\right),$
(21)
$E_{z}(t,\rho,z)=\sum_{n=0}^{N}\sum_{j=1}^{N_{\omega}}\left(\frac{c^{2}}{\omega_{0}}k_{n}^{2}-\omega_{0}\right)\,J_{0}\left(\rho\frac{\omega_{j}}{\omega_{0}}R_{n}\right)\,{\mathcal{I}m}\left[F_{nj}(t,z)\right]+O\left(\frac{1}{N^{4}}\right),$
(22)

with $R_{n}=\sqrt{K_{0}^{2}-k_{n}^{2}}$ and

$F_{nj}(t,z)=\left(\frac{\omega_{j}}{\omega_{0}}\right)^{2}\,a_{n}\exp\left[\mathrm{i}\left(\frac{\omega_{j}}{\omega_{0}}k_{n}z-\omega_{j}\,t\right)\right]\,S_{nj}.$
(23)

Even though these formulae follow from a discretization and from using an
integration rule that is only approximate for the starting integral, they
provide (without the $O\left(\frac{1}{N^{4}}\right)$ remainder) an exact
solution of the source-free Maxwell equations — see §3.3 in Ref. [12].

The arguments of the Bessel functions $J_{0}$ and $J_{1}$ in Eqs. (19),
(20)–(22), as well as the spatial part of the argument of the complex
exponential in Eqs. (19) and (23), are of the order of $\rho/\lambda$ or
$z/\lambda$. For a galaxy, these are numbers of the order of $10^{25}$,
therefore the numerical implementation of the model needs to use a precision
better than quadruple precision — which leads to long computation times [12].

##  3 Calculation of the spectral energy density

The volumic energy density of an EM field, in SI units, is

$U=\frac{\delta W}{\delta V}=\frac{\epsilon_{0}}{2}\left({\bf E}^{2}+c^{2}{\bf
B}^{2}\right),$ (24)

with $\epsilon_{0}=1/(4\pi\times 9\times 10^{9})$. In the relevant case with a
frequency spectrum, the fields corresponding to the different frequencies do
add, hence the quadratic expression (24) is not additive. The time average of
such a quadratic expression is yet known to be additive in the case of a
Fourier expansion [17]. Let us check if this remains true in the present case
of an arbitrary finite frequency spectrum $(\omega_{j})\
(j=1,...,N_{\omega})$. Let $F(t)$ be some component of the EM field at some
given spatial position ${\bf x}$. We have, with some complex coefficients
$C_{j}=C_{j}({\bf x})$:

$F(t)={\mathcal{R}e}\left(\sum_{j=1}^{N_{\omega}}C_{j}e^{-\mathrm{i}\omega_{j}t}\right),$
(25)

whence

$4F^{2}=\left(\sum_{j=1}^{N_{\omega}}C_{j}e^{-\mathrm{i}\omega_{j}t}+C_{j}^{\ast}e^{\mathrm{i}\omega_{j}t}\right)\left(\sum_{k=1}^{N_{\omega}}C_{k}e^{-\mathrm{i}\omega_{k}t}+C_{k}^{\ast}e^{\mathrm{i}\omega_{k}t}\right).$
(26)

Unless $\omega_{j}+\omega_{k}=0$, the time average
$\overline{e^{-\mathrm{i}\left(\omega_{j}+\omega_{k}\right)t}}$ is zero.
Unless $\omega_{j}-\omega_{k}=0$, the time average
$\overline{e^{-\mathrm{i}\left(\omega_{j}-\omega_{k}\right)t}}$ is zero. But
$\omega_{j}+\omega_{k}=0$ does not occur, because $\omega_{j}>0$, and
$\omega_{j}-\omega_{k}=0$ occurs only when $j=k$. Hence, (26) gives indeed the
additivity for the time average of the square of a component admitting an
expansion (25):

$\overline{F^{2}}=\frac{1}{4}\sum_{j=1}^{N_{\omega}}C_{j}C_{j}^{\ast}+C_{j}^{\ast}C_{j}=\frac{1}{2}\sum_{j=1}^{N_{\omega}}\left|C_{j}\right|^{2}.$
(27)

For the EM field (10), its components $B_{\phi},\,E_{\rho},\,E_{z}$ admit
expansions (25) with coefficients $C^{(1)}_{j}$, $C^{(2)}_{j}$, $C^{(3)}_{j}$,
respectively, which are directly read on Eqs. (20)–(22):

$C^{(1)}_{j}(\rho,z)=\sum_{n=0}^{N}R_{n}\,J_{1}\left(\rho\frac{\omega_{j}}{\omega_{0}}R_{n}\right)G_{nj},$
(28)
$C^{(2)}_{j}(\rho,z)=\sum_{n=0}^{N}\frac{c^{2}}{\omega_{0}}k_{n}R_{n}\,J_{1}\left(\rho\frac{\omega_{j}}{\omega_{0}}R_{n}\right)G_{nj},$
(29)
$C^{(3)}_{j}(\rho,z)=\sum_{n=0}^{N}\left(\frac{c^{2}}{\omega_{0}}k_{n}^{2}-\omega_{0}\right)\,J_{0}\left(\rho\frac{\omega_{j}}{\omega_{0}}R_{n}\right)\left(-\mathrm{i}G_{nj}\right),$
(30)

with $R_{n}=\sqrt{K_{0}^{2}-k_{n}^{2}}$ and

$G_{nj}(z)=\exp\left(\mathrm{i}\omega_{j}\,t\right)F_{nj}(t,z)=\left(\frac{\omega_{j}}{\omega_{0}}\right)^{2}\,a_{n}\exp\left(\mathrm{i}\frac{\omega_{j}}{\omega_{0}}k_{n}z\right)\,S_{nj}.$
(31)

From (24) and (27), we thus get for the field (10) a (time-averaged) energy
density with a discrete spectrum:

$\overline{U}(\rho,z):=\overline{\frac{\delta W}{\delta
V}}=\sum_{j=1}^{N_{\omega}}u_{j},\qquad
u_{j}:=\frac{1}{4}\sum_{q=1}^{3}\alpha_{q}\left|C^{(q)}_{j}\right|^{2},$ (32)

where the coefficients $C^{(q)}_{j}\ (q=1,2,3)$ are given by Eqs. (28)–(30),
$\alpha_{1}=\epsilon_{0}c^{2}$, $\alpha_{2}=\alpha_{3}=\epsilon_{0}$. As to
the field (11), it is the dual (5) of the field (10). Therefore, the energy
density (24) for the field (11) is given by the same Eq. (32). Finally, the
energy density for the field (12) is hence just the double of this. Now, in
the next section, the model will be adjusted so as to exactly describe the
measured local EM spectrum, i.e., the values $u_{j}({\bf x}_{\mathrm{loc}})$
will be imposed to be the ones measured at our local position ${\bf
x}_{\mathrm{loc}}$ in the Galaxy. In this context, the precise choice of the
model, i.e., Eq. (10), or (11), or (12), is therefore totally neutral.

##  4 Results and comparison with a radiation transfer model

We first explain the final adjustment of the model. The scalar wave
$\Psi_{({\bf x}_{i})\ (\omega_{j})\ (w_{j})}$ given by Eq. (8) depends on the
numbers $w_{j}>0$, with $\sum_{j}w_{j}=1$. Each of them represents the
relative weight of the corresponding frequency $\omega_{j}$ — though in some
average sense, since $w_{j}$ is common to all “stars” at the points ${\bf
x}_{i}$. There is no exact way to assign values to these relative weights
$w_{j}$, which, together with the set $({\bf x}_{i})$ of the positions of the
stars, determine the functions $S_{j}$ on the r.h.s. of Eq. (9) — or more
exactly, determine the discretized values $S_{nj}$ in Eqs. (14) and (18). In
turn, the values $S_{nj}\ (n=0,...,N)$, for a given value of the frequency
index $j$, determine the intensity of the $\omega_{j}$-component of the scalar
wave obtained after the fitting (the r.h.s. of (18)), and thus determine the
intensity of the $\omega_{j}$-component of the EM field (20)–(22). However,
since we start only from relative weights $w_{j}$, this determination can be
only up to a multiplying factor, and since the very values of the $w_{j}$ ’s
are assigned somewhat arbitrarily, that determination of the $S_{nj}$ ’s $\
(n=0,...,N)$ is really up to a factor depending on $j$, $\xi_{j}>0$. 333 This
would be even clearer if we would separate the least-squares adjustment (18)
into the different frequencies $\omega_{j}$, which is actually easy and even
would allow us to eliminate the common time dependence
$\exp(-\mathrm{i}\omega_{j}t)$ on both sides, thus eliminating the time
variable in the adjustment. This different procedure will be numerically
tested in a future work.  We determine the values $\xi_{j}$ by imposing that
the values $u_{j}$ given by Eq. (32), calculated at our local position ${\bf
x}_{\mathrm{loc}}$ in the Galaxy, be equal to the ones measured. We take
$\rho_{\mathrm{loc}}=8$ kpc and $z_{\mathrm{loc}}=0.02$ kpc (e.g. [18]). We
took the local values of the SED from Fig. 1 of Ref. [2], using the software
CurveUnscan to digitalize the curves. Those local values were determined from
the Apollo 17 mission [7] and from the COBE mission: the DIRBE experiment [8]
and the FIRAS experiment [9]. Figure 1 shows the corresponding SED. For the
calculations, in order to save computer time, we extracted 76 values of
$\lambda$ from the complete list of values of $\lambda$ in the curve; the
corresponding points are also shown on Fig. 1.

Figure 1: SED as measured by space missions. Between the data of Apollo
(ending at $\mathrm{log}_{10}\lambda\simeq-0.8$) and DIRBE (beginning at
$\mathrm{log}_{10}\lambda\simeq+0.1$), the curve has been interpolated.

Then, at any other spatial position ${\bf x}(\rho,z)$ (in the axially
symmetric case considered), the model gives a prediction for the values
$u_{j}({\bf x})$ of the SED corresponding to the frequencies $\omega_{j}$ (or
the wavelengths $\lambda_{j}$) considered in our calculation. The spatio-
temporal grid used was $N_{t}=5,\quad N_{\rho}=10,\quad N_{z}=21$ (thus 1050
events $(t,\rho,z)$), with $t$ varying from $0$ to $T_{0}=\lambda_{0}/c$ with
$\lambda_{0}=10\,\mu\mathrm{m}$, $\rho$ varying from $0$ to $10\,$kpc, and $z$
from $-1\,$ to $+1\,$kpc. Other grids have also been tested, giving very
similar results. Moreover, unless otherwise mentioned, the discretization
number $N$ was $N=48$. (This is the number of subintervals in which any of the
integration intervals $[-K_{j},+K_{j}]$ is subdivided, see Eqs. (13)–(14).)
Figures 2 to 5 show the comparison of the SEDs, obtained at four different
places in our Galaxy, according to either the thus-adjusted present model, or
the radiation transfer model of Popescu et al. [6]. (Again, the relevant
curves, Fig. 9 in Ref. [6], have been digitalized using CurveUnscan.)

Figure 2: SED at $\rho=8\,$kpc, $z=0\,$kpc, according to this model or a
radiation transfer model.

Figure 3: SED at $\rho=8\,$kpc, $z=1\,$kpc, according to this model or a
radiation transfer model.

Figure 4: SED at $\rho=1\,$kpc, $z=0\,$kpc, according to this model or a
radiation transfer model.

Figure 5: SED at $\rho=1\,$kpc, $z=1\,$kpc, according to this model or a
radiation transfer model.

Clearly, the general agreement between these two calculations is far from
perfect. In particular, our calculation oscillates markedly as function of
$\lambda$: at all wavelengths, for the two positions that are closer to the
galactic centre ($\rho=1\,$kpc, Figs. 4 and 5); and merely at short
wavelengths, for the two positions that are farther from the galactic centre
($\rho=8\,$kpc, Figs. 2 and 3). Nevertheless, this agreement is actually
surprisingly good if one remembers how different are the two models. The SED
predicted by our model is higher than the one predicted by the radiation
transfer model of Ref. [6] at the short wavelengths ($\lambda\lesssim
0.5\mu$m, or $\log_{10}\lambda\lesssim-0.3$) — except for the position
($\rho=8\,$kpc, $z=0\,$kpc), that is, relatively speaking, very close to our
local position ($\rho_{\mathrm{loc}}=8\,$kpc, $z_{\mathrm{loc}}=0.02\,$kpc),
where the adjustment is done. The higher levels found with our model at the
short wavelengths might have a relation with the fact that the absorption by
dust, not taken into account by our model, is strong for these radiations.
However, this relation is not obvious, since our model is adjusted on the
measured SED at our local position, and that local SED is of course affected
by the dust absorption to which the radiations directed towards here have been
subjected along their paths in the Galaxy. Actually, one notes also that, in
the long wavelength domain as well as in the very short wavelength domain, the
SED predicted by this model is, in average, higher than that predicted by the
radiation transfer model, at the two positions which have a higher altitude
($z=1\,$kpc): Figs. 3 and 5. Comparing Figs. 2 and 3, and comparing Figs. 4
and 5, it appears that the present model with its current numerical settings
predicts a slower decrease of the energy density as $z$ (or rather
$\left|z\right|$) increases than does the radiation transfer model. Of course
the two models are very different, also in their modelling of the Galaxy: a
discrete star distribution (for the present model) vs. an analytical,
continuous variation of the stellar emissivity as function of $\rho$ and $z$
(for the radiation transfer model of Popescu et al. [6]). However, the $\rho$
and $z$ distribution of the stars in our model has been built to correspond
(at least approximately) with that for the Milky Way [12], as has been built
the $\rho$ and $z$ dependence of the stellar emissivity function for the
radiation transfer model [6]. Of course also, the SEDs of a radiation transfer
model (like that of Ref. [6]) are a prediction, too, not a measurement.
Nevertheless, a further investigation is necessary to quantify the slower
decrease of the SED with increasing $z$ for the present model and to
understand the reason for it.

Figure 6: Predicted maximum values of the SED: $N_{\omega}=23$, three
different grids.

A surprising prediction of our model is that for the maximum values
$u_{j\,\mathrm{max}}$ of the SED in our model of the Galaxy — the SED being
calculated at each point on a regular $(\rho,z)$ grid, as described above.
Figure 6 shows these maximum values for three different spatial grids, with
the same set of $23$ wavelengths or frequencies. It is seen that the
$u_{j\,\mathrm{max}}$ values are quite similar for the different grids, and
these values are extremely high at the short wavelengths, i.e., at high
energies of the radiation field. They are found on the $\rho=0$ axis, thus on
the central axis of the galaxy. (Actually, the values $u_{j}(\rho=0,z)$ are
fairly independent of the value of $z$ in the domain investigated, i.e.,
$-1\,\mathrm{kpc}\leq z\leq 1\,\mathrm{kpc}$.) However, if one refines the
frequency mesh, going from $N_{\omega}=23$ to $N_{\omega}=46$, and then to
$N_{\omega}=76$, one observes that the very high values of
$u_{j\,\mathrm{max}}$ are obtained only at shorter and shorter wavelengths
(Figs. 7 and 8). A similar phenomenon is observed when one refines the
discretization parameter $N$: for a given frequency mesh, increasing $N$
reduces the domain of the high values of $u_{j\,\mathrm{max}}$ towards smaller
values of $\lambda_{j}$, or, equivalently, decreasing $N$ expands that domain
towards larger values of $\lambda_{j}$ (Figs. 9–11). The likely common reason
for the two effects is that, when one refines either the frequency mesh
($N_{\omega}$ increased) or the discretization of the integration intervals
($N$ increased), then the number of the solved-for parameters $S_{nj}\
(n=0,...,N;\,j=1,...,N_{\omega})$ increases. In most of the settings
investigated in the present work, this number: $N_{\mathrm{para}}=(N+1)\times
N_{\omega}$, is actually larger than the number $N_{\mathrm{data}}$ of “data”,
that is the number of points in the spatio-temporal grid, 444 Indeed, at each
of those points one calculates the function on the l.h.s. of Eq. (18). These
are the data which have to be fitted by using the function on the r.h.s., that
depends on the solved-for parameters $S_{nj}$.  thus
$N_{\mathrm{data}}=N_{t}\times N_{\rho}\times N_{z}$. We have
$N_{\mathrm{data}}=1050$ for the grid used in this work, except for the two
grids used for the first two curves in Fig. 6, which give respectively
$N_{\mathrm{data}}=660$ and $N_{\mathrm{data}}=550$. (We have $N_{t}=5$ also
for those two grids.) Hence, we are often in a situation of “overfitting”. 555
However, the least-squares solver: the Matlab backslash operator mldivide,
puts to zero many of the $S_{nj}$ parameters for reason of insufficient
numerical accuracy, so that their effective number was $150$ to $200$, thus
significantly smaller than $N_{\mathrm{data}}$.  Higher values of the ratio
$R=N_{\mathrm{para}}/N_{\mathrm{data}}$ correspond with a more overfitted
situation, in which one expects the predictions to be less accurate. Now, on
each among Figs. 7–8 and 9–11, we saw that the domain of the low values of
$u_{j\,\mathrm{max}}$ increases when $N_{\mathrm{para}}$ increases, the number
$N_{\mathrm{data}}$ being fixed at $N_{\mathrm{data}}=1050$ for these five
figures — thus we saw that the domain of the low values of
$u_{j\,\mathrm{max}}$ increases when $R$ increases, i.e., when more
overfitting occurs. Therefore, it is the low values of $u_{j\,\mathrm{max}}$
which are likely to be a numerical artefact, not the extremely high values.
Note that this effect can be seen also on Fig. 6: $N$ and $N_{\omega}$, hence
$N_{\mathrm{para}}$, are the same for the three curves. So the third and
finest grid $(N_{\rho}=10)\times(N_{z}=21)$ leads to a lower $R$ than do the
other two grids $(N_{\rho}=12)\times(N_{z}=11)$ and
$(N_{\rho}=10)\times(N_{z}=11)$ — and accordingly the domain of the high
values of $u_{j\,\mathrm{max}}$ is larger for the third curve. On the other
hand, the effect of overfitting seems less important at a distance from the
axis (where the values of the SED are not high at all): e.g., if one changes
$N_{\omega}$ from $76$ to $23$, the SEDs considered in Figs. 2–5 stay similar:
they show a smaller number of oscillations due to the smaller number of
points, but the amplitude of the oscillations is comparable. Tentatively, we
may understand this by finding it natural that the very high numbers are more
sensitive to overfitting. In summary, the status of the very high values of
the maxima of $u_{j}(\rho,z)$: physical prediction of the model or numerical
artefact, is not really clear at the present stage of this work — but the
foregoing discussion, with different examples that all show an effect which
can be attributed to the “amount of overfitting”, seems to indicate that those
high values might be a true prediction of the model.

Figure 7: Predicted maximum values of the SED: $N_{\omega}=23$ vs
$N_{\omega}=46$.

Figure 8: Predicted maximum values of the SED: $N_{\omega}=46$ vs
$N_{\omega}=76$.

Figure 9: Predicted maximum values of the SED for $N_{\omega}=23$: $N=48$ vs
$N=24$.

Figure 10: Predicted maximum values of the SED for $N_{\omega}=46$: $N=48$ vs
$N=24$.

Figure 11: Predicted maximum values of the SED for $N_{\omega}=46$: $N=48$ vs
$N=192$.

##  5 Conclusion

In this work, the Maxwell model of the interstellar radiation field [12] has
been first checked for its predictions of the spatial variation of the
spectral energy distributions (SEDs) in our Galaxy. The model has been
adjusted by asking that the SED predicted at our local position in the Galaxy,
Eq. (32) with $\rho=\rho_{\mathrm{loc}}:=8$ kpc and $z=z_{\mathrm{loc}}:=0.02$
kpc, coincide with the observations collected through different spatial
missions [7, 8, 9]. Then the predictions of the model for other positions in
the Galaxy (for which, of course, we have no real observation) have been
compared with the predictions obtained by using a recent radiation transfer
model [6]. The two predictions do not differ too much in magnitude, even
though the predictions of the present model oscillate rather strongly as
function of the wavelength, especially at short wavelengths and not very far
(1 kpc) from the axis of the Galaxy. Also, the SED decreases more slowly with
increasing altitude for the present model. The values of the SED reach their
maximum on the axis itself of the Galaxy; that maximum is extremely high at
short wavelengths. These surprisingly high values occur in a larger or smaller
domain of wavelengths, depending on the settings of the calculation. However,
this dependence is governed by the “amount of overfitting”: less overfitting
increases the range of the high values. This makes it plausible that the high
values might be a true prediction of the model. The next work will aim at
improving the numerical procedure, in order: (i) to try to establish if the
high values are indeed a true prediction of the model; (ii) to check if a
numerically improved model could lead to oscillations of a smaller amplitude;
(iii) to quantify the slower decrease of the SED with increasing altitude for
the present model and to find the reason for it. Nevertheless, the present
results: (1) show that the prediction, by the Maxwell model of the ISRF, of
the spatial variation of the SED in the Galaxy, is quite comparable in
magnitude to that provided by a radiation transfer model (except perhaps for
high energies when the point is very close to the galaxy’s axis); and
therefore (2) indicate that an order-of-magnitude estimate of the “interaction
energy” predicted [13] by an alternative, scalar theory of gravitation, should
be indeed possible with the present model in a future work.

Acknowledgement. I am grateful to the reviewer for his (or her) careful study
of my manuscript, and for the important remarks and questions. In particular,
the remark on the likely overfitting, and the question about the reason for
the dependence of the maximum of energy density on the frequency mesh, taken
together, led me to improve significantly the discussion of the maximum values
of the SED. Also, the question about the reason why, on Fig. 3, this model’s
curve shifts upward from the radiation transfer model in the long wavelength
regime, led me to comment on the dependence of the SED with the altitude $z$.

##  Appendix A Appendix: Extension of the explicit representation to a finite
spectrum

In order to present this extension of the result [14] (the latter being
recalled in Subsect. 2.2), it is necessary to write this result in terms of
mappings. Consider the vector space $\mathcal{S}^{\omega}$ (respectively
$\mathcal{F}^{\omega}$) made of those solutions $\Psi$ of the scalar wave
equation (respectively, made of those source-free Maxwell fields
${{\boldsymbol{F}}}=({\bf E},{\bf B})$) that are axisymmetric and time-
harmonic with frequency $\omega$. 666 As the domain of definition of the
fields in $\mathcal{S}^{\omega}$ and $\mathcal{F}^{\omega}$, we take the whole
spacetime. More exactly, since we assume axial symmetry, we assume that they
are defined for the whole range of values of the spacetime variables
$t,\rho,z$, thus in the domain $\Omega$: ($t,z\in\mathbb{R}$,
$\rho\in\mathbb{R}_{+}$).  We define two mappings $Z_{1}$ and $Z_{2}$ from
$\mathcal{S}^{\omega}$ into $\mathcal{F}^{\omega}$: given $\Psi$ in
$\mathcal{S}^{\omega}$, $Z_{1}(\Psi)$ is the EM field given by Eqs. (2)–(4);
and $Z_{2}(\Psi)$ is the EM field given by Eqs. (2)–(4), followed by Eq. (5).
Obviously, $Z_{1}$ and $Z_{2}$ are linear. We can see that $Z_{1}(\Psi)$
belongs to $\mathcal{F}_{1}^{\omega}$, the subspace of $\mathcal{F}^{\omega}$
with $E_{\phi}=B_{\rho}=B_{z}=0$, whereas $Z_{2}(\Psi)$ belongs to
$\mathcal{F}_{2}^{\omega}$, the subspace of $\mathcal{F}^{\omega}$ with
$B_{\phi}=E_{\rho}=E_{z}=0$. The uniqueness of the decomposition (1) is just
that of the EM fields
$Z_{1}(\Psi)={{\boldsymbol{F}}}_{1}\in\mathcal{F}_{1}^{\omega}$ and
$Z_{2}(\Psi)={{\boldsymbol{F}}}_{2}\in\mathcal{F}_{2}^{\omega}$ (if they
exist), such that the starting EM field
${{\boldsymbol{F}}}\in\mathcal{F}^{\omega}$ is
${{\boldsymbol{F}}}={{\boldsymbol{F}}}_{1}+{{\boldsymbol{F}}}_{2}$. This
results immediately from the complementarity of the non-zero components, noted
in Subsect. 2.2, and inherent in the definition of $\mathcal{F}_{1}^{\omega}$
and $\mathcal{F}_{2}^{\omega}$: the non-zero components of
${{\boldsymbol{F}}}_{1}$ and ${{\boldsymbol{F}}}_{2}$ must be just the
corresponding components of ${{\boldsymbol{F}}}$, hence are known from the
data of ${{\boldsymbol{F}}}$. Henceforth, for any EM field
${{\boldsymbol{F}}}$, we will note ${{\boldsymbol{F}}}_{1}$ (respectively
${{\boldsymbol{F}}}_{2}$) the field whose components $E_{\phi},B_{\rho},B_{z}$
are zero (respectively the field whose components $B_{\phi},E_{\rho},E_{z}$
are zero), the other components being those of ${{\boldsymbol{F}}}$. In
contrast with the uniqueness, the existence of the decomposition (1) is a
rather strong result, which amounts to state that the mapping $Z$ from
$\mathcal{S}^{\omega}\times\mathcal{S}^{\omega}$ to $\mathcal{F}^{\omega}$
defined by

$Z(\Psi,\Psi^{\prime})=Z_{1}(\Psi)+Z_{2}(\Psi^{\prime})$ (33)

is surjective. This also amounts to state that, given any EM field
${{\boldsymbol{F}}}\in\mathcal{F}^{\omega}$, there exist $\Psi$ and
$\Psi^{\prime}$ in $\mathcal{S}^{\omega}$, such that
${{\boldsymbol{F}}}_{1}=Z_{1}(\Psi)$ and
${{\boldsymbol{F}}}_{2}=Z_{2}(\Psi^{\prime})$.

The general axisymmetric solutions of either the scalar wave equation or the
source-free Maxwell equations are got by summing time-harmonic axisymmetric
solutions of the respective equations. We restrict ourselves to such solutions
having a finite set of frequencies, say $(\omega_{j})\ (j=1,...,N_{\omega})$,
for simplicity [12, 14]. The axisymmetric solutions of the wave equation with
this set of frequencies form the vector space
$\mathcal{S}:=\mathcal{S}^{\omega_{1}}\oplus...\oplus\mathcal{S}^{\omega_{N_{\omega}}}$.
The axisymmetric solutions of the source-free Maxwell equations with this set
of frequencies form the vector space
$\mathcal{F}:=\mathcal{F}^{\omega_{1}}\oplus...\oplus\mathcal{F}^{\omega_{N_{\omega}}}$.
The sums are direct because the decomposition in a finite sum of time-harmonic
fields is of course unique. Therefore, the linear mappings $Z_{1}$ and $Z_{2}$
extend immediately to ones from $\mathcal{S}$ to
$\mathcal{F}_{1}:=\mathcal{F}_{1}^{\omega_{1}}\oplus...\oplus\mathcal{F}_{1}^{\omega_{N_{\omega}}}$
or respectively
$\mathcal{F}_{2}:=\mathcal{F}_{2}^{\omega_{1}}\oplus...\oplus\mathcal{F}_{2}^{\omega_{N_{\omega}}}$,
by linearity, and so does the theorem. That is: for any EM field
${{\boldsymbol{F}}}={{\boldsymbol{F}}}^{(1)}+...+{{\boldsymbol{F}}}^{(N_{\omega})}$
in the direct sum $\mathcal{F}$, the result valid for the time-harmonic case
ensures that there exist $2N_{\omega}$ scalar potentials:
$\Psi^{(1)},\Psi^{\prime(1)}\in\mathcal{S}^{\omega_{1}},...,\Psi^{(N_{\omega})},\Psi^{\prime(N_{\omega})}\in\mathcal{S}^{\omega_{N_{\omega}}}$,
such that

${{\boldsymbol{F}}}_{1}^{(1)}=Z_{1}(\Psi^{(1)}),{{\boldsymbol{F}}}_{2}^{(1)}=Z_{2}(\Psi^{\prime(1)}),...,{{\boldsymbol{F}}}_{1}^{(N_{\omega})}=Z_{1}(\Psi^{(N_{\omega})}),{{\boldsymbol{F}}}_{2}^{(N_{\omega})}=Z_{2}(\Psi^{\prime(N_{\omega})}).$
(34)

Hence, for any EM field ${{\boldsymbol{F}}}\in\mathcal{F}$, there are indeed
two scalar potentials $\Psi=\Psi^{(1)}+...+\Psi^{(N_{\omega})}$ and
$\Psi^{\prime}=\Psi^{\prime(1)}+...+\Psi^{\prime(N_{\omega})}$ in the direct
sum $\mathcal{S}$, such that

${{\boldsymbol{F}}}_{1}=Z_{1}(\Psi):=Z_{1}(\Psi^{(1)})+...+Z_{1}(\Psi^{(N_{\omega})}),\quad{{\boldsymbol{F}}}_{2}=Z_{2}(\Psi^{\prime}):=Z_{2}(\Psi^{\prime(1)})+...+Z_{2}(\Psi^{\prime(N_{\omega})}),$
(35)

and thus, finally,

${{\boldsymbol{F}}}={{\boldsymbol{F}}}_{1}+{{\boldsymbol{F}}}_{2}=Z_{1}(\Psi)+Z_{2}(\Psi^{\prime}):=Z(\Psi,\Psi^{\prime}).$
(36)

In the totally propagating case, we have $\Psi^{(j)}=\psi_{\omega_{j}\ S_{j}}$
and $\Psi^{\prime(j)}=\psi_{\omega_{j}\ S^{\prime}_{j}}$, for some
$2N_{\omega}$ functions $S_{j}(k),S^{\prime}_{j}(k),k\in[-K_{j},+K_{j}]\
(j=1,...,N_{\omega})$, with $K_{j}:=\frac{\omega_{j}}{c}$.

## References

  * [1] K. D. Gordon, K. A. Misselt, A. N. Witt, G. C. Clayton. The DIRTY model. I. Monte Carlo radiative transfer through dust. Astrophys. J., 551, 269–276 (2001).
  * [2] T. A. Porter, A. W. Strong. A new estimate of the galactic interstellar radiation field between $0.1\mu$m and $1000\mu$m. In: Proc. 29th International Cosmic Ray Conference, Pune. Mumbai: Tata Institute of Fundamental Research; 2005; vol. 4, pp. 77–80.
  * [3] C. C. Popescu, R. J. Tuffs, M. A. Dopita, J. Fischera, N. D. Kylafis, B. F. Madore. Modelling the spectral energy distribution of galaxies. V. The dust and PAH emission SEDs of disk galaxies. Astron. Astrophys., 527, A109, 40 pages (2011).
  * [4] T. P. Robitaille. HYPERION: an open-source parallelized three-dimensional dust continuum radiative transfer code. Astron. Astrophys. 536, A79, 17 pages (2011).
  * [5] W. J. Maciel. The interstellar radiation field. In: Astrophysics of the interstellar medium. W. J. Maciel, editor. New York: Springer; 2013; Chapter 2, pp. 17–31.
  * [6] C. C. Popescu, R. Yang, R. J. Tuffs, G. Natale, M. Rushton, F. Aharonian. A radiation transfer model for the Milky Way: I. Radiation fields and application to High Energy Astrophysics. Mon. Not. Roy. Astr. Soc., 470, no. 3, 2539–2558 (2017).
  * [7] R. C. Henry, R. C. Anderson, W. G. Fastie. Far-ultraviolet studies. vii. The spectrum and latitude dependence of the local interstellar radiation field. Astrophys. J., 239, 859–866 (1980).
  * [8] R. G. Arendt. The COBE Diffuse Infrared Background Experiment Search for the Cosmic Infrared Background. III. Separation of Galactic Emission from the Infrared Sky Brightness. Astrophys. J., 508, no. 1, 74–105 (1998).
  * [9] D. P. Finkbeiner, M. Davis, D. J. Schlegel. Extrapolation of Galactic Dust Emission at 100 Microns to Cosmic Microwave Background Radiation Frequencies Using FIRAS. Astrophys. J., 524, no. 2, 867–886 (1999).
  * [10] N. D. Kylafis, J. N. Bahcall. Dust distribution in spiral galaxies. Astrophys. J., 317, 637–645 (1987).
  * [11] A. N. Witt, H. A. Thronson, Jr., J. M. Capuano, Jr. Dust and the transfer of stellar radiation within galaxies. Astrophys. J., 393, 611-630 (1992).
  * [12] M. Arminjon. An analytical model for the Maxwell radiation field in an axially symmetric galaxy. Open Physics, 19, 77–90 (2021).
  * [13] M. Arminjon. On the equations of electrodynamics in a flat or curved spacetime and a possible interaction energy. Open Physics, 16, 488–498 (2018).
  * [14] M. Arminjon. An explicit representation for the axisymmetric solutions of the free Maxwell equations. Open Physics, 18, 255–263 (2020).
  * [15] R. L. Garay-Avendaño, M. Zamboni-Rached. Exact analytic solutions of Maxwell’s equations describing propagating nonparaxial electromagnetic beams. Appl. Opt. 53, 4524–4531 (2014).
  * [16] M. Zamboni-Rached, E. Recami, H. E. Hernández-Figueroa. Structure of nondiffracting waves and some interesting applications. In: H. E. Hernández-Figueroa, M. Zamboni-Rached, E. Recami (Eds.), Localized Waves. Hoboken: John Wiley & Sons; 2008; pp. 43–77.
  * [17] L. D. Landau, E. M. Lifshitz. The classical theory of fields, 3rd English edition. Oxford: Pergamon; 1971; pp. 118–119.
  * [18] D. J. Majaess, D. G. Turner, D. J. Lane, Characteristics of the Galaxy according to Cepheids. Mon. Not. Roy. Astron. Soc., 398, 263–270 (2009).


#  Calculation of the entropy for hard-sphere from integral equation method

Purevkhuu Munkhbaatar1 Banzragch Tsednee1 Tsogbayar Tsednee1,2 Tsookhuu
Khinayat1 1Institute of Physics and Technology, Mongolian Academy of Sciences,
Peace Ave 54B, Ulaanbaatar 13330, Mongolia, 2Chemistry Department, University
of North Dakota, Grand Forks, ND 58202, USA

###### Abstract

The Ornstein-Zernike integral equation method has been employed for a single-
component hard sphere fluid in terms of the Percus-Yevick (PY) and Martynov-
Sarkisov (MS) approximations. Virial equation of state has been computed in
both approximations. An excess chemical potential has been calculated with an
analytical expression based on correlation functions, and the entropy has been
computed with a thermodynamic relation. Calculations have been carried out for
a reduced densities of 0.1 to 0.9. It has been shown that the MS approximation
gives better values than those from the PY approximation, especially for high
densities and presents a reasonable comparison with available data in the
literature.

Ornstein-Zernike equation, closure relation, hard-sphere potential, excess
chemical potential, entropy

###### pacs:

Valid PACS appear here

##  I Introduction

In classical statistical physics the physical systems, such as, a liquid can
be described as a spherical symmetric hard-sphere particle model or Lennard-
Jones potentials Hansen06 ; McQuarried10 . Theoretical investigations for the
systems can be performed with various methods, such as an explicit
approaches–Monte-Carlo or Molecular dynamics simulations Hansen06 , and an
implicit approaches represented as an integral equation or a polarizable
continuum model methods Hirataed03 . An integral equation (IE) approach
mentioned here is a mathematical tool we use in our study. As an implicit
approach, the IE method does not consider a number of particles in the system
which in turn may make a calculation a cheap, and a solution of the IE can
give directly the correlation functions determining a general structure of the
system. A one-component, homogeneous system can be successfully investigated
with the Ornstein-Zernike (OZ) OrnsteinZ14  IE theory combined with an
appropriate auxiliary equation.

In this work our purpose is to obtain an excess entropy for single-component
hard-sphere fluid using the OZ IE approach combined with the Percus-Yevick
PYpr68  and Martynov-Sarkisov MSmp68  approximations. To reach it, we will
first compute a virial equation of state and along with it, we will obtain an
excess chemical potential for the system at equilibrium using an analytical
expression based on the correlation function. Using these two quantities, we
will compute the excess entropy employing a thermodynamic relation. We will
compare our findings for thermodynamic properties with available data in
literature Zhou06 ; Carnahan69 . Note that, to our knowledge, the MS
approximation has not been tested for this calculation, yet, even though this
problem had been considered in the past Zhou06 ; Carnahan69 ; Laird92 .
Therefore, we believe that our findings in this work can be considered as some
contributions to this area.

In Section II we will discuss about the Ornstein-Zernike theory and
thermodynamic properties which we will compute. In Section III we present
numerical results and their discussions. In Section IV the conclusion is
given.

##  II Theory

###  II.1 The Ornstein-Zernike equation

In statistical mechanics an structural properties for the liquids existing at
equilibrium can be obtained in terms of the integral equation formalism
Hirataed03 . For a single-component homogeneous system the Ornstein-Zernike
integral equation can be written in the form

$\displaystyle h(r)=c(r)+\rho\int
c(|\mathbf{r}-\mathbf{r}^{\prime}|)h(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}$
(1)

where $h(r)$ and $c(r)$ are the total and direct correlation functions,
respectively, and $\rho$ is density of the system.

In equation (1), the two correlation functions are unknown, therefore, it
cannot be solved directly. To solve this OZ equation (1), we need another
equation which is solved together with the OZ equation (1) in the self-
consistent way. This required equation is called a closure equation (relation)
which may be written in the form

$\displaystyle h(r)=\exp[-\beta u(r)+\gamma(r)-B(r)]-1,$ (2)

where $u(r)$ is a pair potential for particles in the system;
$\gamma(r)=h(r)-c(r)$ is an indirect correlation function; $B(r)$ is a bridge
function; $\beta=1/k_{\mathrm{B}}T$ where $k_{\mathrm{B}}$ is the Boltzmann’s
constant and $T$ is a temperature of the system. The radial distribution
function $g$ can be defined as $g(r)=h(r)+1$ as well.

An inter-particle interaction potential for the hard-sphere with a diameter
$\sigma$ can be written in the form McQuarried10

$\displaystyle u(r)=\begin{cases}0,&r\leq\sigma\\\
\infty,&r>\sigma.\end{cases}$ (3)

The bridge functions we use in this work can be given in the forms:

$\displaystyle B(r)=\begin{cases}\ln(1+\gamma)-\gamma,&\mbox{Percus-Yevick
(PY)}\,\,[5]\\\ \sqrt{1+2\gamma}-\gamma-1,&\mbox{Martynov-Sarkisov
(MS)}\,\,[6].\end{cases}$

###  II.2 Thermodynamic quantities

Once we solve the OZ equation (1), we can have the correlation functions with
which we can compute thermodynamic properties for the system.

####  II.2.1 Virial equation of state

For the hard-sphere system, the virial equation of state can be written in the
form McQuarried10

$\displaystyle\frac{\beta p}{\rho}=1+\frac{2\pi}{3}\rho\sigma^{3}g(\sigma),$
(4)

where $g(\sigma)$ is the contact value of the radial distribution function at
$\sigma$.

####  II.2.2 An excess chemical potential

The excess (e) chemical potential $\beta\mu^{e}$ can be computed with a
following approximated analytical expression

$\displaystyle\beta\mu^{e}$ $\displaystyle\approx$
$\displaystyle\rho\int\Big{(}\frac{1}{2}h(r)^{2}-c(r)-\frac{1}{2}h(r)c(r)\Big{)}d\mathbf{r}$
$\displaystyle+\rho\int\Big{(}B(r)+\frac{2}{3}h(r)B(r)\Big{)}d\mathbf{r}.$

Note that this expression does not depend an explicit form of the bridge
function. Therefore, we can use it for both approximation in this work. A
derivation of this expression can be found in Ref. Tsog_pre19 .

####  II.2.3 An excess entropy

For the hard-sphere system in which an internal energy is zero, the excess
entropy $S^{e}/Nk_{B}$ can be computed as a following thermodynamic relation
McQuarried10

$\displaystyle\frac{S^{e}}{Nk_{\mathrm{B}}}=\frac{\beta
p}{\rho}-\beta\mu^{e}-1.$ (6)

In evaluating expression (6), we use previously obtained values for the
pressure and excess chemical potential.

##  III Results and discussion

In this work we have done calculations for a reduced densities:
$\rho\sigma^{3}=0.1$ to $0.9$. For high density, such as, when
$\rho\sigma^{3}\sim 0.9$, the hard-sphere system behaves like a fluid Hansen06
. We use the Picard iteration method to solve the OZ equation, in which the OZ
equation (1) is solved in the Fourier space while the closure equation (2) is
solved in a coordinate space. A number of grid points in a length interval of
$[0,\,16\sigma]$ is $2^{15}$. The numerical calculation has been done with an
in-house Matlab Matlab15  code. The numerical tolerance for a root-mean-
squared residual of the indirect correlation functions of a successive
iterations was set at $10^{-8}$.

In table I we have compared the values of the equation of states (EOS) and
excess chemical potential obtained from both the PY and MS approximations with
available data in Ref. Zhou06 , which was obtained in terms of the Carnahan
and Starling (CS) EOS calculation Carnahan69 . We note that the CS data
Carnahan69  can be considered as ‘exact’ values Zhou06 .

Table 1: The $\beta p/\rho$ and $\beta\mu^{e}$ as a function of $\rho\sigma^{3}$.  | $\beta p/\rho$ | $\beta\mu^{e}$
---|---|---
$\rho\sigma^{3}$ | PY | MS | Ref. Carnahan69  |  | PY | MS | Ref. Carnahan69
0.1 | $1.24$ | $1.24$ | $1.2397$ |  | $0.46$ | $0.46$ | $0.4637$
0.2 | $1.55$ | $1.55$ |  |  | $1.03$ | $1.03$ |
0.3 | $1.95$ | $1.96$ | $1.9667$ |  | $1.72$ | $1.73$ | $1.7468$
0.4 | $2.48$ | $2.50$ |  |  | $2.58$ | $2.62$ |
0.5 | $3.17$ | $3.23$ | $3.2624$ |  | $3.67$ | $3.77$ | $3.8068$
0.6 | $4.09$ | $4.23$ |  |  | $5.09$ | $5.27$ |
0.7 | $5.32$ | $5.61$ | $5.7102$ |  | $7.00$ | $7.32$ | $7.3593$
0.8 | $7.00$ | $7.58$ | $7.7497$ |  | $9.70$ | $10.2$ | $10.1525$
0.9 | $9.33$ | $10.4$ | $10.7461$ |  | $13.7$ | $14.3$ | $14.1052$

For higher densities, values for both quantities from PY approximation are
lower than those from the MS approximation. The MS values are quite close to
those obtained in Ref. Carnahan69 .

Table II presents our findings for an excess entropy $S^{e}/Nk_{\mathrm{B}}$
obtained in both approximations and their comparison with those of Refs.
Zhou06  and Carnahan69 . As expected, the largest deviations come from the PY
results here, especially for high density, since both the PY pressure and
excess chemical potential values are lower (previous table). Values from Ref.
Zhou06  are obtained with the help of hybrid-bridge function, which enables
better values than ours.

Table 2: An excess entropy $S^{e}/Nk_{\mathrm{B}}$. | $S^{e}/Nk_{\mathrm{B}}$ |
---|---|---
$\rho\sigma^{3}$ | PY | MS | Ref. Zhou06  | Ref. Carnahan69
0.1 | $-0.22$ | $-0.22$ | $-0.2240$ | $-0.2241$
0.2 | $-0.48$ | $-0.48$ |  |
0.3 | $-0.77$ | $-0.77$ | $-0.7790$ | $-0.7801$
0.4 | $-1.10$ | $-1.12$ |  |
0.5 | $-1.50$ | $-1.53$ | $-1.5414$ | $-1.5444$
0.6 | $-2.00$ | $-2.04$ |  |
0.7 | $-2.68$ | $-2.71$ | $-2.6851$ | $-2.6491$
0.8 | $-3.70$ | $-3.60$ | $-3.3824$ | $-3.4028$
0.9 | $-5.39$ | $-4.87$ | $-4.3555$ | $-4.3591$

Based on results shown in tables I and II we note that the MS approximation
makes better calculation than the PY does. A reason why this happens might be
related to the fact that an correlation function from the MS approximation are
better than those from the PY approximation Banz22 ; Banz22_num . Finally,
note that our results obtained are independent on number of grid points and a
length interval employed.

##  IV Conclusion

In this work we have implemented the Orsntein-Zernike integral equation theory
with the Percus-Yevick and Martynov-Sarkisov bridge functions for single-
component hard-sphere system. Analytical expressions based on correlation
functions which have been obtained as a solution of the integral equation for
a density, and an excess chemical potential and entropy have been tested. Our
findings for these thermodynamic quantities from both approximations have been
compared with available accurate data. It has been shown that Matrynov-
Sarkisov approximation does work better than the Percus-Yevick approximation
in which the integral equation can have a closed form solution Laird92 . The
MS values are close to the accurate data. Note that it is better to use a
well-approximated bridge functions, such as hybrid or pressure consistent
bridge approximations Zhou06  to obtain well-consistent data from the integral
equation method.

##  V ACKNOWLEDGEMENT(S)

This research work has been supported by the Mongolian Foundation for Science
and Technology (Project No. ShUTBIKhKhZG-2022/167).

## References

  * (1) J. P. Hansen and I. R. McDonald, Theory of Simple Liquid, Academic Press, New York (2006).
  * (2) D. A. McQuarrie, Statistical Mechanics, University Press Books (2010).
  * (3) F. Hirata, Molecular Theory of Solvation, Kluwer Academic Publishers Dordretch (2003).
  * (4) L. S. Ornstein and F. Zernike, Proc. Natl. Acad. Sci. USA 7, 793 (1914).
  * (5) J. K. Percus and G. J. Yevick, Phys. Rev. 110, 1 (1968).
  * (6) G. A. Martynov and G. N. Sarkisov, Mol. Phys. 49, 1495 (1983).
  * (7) S. Zhou, Chem. Phys. 330, 478 (2006).
  * (8) N. F. Carnahan and K. E. Starling, Chem. Phys. 51, 635 (1969).
  * (9) B. B. Laird and A. D. J. Haymet, Phys. Rev. A 45, 5680 (1992).
  * (10) Ts. Tsednee and T. Luchko, Phys. Rev. E, 99, 032130 (2019).
  * (11) MATLAB R2015a, The MathWorks Inc., Natick, Massachusetts (2015).
  * (12) B. Tsednee, Ts. Tsednee and Ts. Khinayat, Defect and Diffusion Forum 423, 17, (2023).
  * (13) B. Tsednee, Ts. Tsednee and Ts. Khinayat, Sci. trans.$-$Phys., Natl. Univ. Mongolia, 33, 30 (2022).


# On the time fractional heat equation with obstacle

C. Alberini  Dipartimento SBAI, Sapienza Università di Roma,
<EMAIL_ADDRESS>R. Capitanelli  Dipartimento SBAI, Sapienza
Università di Roma<EMAIL_ADDRESS>M. D’Ovidio  Dipartimento
SBAI, Sapienza Università di Roma<EMAIL_ADDRESS>S. Finzi Vita
Dipartimento di Matematica, Sapienza Università di Roma,
<EMAIL_ADDRESS>
Abstract. We study a Caputo time fractional degenerate diffusion equation
which we prove to be equivalent to the fractional parabolic obstacle problem,
showing that its solution evolves for any $\alpha\in(0,1)$ to the same
stationary state, the solution of the classic elliptic obstacle problem. The
only thing which changes with $\alpha$ is the convergence speed.
We also study the problem from the numerical point of view, comparing some
finite different approaches, and showing the results of some tests. These
results extend what recently proved in [1] for the case $\alpha=1$.

##  1 Introduction

In the works of the last decades, the use of fractional calculus has known an
ever increasing use in describing phenomena from the most disparate scientific
fields, from biology to mechanics, to superslow diffusion in porous media till
financial interests.

Aim of this paper is to study the following problem, which generalizes the one
addressed in [1], for $0<\alpha<1$:

$\displaystyle\left\\{\begin{array}[]{ll}\partial^{\alpha}_{t}u-H\left(u-\psi\right)\Delta
u=0&\text{ a.e. in }\Omega,\text{ for all }t\in(0,T),\\\
u\left(0\right)=u^{0}&\text{ in }\Omega,\\\ u=0&\text{ on
}\partial\Omega,\text{ for all }t\in(0,T),\end{array}\right.$ (4)

with $T>0$, where $H$ is the extended Heaviside function such that $H(0)=0$,
that is

$H(r)=\left\\{\begin{array}[]{ll}1&\text{ for }r>0\cr 0&\text{ for }r\leq
0\end{array}\right.$ (5)

and $\Omega$ is a bounded domain in $\mathbb{R}^{n}$, $n\in\mathbb{N}$, with
smooth boundary. For simplicity, we omitted the presence of a forcing term in
the problem.

Here $\partial^{\alpha}_{t}u$ denotes the Caputo fractional derivative, that
is,

$\partial^{\alpha}_{t}u(x,t):=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-s)^{-\alpha}\partial_{s}u(x,s)\
ds,$ (6)

where $\Gamma$ is the Gamma function.

The Riemann-Liouville derivative is defined as

$\displaystyle^{R}\partial^{\alpha}_{t}u(x,t):=\frac{\partial}{\partial
t}\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-s)^{-\alpha}\,u(x,s)\ ds.$ (7)

We observe that

$\,{}^{R}\partial^{\alpha}_{t}C=\frac{1}{\Gamma(1-\alpha)}\partial_{t}\int_{0}^{t}\frac{C}{(t-s)^{\alpha}}ds=C\frac{t^{-\alpha}}{\Gamma(1-\alpha)}$

and the following relation between derivatives holds

$\displaystyle\partial^{\alpha}_{t}u(x,t)=\,^{R}\partial^{\alpha}_{t}\big{(}u(x,t)-u(x,0)\big{)}=\,^{R}\partial^{\alpha}_{t}u(x,t)-u(x,0)\frac{t^{-\alpha}}{\Gamma(1-\alpha)}.$
(8)

We refer to the book [5] for further details on fractional derivatives.

Let us consider the complementarity system, for all $t>0$,

$\displaystyle\left\\{\begin{array}[]{ll}w(x,t)\geq\psi(x)&\ \text{in
}\Omega\\\ \partial^{\alpha}_{t}w\geq\Delta w&\ \text{a.e. in }\Omega\\\
\displaystyle(w-\psi)(\partial^{\alpha}_{t}w-\Delta w)=0&\ \text{a.e. in
}\Omega\\\ w\left(x,0\right)=u^{0}&\ \text{in }\Omega\\\ w(x,t)=0&\ \text{on
}\partial\Omega.\end{array}\right.$ (14)

In Section 2 we will prove that

-

by analogy with the classical case, (4) is equivalent to the complementarity
system (14);

-

asymptotically, the solution evolves for each $\alpha$ towards the same steady
state of the classical problem, that is

$\bar{u}\geq\psi,\quad-\Delta\bar{u}\geq
0,\quad(\bar{u}-\psi)\Delta\bar{u}=0\quad\hbox{in }\Omega\times(0,T)$ (15)

with different convergence speed (for $\alpha=1$ exponential, for
$\alpha\in(0,1)$ polynomial).

In Section 3 we present three possible finite difference schemes for the
numerical approximation of problem (4) or its equivalent form given by the
complementary system (14). The Caputo derivative has been discretized using
standard methods found in the literature, the so-called L1 or Convolution
Quadrature (CQ) approaches (see e. g. [9] or [7]). The space discretization
has been carried out through the semi implicit finite differences scheme
introduced in [1] for problem (4), or through the implicit scheme of [3] for
the evolutive obstacle problem (14). Note that for $\alpha\to 1^{-}$ all these
schemes give back the known results of the classic heat equation with
obstacle.

Our aim was not to find an optimal strategy of approximation for the problem,
but only to derive working schemes in orders to confirm through explicit
simulations the behavior of the solution as characterized by the results of
Section 2. For that, in Section 4, we have tested the previous schemes on a
couple of one-dimensional examples. Such simulations also allow to compare
their reliability and computational cost. The semi implicit approach requires
time step restrictions (strongly increasing as $\alpha\to 0^{+}$) in order to
detect the correct contact set with the obstacle, restrictions unnecessary for
the implicit approach. On the other hand, in the first case each time
iteration is much less expensive, since it reduces to a single linear system
solution. So, after all, all the proposed schemes appear to be competitive.

For the following, we assume that the initial datum $u^{0}$ and the
(independent of time) obstacle $\psi$ satisfy the following conditions

$u^{0}\in H^{1}_{0}\left(\Omega\right),\,\psi\in H^{2}\left(\Omega\right),\
\psi\leq 0\text{ on }\partial\Omega.$ (16)

We define as solution of problem (4) a function $u\in
L^{2}\left(0,T;H^{1}_{0}\left(\Omega\right)\cap
H^{2}\left(\Omega\right)\right)$ with $\partial_{t}^{\alpha}u(x,t)\in
L^{2}\left(0,T;L^{2}(\Omega)\right)$ which solves problem (4).

##  2 Equivalence of the problems and asymptotic behaviour

In the present section we analyze problem (4), showing that, under suitable
conditions, it is equivalent to complementarity system (14). Let us introduce
the following hypoteses:

H1:

$u^{0}>\psi$ a.e. in $\Omega$;

H2:

$\Delta\psi\leq 0$ a.e. in $\Omega$.

######  Proposition 2.1.

Assume that conditions (16), H1 and H2 hold. Then problems (4) and (14) are
equivalent.

Proof.

First we note that if $w$ solves problem (14), then $w$ solves (4). In fact,
initial and boundary conditions are the same in (4) and (14). Since
$(w-\psi)(\partial^{\alpha}_{t}w-\Delta w)=0$ a.e. in $\Omega$, where
$w>\psi$, then we obtain $\partial^{\alpha}_{t}w-\Delta w=0$ a.e $\text{in
}\Omega$ whereas, where $w(x,t)=\psi(x)$ then $\partial_{t}w=0$ which implies
$\partial_{t}^{\alpha}w=0$.

Now we prove that if $u$ solves problem (4) then $u$ coincides with solution
$w$ of (14). Initial and boundary conditions are the same in (4) and (14). Let
us now prove that if $u(x,t)$ is a solution of (4) then necessarily
$u(.,t)\geq\psi(.)$ in $\Omega$ for any time $t$. Assume that $u<\psi$. Then
$\partial^{\alpha}_{t}u=0$ from (4). From ${\bf H}_{1}$, $u^{0}>\psi$.
Moreover, due to the fact that $u<\psi<u^{0}$, we obtain from (8)

$\displaystyle\partial^{\alpha}_{t}u=\,^{R}\partial^{\alpha}_{t}(u-u^{0})<0,$

which does not agree with $\partial^{\alpha}_{t}u=0$. Then $u\geq\psi$ by
contradiction and the first inequality of (14) holds. The equation in the
third line of (14) is trivially satisfied where $u(x,t)=\psi(x)$; where
$u(x,t)>\psi(x),$ from (4) we obtain $\partial^{\alpha}_{t}-\Delta u=0$, so it
is always true. Concerning the second inequality of (14), we have already seen
that it is satisfied (with the equal sign) when $u>\psi$. But when
$u(x,t)=\psi(x)$, (4) and assumption H2 imply that

$\partial^{\alpha}_{t}u-\Delta u=-\Delta\psi\geq 0.$

Then $u$ solves problem (14).∎

We recall that the same problem for $\alpha=1$ has been faced in [1], proving
in particular the equivalence of problem (4) with a parabolic obstacle
problem. Let $v$ be the unique solution of problem (14) with $\alpha=1$, that
is

$\displaystyle\left\\{\begin{array}[]{ll}v(x,t)\geq\psi(x)&\ \text{in
}\Omega\\\ \partial_{t}v\geq\Delta v&\ \text{a.e. in }\Omega\\\
\displaystyle(v-\psi)(\partial_{t}v-\Delta v)=0&\ \text{a.e. in }\Omega\\\
v\left(x,0\right)=u^{0}&\ \text{in }\Omega\\\ v(x,t)=0.&\ \text{on
}\partial\Omega\end{array}\right.$ (22)

Let us introduce

$u(x,t)=\int_{0}^{\infty}v(x,s)\,\ell(s,t)\,ds$ (23)

where the probability density $\ell:[0,\infty)\times[0,\infty)\to[0,\infty)$
satisfies (see [6])

$\begin{cases}\displaystyle\partial^{\alpha}_{t}\ell=-\frac{\partial\ell}{\partial
s},\quad(s,t)\in(0,\infty)\times(0,\infty)\\\
\displaystyle\ell(s,0)=\delta(s),\quad s>0,\\\ \displaystyle\ell(s,t)=0,\quad
s<0.\end{cases}\quad\textrm{and}\quad\begin{cases}{}^{R}\partial^{\alpha}_{t}\ell=-\frac{\partial\ell}{\partial
s},\quad(s,t)\in(0,\infty)\times(0,\infty),\\\
\displaystyle\ell(s,0)=\delta(s),\quad s>0,\\\
\displaystyle\ell(0,t)=\frac{t^{-\alpha}}{\Gamma(1-\alpha)},\quad t>0,\\\
\displaystyle\ell(s,t)=0,\quad s<0.\end{cases}$

######  Proposition 2.2.

The function $u$ defined in (23) solves (14).

Proof.

As $v\geq\psi$, then

$\displaystyle\int_{0}^{\infty}v(x,s)\,\ell(s,t)\,ds\geq\int_{0}^{\infty}\psi(x)\,\ell(s,t)\,ds=\psi(x),$

and we obtain that $u\geq\psi.$

Now we prove that $(u-\psi)(\partial^{\alpha}_{t}u-\Delta u)=0.$ For $u>\psi$,
it holds the classical theory about fractional Cauchy problems (see [2], [10]
and [4, Theorem 5.2]) and therefore we have that, in $L^{2}(\Omega)$,

$\displaystyle\partial_{t}^{\alpha}u(x,t)=\Delta u(x,t),\quad
u(x,0)=u^{0}(x).$

Concerning the second inequality of (14), we observe that if
$\partial_{t}v-\Delta v\geq 0$, then

$\displaystyle\Delta\int_{0}^{\infty}v(x,s)\,\ell(s,t)\,ds\leq$
$\displaystyle\int_{0}^{\infty}\partial_{s}v(x,s)\,\ell(s,t)\,ds$
$\displaystyle=$
$\displaystyle-v(x,0)\frac{t^{-\alpha}}{\Gamma(1-\alpha)}-\int_{0}^{\infty}v(x,s)\,\partial_{s}\ell(s,t)\,ds$
$\displaystyle=$
$\displaystyle-v(x,0)\frac{t^{-\alpha}}{\Gamma(1-\alpha)}+\,^{R}\partial^{\alpha}_{t}u(x,t).$

Since $u(x,0)=v(x,0)$ from the relation (8) we obtain
$\partial^{\alpha}_{t}u-\Delta u\geq 0$.

The initial condition and the boundary condition for the function $u$ follow
trivially from the initial condition and the boundary condition for the
function $v$. ∎

In [1] the asymptotic solution of (22) has been characterized as the solution
of the corresponding stationary (elliptic) obstacle problem. In particular, it
has been proved that under conditions (16), H1 and H${}_{2},$ the solution
$v(t)$ to the parabolic obstacle problem (22) converges strongly in
$H^{1}_{0}(\Omega),$ for $t\to\infty,$ to the unique solution $\overline{u}$
of the corresponding stationary obstacle problem

$\overline{u}\in
H^{1}_{0}(\Omega),\quad\overline{u}\geq\psi,\quad-\Delta\overline{u}\geq
0,\quad(\overline{u}-\psi)(\Delta\overline{u})=0\ \hbox{ a.e. in }\Omega.$
(24)

Moreover, it has been proved that there is a constant $C>0$ such that, for
every $t\geq 1$,

$\left\|v(t)-\overline{u}\right\|_{H^{1}\left(\Omega\right)}\leq e^{-Ct}\ .$
(25)

In the next Theorem we prove that a similar result holds for any
$\alpha\in(0,1).$

######  Theorem 2.1.

Assume conditions (16), H1 and H2 hold. Let $u(t)$ be the function defined in
(23) solution of (14). Let $\overline{u}$ be the unique solution of the
obstacle problem (24). Then, $u(t)$ converges to $\overline{u}$ for
$t\to\infty,$ and there is a constant $C>0$ such that, for every $t\geq 1$,

$\left\|u(t)-\overline{u}\right\|_{L^{1}\left(\Omega\right)}\leq\kappa_{\Omega}\,E_{\alpha}(-Ct^{\alpha})\
$ (26)

for any $\alpha\in(0,1),$ where
$\displaystyle\kappa_{\Omega}=\left(\int_{\Omega}dx\right)^{1/2}$ and

$\displaystyle\frac{E_{\alpha}(-Ct^{\alpha})}{J(t)}\to 1,\quad\textrm{as}\quad
t\to\infty,$ (27)

with

$\displaystyle J(t)=\frac{1}{C}\frac{t^{-\alpha}}{\Gamma(1-\alpha)}.$ (28)

Proof. Since

$\displaystyle\|v-\bar{u}\|_{L^{1}(\Omega)}\leq\kappa_{\Omega}\;\|v-\bar{u}\|_{L^{2}(\Omega)}$

we have

$\displaystyle\|u(t)-\bar{u}\|_{L^{1}(\Omega)}=$
$\displaystyle\int_{\Omega}\bigg{|}\int_{0}^{\infty}v(x,s)\,\ell(s,t)\,ds-\bar{u}(x)\bigg{|}dx$
$\displaystyle=$
$\displaystyle\int_{\Omega}\bigg{|}\int_{0}^{\infty}\left(v(x,s)-\bar{u}(x)\right)\ell(s,t)\,ds\bigg{|}dx$
$\displaystyle\leq$
$\displaystyle\int_{\Omega}\int_{0}^{\infty}\bigg{|}v(x,s)-\bar{u}(x)\bigg{|}\ell(s,t)\,dsdx$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\|v(x,s)-\bar{u}(x)\|_{L^{1}(\Omega)}\,\ell(s,t)\,ds$
$\displaystyle\leq$
$\displaystyle\kappa_{\Omega}\int_{0}^{\infty}e^{-Cs}\,\ell(s,t)\,ds$
$\displaystyle=$
$\displaystyle\kappa_{\Omega}\,E_{\alpha}(-Ct^{\alpha}),\quad\forall\,\alpha\in(0,1)$

where, in the last step, we used the fact that $E_{\alpha}$ is the Laplace
transform of the density $\ell$. For the asymptotic behaviour of the Mittag-
Leffler, consult the book [8, formula (4.4.17)].

∎

##  3 Numerical approximation

The starting idea for a numerical approximation of the problem has been to
combine classical time-stepping discretization schemes for the Caputo
derivative, such as the convolution quadrature (CQ) or the finite difference
(L1) scheme, see e.g. [9], with schemes usually working for the parabolic
obstacle problem (14), such as the one proposed in [3], but also the semi-
implicit f.d. scheme tested in [1] for the equivalent Heaviside function
formulation (4) of the problem.

For sake of simplicity we work in the one-dimensional case, with
$\Omega=(a,b)$. Let us call $h>0$ the space discretization step ($h=(b-a)/N$,
so that we have $(N-1)$ internal nodes in $\Omega$, $x_{i}=a+ih$, for
$i=1,...,N-1$) and $\tau=T/M$ the time discretization step (with $M$ time
instants $t^{m}=m\tau$, for $m=1,...,M$); with
$\gamma_{\alpha}=\tau^{\alpha}/h^{2}$ we denote the parabolic ratio between
the steps related to a specific $\alpha$. Note that for fixed steps $h$ and
$\tau$, if $\alpha$ decreases to zero then $\gamma_{\alpha}$ quickly grows: in
other words for small $\alpha$, in order to keep $\gamma_{\alpha}$ small on a
fixed mesh, the step $\tau$ has to be considerably reduced, with a significant
increase of computational costs.

Here we analyze three possible approaches:

  1. 1.

Scheme S1: solves problem (4) with L1 for the Caputo derivative and the semi
implicit f.d. scheme of [1] in space.

The time discretization of the Caputo derivative by the L1 scheme leeds to the
formula (see [9]):

$\partial^{\alpha}_{t}u(x,t^{m})\simeq\frac{1}{\Gamma(2-\alpha)\tau^{\alpha}}\left\\{u(x,t^{m})-\sum_{k=0}^{m-1}C_{m,k}u(x,t^{k})\right\\}\
,$

where

$C_{m,0}:=f(m),\quad C_{m,k}:=f(m-k)-f(m-(k-1))\ \ \hbox{ for }k=1,...,m-1,$

and

$f(r):=r^{1-\alpha}-(r-1)^{1-\alpha},\hbox{ for }r\geq 1.$

Then, after semidiscretization in time, we need to solve for any instant
$t^{m}$ the equation

$\frac{1}{\Gamma(2-\alpha)\tau^{\alpha}}\left\\{u(x,t^{m})-\sum_{k=0}^{m-1}C_{m,k}u(x,t^{k})\right\\}-H(u-\psi)u_{xx}=0\quad\hbox{
in }\Omega,$ (29)

with the same boundary conditions of (4). Since $\Omega$ was splitted in $N$
subintervals through the nodes $x=\\{x_{i}\\}_{i}$, the initial data will be
the vector $u^{0}\in\mathbb{R}^{N-1}$, with $u^{0}_{i}=u^{0}(x_{i})$; applying
the semi implicit finite difference scheme of [1] in space, the solution
$u^{1}$ at the first discrete time $t^{1}=\tau$ will solve at any node the
relation:

$\frac{1}{g\tau^{\alpha}}\left(u^{1}_{i}-C_{1,0}u^{0}_{i}\right)=H(u^{0}_{i}-\psi_{i})\delta^{2}u^{1}_{i}:=H(u^{0}_{i}-\psi_{i})\frac{u^{1}_{i-1}-2u^{1}_{i}+u^{1}_{i+1}}{h^{2}},$

with $g=\Gamma(2-\alpha)$ (note that $0.8862\leq g\leq
1,\forall\alpha\in[0,1]$, so that this term will be negligible with respect to
$\gamma_{\alpha}$). If we set $v^{k}_{i}=u^{k}_{i}-\psi_{i}$, redistributing
all the terms between the two members, it is equivalent to solve

$u^{1}_{i}-g\tau^{\alpha}H(v^{0}_{i})\delta^{2}u^{1}_{i}=C_{1,0}u^{0}_{i}\quad\hbox{
for every }i;$

with vector notations it means that $u^{1}$ solves the linear system:

$B^{0}u^{1}:=(I+g\gamma_{\alpha}H(v^{0})*A)u^{1}=C_{1,0}u^{0},$

where $A$ is the usual tridiagonal matrix $(N-1)\times(N-1)$ with values $2$
on the main diagonal and $-1$ on the two adjacent diagonals, and we denoted

$\\{(H(v)*A)u\\}_{i}:=H(v_{i})(Au)_{i}=H(v_{i})\sum_{j=1}^{N-1}A_{i,j}u_{j}\
.$

Since the discrete solution could overstep the obstacle at some nodes, in
particular when a large value of $\gamma_{\alpha}$ is used, the following
correction is needed at any iteration:

$u^{1}_{i}=\max(u^{1}_{i},\psi(x_{i}))\ .$

In the same way we see that $u^{2}$ solves for any $i$

$\frac{1}{g\tau^{\alpha}}\left(u^{2}_{i}-C_{2,0}u^{0}_{i}-C_{2,1}u^{1}_{i}\right)=H(u^{1}_{i}-\psi_{i})\delta^{2}u^{2}_{i},$

that is the system

$B^{1}u^{2}=(I+g\gamma_{\alpha}H(v^{1})*A)u^{2}=C_{2,0}u^{0}+C_{2,1}u^{1},$

with the same matrix $A$ and the subsequent correction. In general, at any
time step $u^{m}$ solves the linear system

$B^{m-1}u^{m}=b^{m},$ (30)

where we have set

$B^{m-1}=I+g\gamma_{\alpha}H(v^{m-1})*A,\qquad
b^{m}=\sum_{k=0}^{m-1}C_{m,k}u^{k}\ ,$ (31)

followed by the correction

$u^{m}_{i}=\max(u^{m}_{i},\psi(x_{i})).$ (32)

Note that all the matrices $B^{m}$ are symmetric positive definite (and
M-matrices), since so it is $A$, while $H(v)\geq 0$. Then all the previous
linear systems are well posed.

With respect to the classic parabolic obstacle problem approach discussed in
[1] there is here an important difference. In that case (which corresponds to
the case $\alpha=1$), when the solution at time $t^{m-1}$ touches the obstacle
at node $x_{i}$, then $v^{m-1}_{i}=0$, $H(v^{m-1}_{i})=0$, and (30) trivially
yields:

$u^{m}_{i}=u^{m-1}_{i}\ ;$

in other words, once touched the obstacle at a particular node the solution
does not change anymore there, but only at the remaining free nodes. In the
general case of $\alpha\in(0,1)$ on the contrary, at the contact time system
(30) immediately yields for the $i$-th component:

$u^{m}_{i}=\sum_{k=0}^{m-1}C_{m,k}u^{k}_{i}=C_{m,0}u^{0}_{i}+...+C_{m,m-1}\psi_{i}>\psi_{i}\sum_{k=0}^{m-1}C_{m,k}=\psi_{i},$

since at least $u^{0}_{i}>\psi_{i}$ and $\sum_{k=0}^{m-1}C_{m,k}=1$. It
follows that the solution has a little rebound at $x_{i}$ which detaches it
again from the obstacle, and produces an (innatural) oscillating evolution
from that time on. The width of such rebound depends on the size of
$\gamma_{\alpha}$. Then, even if in principle the semi implicit scheme does
not not require stability restrictions on the discretization steps, a small
value of $\gamma_{\alpha}$ will be necessary to reduce the oscillations (they
would vanish for $\tau\to 0$). The way to solve this difficulty is to remove
the memory effect at the contact nodes of the solution with the obstacle, that
is where the obstacle retains the solution. This suggests to modify the scheme
replacing the vector $b^{m}$ of (31) by the vector $\hat{b}^{m}$ defined by

$\hat{b}^{m}=\max(H(v^{m-1})*b^{m},u^{m-1});$ (33)

if $H(v^{m-1}_{i})=0$, that is $u^{m-1}_{i}=\psi_{i}$, then
$\hat{b}_{i}^{m}=u^{m-1}_{i}$ and $u^{m}_{i}=u^{m-1}_{i}$; otherwise
$\hat{b}_{i}^{m}=b_{i}^{m}$ , since $b_{i}^{m}\geq u^{m-1}_{i}$, and all
remains as before. No rebound is still possible after a contact.

  2. 2.

Scheme S2: solves problem (4) with CQ for the Caputo derivative and the semi
implicit f.d. scheme of [1] in space.

In this case the Caputo derivative is approximated through the so-called
convolution quadrature (CQ) method, proposed by Lubich for the discretization
of Volterra integral equations. In particular, for a function $\phi(t)$ (with
$\phi(0)=0$), the Riemann-Liouville derivative
${}^{R}\partial^{\alpha}_{t}\phi$ defined in (7) can be approximated by the
discrete convolution:

${}^{R}\partial^{\alpha}_{\tau}\phi^{m}:=\frac{1}{\tau^{\alpha}}\sum_{j=0}^{m}c_{j}\phi^{m-j},$

where $\phi^{m}=\phi(t_{m})$, and the coefficients $\\{c_{j}\\}$ are obtained
from a suitable power series expansion, connected to a specific approximation
method for the ODE (see [9]). In the case of the Euler backward method, it is
known as the Grunwald-Letnikov approximation, and provides the following
recursive formula for the coefficients:

$c_{0}=1,\quad c_{j}=-\frac{\alpha-j+1}{j}c_{j-1}.$

Then, using the relation (8) between the Caputo and the Riemann-Liouville
derivatives we can rewrite the initial problem as

${}^{R}\partial^{\alpha}_{t}(u-u^{0})-H(u-\psi)\Delta u=0,$

which discretized in time and space (with the same notations of S1) becomes:

$\frac{1}{\tau^{\alpha}}\sum_{j=0}^{m}c_{j}(u^{m-j}-u^{0})+\frac{1}{h^{2}}H(u^{m-1}-\psi)Au^{m}=0\
,$

equivalent to the solution at any iteration of the linear system

$B^{m-1}u^{m}=b^{m},$ (34)

where this time we have set

$B^{m-1}=I+\gamma_{\alpha}H(u^{m-1}-\psi)*A,\qquad
b^{m}=u^{0}-\sum_{j=1}^{m-1}c_{j}(u^{m-j}-u^{0})\ ,$ (35)

followed again by the correction

$u^{m}_{i}=\max(u^{m}_{i},\psi(x_{i}))\ .$ (36)

Even in this case the vector $b^{m}$ has to be modified inside the contact set
in order to remove the memory effect and prevent rebounds, as done in (33). In
fact when $u^{m-1}_{i}=\psi_{i}$ then again $b^{m}_{i}=0$, and from (35) we
get

$u^{m}_{i}=u^{0}_{i}-\sum_{j=1}^{m-1}c_{j}(u^{m-j}_{i}-u^{0}_{i})>u^{0}_{i}-(\psi_{i}-u^{0}_{i})\sum_{j=1}^{m-1}c_{j}>\psi_{i}\
,$

since $u^{m-j}_{i}\geq\psi_{i}$ for any $j$, and
$\sum_{j=1}^{m-1}c_{j}\geq-1$.

  3. 3.

Scheme S3: solves problem (14) with L1 for the Caputo derivative and the
scheme of [3] for the evolutive obstacle problem.

If we discretize the equation of system (14) through finite differences, using
the L1 scheme for the Caputo derivative, we get the equation

$(u^{m}-\psi)^{T}\left(\frac{1}{g\tau^{\alpha}}(u^{m}-\sum_{k=0}^{m-1}C_{m,k}u^{k})+\frac{1}{h^{2}}Au^{m}\right)=0\
.$

Setting $y^{m}=u^{m}-\psi$, and remembering that $\sum_{k=0}^{m-1}C_{m,k}=1$,
it is equivalent to

$y^{m}\left(y^{m}+g\gamma_{\alpha}A(y^{m}+\psi)-\sum_{k=0}^{m-1}C_{m,k}y^{k}\right)=0\
.$

Then $y^{m}=\max(0,x^{m})$ is solution of the previous equation if $x^{m}$
solves

$(I+g\gamma_{\alpha}AP(x))x=b^{m}\ ,$ (37)

where now

$b^{m}=\sum_{k=0}^{m-1}C_{m,k}y^{k}-g\gamma_{\alpha}A\psi=\sum_{k=0}^{m-1}C_{m,k}u^{k}-\psi-g\gamma_{\alpha}A\psi\
,$ (38)

while $P(x)=diag(p(x_{i}))$ denotes the diagonal matrix with $p(x_{i})=1$ if
$x_{i}>0$ and $p(x_{i})=0$ otherwise. As seen in [3], (37) can be solved by
the so-called Picard iterations:

$P^{0}=O,\quad(I+g\gamma_{\alpha}AP^{n})x^{n+1}=b^{m},\quad
P^{n+1}=diag(p(x^{n+1}))\quad\hbox{ per }n=0,1,...$

until $P^{n+1}=P^{n}$ ($O$ is the null matrix); at that point $x^{n+1}$ is the
sought solution.

Of course other schemes could be obtained by different combinations of
specific numerical approaches, but for our purposes the three previous schemes
were sufficient to perform in the next section explicit simulations of the
problem.

##  4 Numerical tests

We have applied the schemes described in Section 3 to some specific examples,
for different values of $\alpha$. We have choosen a sufficiently large final
time $T$, but also added a stopping time criterium in order to put in evidence
the convergence towards the asymptotic solution. Since this convergence
corresponds to the stabilization of the solution vector and to the
satisfaction of the asymptotic complementarity relation

$(u-\psi)\Delta u=0,$

which means $u$ harmonic (linear in 1D) outside the contact set, we have used
the criterium

$STOP\quad
when\quad\max(\|u^{m}-u^{m-1}\|_{\infty},\|\left(u^{m}-\psi)Au^{m}\right\|_{\infty})<tol\
,$ (39)

for a given tolerance parameter $tol$.

Example 1. $\Omega=(-1,1)$, $u^{0}(x)=0.7-0.7x^{2}$, $\psi(x)=0.5-2x^{2}$ (see
Fig.1).

Figure 1: Data of Example 1.

In the following Table 1 we reported some results obtained by the simulations
with schemes S1, S2 and S3, for different values of $\alpha$, $N$ and
$\gamma_{\alpha}$, with $tol=10^{-4}$. We adopted the following notations:

  * •

FC time = full contact time, that is the first time at which the solution has
reached the whole contact set (no more changing in the successive iterations);

  * •

STOP time = the exit time according to criterium (39);

  * •

# iter. = final number of time iterations;

  * •

# Pic. = average number of Picard iterations for each time step in scheme S3;

  * •

# LS = approximate number of linear systems to be solved (essentially the
product of the previous two quantities)

  * •

S = working schemes (the ones detecting the correct contact set).

Looking at the table some remarks and comments are possible:

Table 1:  $\alpha$ | $N$ | $\gamma_{\alpha}$ | $\tau$ | FC time | STOP time | # iter. | # Pic. | # LS | S
---|---|---|---|---|---|---|---|---|---
0 | 32 | 256 | any | first it. | no conv. | 1 | 14 |  | S3
| 64 | 1024 | any | first it. | no conv. | 1 | 29 |  | S3
| 128 | 4096 | any | first it. | second it. | 1 | 57 |  | S3
0.3 | 32 | 60 | 0.0079 | 0.10 | $93.16$ | 11739 | 12 | 140868 | all
| 32 | 75 | 0.016 | 0.11 | $93.15$ | 5580 | 13 | 72540 | S3
| 64 | 220 | 0.0059 | 0.49 | 1.33 | 226 | 21 | 4746 | all
| 128 | 850 | 0.005 | 0.2 | 1.66 | 315 | 42 | 13230 | all
0.5 | 32 | 25 | 0.009 | 0.16 | 8.38 | 880 | 8 | 7040 | all
| 32 | 50 | 0.038 | 0.19 | 8.39 | 221 | 11 | 2431 | S3
| 64 | 100 | 0.009 | 0.39 | 2.13 | 225 | 15 | 3375 | all
| 128 | 400 | 0.009 | 0.8 | 2.67 | 282 | 28 | 7896 | all
0.7 | 32 | 50 | 0.097 | 0.29 | 5.81 | 61 | 11 | 671 | all
| 32 | 100 | 0.026 | 0.52 | 10.44 | 41 | 14 | 574 | S3
| 64 | 200 | 0.097 | 0.38 | 7.08 | 74 | 21 | 1554 | all
| 128 | 400 | 0.036 | 0.5 | 4.82 | 135 | 29 | 3915 | all
1 | 32 | 15 | 0.058 | 0.27 | 0.99 | 18 | 7 | 126 | all
| 32 | 20 | 0.078 | 0.31 | 1.09 | 15 | 8 | 120 | S3
| 64 | 60 | 0.058 | 0.35 | 1.05 | 19 | 13 | 247 | all
| 128 | 240 | 0.058 | 0.64 | 1.05 | 19 | 23 | 437 | all

  * •

The stationary solution is the same for any $\alpha$, as stated by Theorem
2.1, and corresponds to the one of the stationary problem (15). Different is
only the speed of convergence. The detected right extremum of the symmetric
contact set $C$ on the used meshes is 0.125 (the continuous value should be
approximately $0.132$).

  * •

Schemes S1 and S2 have a very similar behaviour; in both cases for large
values of $\gamma_{\alpha}$ the contact set can be overestimated, a problem
not present for S3, due to the implicit nature of the quasi-Newton scheme of
[3]. The semi implicit schemes S1 and S2 pay for the delay with which the
contact information with the obstacle is achieved, allowing an uncorrect
evolution of the solution. To avoid that, strong $\tau$-step restrictions are
necessary: experiments show that all the schemes work correctly for example if
$\tau^{\alpha}<0.1$, a bound which becomes particularly heavy when $\alpha$ is
small. In the Table we reported approximately for any number of nodes the
largest values of $\gamma_{\alpha}$ for which all the three schemes give the
same correct solution.

  * •

Computational costs: the previous remark suggests that S3 is the more reliable
and even the less expensive of the three schemes, allowing larger time steps
and hence less iterations. Anyway, any single time iteration of S3 is much
more expensive, many Picard iterations (growing with $N$) with respect to a
single linear system necessary to be solved for S1 and S2. Then, comparing the
computational cost of the schemes, even these two schemes reveal competitive
in terms of the total number of linear system solved.

  * •

For a fixed $\alpha$ the full contact time grows with the number of nodes, and
does not seem to depend from $\gamma_{\alpha}$. On the contrary the
stabilization time grows with $\gamma_{\alpha}$ but decreases with the number
of nodes, since a finer mesh reduces the error at the boundary of the contact
set.

  * •

All the schemes correctly work also for the case $\alpha=1$: it easy to see
that in that case both the used approximations of the Caputo derivative reduce
to the standard incremental ratio in time.

  * •

The case $\alpha=0$ is a sort of control test: since in that case
$\partial^{\alpha}_{t}u=u-u^{0}$, in absence of an obstacle the equation (4)
would reduce to the stationary equation

$-\Delta\overline{u}+\overline{u}=u^{0}.$ (40)

In presence of an obstacle we expect the discrete solution to satisfy (40)
only outside the contact set, and from the first iteration. It is in fact
clear from (31) and (35) that since $H(u^{0}-\psi)=1$ and
$\gamma_{\alpha}=1/h^{2}$, the first iteration of all the schemes becomes

$B^{0}u^{1}=\left(I+\frac{1}{h^{2}}A\right)u^{1}=u^{0}\ ,$

which is essentially the discrete version of (40). If $u^{1}$ goes over the
obstacle, such identity will be satisfied only where $u^{1}>\psi$. In Figure 2
it can be seen what happens in our example with $N=128$ nodes and scheme S3:
the solution $u$ is plotted in blue, the obstacle in red, the initial datum
$u^{0}$ in black and the quantity $-\Delta u+u$ through asterisks: the last
two quantities coincide in the detachment set, with a natural discontinuity at
the boundary of such a region. Since the time step $\tau$ has no effect on the
solution, the semi implicit schemes $S1$ and $S2$ for this example always
overestimate the contact set.

Figure 2: Solution for $\alpha=0$ and 128 nodes, scheme S3.

Example 2. $\Omega=(-1,1)$, $u^{0}(x)=1-x^{2}$, $\psi(x)=0.5-(2x^{2}-0.5)^{2}$
(see Fig.3).

Figure 3: Data of Example 2.

On this example we tested numerically the estimate (26) of Theorem 2.1, using
scheme S3. The $L^{1}$ norm of the error at time $t^{m}$ on the given mesh was
approximated by a natural quadrature formula, that is

$\left\|u(t^{m})-\overline{u}\right\|_{L^{1}\left(\Omega\right)}\simeq
h\sum_{i}|u^{m}_{i}-\bar{u}_{i}|\ $

(the vector $\bar{u}$ on the mesh was computed in advance with a sufficiently
high precision). Such an error was computed for different values of $\alpha$
and of the time $t^{m}$. In Table 2 we reported (in the third column) the
discrete $L^{1}$ error at $T=10$ with $N=32$ nodes and the same
$\gamma_{\alpha}=50$ for any $\alpha$. In the fourth column the corresponding
values of quantity $J(T)=\frac{1}{C}\frac{T^{-\alpha}}{\Gamma(1-\alpha)}$ of
(28) are shown; for the constant $C$ we adopted a computed estimate ($C=26$).
It is evident that the two quantities decay at the same rate. In the second
column it is also possible to see the number of instant times $M$ needed for
each $\alpha$ to keep the same $\gamma_{\alpha}$, and the consequent
increasing complexity of the computation. In order to confirm the order of
decay of the error, polynomial for $\alpha\in(0,1)$ and exponential (up to the
machine precision) for $\alpha=1$, we illustrate in Figure 4 such behavior
plotting the continuous curves $J(t)$ in the time interval $(0,20)$ for
$\alpha=0.3,0.6$ and $0.9$, and the computed errors on ten different discrete
times (marked by asterisks). These values correctly lie on the curves.

Figure 4: Example 2: polynomial speed of stabilization for $\alpha=0.3,0.6$ and $0.9$. Table 2: Example 2: $N=32$, $\gamma_{\alpha}=50$, scheme $S3$. $\alpha$ | $M$ |  $L^{1}$ error at $T=10$ | ${T^{-\alpha}}/{(C\Gamma(1-\alpha))}$
---|---|---|---
1 | 51 | $8.07\ 10^{-6}$ | $\simeq 0$
0.9 | 61 | $5.27\ 10^{-4}$ | $5.09\ 10^{-4}$
0.8 | 77 | $1.37\ 10^{-3}$ | $1.32\ 10^{-3}$
0.7 | 103 | $2.62\ 10^{-3}$ | $2.56\ 10^{-3}$
0.6 | 152 | $4.41\ 10^{-3}$ | $4.35\ 10^{-3}$
0.5 | 262 | $6.87\ 10^{-3}$ | $6.86\ 10^{-3}$
0.4 | 593 | $1.01\ 10^{-2}$ | $1.02\ 10^{-2}$
0.3 | 2313 | $1.44\ 10^{-2}$ | $1.48\ 10^{-2}$
0.2 | 35184 | $1.98\ 10^{-2}$ | $2.08\ 10^{-2}$
0.1 | 123000000 |  $2.7\ 10^{-2}$ (estimate) | $2.85\ 10^{-2}$

## References

  * [1] Alberini C., Capitanelli R. and Finzi Vita S., A numerical study of an Heaviside function driven degenerate diffusion equation, preprint 2020, submitted.
  * [2] Bazhlekova E.G., Subordination principle for fractional evolution equations, Fractional Calculus and Applied Analysis. 3, No. 3, (2000), 213–230.
  * [3] Brugnano L. and Sestini A., Iterative solution of piecewise linear systems for the numerical solution of obstacle problems, J. Num. Anal. Ind. Appl. Math. (JNAIAM), 6, No. 3-4, (2011), 67–82.
  * [4] Capitanelli R. and D’Ovidio M., _Fractional equations via convergence of forms_ , Fractional Calculus and Applied Analysis, 22, (2019), 844–870.
  * [5] Diethelm K., _The analysis of fractional differential equations. An application-oriented exposition using differential operators of Caputo type_ , Lecture Notes in Mathematics, 2004. Springer-Verlag, Berlin, (2010).
  * [6] D’Ovidio M., _On the fractional counterpart of the higher-order equations_ , Statistics and Probability Letters, 81, (2011), 1929–1939.
  * [7] Giga Y., Liu Q. and Mitake H., On a discrete scheme for time fractional fully nonlinear evolution equations, Asymptotic Analysis, 120, No. 1-2, (2020), 151–162.
  * [8] Gorenflo R., Kilbas A.A., Mainardi F. and Rogosin S.V., _Mittag-Leffler Functions, Related Topics and Applications_ , Springer Monographs in Mathematics, Springer-Verlag Berlin Heidelberg (2014).
  * [9] Jin B, Lazarov R. and Zhou Z., Numerical methods for time-fractional evolution equations with nonsmooth data: a concise overview, Comput. Methods Appl. Mech. Engrg., 346, (2019), 332–358.
  * [10] Kochubei A.N., The Cauchy problem for evolution equations of fractional order, Differential Equations, 25, (1989), 967–974.