††thanks: Corresponding author

# Carathéodory’s thermodynamics of the Schwarzschild black hole surrounded by
quintessence

Mohsen Fathi  Departamento de Física, Universidad de Santiago de Chile,
Avenida Víctor Jara 3493, Santiago, Chile  Martín Molina  Instituto de Física,
Pontificia Universidad Católica de Valparaíso, Avenida Brasil 2950,
Valparaíso, Chile  J.R. Villanueva<EMAIL_ADDRESS>Instituto de Física
y Astronomía, Universidad de Valparaíso, Avenida Gran Bretaña 1111,
Valparaíso, Chile

###### Abstract

In this paper, we apply the Carathéodory’s method of geometrothermodynamics to
investigate the behavior of the main thermodynamic parameters associated with
a Schwarzschild black hole surrounded by quintessence. The corresponding
Pfaffian form is constructed by means of the Schwarzschild radius $r_{s}$, and
the quintessential radius $r_{\gamma}$, as independent variables. This form is
then used to characterize the thermodynamic manifold. The homogeneity of the
system allows for the recognition of the empirical temperature and entropy,
and thus, connects with the usual laws of thermodynamics. In particular, we
show that the Helmholtz and Gibbs free energies lead to the same value for the
Schwarzschild black hole, in the case of the vanishing cosmological term.

keywords: Adiabatic processes, black hole thermodynamics, quintessence

PACS numbers: 04.20.Fy, 04.20.Jb, 04.25.-g

###### Contents

  1. I Introduction and Motivation
  2. II The black hole solution in the dark background
  3. III The Carathéodory thermodynamics applied to the black hole
  4. IV The adiabatic-isoareal processes and the extremal limit
    1. IV.1 Thermodynamic limit
  5. V The heat capacity and the free energies
  6. VI Conclusions
  7. A Derivation of the solution to the Cauchy problem

##  I Introduction and Motivation

Attributing the laws of thermodynamics to black holes, as proposed in the
1970’s in a series of seminal papers Bardeen:1973gs ; Bekenstein:1972 ;
Bekenstein:1973 ; Bekenstein:1974 ; Bekenstein:1975 ; Hawking:1974sw , has
opened new gates in the study of these mysterious objects, from both
theoretical and astrophysical points of view. As an example, and from the
theoretical side, the concept of entropy is applied to black holes by means of
the famous Bekenstein-Hawking (B-H) entropy formula, although its direct
application to the extremal black holes (EBHs) is controversial because the
zero entropy conjecture for EBHs Teitelboim:1994 ; Carroll:2009 , neglects the
direct relationship between the entropy and the event horizon’s area.
Regarding the experimental and observational evidences, the technological
advancements have then provided facilities to do some tests on analogue
models, in order to unravel the links between thermodynamic entities and black
hole dynamics Unruh:1981 ; Novello:2002 ; Schutzhold:2005 ; Carusotto:2008 ;
Belgiorno:2010 ; Weinfurtner:2011 ; Castelvecchi:2016 ; Steinhauer:2016 ;
Lima:2019 ; Kolobov:2021 .

It is, however, worth mentioning that the axiomatic settings of the laws of
classical thermodynamics, themselves, are not always given in the conventional
ways. In fact, modern thermodynamics is the result of a long-term process, in
which, the laws have been given rigorous tests and undergone different
experiments (see the reviews in Refs. Hutter:1977 ; Neumaier:2007 ). Among
all, the geometric approach proposed by Carathéodory in Ref. caratheodory09 ,
is significant (see the reviews in Refs. Landsberg:1956 ; Antoniou:2002 ).
Together with the approach given by Gibbs Gibbs:1949 , the Carathéodory’s
method and its further developments by Born Born:1949 , form the foundations
of the so-called Geometrothermodynamics Quevedo:2007 . The link between the
methods of Carathéodory and Gibbs, is however, argued to be established in
terms of the homogeneity of the Pfaffian form $\delta Q_{\mathrm{rev}}$, as
the infinitesimal heat exchange reversibly Belgiorno:2002 , and this method
has been applied to the laws of black hole thermodynamics in Refs.
Belgiorno:2002iv ; Belgiorno:2002iw ; Belgiorno:2003a ; Belgiorno:2003b ;
Belgiorno:2004 , and recently, in Refs. Molina:2021 ; Fathi:2021EPJC ;
Fathi:2021PLB  regarding the adiabatic (isoareal) processes of Hayward and BTZ
BHs. The geometric formulation of the Carathéodory’s method, makes it possible
to have a self-contained study of the black hole thermodynamics, by using only
the respected spacetime structure.

If the spacetime is coupled with cosmological parameters, its thermodynamics
also reveals the evolutionary structure of the universe, in which the black
hole resides. This problem is of our interest in this paper. As it is
elaborated in the next section, we take into account a black hole spacetime
that is coupled with a quintessential dark field and we apply the Carathéodory
thermodynamics to explore the possible adiabatic processes, based on the
solutions to the Pfaffian in the context of the black hole geometry. To
elaborate this, we calculate the analytical solutions to the corresponding
Cauchy problem, and accordingly, we determine the allowed physical paths on
the thermodynamic manifold.

The paper is organized as follows: In Sect. II, we introduce the spacetime and
its causal structure. In Sect. III, we introduce the Carathéodory
geometrothermodynamics and express the black hole thermodynamic parameters as
functions of its spacetime components. In this section, we calculate the
geometric entropy and temperature, which lead to the interpretation of the
Pfaffian in the context of the first law. These parameters are then applied in
order to demonstrate the entropy-temperature behavior within the thermodynamic
foliations. In Sect. IV, we present the Cauchy problem in the context of
isoareal processes, to find the permissible trajectories for adiabatic
processes. In Sect. V, the previously calculated thermodynamic parameters are
applied to calculated the heat capacity and the free energies of the black
hole, and the appropriate limits are discussed. We conclude in Sect. VI.
Throughout this work, we apply a geometrized system of units, in which
$G=c=1$.

##  II The black hole solution in the dark background

The current study is aimed at the continuation of applying the Carathéodory’s
method to black hole thermodynamics. To elaborate this purpose, we however,
choose to include the cosmological dynamics that feature as an evolutionary
characteristic of black hole geometries. This way, one needs to take into
account the effects from the dark side of the universe. Such features have
been discussed extensively in general relativity and alternative gravity
theories (see for example Ref. JimenezMadrid:2005 ; Jamil:2009 ; Li:2019 ;
Roy:2020 ). These may include the presence of a dark matter halo Xu:2018 ;
Das:2021 , or coupling with a quintessential field Kiselev:2003 ; Saadati:2019
; AliKhan:2020 ; CFOV21 . The standard thermodynamics of static black holes in
quintessence has been studied extensively, for example in Refs. Tharanath:2013
; Ghaderi2016 ; Ma:2017 ; Ghosh:2018 ; Rodrigue:2018 ; Nam:2020 ;
Lutfuoglu:2021 . We are, however, interested in investigating the
geometrothermodynamics of such a black hole (a Schwarzschild black hole)
through the Carathéodory’s method, in order to find the limits imposed on the
corresponding thermodynamic manifold, that are not accessible by adiabatic
processes.

The static, spherically symmetric black hole solution associated with
quintessence is described by the line element

$\mathrm{d}s^{2}=-B(r)\mathrm{d}t^{2}+B^{-1}(r)\mathrm{d}r^{2}+r^{2}\mathrm{d}\theta^{2}+r^{2}\sin^{2}\theta\mathrm{d}\phi^{2}$
(1)

in the $x^{\mu}=(t,r,\theta,\phi)$ coordinates, where the lapse function is
given by Kiselev:2003

$B(r)=1-\frac{r_{s}}{r}-\frac{\gamma}{r^{3w_{q}+1}},$ (2)

with $r_{s}=2M$, $\gamma$ and $w_{q}$, representing the parameters of
quintessence and the equation of state (EoS), and $M$ is the black hole’s
mass. For an accelerating universe, the EoS parameter respects the range
$-1<w_{q}<-\frac{1}{3}$, and the particular case of $w_{q}=-1$ recovers the
cosmological constant. In this paper, we confine ourselves to the case of
$w_{q}=-\frac{2}{3}$, that recovers

$B(r)=1-\frac{r_{s}}{r}-\gamma r,$ (3)

and accordingly, $\gamma$ has dimension of $L^{-1}$. Now let us recast the
line element as

$B(r)=1-\frac{r_{s}}{r}-\frac{r}{r_{\gamma}},$ (4)

by defining $r_{\gamma}\doteq\frac{1}{\gamma}$ as the quintessential radius.
The black hole horizons that are located at the radial distances $r_{h}$ at
which $B(r_{h})=0$, are therefore given by

$\displaystyle
r_{++}=\frac{r_{\gamma}}{2}\left(1+\sqrt{1-\frac{4r_{s}}{r_{\gamma}}}\right),$
(5) $\displaystyle
r_{+}=\frac{r_{\gamma}}{2}\left(1-\sqrt{1-\frac{4r_{s}}{r_{\gamma}}}\right),$
(6)

which are, respectively, the (quintessential) cosmological and the event
horizons. Hence, the extremal black hole, for which
$r_{+}=r_{++}=r_{e}=2r_{s}$, corresponds to the case of $r_{\gamma}=r_{\gamma
e}=4r_{s}$, and a naked singularity is occurred when $r_{\gamma}<r_{\gamma
e}$.

##  III The Carathéodory thermodynamics applied to the black hole

The Carathéodory’s framework of thermodynamics is based on the Carathéodory’s
principle which reads: in the neighbourhood of any arbitrary state $J$ of a
thermally isolated system $\Sigma$, there are states $J^{\prime}$ which are
inaccessible from $J$ chandrasekhar39 ; Buchdahl:1966 ; adkins_1983 . This
inaccessibility may be established on the integrability of the appropriate
Pfaffian form $\delta Q_{\mathrm{rev}}$ for the system, which implies that it
can be written as

$\delta Q_{{\rm rev}}=\tau\mathrm{d}\sigma,$ (7)

where $\tau$ is an integrating factor which is considered to be the empirical
temperature, and $\sigma$ is the empirical entropy chandrasekhar39 ;
Buchdahl:1966 ; adkins_1983 . The existence of an integrating factor ensures
the existence of an infinite number of them. So if the integrating factor is
considered to be the absolute temperature $T$, then the Pfaffian form reads

$\delta Q_{{\rm rev}}=T{\rm d}S,$ (8)

where $S$ is the metric entropy, that is related to the second law for the
irreversible processes. The integrating factor can be calculated, once $\delta
Q_{\mathrm{rev}}$ represents a symmetry. In this sense, the thermodynamics
manifold is constructed by the foliation of adiabatic hyper-surfaces on which,
$\delta Q_{\mathrm{rev}}$$=0$. Therefore, the homogeneity properties of the
Pfaffian form allows to connect the Carathéodory’s framework with that of
Gibbs Belgiorno:2002 ; Belgiorno:2002iv , and then, all thermodynamics can be
applied to obtain the relevant macroscopic magnitudes. Thus, the black hole
thermodynamics can be studied using this approach which leads to interesting
properties.

In particular, our interest is in the study of the system characterized by the
variables $(r_{s},r_{\gamma})$, and hence, are chosen as the independent
extensive variables of the thermodynamic manifold with the constraint
$1-\frac{4r_{s}}{r_{\gamma}}>0$, as inferred from Eqs. (5) and (6).
Accordingly, this manifold is bounded by the extremal sub-manifold,
corresponding to $r_{\gamma}>4r_{s}$.

With this in mind, we postulate that the Pfaffian form $\delta Q_{{\rm rev}}$
can be written as

$\delta Q_{\mathrm{rev}}=\mathrm{d}r_{s}-\Gamma\mathrm{d}r_{\gamma},$ (9)

in terms of the system variables $(r_{s},r_{\gamma})$, where
$\Gamma\equiv\Gamma(r_{s},r_{\gamma})$ is regarded as the ”generalized force”
associated with the quintessence contribution. This coefficient is supposed to
be non-zero everywhere on the thermodynamic domain, so that the Pfaffian is
always non-singular, and its integrability is guaranteed by the condition
$\delta Q_{\mathrm{rev}}$$\wedge\mathrm{d}$($\delta Q_{\mathrm{rev}}$)$\,=0$
Belgiorno:2002iv . In this way, the $\mathrm{d}r_{s}$ term in Eq. (9) plays
the role of the infinitesimal changes in the black hole’s internal energy,
whereas $-\Gamma\mathrm{d}r_{\gamma}$ is a work term which is identified,
completely, in geometrical contexts.

The determination of these coefficients is based on the contribution of the
metric entropy in the Pfaffian. In fact, the B-H entropy relation
Hawking:1974sw ; Majumdar:1998

$S=\frac{k_{\mathrm{B}}\mathcal{A}_{+}}{4\ell_{\mathrm{p}}^{2}},$ (10)

with $k_{\mathrm{B}}$, $\mathcal{A}_{+}=4\pi r_{+}^{2}$ and
$\ell_{\mathrm{p}}$, being respectively the Boltzmann constant, the event
horizon area, and the Planck length, implies that $S\equiv
S(r_{s},r_{\gamma})$. Introducing $\tilde{a}\equiv\frac{\pi
k_{B}}{\ell_{p}^{2}}$, one can define an entropy function as

$\mathcal{S}(r_{s},r_{\gamma})\equiv\frac{S}{\tilde{a}}=r_{s}^{2}\,R_{+}^{2}(r_{s},r_{\gamma}),$
(11)

where

$R_{+}(r_{s},r_{\gamma})=\frac{r_{\gamma}}{2r_{s}}\left(1-\sqrt{1-\frac{4r_{s}}{r_{\gamma}}}\right).$
(12)

It can be verified that, for any real-valued constant $\lambda$, we have
$R_{+}(\lambda r_{s},\lambda r_{\gamma})=R_{+}(r_{s},r_{\gamma})$. So $R_{+}$
is homogeneous of degree zero and, therefore, is an ”intensive” parameter.
Since the temperature function $\mathcal{T}$ is an integration factor for the
Pfaffian form $\delta Q_{{\rm rev}}=\mathcal{T}\mathrm{d}\mathcal{S}$, it is
obtained from the relation

$\mathcal{T}=\left(\frac{\partial\mathcal{S}}{\partial
r_{s}}\right)^{-1}_{r_{\gamma}}=\frac{1}{r_{\gamma}}\frac{\sqrt{1-\frac{4r_{s}}{r_{\gamma}}}}{1-\sqrt{1-\frac{4r_{s}}{r_{\gamma}}}},$
(13)

which is homogeneous of degree $-1$. It is informative to demonstrate the
mutual behavior of the above thermodynamic parameters in a
$\mathcal{S}$-$\mathcal{T}$ diagram (see Fig. 1). As observed from the figure,
for a fixed $r_{s}$, the $\Delta\mathcal{S}>0$ condition corresponds to
$\Delta r_{\gamma}<0$. Therefore, by varying $r_{\gamma}$ in a particular
$r_{s}$-constant foliation, the system transits towards the Schwarzschild
black hole (SBH). It is also straightforward to verify that, going from the
state (1) to the state (2) (for which $\mathcal{T}_{1}<\mathcal{T}_{2}$), the
variable $r_{s}$ increases. Hence, one can infer that in an adiabatic process,
both of the variables $(r_{s},r_{\gamma})$ increase.

Figure 1: The $\mathcal{S}$-$\mathcal{T}$ diagram indicating both the EBH and
the SBH limits. From bottom to the top, the solid curves indicate the
thermodynamic foliations for $r_{s}=0.4,0.5,0.625,0.75,0.875,1$. The green
arrow indicates an adiabatic process, going from state (1) to state (2).

The generalized force is obtained from

$\Gamma=\mathcal{T}\left(\frac{\partial\mathcal{S}}{\partial
r_{\gamma}}\right)_{r_{s}}=\frac{1}{2}\left(1-\sqrt{1-\frac{4r_{s}}{r_{\gamma}}}-\frac{2r_{s}}{r_{\gamma}}\right),$
(14)

which is an intensive function. Note that, the extremal case corresponds to
$\Gamma_{e}\equiv\Gamma(r_{s},r_{\gamma e})=\frac{1}{4}$, which according to
Eq. (13), is the generalized force for a black hole of zero temperature
($\mathcal{T}=0$). Thus, performing $(r_{s},r_{\gamma})\mapsto(\lambda
r_{s},\lambda r_{\gamma})$, we get $\delta Q_{{\rm rev}}\mapsto\lambda\delta
Q_{{\rm rev}}$, which means that the Pfaffian form is homogeneous of degree
one. In this way, we have an Euler vectorial field, or a Liouville operator,
as the infinitesimal generator of the homogeneous transformations

$D\equiv r_{s}\frac{\partial}{\partial
r_{s}}+r_{\gamma}\frac{\partial}{\partial r_{\gamma}}.$ (15)

In Fig. 2, the $\Gamma$-$\mathcal{T}$ diagram has been plotted for several
fixed values of $r_{\gamma}$, indicating that the raise in $r_{\gamma}$ (i.e.
the decrease in $\gamma$), increases the negative slope of the curves towards
the asymptot, while $\mathcal{T}$ is increasing.

Figure 2: The adiabatic foliations given in terms of the mutual behavior of
$\Gamma$ and $\mathcal{T}$, for some values of $r_{\gamma}$ and varying
$r_{s}$.

##  IV The adiabatic-isoareal processes and the extremal limit

The correct construction of the thermodynamic manifold leads us to the study
of the adiabatic processes that represent the foliations of this manifold.
Thus, we must ensure that this construction is consistent with the usual laws
of black hole thermodynamics, and in particular, the status of the third law
and the connection between extremal and non-extremal states.

###  IV.1 Thermodynamic limit

We start by analyzing the surface generated by the extremal states and its
implication to the construction of the entire thermodynamic manifold. In fact,
if $\mathbf{r}_{s}\equiv r_{s}^{e}=\Gamma_{e}r_{\gamma}$ is the extremal value
of $r_{s}$, then the area of the extremal states can be written as

$\mathcal{A}_{e}=4\pi(r_{+}^{e})^{2}=4\pi(r_{s}^{e}R_{+}^{e})^{2}=16\pi\mathbf{r}_{s}^{2},$
(16)

which implies that

${\rm d}\mathcal{A}_{e}=\pi r_{\gamma}{\rm d}r_{\gamma},$ (17)

and thus, the isoareal condition d$\mathcal{A}_{e}=0$ is only satisfied by
$r_{\gamma}=$ const. Therefore, although the transformations between extremal
states are adiabatic, they are not isoareal. This indicates that the area-
entropy law is not valid for extremal states. We can fix this as it has been
established previously by letting $\mathcal{S}=0$ on this manifold
Teitelboim:1994 ; Carroll:2009 . Despite this, we can also try using a
different criterion, as proposed by Lemos in lower-dimensional gravity
Lemos:2017aol . For now, we consider that the extremal sub-manifold
constitutes the boundary of the thermodynamic manifold, which will be studied
later.

As discussed in the previous section, the Pfaffian form is responsible for
constructing the physically accepted thermodynamic manifold. In this sense,
the Cauchy problem $\delta Q_{\mathrm{rev}}$$=0$ generates the non-extremal
isentropic (i.e. adiabatic and reversible) sub-manifolds of the thermodynamic
foliations Belgiorno:2002iv . Let us rewrite the Pfaffian (9) as

$\delta
Q_{\mathrm{rev}}=\frac{\mathrm{d}x}{2\sqrt{x}}-\frac{\Gamma(x,y)}{\Gamma_{e}}\frac{\mathrm{d}y}{2\sqrt{y}},$
(18)

by defining $x\equiv r_{s}^{2}$ and $y\equiv\Gamma_{e}^{2}r_{\gamma}^{2}$.
After some arrangements, these changes of variables yield

$\frac{\Gamma(x,y)}{\Gamma_{e}}=\left(1-\sqrt{1-\sqrt{\frac{x}{y}}}\,\right)^{2}.$
(19)

This way, the Cauchy problem leads to the differential equation

$\frac{\mathrm{d}y}{\mathrm{d}x}=\frac{\sqrt{\frac{y}{x}}}{\left(1-\sqrt{1-\sqrt{\frac{x}{y}}}\,\right)^{2}}=F(x,y),$
(20)

that mandates the condition $x<y$ on the thermodynamic manifold. Note that,
$y(x)=x$ is the solution to the above equation which means that extremal
states are adiabatically interconnected. According to our discussions, this is
in conflict with the statement of the second law and its connection with the
area (i.e. the B-H) formula. But they can still be reconciled if we consider
that both varieties are disconnected, and thus, the third law is also
preserved.

The thermodynamic manifold is, therefore, composed by the two mutually
inaccessible sub-manifolds, inferred from the following problems for
$(x_{0},y_{0})$ being an initial thermodynamic state:

  * •

On the $\mathcal{T}\neq 0$ sub-manifold,

$\displaystyle\frac{\mathrm{d}y}{\mathrm{d}x}=F(x,y),$ (21a) $\displaystyle
y(x_{0})=y_{0}>x_{0}.$ (21b)

  * •

On the $\mathcal{T}=0$ sub-manifold (extremal limit),

$\displaystyle\frac{\mathrm{d}y}{\mathrm{d}x}=F(x,y),$ (22a) $\displaystyle
y(x_{0})=y_{0}=x_{0}.$ (22b)

The Cauchy problem (20) has the general solutions (see appendix A)

$y_{i}(x;x_{0},y_{0})=\frac{x}{\left[1-\left(\alpha_{i}\sqrt{\frac{x}{x_{0}}}-1\right)^{2}\right]^{2}},$
(23)

where the $\alpha_{i}$’s depend on the initial condition by

$\displaystyle\alpha_{1}(x_{0},y_{0})\equiv\alpha_{1}=1+\sqrt{1+\sqrt{\frac{x_{0}}{y_{0}}}},$
(24a)
$\displaystyle\alpha_{2}(x_{0},y_{0})\equiv\alpha_{2}=1+\sqrt{1-\sqrt{\frac{x_{0}}{y_{0}}}},$
(24b)
$\displaystyle\alpha_{3}(x_{0},y_{0})\equiv\alpha_{3}=1-\sqrt{1-\sqrt{\frac{x_{0}}{y_{0}}}},$
(24c)
$\displaystyle\alpha_{4}(x_{0},y_{0})\equiv\alpha_{4}=1-\sqrt{1+\sqrt{\frac{x_{0}}{y_{0}}}}.$
(24d)

Since $y_{0}>x_{0}$, the hierarchy order for this parameters is
$\alpha_{1}>\alpha_{2}>\alpha_{3}>\alpha_{4}$. It is interesting to see the
particular cases obtained from some initial states. Firstly, if the initial
state is a pure quintessence (i.e. $x_{0}=0$), then $\alpha_{1}=\alpha_{2}=2$
and $\alpha_{3}=\alpha_{4}=0$, which implies that $y_{1}(x)=y_{2}(x)=0$ for
all $y_{0}$, and remains in that state indefinitely. This is while $y_{3}(x)$
and $y_{4}(x)$ diverge. Furthermore, if the initial state corresponds to the
extremal case, then this functions take the forms

$\displaystyle
y_{1}^{e}(x)=\frac{x}{\left[1-\left(\sqrt{\frac{x}{x_{0}}}\left(\sqrt{2}+1\right)-1\right)^{2}\right]^{2}},$
(25a) $\displaystyle
y_{2}^{e}(x)=\frac{x}{\left[1-\left(\sqrt{\frac{x}{x_{0}}}-1\right)^{2}\right]^{2}},$
(25b) $\displaystyle y_{3}^{e}(x)=y_{2}^{e}(x),$ (25c) $\displaystyle
y_{4}^{e}(x)=\frac{x}{\left[1-\left(\sqrt{\frac{x}{x_{0}}}\left(\sqrt{2}-1\right)+1\right)^{2}\right]^{2}},$
(25d)

which, obviously, must be removed from the set of permissible (physical)
solutions.

We can now summarize the properties of these functions considering the
following statements:

Proposition 1: Each function is finite at $x=0$ with a value

$y_{i0}\equiv\lim_{x\rightarrow 0^{+}}y_{i}(x)=\frac{x_{0}}{4\alpha_{i}^{2}},$
(26)

and thus, $0<y_{10}<y_{20}<y_{30}<y_{40}.$

Proposition 2: Each function diverges at

$x_{\infty}^{(i)}=\frac{4x_{0}}{\alpha_{i}^{2}},$ (27)

where $0<x_{\infty}^{(1)}<x_{\infty}^{(2)}<x_{\infty}^{(3)}<x_{\infty}^{(4)}$.

Proposition 3: For all $x_{0}\leq y_{0}$, the condition

$x_{\infty}^{(1)}<x_{0}<x_{\infty}^{(2)}$ (28)

holds.

Proposition 4: Each function intersects the straight line of the extremal
states (i.e. $y=x$), at

$\displaystyle x_{e}^{i}=\frac{x_{0}(\sqrt{2}+1)^{2}}{\alpha_{i}^{2}},$ (29a)
$\displaystyle x_{e^{\prime}}^{i}=\frac{x_{0}}{\alpha_{i}^{2}}.$ (29b)

In Fig. 3, the solutions in Eqs. (24) have been plotted for a particular
initial condition, where the extremal limit $y(x)=x$ has been distinctively
shown.

(a)  (b)

Figure 3: (a) The branches of $y(x)$ for $(x_{0},y_{0})=(2,3)$, shown as a
single blank point of intersection, plotted together with the extremal limit
(the straight line). The physically accepted paths are shown with solid
curves, whereas the dashed curves are not allowed. (b) The physically accepted
paths shown in a single diagram.

The physically acceptable branches, are however, those that respect the
correct changes of the thermodynamics coordinates. Accordingly, the permitted
paths are those that allow for simultaneous raise of $y(x)$ and $x$ for the
values $x<x_{0}$. This is while for $x>x_{0}$, the adiabatic paths should
respect a raise in $y(x)$ with a reduction in $x$. The corresponding adiabatic
surface that corresponds to the Cauchy problem has been also plotted in Fig.
4.

Figure 4: The adiabatic surface constructed by the solution (24). The extremal
(i.e. $\mathcal{T}=0$) and the SBH limits have been also indicated.

##  V The heat capacity and the free energies

The behavior of a thermodynamic system depends strongly on the extent, to
which, it can absorb heat. Essentially, we study this heat capacity by keeping
$r_{\gamma}=\mathrm{const.}$, and thus

$C_{\gamma}=\mathcal{T}\left(\frac{\partial\mathcal{S}}{\partial\mathcal{T}}\right)_{r_{\gamma}}=\left(\frac{\partial
r_{s}}{\partial\mathcal{T}}\right)_{r_{\gamma}}=\left(\frac{\partial\mathcal{T}}{\partial
r_{s}}\right)^{-1}_{r_{\gamma}}.$ (30)

Assuming the thermal equilibrium between the black hole and its environment,
then by applying Eq. (13) in Eq. (30), we get

$C_{\gamma}=-\frac{r_{\gamma}^{2}}{2}\left(1-\sqrt{1-\frac{4r_{s}}{r_{\gamma}}}\right)^{2}\sqrt{1-\frac{4r_{s}}{r_{\gamma}}},$
(31)

that reaches its minimum

$C_{\gamma}^{\mathrm{m}}(r_{s}^{\mathrm{m}},r_{\gamma})\equiv
C_{\gamma}^{\mathrm{m}}=-\frac{3\left(r_{s}^{\mathrm{m}}\right)^{2}}{2},$ (32)

at $r_{s}^{\mathrm{m}}=\frac{2r_{\gamma}}{9}$. It is of worth mentioning that,
the heat capacity of this black hole in the context of the generalized
uncertainty principle (GUP) has been calculated in Ref. Lutfuoglu:2021 , which
reads as

$\mathfrak{C}_{\gamma}=-\frac{\beta}{4}\frac{(1-\frac{2r_{+}}{r_{\gamma}})\sqrt{1-\frac{\beta}{4r_{+}^{2}}}}{(1-\frac{4r_{+}}{r_{\gamma}})\left(1-\sqrt{1-\frac{\beta}{4r_{+}^{2}}}\right)+\frac{\beta}{2r_{+}r_{\gamma}}},$
(33)

where $0<\beta<1$ is the deformation parameter proportional to the Planck
length, and accounts for the generalization of the Heisenberg uncertainty
principle (HUP)111Note that, there is also a $\pi$ factor included in Eq.
(33), which comes from the authors’ version of definition of the HUP.. Taking
into account the expression of $r_{+}$ in Eq. (6), it can be seen by direct
calculation, that $C_{\gamma}\rightarrow\mathfrak{C}_{\gamma}$ in the case of
$\beta\rightarrow 0$. In Fig. 5, the behavior of $C_{\gamma}$ has been plotted
for different values of $r_{\gamma}$. Each of the resultant branches
corresponding to fixed values of $r_{\gamma}$, has a minimum that satisfy the
condition $\frac{r_{s}^{\mathrm{m}}}{r_{\gamma}}=\frac{2}{9}$.

Figure 5: The behavior of $C_{\gamma}$ plotted for $r_{\gamma}=30,50,70,90$
and $110$. The dashed curve connects the minimum values
$C_{\gamma}^{\mathrm{m}}$, whereas the dotted-dashed curve corresponds to
$C_{\gamma}^{\mathrm{Sch}}$.

In the case of $\gamma\rightarrow 0$ (or $r_{\gamma}\rightarrow\infty)$, we
have that

$C_{\gamma}\approx-2r_{s}^{2}\left(1-\frac{2r_{s}}{r_{\gamma}}\right)\cong-2r_{s}^{2}=C_{\gamma}^{\mathrm{Sch}},$
(34)

which is the heat capacity for the SBH. Note that, the negativity of
$C_{\gamma}$, as seen in Fig. 5, implies that the black hole gets hotter as it
radiates. Furthermore, one can observe that $C_{\gamma}$ is continuous in the
range $0<r_{s}<\frac{r_{\gamma}}{4}$, which indicates the absence of any kind
of phase transition. It is also fruitful to calculate the free energies
associated with this heat capacity. According to the notions applied here and
in Ref. Molina:2021 , the Helmholtz free energy is give by
$\mathcal{F}=r_{s}-\mathcal{T}\mathcal{S}$, that using Eqs. (11) and (13),
yields

$\mathcal{F}_{\gamma}=\frac{r_{\gamma}}{4}\left(1-\sqrt{1-\frac{4r_{s}}{r_{\gamma}}}\right),$
(35)

which is presented for the case of $r_{\gamma}=$ const. In Fig. 6, the
behavior of the $\mathcal{F}_{\gamma}$ has been ramified for the same values
of $r_{\gamma}$ as those used in Fig. 5.

Figure 6: The behavior of $\mathcal{F}_{\gamma}$ plotted for
$r_{\gamma}=30,50,70,90$ and $110$. The dashed line connects the maximum
values $\mathcal{F}_{\gamma}^{\mathrm{M}}$, whereas the dotted-dashed line
corresponds to $\mathcal{F}_{\gamma}^{\mathrm{Sch}}$.

As it is inferred from the figure, the Helmholtz free energy exhibits the
upper limit $\mathcal{F}_{\gamma}^{\mathrm{M}}=r_{s}$, for each of the
branches. Furthermore, it respects the limit
$\mathcal{F}_{\gamma}\rightarrow\frac{r_{s}}{2}=\mathcal{F}_{\gamma}^{\mathrm{Sch}}$
for very large $r_{\gamma}$, indicating its value for the SBH. This way, and
as inferred from Fig. 6, the Helmholtz free energy is confined within the
domain
$\mathcal{F}_{\gamma}^{\mathrm{Sch}}\leq\mathcal{F}_{\gamma}\leq\mathcal{F}_{\gamma}^{\mathrm{M}}$.
Moreover, the Gibbs free energy can be calculated by means of the relation
$\mathcal{G}=r_{s}-\mathcal{T}\mathcal{S}-\Gamma r_{\gamma}$ Molina:2021 ,
which by means of Eqs. (11), (13) and (14), provides

$\mathcal{G}_{\gamma}=r_{s}-\frac{r_{\gamma}}{4}\left(1-\sqrt{1-\frac{4r_{s}}{r_{\gamma}}}\right),$
(36)

for constant values of $r_{\gamma}$. This quantity has a maximum of

$\mathcal{G}_{\gamma}^{\mathrm{M}}(r_{s}^{\mathrm{M}},r_{\gamma})\equiv\mathcal{G}_{\gamma}^{\mathrm{M}}=\frac{r_{s}^{\mathrm{M}}}{3},$
(37)

for $r_{s}^{\mathrm{M}}=\frac{3r_{\gamma}}{16}$. The behavior of
$\mathcal{G}_{\gamma}$ has been shown in Fig. 7. Either of the branches in the
diagram, possesses a maximum satisfying the relation
$\frac{r_{s}^{M}}{r_{\gamma}}=\frac{3}{16}$. Note that, as expected, for very
large values of $r_{\gamma}$, we have
$\mathcal{G}_{\gamma}\rightarrow\frac{r_{s}}{2}=\mathcal{G}_{\gamma}^{\mathrm{Sch}}$,
that corresponds to its value for the SBH. So, in the absence of quintessence,
both of the Helmholtz and Gibbs free energies will lead to the same value,
which is that of the SBH.

Figure 7: The behavior of $\mathcal{G}_{\gamma}$ plotted for
$r_{\gamma}=30,50,70,90$ and $110$. The dashed line connects the maximum
values $\mathcal{G}_{\gamma}^{\mathrm{M}}$, whereas the dotted-dashed line
corresponds to $\mathcal{G}_{\gamma}^{\mathrm{Sch}}$.

##  VI Conclusions

The thermodynamic features of black holes, beside being of interest on their
own as the physical tools for attributing sensible physical phenomena to
extremally gravitating systems, can also be important regarding modern
cosmology, when these systems are associated with some dark fields. Basically,
after redefining the causality of the spacetime by introducing new horizons,
such cosmological components also affect the thermodynamics of the black hole
through the extra terms they induce to the spacetime metric. In this paper, we
showed that the quintessential dark field surrounding a SBH, can be regarded
as a thermodynamic coordinate, and be used as a parameter together with the
black hole mass. Applying the Pfaffian form of the infinitesimal heat exchange
reversibly, we calculated the geometric entropy $\mathcal{S}$ and temperature
$\mathcal{T}$ in terms of these coordinates. Accordingly, we foliated the
$\mathcal{S}$-$\mathcal{T}$ curves for different black hole masses, and showed
that during the adiabatic processes, both of the thermodynamic coordinates
increase. This inference is crucial for this particular black hole, because we
saw further that the solutions to the Cauchy problem which are given in the
context of isoareal (adiabatic) processes, lead to several different paths on
the thermodynamic manifold, that not all of them are physically meaningful.
Accordingly, one needs to choose, among these paths, those that rely on the
increase of both of the parameters. We also showed the acceptable paths in
two-dimensional and three-dimensional plots. Furthermore, we used the formerly
derived thermodynamic quantities, to calculate the free energies associated
with the black hole horizon. We found that for constant values of the
cosmological component, the free energies are limited, and in the absence of
the cosmological term, both of them lead correctly to the same value for the
SBH. The discussion presented in this paper, can be generalized to the
rotating counterpart of the black hole which may be accounted for some
astrophysical support. This investigation is left for future studies.

## Acknowledgements

M. Fathi acknowledges Universidad de Santiago de Chile for financial support
through the Proyecto POSTDOCDICYT, Código 042331 CM-Postdoc. J.R. Villanueva
was partially supported by the Centro de Astrofísica de Valparaíso (CAV).

##  Appendix A Derivation of the solution to the Cauchy problem

Applying the change of variable $z\doteq\sqrt{\frac{x}{y}}$, the differential
equation (20) can be rewritten as

$\frac{1}{z^{2}}-\frac{2x}{z^{3}}\frac{\mathrm{d}z}{\mathrm{d}x}=\frac{1}{z(1-\sqrt{1-z})^{2}},$
(38)

which can be recast as

$\frac{\mathrm{d}w}{\mathrm{d}x}=\frac{1+w}{2x},$ (39)

using the second change of variable $w\doteq\sqrt{1-z}$. Considering the
initial states $(x_{0},w_{0})$, the above equation results in the solution

$\frac{\sqrt{x}}{1+w}=\bar{c}_{1},$ (40)

where $\bar{c}_{1}=\frac{\sqrt{x_{0}}}{1+w_{0}}$. After doing some algebraic
arrangements, and taking into account the applied changes of available, we
finally get to the solutions in Eqs. (24).

## References

  * (1) J. M. Bardeen, B. Carter, and S. Hawking, “The Four laws of black hole mechanics,” Commun. Math. Phys., vol. 31, pp. 161–170, 1973.
  * (2) J. D. Bekenstein, “Black holes and the second law,” Lett. Nuovo Cimento, vol. 4, pp. 737–740, Aug 1972.
  * (3) J. D. Bekenstein, “Black holes and entropy,” Phys. Rev. D, vol. 7, pp. 2333–2346, Apr 1973.
  * (4) J. D. Bekenstein, “Generalized second law of thermodynamics in black hole physics,” Phys. Rev. D, vol. 9, pp. 3292–3300, Jun 1974.
  * (5) J. D. Bekenstein, “Statistical black hole thermodynamics,” Phys. Rev. D, vol. 12, pp. 3077–3085, Nov 1975.
  * (6) S. Hawking, “Particle creation by black holes,” Commun. Math. Phys., vol. 43, pp. 199–220, 1975.  [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * (7) C. Teitelboim, “Action and entropy of extreme and nonextreme black holes,” Phys. Rev. D, vol. 51, p. 4315, 1995.  [Erratum: Phys.Rev.D 52, 6201 (1995)].
  * (8) S. M. Carroll, M. C. Johnson, and L. Randall, “Extremal limits and black hole entropy,” J. High Energy Phys., vol. 2009, pp. 109–109, Nov. 2009.
  * (9) W. G. Unruh, “Experimental Black-Hole Evaporation?,” Physical Review Letters, vol. 46, pp. 1351–1353, May 1981.
  * (10) M. Novello, M. Visser, and G. Volovik, Artificial Black Holes.  World Scientific, 2002.
  * (11) R. Schützhold and W. G. Unruh, “Hawking Radiation in an Electromagnetic Waveguide?,” Physical Review Letters, vol. 95, p. 031301, July 2005.
  * (12) I. Carusotto, S. Fagnocchi, A. Recati, R. Balbinot, and A. Fabbri, “Numerical observation of Hawking radiation from acoustic black holes in atomic Bose-Einstein condensates,” New Journal of Physics, vol. 10, p. 103001, Oct. 2008.
  * (13) F. Belgiorno, S. L. Cacciatori, M. Clerici, V. Gorini, G. Ortenzi, L. Rizzi, E. Rubino, V. G. Sala, and D. Faccio, “Hawking radiation from ultrashort laser pulse filaments,” Phys. Rev. Lett., vol. 105, p. 203901, Nov 2010\.
  * (14) S. Weinfurtner, E. W. Tedford, M. C. J. Penrice, W. G. Unruh, and G. A. Lawrence, “Measurement of Stimulated Hawking Emission in an Analogue System,” Phys. Rev. Lett., vol. 106, p. 021302, Jan. 2011.
  * (15) D. Castelvecchi, “Artificial black hole creates its own version of Hawking radiation,” Nature, vol. 536, pp. 258–259, Aug. 2016.
  * (16) J. Steinhauer, “Observation of quantum Hawking radiation and its entanglement in an analogue black hole,” Nature Physics, vol. 12, pp. 959–965, Oct. 2016.
  * (17) C. A. U. Lima, F. Brito, J. A. Hoyos, and D. A. T. Vanzella, “Probing the Unruh effect with an accelerated extended system,” Nature Communications, vol. 10, p. 3030, Dec. 2019.
  * (18) V. I. Kolobov, K. Golubkov, J. R. Muñoz de Nova, and J. Steinhauer, “Observation of stationary spontaneous Hawking radiation and the time evolution of an analogue black hole,” Nature Physics, vol. 17, pp. 362–367, Mar. 2021.
  * (19) K. Hutter, “The foundations of thermodynamics, its basic postulates and implications. A review of modern thermodynamics,” Acta Mechanica, vol. 27, pp. 1–54, Mar. 1977.
  * (20) A. Neumaier, “On the foundations of thermodynamics,” arXiv e-prints, p. arXiv:0705.3790, May 2007.
  * (21) C. Carathéodory, “Untersuchungen über die Grundlagen der Thermodynamik.,” Math. Ann., vol. 67, pp. 355–386, 1909.
  * (22) P. T. Landsberg, “Foundations of thermodynamics,” Rev. Mod. Phys., vol. 28, pp. 363–392, Oct 1956.
  * (23) I. E. Antoniou, “[No title found],” Foundations of Physics, vol. 32, no. 4, pp. 627–641, 2002.
  * (24) J. Gibbs, Thermodynamics, Vol. 1.  Yale University Press, New Haven, CT, 1948.
  * (25) M. Born, Natural Philosophy of Cause and Chance.  Clarendon Press, Oxford, 1949.
  * (26) H. Quevedo, “Geometrothermodynamics,” Journal of Mathematical Physics, vol. 48, no. 1, p. 013506, 2007.
  * (27) F. Belgiorno, “Homogeneity as a bridge between Carathéodory and Gibbs,” arXiv e-prints, pp. math–ph/0210011, Oct. 2002.
  * (28) F. Belgiorno, “Black hole thermodynamics in Caratheodory’s approach,” Phys. Lett. A, vol. 312, pp. 324–330, 2003.
  * (29) F. Belgiorno, “Quasihomogeneous thermodynamics and black holes,” J. Math. Phys., vol. 44, pp. 1089–1128, 2003.
  * (30) F. Belgiorno, “Notes on the third law of thermodynamics: I,” J. Phys. A: Math. Gen., vol. 36, pp. 8165–8193, Aug. 2003.
  * (31) F. Belgiorno, “Notes on the third law of thermodynamics: II,” J. Phys. A: Math. Gen., vol. 36, pp. 8195–8221, jul 2003.
  * (32) F. Belgiorno and M. Martellini, “Black holes and the third law of thermodynamics,” Int. J. Mod. Phys., vol. D13, pp. 739–770, 2004.
  * (33) M. Molina and J.R. Villanueva, “On the thermodynamics of the Hayward black hole,” Class. Quant. Grav., vol. 38, no. 10, p. 105002, 2021.
  * (34) M. Fathi, S. Lepe, and J. R. Villanueva, “Adiabatic analysis of the rotating BTZ black hole,” The European Physical Journal C, vol. 81, p. 499, June 2021.
  * (35) M. Fathi, M. Molina, and J. Villanueva, “Adiabatic evolution of hayward black hole,” Physics Letters B, vol. 820, p. 136548, 2021.
  * (36) J. A. Jimenez Madrid and P. F. Gonzalez-Diaz, “Evolution of a kerr-newman black hole in a dark energy universe,” Grav. Cosmol., vol. 14, pp. 213–225, 2008.
  * (37) M. Jamil, “Evolution of a Schwarzschild black hole in phantom-like Chaplygin gas cosmologies,” Eur. Phys. J. C, vol. 62, pp. 609–614, 2009.
  * (38) X.-Q. Li, B. Chen, and L.-l. Xing, “Charged Lovelock black holes in the presence of dark fluid with a nonlinear equation of state,” Eur. Phys. J. Plus, vol. 135, no. 2, p. 175, 2020.
  * (39) R. Roy and U. A. Yajnik, “Evolution of black hole shadow in the presence of ultralight bosons,” Physics Letters B, vol. 803, p. 135284, 2020.
  * (40) Z. Xu, X. Hou, X. Gong, and J. Wang, “Black hole space-time in dark matter halo,” Journal of Cosmology and Astroparticle Physics, vol. 2018, pp. 038–038, sep 2018.
  * (41) A. Das, A. Saha, and S. Gangopadhyay, “Investigation of circular geodesics in a rotating charged black hole in the presence of perfect fluid dark matter,” Classical and Quantum Gravity, vol. 38, p. 065015, feb 2021.
  * (42) V. V. Kiselev, “Quintessence and black holes,” Classical and Quantum Gravity, vol. 20, pp. 1187–1197, Mar. 2003.
  * (43) R. Saadati and F. Shojai, “Bending of light in a universe filled with quintessential dark energy,” Phys. Rev. D, vol. 100, p. 104041, Nov 2019\.
  * (44) I. Ali Khan, A. Sultan Khan, and S. Islam, “Dynamics of the particle around de sitter–schwarzschild black hole surrounded by quintessence,” International Journal of Modern Physics A, vol. 35, no. 23, p. 2050130, 2020\.
  * (45) V. H. Cárdenas, M. Fathi, M. Olivares, and J. R. Villanueva, “Probing the parameters of a Schwarzschild black hole surrounded by quintessence and cloud of strings through four standard astrophysical tests,” Eur. Phys. J. C, vol. 81, no. 10, p. 866, 2021.
  * (46) R. Tharanath and V. C. Kuriakose, “Thermodynamics and spectroscopy of Schwarzschild black hole surrounded by quintessence,” Modern Physics Letters A, vol. 28, p. 1350003, Feb. 2013.
  * (47) K. Ghaderi and B. Malakolkalami, “Thermodynamics of the schwarzschild and the reissner–nordström black holes with quintessence,” Nuclear Physics B, vol. 903, pp. 10–18, 2016.
  * (48) M.-S. Ma, R. Zhao, and Y.-Q. Ma, “Thermodynamic stability of black holes surrounded by quintessence,” General Relativity and Gravitation, vol. 49, p. 79, June 2017.
  * (49) S. G. Ghosh, S. D. Maharaj, D. Baboolal, and T.-H. Lee, “Lovelock black holes surrounded by quintessence,” The European Physical Journal C, vol. 78, p. 90, Feb. 2018.
  * (50) K. K. J. Rodrigue, M. Saleh, B. B. Thomas, and T. C. Kofane, “Thermodynamics phase transition and Hawking radiation of the Schwarzschild black hole with quintessence-like matter and a deficit solid angle,” General Relativity and Gravitation, vol. 50, p. 52, May 2018.
  * (51) C. H. Nam, “Higher dimensional charged black hole surrounded by quintessence in massive gravity,” General Relativity and Gravitation, vol. 52, p. 1, Jan. 2020.
  * (52) B. C. Lütfüoğlu, B. Hamil, and L. Dahbi, “Thermodynamics of Schwarzschild black hole surrounded by quintessence with generalized uncertainty principle,” The European Physical Journal Plus, vol. 136, p. 976, Sept. 2021.
  * (53) S. Chandrasekhar, An Introduction to the Study of Stellar Structure.  Astrophysical monographs, University of Chicago Press, 1939.
  * (54) H. A. Buchdahl, The Concept of Classical Thermodynamics.  Cambridge University Press, Cambridge, 1966.
  * (55) C. J. Adkins, Equilibrium Thermodynamics.  Cambridge University Press, 1983.
  * (56) P. Majumdar, “Black hole entropy and quantum gravity,” Indian J. Phys. B, vol. 73, p. 147, 1999.
  * (57) J. P. S. Lemos, M. Minamitsuji, and O. B. Zaslavskii, “Unified approach to the entropy of an extremal rotating BTZ black hole: Thin shells and horizon limits,” Phys. Rev. D, vol. 96, no. 8, p. 084068, 2017.


# Optimizing the optimizer for data driven deep neural networks and physics
informed neural networks

John Taylor 1,2,3,∗, Wenyi Wang 3, Biswajit Bala 3, Tomasz Bednarz 1

1CSIRO Data61, Canberra, Australian Capital Territory, Australia
2College of Engineering and Computer Science, The Australian National
University, Canberra, Australia
3Defence Science and Technology Group, Department of Defence, Canberra
Australia

*Corresponding Author Address: John Ashley Taylor, CSIRO Data61, CSIRO Synergy Building, Clunies Ross Street, Black Mountain, Canberra ACT 2601.
E-mail<EMAIL_ADDRESS>
## Abstract

We investigate the role of the optimizer in determining the quality of the
model fit for neural networks with a small to medium number of parameters. We
study the performance of Adam, an algorithm for first-order gradient-based
optimization that uses adaptive momentum, the Levenberg-Marquardt (LM)
algorithm a second order method, Broyden–Fletcher–Goldfarb–Shanno algorithm
(BFGS) a second order method and L-BFGS, a low memory version of BFGS. Using
these optimizers we fit the function $y=sinc(10x)$ using a neural network with
a few parameters. This function has a variable amplitude and a constant
frequency. We observe that the higher amplitude components of the function are
fitted first and the Adam, BFGS and L-BFGS struggle to fit the lower amplitude
components of the function. We also solve the Burgers equation using a physics
informed neural network(PINN) with the BFGS and LM optimizers. For our example
problems with a small to medium number of weights, we find that the LM
algorithm is able to rapidly converge to machine precision offering
significant benefits over other optimizers. We further investigated the Adam
optimizer with a range of models and found that Adam optimiser requires much
deeper models with large numbers of hidden units containing up to $\approx$26x
more parameters, in order to achieve a model fit close that achieved by the LM
optimizer. The LM optimizer results illustrate that it may be possible build
models with far fewer parameters. We have implemented all our methods in Keras
and TensorFlow 2.

##  1 Introduction

Neural networks have been successfully applied to a broad range of
applications and are now increasingly employed to gain new insights into
complex and demanding science problems [Lecun, Bengio, and HintonLecun et
al.2015]. Many science problems, in the context of deep learning, are
considered as regression problems where obtaining the best model fit is a key
measure of success. The best model fit is typically the result of applying an
optimisation method that minimises a loss function. Improvements in the model
fit are often a key criteria used to establish that a new model delivers
improved performance over prior modelling approaches.

Training neural networks usually continues until the optimizer achieves a
minimum value where the loss function no longer decreases or where a
comparison of the loss achieved on the training and test data indicates that
over-fitting is taking place. The optimizer may not succeed and cannot
guarantee that the global minimum has been achieved. The selection and
application of the optimizer therefore warrants careful consideration.

Neural networks are considered to be universal approximators [Hornik,
Stinchcombe, and WhiteHornik et al.1989]. While this property gives neural
networks incredible flexibility, as evidenced by their application to a very
wide range of problems, is not straightforward to determine an appropriate
model neural network architecture and the corresponding set of model weights
required to achieve a useful approximator. In this paper we examine the
challenge of finding the optimal weights for a given neural network model
architecture by comparing the performance of several well known optimisation
algorithms, Adam [Kingma and BaKingma and Ba2015], LM [LevenbergLevenberg1944]
[MarquardtMarquardt1963], BFGS [Nocedal and WrightNocedal and Wright2006] and
BFGS-L [Liu and NocedalLiu and Nocedal1989]. We apply these algorithms to a
well known but challenging problem where the frequency is fixed and the
amplitude varies. We also investigate the process of model optimisation by
studying the fit achieved in incremental steps.

Finally, the results presented here have significant practical implications.
This study was motivated by earlier studies where we examined applying machine
learning (ML) to detecting faults of fielded machinery using unsupervised
learning [Wang, Taylor, and ReesWang et al.2021a] [Wang, Taylor, and ReesWang
et al.2021b] where the function approximation and the choice of the optimizer
was a key task in delivering breakthrough results.

##  2 Methods

There are many quasi-Newton based optimisation algorithms available for
training neural networks. All algorithms are a compromise between accuracy,
robustness, computational efficiency, memory usage, scalability and
flexibility. We have selected four very well known optimisation algorithms
that have seen widespread application for our study, briefly summarized as
follows.

###  2.1 Adam

The Adam algorithm [Kingma and BaKingma and Ba2015] is a first-order gradient
based optimization of stochastic loss functions and has enjoyed widespread
application to fitting large and complex deep neural networks. The method is
robust, computationally efficient and works well on deep neural networks with
millions of parameters. Tuning of the learning rate, learning rate warmup,
learning rate schedules, momentum options, selection of regularizers and their
parameters is needed to achieve the best optimisation result. We use the Adam
optimiser algorithm as implemented in TensorFlow 2.

###  2.2 BFGS

The BFGS algorithm [Nocedal and WrightNocedal and Wright2006] is a second-
order iterative gradient method that computes the Hessian matrix of the loss
function using gradient evaluations, ${\displaystyle{\mathcal{O}}(n^{2})}$,
where $n$ is the number of parameters. Updates to the BFGS curvature matrix do
not need a matrix inversion which reduces the computational cost
significantly. As the Hessian matrix is the basis of the BFGS algorithm,
memory usage grows as the square of the number of parameters which results in
rapid growth in memory usage, making the approach unsuitable for neural
networks with large numbers of parameters. The selection of BFGS algorithm
parameters is straightforward. We use the BFGS optimiser as implemented in
TensorFlow Probability [Abadi, Agarwal, Barham, Brevdo, Chen, Citro, Corrado,
Davis, Dean, Devin, Ghemawat, Goodfellow, Harp, Irving, Isard, Jozefowicz,
Jia, Kaiser, Kudlur, Levenberg, Mané, Schuster, Monga, Moore, Murray, Olah,
Shlens, Steiner, Sutskever, Talwar, Tucker, Vanhoucke, Vasudevan, Viégas,
Vinyals, Warden, Wattenberg, Wicke, Yu, and ZhengAbadi et al.2015]
[TensorFlowTensorFlow2021].

###  2.3 L-BFGS

L-BFGS [Liu and NocedalLiu and Nocedal1989] addresses the problem of large
memory usage by the BFGS algorithm by saving only a few vectors that represent
an estimate of the full Hessian matrix. L-BFGS is computationally more
efficient than BFGS, uses less memory and can be applied to problems with
larger numbers of parameters than BFGS. We use the L-BFGS optimiser as
implemented in TensorFlow Probability [Abadi, Agarwal, Barham, Brevdo, Chen,
Citro, Corrado, Davis, Dean, Devin, Ghemawat, Goodfellow, Harp, Irving, Isard,
Jozefowicz, Jia, Kaiser, Kudlur, Levenberg, Mané, Schuster, Monga, Moore,
Murray, Olah, Shlens, Steiner, Sutskever, Talwar, Tucker, Vanhoucke,
Vasudevan, Viégas, Vinyals, Warden, Wattenberg, Wicke, Yu, and ZhengAbadi et
al.2015] [TensorFlowTensorFlow2021].

Both BFGS and L-BFGS can fall into local minima. By combining Adam as an
initial optimisation algorithm, run over many hundreds of epochs, then
continuing with BFGS and/or L-BFGS using the model weights delivered by Adam,
the risk of falling into local minima can be reduced. This approach has been
used in physics-Informed machine learning in order to solve systems of partial
differential equations(PDEs) [Karniadakis, Kevrekidis, Lu, Perdikaris, Wang,
and YangKarniadakis et al.2021], for example, and is also used in this study.

###  2.4 Levenberg-Marquardt

Levenberg [LevenbergLevenberg1944] was the original developer of this method
for non-linear least squares based on the finding that simple gradient descent
and Gauss-Newton iteration are complementary methods. Marquardt
[MarquardtMarquardt1963] extended the original method by scaling each
component of the gradient according to the curvature by including the diagonal
of the Hessian. Ananth [RanganathanRanganathan2004] provides an intuitive
summary of the method. The LM method does require a matrix inversion which is
computationally expensive and consumes significant amounts of memory, limiting
the application of LM to moderately sized models with thousands of parameters.
It is worth noting that modern computing hardware, with substantial
computational capability and memory, has expanded the size of models to which
LM can usefully be applied. As there is no built-in implementation of LM in
Tensorflow and Keras, we imported a user contributed code for Levenberg-
Marquardt (LM) algorithm by Fabio Di Marco [Fabio Di MarcoFabio Di Marco2021].

###  2.5 Test model

Our test model is the function $y=sinc(10x)$ over the domain [-1.5, 1.5]. We
build a neural network model consisting of an input layer, 2 Dense layers each
with 20 hidden units, and an output layer with a total of 481 trainable
parameters. We use the $tanh$ activation function in the two dense layers. In
all cases we minimize the mean square error loss. We generate a training data
set consisting of 20,000 evenly spaced points. We implemented the model in
Python using Keras [CholletChollet2015] and TensorFlow 2 [Abadi, Agarwal,
Barham, Brevdo, Chen, Citro, Corrado, Davis, Dean, Devin, Ghemawat,
Goodfellow, Harp, Irving, Isard, Jozefowicz, Jia, Kaiser, Kudlur, Levenberg,
Mané, Schuster, Monga, Moore, Murray, Olah, Shlens, Steiner, Sutskever,
Talwar, Tucker, Vanhoucke, Vasudevan, Viégas, Vinyals, Warden, Wattenberg,
Wicke, Yu, and ZhengAbadi et al.2015].

###  2.6 PINN model

Physics Informed Neural Networks (PINNs) are neural networks designed to solve
a variety of computational problems while accounting for the physical
equations which govern their respective natural phenomena. PINN is an emerging
deep learning technique that includes linear/nonlinear Ordinary Differential
Equations (ODEs), Partial Differential Equations (PDEs) into the loss function
of the NNs via automatic differentiation, in addition to any collected data,
essentially acting as an approximate differential equation solver which may
function well even with little to no training data. The PINN approach was
originally developed by [Raissi, Perdikaris, and KarniadakisRaissi et
al.2017], with two main applications in mind: data-driven solution (solving
differential equations given data/initial conditions/equation parameters) and
data-driven discovery of PDEs (find the parameters of the governing PDEs given
collected data), by using the fact that deep neural networks can act as
universal function approximators [Hornik, Stinchcombe, and WhiteHornik et
al.1989].

PINN’s initial design philosophy addressed the issue that for certain complex
physical systems the collection of data may be restrictive. PINNs shine as
they do not require huge data sets to produce a sufficiently good result. In
most machine learning techniques, more training data means more accurate and
robust results, and a sufficient amount of training data is needed for the
models to produce anything of quality. The creators of PINNs took advantage of
the abundance of domain specific information surrounding physical systems
which may not be included in standard neural network methods of modelling.
PINNs maximise the available knowledge of a problem by including already known
mathematical/physical descriptions of the problems, in addition to the
available training data, into the learning algorithms. Doing so may induce
faster convergence, and result in a much more generalisable solution when
compared to traditional machine learning methods, even with little to no
training data [Raissi, Perdikaris, and KarniadakisRaissi et al.2019].

This module implements a PINN model for Burgers’ equation [Raissi, Perdikaris,
and KarniadakisRaissi et al.2019]. Burgers’ equation is given by:

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}=\nu\frac{\partial^{2}u}{\partial x^{2}}$ (1)

where $\nu$ is the kinematic viscosity, $x$ is the spatial coordinate, $t$ is
the temporal coordinate, and $u(x,t)$ is the velocity of the fluid at point
$(x,t)$. This is a simplification of Navier-Stokes equation where the velocity
is in one spatial dimension and where external forces and pressure gradient
terms are ignored. We apply the initial conditions $u(x,t=0)=-sin(x)$ and the
boundary conditions $u(x=-1,+1,t)=0.$ The PINN model predicts $u(x,t)$ at the
point $(x,t)$. $u(x,t)$ represents a typical engineering problem where abrupt
changes are expected, which can be challenging for optimizers

We use an implementation of the solution for Burgers’ equation [Raissi,
Perdikaris, and KarniadakisRaissi et al.2019] written in TensorFlow provided
by [Okada39Okada392022]. We compare solutions of the PINN model obtained using
the L-BFGS-B and LM optimizers.

##  3 Results and Discussion

We first explored the performance of the Adam optimiser by allowing the model
to train for a fixed number of epochs, at which point training was paused and
the solution evaluated. The model training was then restarted using the model
weights obtained when training was paused. In this way we can inspect the
performance of the optimiser as training proceeds.

Figure 1 shows the fit achieved to the function $y=sinc(10x)$ of the Adam
optimizer after 250, 500, 750,1000, 1250 and 1500 epochs. As expected the
model fit improves as we increase the number of epochs. The Adam optimizer
fits the higher amplitude oscillations first, then fits each subsequent lower
amplitude oscillation in the function. Note that the frequency of the
oscillations is constant i.e. it is the amplitude only that varies. After 1500
epochs we found that the model fit does not improve and from inspection of the
graph we can see that the model has not achieved a highly accurate fit. Simply
relying on having obtained a minimum value for the MSE loss function cannot be
relied on to ensure that we have obtained a satisfactory fit. In the following
experiment we compare the fit obtained from multiple loss functions.

Figure 2 provides a comparison of the performance of the Adam optimizer, the
LM, BFGS and BFGS-L optimizers, where we fit the function $y=sinc(10x)$ using
the identical neural network, settings and input data for all training runs.
BFGS and L-BFGS use the weights generated by the Adam optimizer as their
starting point as suggested in previous studies [Raissi, Perdikaris, and
KarniadakisRaissi et al.2019]. All model runs show the results of the best fit
achieved for each optimiser. The Adam, L-BFGS and BFGS optimizers have
difficulty fitting the lower amplitude tails of the function with the LM model
achieving the best fit. For Adam, the minimum value of MSE loss function was
0.0002018958 and occurred at optimisation step 1,464, for BFGS-L, 8.655213e-05
at step 14,558, for BFGS, 4.2915726e-06 at step 12,477 and for LM,
1.642267e-07 at step 150, the maximum number of epochs set for LM model
training. The LM MSE value is $\approx$26x lower than that obtained by BFGS,
$\approx$527x lower than L-BFGS and $\approx$1,230x lower than Adam. The
compute time for each optimisation method, measured using the Python time
module, was for Adam 76s, BFGS 58s, L-BFGS 172 s and LM 116s, noting that the
total time for BFGS and L-BFGS optimisation includes Adam as the starting
point, results in a total time for BFGS of 134s and for L-BFGS of 248s.

Figure 3 shows MSE loss value, plotted on a linear scale, using the Adam
optimizer for each training epoch where we fit the function $y=sinc(10x)$
using a neural network. The upper panel plots the MSE loss over the entire
training run with the bottom panel focussing on the end of the model run. The
Adam optimiser exhibits significant variability between epochs before reaching
a minimum value after $\approx$1300 epochs.

Figure 4 and Figure 5 compare the MSE loss value, plotted on a linear scale,
using the L-BFGS and BFGS optimizers. Note that the L-BFGS and BFGS optimizers
were initialised with the set of weights from the conclusion of the Adam
optimzer run. Figure 4 shows that the L-BFGS optimiser takes many more steps
before converging at a higher value of the MSE than BFGS. This would be
expected given that L-BFGS does not calculate the full Hessian matrix. Where
parameter numbers are in the low to medium range BFGS is likely the better
choice of optimizer than L-BFGS. Only where memory usage or computational time
becomes too large would L-BFGS be recommended.

Figure 6 presents the MSE loss value, plotted on a linear scale, using the LM
optimiser. The LM optimizer reduces the MSE smoothly and in only a few epochs.
The MSE loss value rapidly falls below the Adam, L-BFGS and BFGS optimizer MSE
values. The benefits of the high quality model fit achieved by the LM
optimizer are clearly demonstrated in Figure 2 where, in contrast with the
other optimizers, the low amplitude (not frequency) components of the function
are well fitted.

Figure 7 compares the results obtained using all optimizers plotted on a log
scale. The log scale makes it possible to more easily compare optimiser
performance. Figure 7 clearly demonstrates the superior performance of the LM
method with a rapid minimization of the MSE to levels several orders of
magnitude lower than that achieved by other methods. The LM optimizer requires
no tuning to produce these results. While these results were obtained using a
simple test function, earlier studies where we examined applying machine
learning (ML) to detecting faults of fielded machinery using unsupervised
learning [Wang, Taylor, and ReesWang et al.2021a] [Wang, Taylor, and ReesWang
et al.2021b] [Wang, Taylor, and BalaWang et al.2021]with noisy data collected
from aircraft sensors, support the conclusion that the LM optimizer can
deliver superior performance in a real-world setting.

Having determined that the LM optimizer can deliver improved performance on
problems with a small to medium number number of weights we have also
considered the application of the LM optimizer to PINNs, an important new
class of NN applications. Many PINN models require a small to medium number of
weights and would benefit from the high accuracy solutions, often required of
science and engineering problems, that the LM model is able to deliver.

Figure 8 presents the results of fitting a PINN model, as described above, to
Burgers’ equation [Raissi, Perdikaris, and KarniadakisRaissi et al.2019] using
the L-BFGS-B optimizer. This version of the BFGS model extends L-BFGS to take
into account the box constraints needed to solve Burgers’ equation. The
results obtained are identical to that obtained previously [Raissi,
Perdikaris, and KarniadakisRaissi et al.2019]. Figure 9 presents the results
obtained using the LM optimizser and are identical to the results obtained
using L-BFGS-B. Finally, Figure 10 shows the MSE loss curve for the LM
optimiser over 50 epochs, with the MSE value approaching
$1\text{\times}{10}^{-7}$, several orders of magnitude lower than previously
reported MSE values, with the potential for further improvement. The
application of the LM model required no hyperparameter tuning.

The previous figures all show results with the model architecture fixed for
all optimizers. As the Adam optimizer is widely used when fitting deep
learning models, we have also investigated how to obtain the best fit to the
test function using only the Adam optimizer by modifying the model
architecture. As previously shown, the Adam optimizer delivered an MSE value
approximately 3 orders of magnitude higher than the LM optimizer. To achieve
an improved result for the Adam optimizer we investigated a range of hidden
units (HU) using 16,32,48,64 and 80 HU and varied the number of dense layers
in the model from 1-4. We use the same number of hidden units for each layer.
Using this basic grid search, we have a total of 20 model fit runs using a
combination of 5 hidden units and 4 layer combinations. We fit all models over
5000 epochs.

Figure 11 shows the results of this grid search organised as four panels with
an increasing number of model layers with each panel including the model fit
for all five HU values. Using only a single layer the model fit is very
similar across all HU values and does not improve on the MSE of the original
model. The maximum number of parameters using a 1 layer model with 80 HUs is
241 compared with the 481 parameters in the original test function with 2
layers and 20 HUs. With the 2 layer model all models achieve improved
performance compared with the single layer models. The model with 80 HUs and 2
layers (6,721 parameters) achieves the best result with the MSE value falling
below $1\text{\times}{10}^{-4}$, an improvement on the original model. The
model with 3 layers further improves on these results, however the model with
64 HUs (8,513 parameters), rather than 80 HUs (13,201 parameters), produces
the best results with the MSE value dropping below $1\text{\times}{10}^{-6}$.
With 4 layers the model is able to achieve good model results 48, 64, and 80
HUs with the best results at 48 HUs (7,201 parameters) and 64 HUs (12,673
parameters). The models with the larger number of parameters were approaching
the total number of training samples (20,000). Figure 12 shows the
corresponding model fit achieved using the models in Figure 11. The model fit
is consistent with the MSE values achieved during model fitting.

Figure 11 and 12 illustrate that the Adam optimizer performs best on this
problem primarily by increasing the number of layers and then by increasing
the number of HUs. Note that increasing the number of HUs to 80 did not
produce the best model for models with 3 and 4 layers. The best Adam optimizer
results required $\approx$15x more parameters with a 3 layer model and
$\approx$26x more parameters with 4 layer model as the LM optimized model.
While the Adam optimiser results can be improved, even with this much greater
number of model parameters, the Adam model cannot outperform the LM optimiser.
This result implies that Adam optimiser may require deeper models with large
numbers of hidden units in order to achieve the best model fit. The LM
optimiser results illustrate that it may be possible build models with far
fewer parameters.

##  4 Conclusions

We compared a range of well known optimizers applied to fitting neural
networks with a small to medium number of weights. We studied the performance
of Adam, the Levenberg-Marquardt (LM) algorithm, BFGS and L-BFGS and concluded
that the LM optimiser delivered significant advantages over the other methods,
including orders of magnitude improvement in the MSE values, rapid
optimisation and straightforward application as no hyperparameter tuning was
required.

Using these Optimizers we fit the function $y=sinc(10x)$ using a neural
network with a few parameters. This function has a variable amplitude and a
constant frequency. We observed that as the model fit progressed, the higher
amplitude components of the function were fit first. We found that only the LM
optimiser could fit this function well across the full range of amplitudes.

We have also demonstrated the usefulness of the LM optimiser to PINNs by
solving the Burgers equation to a higher accuracy than has previously been
achieved. Combined with the observation that the LM method can achieve an
accurate fit where other optimizers cannot, the LM method is likely to be very
useful for PINNs and to the broader class of NN models with a small to medium
number of parameters.

It would likely benefit the broader ML community if the LM and BFGS optimizers
included in could be made available without additional wrapping in TensorFlow.
This is recommended as there are a large number of potential applications that
would benefit significantly from the ready availability of these optimizers.

## References

  * [Abadi, Agarwal, Barham, Brevdo, Chen, Citro, Corrado, Davis, Dean, Devin, Ghemawat, Goodfellow, Harp, Irving, Isard, Jozefowicz, Jia, Kaiser, Kudlur, Levenberg, Mané, Schuster, Monga, Moore, Murray, Olah, Shlens, Steiner, Sutskever, Talwar, Tucker, Vanhoucke, Vasudevan, Viégas, Vinyals, Warden, Wattenberg, Wicke, Yu, and ZhengAbadi et al.2015] Abadi, M., A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, R. Jozefowicz, Y. Jia, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, M. Schuster, R. Monga, S. Moore, D. Murray, C. Olah, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng (2015).  TensorFlow: Large-scale machine learning on heterogeneous systems.
  * [CholletChollet2015] Chollet, F. (2015).  Keras.
  * [Fabio Di MarcoFabio Di Marco2021] Fabio Di Marco (2021).  Setup Levenberg-Marquardt optimization algorithm on TensorFlow.
  * [Hornik, Stinchcombe, and WhiteHornik et al.1989] Hornik, K., M. Stinchcombe, and H. White (1989).  Multilayer feedforward networks are universal approximators.  Neural Networks 2(5), 359–366.
  * [Karniadakis, Kevrekidis, Lu, Perdikaris, Wang, and YangKarniadakis et al.2021] Karniadakis, G. E., I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang (2021).  Physics-informed machine learning.  Nature Reviews Physics 3(6), 422–440.
  * [Kingma and BaKingma and Ba2015] Kingma, D. P. and J. L. Ba (2015).  Adam: A method for stochastic optimization.  In 3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings. ArxIV.
  * [Lecun, Bengio, and HintonLecun et al.2015] Lecun, Y., Y. Bengio, and G. Hinton (2015).  Deep learning.
  * [LevenbergLevenberg1944] Levenberg, K. (1944).  A method for the solution of certain non-linear problems in least squares.  Quarterly of Applied Mathematics 2(2), 164–168.
  * [Liu and NocedalLiu and Nocedal1989] Liu, D. C. and J. Nocedal (1989).  On the limited memory BFGS method for large scale optimization.  Mathematical Programming 45(1-3), 503–528.
  * [MarquardtMarquardt1963] Marquardt, D. W. (1963).  An Algorithm for Least-Squares Estimation of Nonlinear Parameters.  Journal of the Society for Industrial and Applied Mathematics 11(2), 431–441.
  * [Nocedal and WrightNocedal and Wright2006] Nocedal, J. and S. J. Wright (2006).  Numerical optimization.  In Springer Series in Operations Research and Financial Engineering, pp. 1–664.
  * [Okada39Okada392022] Okada39 (2022).  Physics Informed Neural Network (PINN) for Burgers’ equation.
  * [Raissi, Perdikaris, and KarniadakisRaissi et al.2019] Raissi, M., P. Perdikaris, and G. Karniadakis (2019, feb).  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics 378, 686–707.
  * [Raissi, Perdikaris, and KarniadakisRaissi et al.2017] Raissi, M., P. Perdikaris, and G. E. Karniadakis (2017).  Physics Informed Deep Learning (Part II): Data-driven Discovery of Nonlinear Partial Differential Equations.  arXiv preprint arXiv:1711.10561.
  * [RanganathanRanganathan2004] Ranganathan, A. (2004).  The Levenberg-Marquardt Algorithm.  Internet httpexcelsior cs ucsb educoursescs290ipdfL MA pdf 142(June).
  * [TensorFlowTensorFlow2021] TensorFlow (2021).  TensorFlow Probability.
  * [Wang, Taylor, and BalaWang et al.2021] Wang, W., J. Taylor, and B. Bala (2021).  Exploiting the Power of Levenberg-Marquardt Optimizer with Anomaly Detection in Time Series.
  * [Wang, Taylor, and ReesWang et al.2021a] Wang, W., J. Taylor, and R. J. Rees (2021a).  Recent Advancement of Deep Learning Applications to Machine Condition Monitoring Part 1: A Critical Review.
  * [Wang, Taylor, and ReesWang et al.2021b] Wang, W., J. Taylor, and R. J. Rees (2021b).  Recent Advancement of Deep Learning Applications to Machine Condition Monitoring Part 2: Supplement Views and a Case Study.  Acoustics Australia 49(2), 221–228.

## Figures

Figure 1: The fit achieved to the function $y=sinc(10x)$ of the Adam optimizer
after 250, 500, 750, 1000, 1250 and 1500 epochs. As expected the model fit
improves as we increase the number of epochs. The Adam optimizer fits the
higher amplitude oscillations first, then fits each subsequent lower amplitude
oscillation in the function. Note that the frequency of the oscillations is
constant and after 1500 epochs the model fit does not improve. Figure 2: A
comparison of the performance of the Adam optimizer, an algorithm for first-
order gradient-based optimisation that uses adaptive momentum, the Levenberg-
Marquardt (LM) algorithm a second order method, BFGS a second order method and
BFGS-L, a low memory version of BFGS where we fit the function $y=sinc(10x)$
using the same neural network. All model runs show the results of the best fit
achieved. The Adam, L-BFGS and BFGS Optimizers have difficulty fitting the
lower amplitude tails of the function. Figure 3:  The MSE loss value, plotted
on a linear scale, using the Adam optimizer for each training epoch where we
fit the function $y=sinc(10x)$ using a neural network. The upper panel plots
the MSE loss over the entire training run with the bottom panel focussing on
the end of the model run. Note that the scale of the two plots is not the
same.  Figure 4: The MSE loss value, plotted on a linear scale, using the
L-BFGS optimizer for each training epoch where we fit the function
$y=sinc(10x)$ using a neural network. The upper panel plots the MSE loss over
the entire training run with the bottom panel focussing on the end of the
model run. Note that the scale of the two plots is not the same. Figure 5:
The MSE loss value, plotted on a linear scale, using the BFGS optimizer for
each training epoch where we fit the function $y=sinc(10x)$ using a neural
network. The upper panel plots the MSE loss over the entire training run with
the bottom panel focussing on the end of the model run. Note that the scale of
the two plots is not the same.  Figure 6:  The MSE loss value, plotted on a
linear scale, using the LM optimizer for each training epoch where we fit the
function $y=sinc(10x)$ using a neural network. The upper panel plots the MSE
loss over the entire training run with the bottom panel focussing on the end
of the model run. Note that the scale of the two plots is not the same. Figure
7: A comparison of the performance of the Adam optimizer, an algorithm for
first-order gradient-based optimisation that uses adaptive momentum, the
Levenberg-Marquardt (LM) algorithm a second order method, BFGS a second order
method and BFGS-L, a low memory version of BFGS, all plotted on a log scale.
The BFGS and L-BFGS algorithms commence fitting using the final state of the
Adam optimizer after the Adam optimizer has completed 1500 time steps. The
Levenberg-Marquardt (LM) algorithm achieves an improved MSE loss value
compared to the Adam optimizer after 10 epochs and the BFGS algorithm after 20
epochs. Figure 8: Solution of the Burgers equation solved using a PINN model
using the L-BFGS-B optimizer. Figure 9: Solution of the Burgers equation
solved using a PINN model using the Lavenburg-Marquardt optimizer. The
solution matches that obtained using the L-BFGS-B optimizer. Figure 10: The
MSE loss as a function of epoch using the LM optimizer for the Burgers
equation solved using a PINN model. The upper panel shows the loss on a linear
scale with the lower panel showing the loss on a log scale. The final loss
value is lower than that obtained in previous studies.

Figure 11: The MSE loss as a function of epoch using the Adam optimizer for
the test function over 5000 epochs. Each panel shows the results with a model
with 1, 2, 3 or 4 layers with the number of hidden units in all layers of the
model equal to 16, 32, 48, 64 or 80. Note that some curves are hidden beneath
other results.

.

Figure 12: The model fit achieved using the Adam optimiser using a model with
1, 2, 3 or 4 layers with the number of hidden units in all layers of the model
equal to 16, 32, 48, 64 or 80. Note that some curves are hidden beneath other
results.