Modeling the Complexity of Elliptic Black Hole Solution In 4D Using Hamiltonian Monte Carlo with Stacked Neural Networks

1. Introduction

Black holes are known by their mass, angular momentum and their charge. Choptuik [1] illustrates that there may be another parameter explaining the gravitational collapse solutions. Christodolou [2] elaborated on the spherically symmetric collapse of the scalar field. Choptuik [1] empirically showed the critical behaviour of the solutions that illustrated the discrete self-similarity for the real scalar field gravitational collapse. The gravitational collapse solution shows space-time self-similarity so that dilations occur. The critical solution provides the scaling law. We can illustrate the initial condition of the real scalar field by parameter p, which is related to the field amplitude. Let $p=p_{\mathrm{crit}}$ define the critical solution. The black hole is then formed when p gets bigger than $p_{\mathrm{crit}}$. If $p>p_{\mathrm{crit}}$, then one finds out the mass of the black hole or the Schwarzschild radius by a scaling law as follows

$$r_S\left(p\right)\propto M_{\mathrm{bh}}\left(p\right)\propto {(p-p_{\mathrm{crit}})}^{\gamma }.$$

(1)

The critical exponent for a real scalar field in four dimensions is found to be $\gamma \simeq 0.37$ [1,3,4], while for the other dimensions ($d\ge 4$), the mass of black hole is given by [5,6] as

$$r_S\left(p\right)\propto {(p-p_{\mathrm{crit}})}^{\gamma },M_{\mathrm{bh}}\left(p\right)\sim {(p-p_{\mathrm{crit}})}^{(D-3)\gamma }.$$

(2)

The numerical optimization for different matter content can be read from [7,8,9,10,11,12]. The collapse solutions of the perfect fluid had been studied in [5,13,14,15] and its critical exponent $\gamma \simeq 0.36$ was also found in [14]. It has also been discussed in [16] that $\gamma$ may have a universal value for all the fields that are coupled to gravity in four dimensions. One method for obtaining the critical exponent is based on perturbations of self-similar solutions. This method was established in [5,15,17]. Let h be any field, its perturbation is given by the following

$$h=h_0+\epsilon h_{-\kappa },$$

(3)

where $h_{-\kappa }$ has the scaling $-\kappa \in \mathbb{C}$ which is related to different modes. The most relevant mode ${\kappa }^{*}$ is related to the biggest value of $\Re \left(\kappa \right)$. The minus sign demonstrates a growing mode near the black hole formation time $t\to 0$. It was proved in [5,15,17] that ${\kappa }^{*}$ can be related to Choptuik exponent by the following relation

$$\gamma =\frac{1}{\Re {\kappa }^{*}}.$$

(4)

The axial symmetry was studied in [18]. On the other hand, Alvarez-Gaume et al. [19] explore other solutions as well as shock waves. The critical exponent $\gamma \simeq 0.2641$ for the axion-dilaton system in four dimensions was first discovered in [20]. The analysis of [3,21] had been done in four dimensions for the elliptic case; however, their techniques can be used in hyperbolic and parabolic cases as well as in other dimensions, for further results see [22,23].

As the first motivation for the critical collapse of the axion-dilaton system, one can consider the AdS/CFT correspondence [24,25,26,27], relating the Choptuik exponent, the imaginary part of quasinormal modes and also the dual conformal field theory [28]. In the context of AdS/CFT, it is desirable to consider the spaces which tend asymptotically to $Ad{S}_5\times S^5$. From the bosonic fields of the theory, one can employ the system of axion-dilaton and the self-dual 5-form to analyze the critical collapse solutions in different dimensions. The other motivation of the axion-dilaton system is actually related to the holographic description of black hole formation [13,29]. The next motivation could potentially be the application of the system to black hole physics [30,31] as well as the important role of S-duality in self similar solutions [32].

The families of continuous self-similar solutions were explored by [33] in four and five dimensions for all classes of SL(2,R). These solutions are the extended results of [34,35]. More specifically, [36] studied the perturbations and reproduced the existing value [20] of $\gamma \sim 0.2641$ in four dimensions. Some other critical exponents for four and five dimensions were proved in [23]. Hatefi and Hatefi [37] proposed the Fourier-based regression models for the nonlinear critical collapse functions. To address the challenges of [37], Hatefi and Hatefi [38] employed truncated power basis, natural spline and penalized B-spline regression models to estimate the non-linear critical collapse functions. Recently, Hatefi et al. [39] applied artificial neural networks to analyze the black hole solutions in parabolic class in higher dimensions in detail.

In this paper, we propose a new methodology, which allows us to approach the problem of the complexity of elliptical black holes in 4 dimensions by means of Bayesian statistical modeling. Technically, we address the problem of modeling the self-similar solutions for the described spherical gravitational collapse [3,21] in four dimensions, using Hamiltonian Monte Carlo with Stacked Neural Networks. The modeling of black holes in the 4-dimensional elliptic case is reduced to a set of differential equations, which have always been approached by numerical methods that try to find a unique solution that satisfies the initial conditions. On the contrary, we argue that it is possible to attack the problem with Bayesian statistical modeling of the differential equations. In order to simplify the equations and the parameters of the equations of motion, researchers often have to use a variety of numerical approaches due to the complicated and highly nonlinear form of the equations of motion in the physics of black holes (e.g. [22,38,39]). For this reason, while calculating the parameters of the equations of motion, they must disregard the measurement errors imposed by the numerical calculations.

To the best of our knowledge, for the first time in the literature on the axion-dilation black hole, we investigate the equations of motion of the critical collapse functions within a Bayesian framework. This Bayesian framework enables us to address the measurement errors imposed by the numerical calculations in the estimation process. The parameters of the equations of motion are handled as random variables due to this statistical modelling. Because of this, we can derive the posterior distributions for these parameters, which we can then utilize as legitimate priors for a collection of stacked neural network-based solvers that will calculate the critical collapse function. In a nutshell, this paper describes a novel pipeline that starts solving the elliptic black hole equations in 4 dimensions using Hamiltonian Monte Carlo statistical modelling. Then, the posterior distributions discovered, for the parameters of the equations, are used with stacked neural networks to estimate the critical collapse functions. Unlike methods in the literature, this probabilistic perspective enables us to recover the deterministic solution in the literature and find all possible solutions that may occur due to numerical measurement errors in the parameter estimation phase. This paper paves the path to view the problem from a stochastic rather than a deterministic perspective.

This paper is organized as follows. Section 2 discusses the effective action of the axion-dilaton system in four dimensions, its equations of motion, and the initial conditions for the axion-dilaton system that come from the requirement of continuous self-similarity. Section 3 describes our novel approach for modelling the complexity of elliptic black hole solution in 4d from Hamiltonian Monte Carlo with stacked neural networks. Through extensive numerical studies, we develop the Bayesian posterior mean and median neural networks and the 95% credible intervals for all the critical functions in the elliptic class 4d in Section 4. Finally, in Section 5, we present the summary and concluding remarks.

2. The Relevant System

The axion-dilaton fields are combined into a single complex scalar field $\tau \equiv a+i{e}^{-\mathbf{OE}}$. The four-dimensional axion-dilaton ($a,\Phi$) action is given by

$$S=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}\left(R-\frac{1}{2}\frac{{\partial }_a\tau {\partial }^a\overline{\tau }}{{\left(\mathrm{Im}\tau \right)}^2}\right),$$

(5)

where R is the scalar curvature. By taking variations from the metric and $\tau$, one finds all the equations of motion as follows

$$R_{ab}=\frac{1}{4{\left(\mathrm{Im}\tau \right)}^2}({\partial }_a\tau {\partial }_b\overline{\tau }+{\partial }_a\overline{\tau }{\partial }_b\tau ),$$

(6)

$${\nabla }^a{\nabla }_a\tau +\frac{i{\nabla }^a\tau {\nabla }_a\tau }{\mathrm{Im}\tau }=0.$$

(7)

Note that the theory is classically invariant under SL(2,R) transformations, that is, if $\tau$ gets exchanged by

$$\tau \to \frac{a\tau +b}{c\tau +d},$$

(8)

where $(a,b,c,d)\in R$, $ad-bc=1$ and $g_{ab}$ as well as the action remain to be invariant. This group can be broken to SL(2,Z) as a consequence of duality in String Theory (for review of duality see for instance [40,41,42,43]).

The spherically symmetric metric is given by [21] as

$$d{s}^2=\left(1+u(t,r)\right)\left(-b{(t,r)}^{2}d{t}^2+d{r}^2\right)+r^{2}d{\Omega }_{d-2}^2.$$

(9)

One can implement time re-definitions for (9), hence as in [21] we set $b(t,0)=1$ and regularity condition implies $u(t,0)=0$. The continuous self-similarity means the existence of a killing vector $\xi$ that produces global scale transformation as follows

$${\mathcal{L}}_{\xi }{g}_{ab}=2g_{ab}.$$

(10)

In spherical coordinates, one chooses $\xi =t\partial /\partial t+r\partial /\partial r$. Now by defining the scale invariant variable $z=-r/t$, self-similarity of the metric reflects the fact that all the functions $u(t,r),b(t,r)$ can be expressed in terms of z.

The action in (5) is SL(2,R)-invariant, so one compensates a scale transformation of $(t,r)$ by an SL(2,R) transformation. Hence, by considering the change of variables to $(t,z)$, one derives a differential condition (as originally pointed out in [34]) for $\tau (t,z)$ as

$$t\frac{\partial }{\partial t}\tau (t,z)=a+b\tau +c{\tau }^2,$$

(11)

where all $a,b,c$ are real numbers. The above equation has two roots which can be either two complex conjugate numbers, two real numbers or a double real root where they are related to compensating the scaling transformation in SL(2,R). In the elliptic class, the general form of the ansatz is given by

$$\tau (t,r)=i\frac{1-{(-t)}^{i\omega }f\left(z\right)}{1+{(-t)}^{i\omega }f\left(z\right)},$$

(12)

where under a scaling transformation $t\to \lambda t$, $\tau (t,r)$ changes by a SL(2,R) rotation, which means that all equations are invariant under the following transformation

$$f\left(z\right)\to e^{i\theta }f\left(z\right),$$

(13)

where $f\left(z\right)$ is a complex function that satisfies |f(z)| < 1 to be determined by solutions to non-linear ordinary differential equations, and $\omega$ is a real constant.

In the following, we describe briefly the derivation of the equations of motion for the elliptic class in four dimensions. Applying continuous self-similarity ansätze (12) to all the equations of motion (6) and (7), one finds the ordinary differential equations for $u\left(z\right)$, $b\left(z\right)$, $f\left(z\right)$. By taking the spherical symmetry, $u\left(z\right)$ can be expressed in terms of $b\left(z\right)$ and $b^{\prime }\left(z\right)$ as

$$u\left(z\right)=\frac{z{b}^{\prime }\left(z\right)}{(d-3)b\left(z\right)}.$$

(14)

Note that the first derivative of $u\left(z\right)$ can be eliminated from all equations of motion by the following constraint

$$\frac{q{u}^{\prime }\left(z\right)}{2(1+u(z\left)\right)}=\frac{4z{f}^{\prime }\left(z\right){\overline{f}}^{\prime }\left(z\right)+2i(w\overline{f}\left(z\right)f^{\prime }\left(z\right)-wf\left(z\right){\overline{f}}^{\prime }\left(z\right))}{2{(f\left(z\right)\overline{f}\left(z\right)-1)}^2}.$$

(15)

Hence, the equations of motion in the elliptic class in four dimensions are then given by

$$\begin{array}{ccc}\hfill 0& =& b^{\prime }+\frac{z(b^2-z^2)}{{b(-1+|f|}^2{)}^2}{f}^{\prime }{\overline{f}}^{\prime }-\frac{i\omega (b^2-z^2)}{{b(-1+|f|}^2{)}^2}(f{\overline{f}}^{\prime }-\overline{f}{f}^{\prime })-\frac{{\omega }^{2}z{\left|f\right|}^2}{{b(-1+|f|}^2{)}^2},\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill 0& =& f^{\prime \prime }-\frac{z(b^2+z^2)}{b^2{(-1+|f|}^2{)}^2}{f}^{\prime 2}{\overline{f}}^{\prime }+\frac{2}{(1-|f{|}^2)}\left(1-\frac{i\omega (b^2+z^2)}{2b^2{(1-|f|}^2)}\right)\overline{f}{f}^{\prime 2}\hfill \\ & & +\frac{i\omega (b^2+2z^2)}{b^2{(-1+|f|}^2{)}^2}f{f}^{\prime }{\overline{f}}^{\prime }+\frac{2}{z}\left(1+\frac{i\omega z^2{(1+|f|}^2)}{(b^2-z^2){(1-|f|}^2)}\right.\hfill \\ & & +\left.\frac{{\omega }^2{z}^4{\left|f\right|}^2}{b^2(b^2-z^2){(1-|f|}^2{)}^2}\right)f^{\prime }+\frac{{\omega }^{2}z}{b^2{(-1+|f|}^2{)}^2}{f}^2{\overline{f}}^{\prime }+\hfill \\ & & \frac{2i\omega }{(b^2-z^2)}\left(\frac{1}{2}-\frac{i\omega (1+|f{|}^2)}{2(1-|f{|}^2)}\right.-\left.\frac{{\omega }^2{z}^2{\left|f\right|}^2}{2b^2{(-1+|f|}^2{)}^2}\right)f.\hfill \end{array}$$

(17)

The first equation of motion includes $b\left(z\right),f\left(z\right),f^{\prime }\left(z\right)$ where b is a first-order linear in-homogeneous equation where its initial condition is given by $b\left(0\right)=1$. The initial conditions for $f\left(z\right),f^{\prime }\left(z\right)$ are known by applying the smoothness of the critical solution. Applying polar coordinate representation $f\left(z\right)=f_m\left(z\right)e^{i{f}_a\left(z\right)}$, one can easily show that all equations are invariant under a global redefinition of the phase of $f\left(z\right)$. This implies that $f_a\left(0\right)=0$. From the regularity condition at $z=0$ and the residual symmetries in the equations of motions (17), one can derive the initial boundary conditions of the equations of motion by

$$\begin{array}{c}\hfill b\left(0\right)=1,f_m\left(0\right)=x_0,f_m^{\prime }\left(0\right)=f_a^{\prime }\left(0\right)=f_a\left(0\right)=0\end{array}$$

(18)

3. Statistical Methods

Let $\mathbf{x}\left(t\right)=(x_1\left(t\right),\dots ,x_H\left(t\right))$ denote the solutions to a system of H differential equations (DEs), where $x_i\left(t\right)$ represents the solution to the i-th equation; henceforth, $x_i\left(t\right)$ is called the DE variable. Combining all the DE variables, we can model the system by

$$\begin{array}{c}\hfill \frac{d}{dt}{x}_i\left(t\right)=g_i\left(\mathbf{x}\left(t\right)\right|\theta ),\end{array}$$

(19)

where t denotes the space-time argument, $\theta$ represents the vector of all unknown parameters in the DE system with $x_i\left(0\right)=x_{i0}$ being the initial condition for the i-th DE variable.

3.1. The Likelihood Function

According to the complex and nonlinear pattern of the equations of motion, the true trajectory of the DE variables is not observable. Instead, one obtains the trajectory of the DE variables with measurement error implied by numerical experiments. To take this uncertainty into our statistical models, let ${\mathbf{y}}_i=(y_{i1},\dots ,y_{i{n}_i})$ denote the observed trajectory of DE variable $x_i\left(t\right)$ at $n_i$ points where $y_{ij}$ represents the observed value of $x_i\left(t\right)$ with measurement error at space-time point $j,j=1,\dots ,n_i$. To take into account the errors involved in measuring the DE variables, we assume that $y_{ij}$ follows Gaussian distribution with mean $x_i\left(t_j\right|\theta )$ and variance ${\sigma }_i^2$ for $i=1,\dots ,H$; that is

$$y_{ij}\sim N\left(x_i\left(t_j\right|\theta ),{\sigma }_i^2\right),j=1,\dots ,n_i,$$

where $x_i\left(t_j\right|\theta )$ denotes the true value of the DE variable given the unknown parameters $\theta$ and $t_j$ denotes the space-time point where we observed $y_{ij}$ as the j-th value of the DE variable $x_i\left(t_j\right|\theta )$. Let $\Omega =(\theta ,\sigma )$ be the set of all unknown parameters of the Gaussian process, where $\sigma =({\sigma }_1,\dots ,{\sigma }_H)$. Given $\mathbf{y}=({\mathbf{y}}_1,\dots ,{\mathbf{y}}_H)$, the likelihood function of $\Omega$ is given by

$$\begin{array}{c}\hfill L(\Omega |\mathbf{y})=\prod _{i=1}^H\prod _{j=1}^{n_i}{\left(\frac{1}{{\sigma }_i^2}\right)}^{-1/2}exp\left\{-\frac{{(y_{ij}-x_i\left(t_j\right|\theta ))}^2}{2{\sigma }_i^2}\right\}.\end{array}$$

(20)

The likelihood function (20) translates the DE system into a Gaussian process. Thus, one can estimate the unknown parameter $\Omega$ by maximizing the likelihood function (20) and obtain the maximum likelihood estimates of $\Omega$. According to the complex and nonlinear forms of the equations of motion, the likelihood function surface may lead to multiple isolated optimum points. Therefore, the likelihood function may appear very sensitive regarding parameter changes.

3.2. Hamiltonian Monte Carlo

We develop a Bayesian framework to estimate the unknown parameters of DE system (19). One must incorporate prior distribution for the model parameters $\Omega$ in the Bayesian framework. Suppose $\pi (\Omega ;\alpha )$ denotes the prior distribution for $\Omega$, where $\alpha$ denotes the vector of all hyper-parameters of the prior distribution. From (20) and the prior distribution ${\pi }_0(\Omega ;\alpha )$, the posterior distribution $\Omega$ is given by

$$\pi (\Omega |\mathbf{y})=\frac{L(\Omega |\mathbf{y}){\pi }_0(\Omega ;\alpha )}{{\int }_{\Omega }d\Omega L(\Omega |\mathbf{y}){\pi }_0(\Omega ;\alpha )}.$$

Due to the complex uncertainty involved in the multidimensional integrals of the marginal distributions, there is no closed form for the posterior distribution.

Markov Chain Monte Carlo (MCMC) is a well-established method in the Bayesian framework to make inferences about the posterior distribution of differential equations [44]. The MCMC approach translates the estimation problem into the sampling from the target non-normalized posterior distribution $\pi (\Omega |\mathbf{y})\propto L(\Omega |\mathbf{y}){\pi }_0(\Omega ;\alpha )$. There is a rich literature on the properties of the MCMC methods [45,46]. MCMC employs various stochastic steps to construct a Markov chain from the target posterior distribution numerically. If the chain is run long enough, it is guaranteed that the MCMC converges to the true posterior distribution. It reaches a stationary state where the sampler enables us to generate samples from the posterior distribution.

Metropolis-Hastings (MH) is one of the most common MCMC approaches [45,47]. The MH uses a transition distribution $q\left({\Omega }^{*}\right|\Omega )$ computing the probability of making transition from $\Omega$ to ${\Omega }^{*}$ in the domain of the parameter space. MH then employs a stochastic step to accept or reject the transition. Let ${\Omega }^{\left(t\right)}$ denote the t-th state of the MCMC. The probability of accepting the transition to candidate state ${\Omega }^{*}$ is then calculated by

$$\begin{array}{c}\hfill min\left\{1,\frac{\pi \left({\Omega }^{*}\right|\mathbf{y})q(\Omega |{\Omega }^{*})}{\pi (\Omega |\mathbf{y})q\left({\Omega }^{*}\right|\Omega )}\right\}.\end{array}$$

(21)

Random Walk (RW) is a common MH approach using a guess-and-check strategy to generate proposals from $\pi (\Omega |\mathbf{y})$ in a neighbourhood of the current state of the MCMC. The RW is a popular approach because it is simple and intuitive to generate samples from a neighbourhood of the current state of the chain. Despite these advantages, the RW may appear very inefficient in high dimensions and complex systems where the posterior distribution may have multiple isolated modes. If the RW is tuned to take high random jumps, the MCMC chain may randomly jump from one low-density area in the parameter space to another. Consequently, this leads to a high probability of rejection and thus results in a biased estimate. On the other side, if RW is restricted to small jumps, the MCMC keeps accepting all the proposals and will be unable to visit other high-density posterior areas. The MCMC again leads to a biased estimate.

Hamiltonian Monte Carlo (HMC) can remedy the above challenges and find the posterior distribution in the DE systems [48,49]. The HMC augments the MH sampler with a momentum variable $\mathbf{OE}$ to make the transition much more rapidly across the parameter space. This enables the MCMC to sample from almost all high-density areas and their neighbourhoods. Consequently, the MCMC estimates more effectively the posterior distribution. The HMC, as an iterative technique, consists of three steps in each iteration. Let $({\Omega }^{(t-1)},{\mathbf{OE}}^{(t-1)})$ be the current state of the MCMC in the t-th iteration. Thus, the $\left(t\right)$-th iteration of the HMC is carried out as follows:

• Momentum Initialization Step: We initialize the candidate ${\mathbf{OE}}^{*}$ from the prior distribution $p\left(\mathbf{OE}\right)$.
• Leapfrog Step: This is the main step of the HMC consisting of three steps to upload simultaneously the DE parameters $\Omega$ and momentum parameters $\mathbf{OE}$. The leapfrog step is then alternated for L iterations where each iteration involves the (i)-(iii) steps:

  (i)

  We update the half-step of the momentum parameter $\mathbf{OE}$ using the gradient of the log-posterior distribution $\pi (\Omega |\mathbf{y})$ scaled by tuning factor $ϵ$ as

  $$\begin{array}{c}\hfill {\mathbf{OE}}^{*}=\mathbf{OE}+\frac{1}{2}ϵ\frac{d}{d\Omega }log\pi (\Omega |\mathbf{y}),\end{array}$$

  (22)

  (ii)

  The half-step update (22) is used to full-step update of the DE parameters $\Omega$ by

  $$\begin{array}{c}\hfill {\Omega }^{*}=\Omega +ϵM^{-1}{\mathbf{OE}}^{*},\end{array}$$

  (23)

  where M is the variance-covariance matrix of the momentum prior distribution $p\left(\mathbf{OE}\right)$.

  (iii)

  The updated ${\Omega }^{*}$ from (23) is then used to update the the second half-step of $\mathbf{OE}$ by

  $$\begin{array}{c}\hfill {\mathbf{OE}}^{*}={\mathbf{OE}}^{*}+\frac{1}{2}ϵ\frac{d}{d\Omega }log\pi (\Omega \left|{\mathbf{y}\left)\right|}_{\Omega ={\Omega }^{*}}\right.,\end{array}$$

  (24)

• Accept-Reject Step: HMC accepts the candidate ${\Omega }^{*}$ as the next state of the MCMC chain, namely ${\Omega }^{\left(t\right)}$, with probability

  $$\begin{array}{c}\hfill min\left\{1,\frac{\pi \left({\Omega }^{*}\right|\mathbf{y})p\left({\mathbf{OE}}^{*}\right)}{\pi \left({\Omega }^{(t-1)}\right|\mathbf{y})p\left({\mathbf{OE}}^{(t-1)}\right)}\right\}.\end{array}$$

  (25)
  Note that we do not need to accept or reject the candidate ${\mathbf{OE}}^{*}$, because $\mathbf{OE}$ is a latent variable introduced to facilitate the transition of the HMC in the parameter space. Moreover, $\mathbf{OE}$ will be updated by the momentum initialization step in the next run of the HMC. The HMC finally employs the No-U-Turn strategy [50] to automatically tune parameters of $ϵ$ and L in the leapfrog step.

3.3. Stacked Neural Networks

In this subsection, we employ the properties of artificial neural networks (NNs) to incorporate information contained in the posterior distribution of the parameters in estimating the critical collapse functions. As a state-of-the-art predictive method, NNs find the solutions to the DE system by translating the problem into optimizing a loss function in estimating the DE variables [39]. NNs fit a multi-layer perceptron to estimate the DE variables of the system. Each layer of the NNs consists of several nodes learning the curvature and form of the DE variables over the system domain. The nodes in the layers are linked by

$$\begin{array}{c}\hfill N^m(\mathbf{x},t,\psi )=\eta \left({\mathbf{A}}^m{N}^{m-1}(\mathbf{x},t,\psi )+{\mathbf{b}}^m\right),\end{array}$$

(26)

where $N^{m-1}(·)$ denotes the observed response from the m-th layer of the NNs after applying the activation function $\eta (·)$ with $\mathbf{x}$ the input of the NNs, ${\mathbf{A}}^m$ a weight matrix of the m-th layer, ${\mathbf{b}}^m$ the bias of the m-th layer. Also, let $\psi$ encode the set of all unknown parameters of (26). Let $\mathcal{N}(\mathbf{x},t,\psi )$ represent the output of the NNs in estimating the DE variables $\mathbf{g}=(g_1,\dots ,g_H)$ of (19). The solution to DE system (19) is finally derived by optimizing the loss function

$$\widehat{\mathcal{N}}(\mathbf{x},t,\psi )=arg\underset{\psi }{min}{\left(\mathcal{N}(\mathbf{x},t,\psi )-\mathbf{g}\right)}^2.$$

For more details about NNs and how they are used to estimate the DE variables, readers are referred to [39,51] and references therein.

In this paper, we first use the HMC sampler to find the posterior distribution for the parameter of the DE system. Let ${\theta }_l^{*},l=1,\dots ,L$ denote L samples from the domain of the posterior distribution $\pi \left(\theta \right|\mathbf{y})$. Applying these samples to the NN solver, let $\widehat{\mathcal{N}}(\mathbf{x}\left(t\right|{\theta }_l^{*}),t,\psi )$ denote the NN-based estimator given by ${\theta }_l^{*},l=1,\dots ,L$ from the posterior distribution. One can view the $\widehat{\mathcal{N}}(\mathbf{x}\left(t\right|{\theta }_l^{*}),t,\psi )$ from a probabilistic viewpoint, as an estimator for the DE variables $\mathbf{x}\left(t\right|\theta )$ whose probability distribution corresponds to the posterior distribution $\pi \left(\theta \right|\mathbf{y})$.

Accordingly, we propose a Bayesian model averaging strategy to stack the NNs in estimating the DE variables $\mathbf{x}\left(t\right|\theta )$. Let ${\theta }_l,l=1,\dots ,L$ denote a set of candidate models from the training set $\mathbf{y}$. One can compute the distribution of $\widehat{\mathcal{N}}(\mathbf{x}\left(\mathbf{t}\right|{\theta }^{*}),t,\psi )$ at some fixed space-time point t by

$$\begin{array}{c}\hfill P\left(\widehat{\mathcal{N}}(\mathbf{x}\left(t\right|\theta ),t,\psi )\right|\mathbf{y})=\sum _{l=1}^{L}P\left(\widehat{\mathcal{N}}(\mathbf{x}\left(t\right|{\theta }^l),t,\psi )\right|\mathbf{y})\pi \left({\theta }^l\right|\mathbf{y}).\end{array}$$

(27)

Thus, from (27), the posterior mean of $\widehat{\mathcal{N}}(\mathbf{x}\left(\mathbf{t}\right|{\theta }^{*}),t,\psi )$ at some fixed space-time point t is given by

$$\begin{array}{c}\hfill E\left(\widehat{\mathcal{N}}(\mathbf{x}\left(t\right|\theta ),t,\psi )\right|\mathbf{y})=\sum _{l=1}^{L}E\left(\widehat{\mathcal{N}}(\mathbf{x}\left(t\right|{\theta }^l),t,\psi )\right|\mathbf{y})\pi \left({\theta }^l\right|\mathbf{y}).\end{array}$$

(28)

This Bayesian estimate can be viewed simply as the weighted average of the individual NN-based estimates using the weights regulated by the posterior distribution $\pi \left(\theta \right|\mathbf{y})$.

The Bayesian model averaging (28) leads to multiple proposals to stack the NN estimators. The first stacked NN-based method is to take the unweighted average of the NNs using the samples ${\theta }^l,l=1,\dots ,L$ from $\pi \left(\theta \right|\mathbf{y})$ at each space-time. Although the mean stacked NNs are computed based on the unweighted average of the NNs, the random sampling of the candidates from the posterior distribution indirectly accommodates $\pi \left(\theta \right|\mathbf{y})$ in (28). Thus, the mean stacked NN estimate can be considered the posterior mean estimate. In the second approach, we propose the posterior median estimate. Here, we first stack the NNs estimates $\widehat{\mathcal{N}}(\mathbf{x}\left(t\right|{\theta }^l),t,\psi ),l=1,\dots ,L$ and then compute the median of them at each space-time point in the domain of the DE system.

4. Numerical Studies

In this section, we estimate the critical collapse functions of the equations of motion in the elliptic class of 4 dimensions in a Bayesian framework. In this framework, the parameters of the equations are treated as random variables to account for the measurement error involved in numerically estimating the parameters of the equations. We also explore the effect of this complexity on the estimation of the critical collapse functions.

Various research articles in the literature [36,39] investigated the non-linear equations of motion of the axion-dilaton system in different dimensions and for different ansatz. Hatefi et al. [37,38] used the properties of the Fourier-based and spline nonlinear smothers in predicting the critical collapse functions, respectively. Recently, Hatefi et al. [39] developed an NN-based solver to show that there is no solution in the parabolic class of black holes and [36] used a grid search root-finding method to estimate the parameters of the equations of motion numerically. According to the complex and highly nonlinear pattern of the equations, all these methods in the literature use almost a series of numerical techniques to make the equations tractable and approximate the parameters of the models. These techniques include, for instance, imposing some constraints on the equations of motions such as finiteness of $f^{\prime \prime }\left(z\right)$ as $z\to z_{+}$. Another constraint that one employs is the vanishing of the divergent part of $f^{\prime \prime }\left(z\right)$ that produces a complex-valued constraint at $z_{+}$. As the equations are not tractable analytically in elliptic equations of motion, for example, [33] used a numerical profile-root finding method and applied numerically a grid search discrete optimization on the coordinates of the extended parameter space of the equations of motion. To do so, the equations are simplified by approximating the target functions by the first two orders of the Taylor expansions. Applying a grid search optimization on the coordinates of the extended parameter space in the equation of motion, [33] could estimate the required parameters of the equations of motion. Removing the spurious roots of the equations, they used parameter estimates to find the solution to the equations of motion and estimate the critical collapse functions.

Unlike [36,37,38,39], to our best knowledge, for the first time in the literature, through a probabilistic viewpoint, we develop a Bayesian method to model the measurement errors into the estimation of the critical collapse functions. Hence the framework provides a robust approach against measurement errors in parameter estimation. These Bayesian estimates incorporate information from all possible parameters’ outcomes, that may occur in the numerical experiments, from the posterior distribution and accommodate this complexity in estimating critical collapse functions.

Here all the critical collapse functions are treated as DE variables of DE equations of motion 16 and 17. Since all the DE variables are numerically solved simultaneously, it is reasonable to assume that the observed DE variables have the same standard deviation $\sigma$ parameter representing the variability in the numerical experiments. Therefore $\Omega =(\omega ,\sigma )$ will be the set of all unknown parameters of the likelihood function under equations of motion. We first use the HMC sampler, described in SubSection 3.2, to find the posterior distribution of the likelihood parameter $\Omega$. We used the Python package Pymc3 [52] to implement the HMC for the equations of motion. To evaluate the effect of the prior information of likelihood parameters in estimating the critical functions, we considered two different prior distributions for the parameter $\omega$. We set the first prior distribution as Gaussian distribution with a mean $1.15$ and standard deviation $0.02$; that is ${\pi }_1\left(\omega \right)\sim N(\mu =1.15,\sigma =0.02)$. In contrast to the unimodal Gaussian, we assumed the second prior followed non-informative Uniform distribution between (1,1.5); that ${\pi }_2\left(\omega \right)\sim U(a=1,b=1.5)$. We also used the half-Cauchy distribution with scale parameter $\beta =0.1$ as the prior distribution for $\sigma$ to incorporate the uncertainty of the observed DE variables in the likelihood function.

Figure 1 and Figure 6 show the posterior distributions $\omega$ and $\sigma$ under the Gaussian and Uniform prior distributions, respectively. Comparing Figure 1 and Figure 6, under both prior distributions, we observe that the posterior distribution of $\omega$ converges to a roughly uni-modal, bell-shaped distribution where the 94% highest density interval (HDI) roughly is (1.165,1.205). In other words, considering the measurement errors from numerical analysis of the parameter estimation, the true value of $\omega$ in the equations of motion, with probability 94%, will vary between (1.165,1.205). Interestingly this is compatible with the solution for the elliptic case for four dimensions, obtained using various numerical experiments [21,34], which as reported as $\omega =1.176$. We see that the 94% HDI contains (under both prior distributions) $w=1.176$ as one possible outcome of the numerical analysis. It can be concluded that given the numerical measurement errors involved in the DE variables, the probability that the true value of $\omega$ appears $1.176$ is roughly 13%; that is $\pi (w\in 1.176\pm {10}^{-3}|\mathbf{y})\propto 0.13$.

In the second part of the analysis, we incorporated the posterior information of the parameters into the NNs, described in SubSection 3.3, to obtain the NN-based posterior mean and posterior median estimators of the critical collapse functions $b_0\left(z\right),\mathrm{Im}\left(f\left(z\right)\right)$ and $\mathrm{Re}\left(f\right(z\left)\right)$. To do so, we randomly took $L=200$ realizations ${\omega }_l^{*},l=1,\dots ,L$ from the posterior distribution $\pi \left(\omega \right|\mathbf{y})$. We then applied these 200 realizations iteratively and obtained the NN-based estimates of $b_0\left(z\right|{\omega }_l^{*}),\mathrm{Re}\left(f\left(z\right|{\omega }_l^{*})\right)$ and $\mathrm{Im}\left(f\left(z\right|{\omega }_l^{*})\right)$ given that $\omega ={\omega }_l^{*},l=1,\dots ,L$ at $n=1000$ equally spaced $z\in [0,5]$ in the domain of the equations of motion. These NN-based estimators are denoted by $\widehat{\mathcal{N}}\left(b_0\left(z_i\right|{\omega }_l^{*}),z,\psi \right)$, $\widehat{\mathcal{N}}\left(\mathrm{Re}(f\left(z_i\right|{\omega }_l^{*}),z,\psi \right)$ and $\widehat{\mathcal{N}}\left(\mathrm{Im}(f\left(z_i\right|{\omega }_l^{*}),z,\psi \right)$ for $l=1,\dots ,L$ and $i=1,\dots ,n$. In each iteration, we implement the NNs solver using the Python Package NeuroDiffEq [53] configured with four hidden layers, each of 16 nodes, following [39]. We ran the solvers for 1000 epochs using the initial boundary conditions in Eq. 18.

Before establishing the Bayesian stacked NN estimators, we would like to explain the probabilistic perspective in sampling from the posterior distribution obtained. The posterior distribution $\pi \left(\omega \right|\mathbf{y})$ lists all possible values of $\omega$ and how often they are observed as the posterior estimate of the parameter of the equations of motion. The posterior candidate ${\omega }^{*}$, as a sample from the posterior distribution, changes from one sample to another. Accordingly, given the posterior candidate $\omega ={\omega }^{*}$, the posterior NN candidates, as estimators for the critical collapse functions, will change from one sample to another. The probability that each of these posterior NN candidates is observed in estimating the true critical collapse functions corresponds to the posterior distribution $\pi \left(\omega \right|\mathbf{y})$. To show this probabilistic perspective, Figure 2 and Figure 7 represent 10 posterior NN candidates of the critical collapse functions based on 10 randomly selected ${\omega }_l^{*},l=1,\dots ,10$ from $L=200$ realizations of $\pi \left(\omega \right|\mathbf{y})$ using the Gaussian and Uniform prior distributions, respectively.

Because we have L Bayesian NNs candidates corresponding to L realizations from the posterior distribution, we applied the Bayesian model averaging. We obtained the two mean and median stacked NNs for all the critical functions at each space-time $z_i,i=1,\dots ,n$. In addition to posterior mean and median estimates, we also computed the 95% Bayesian credible intervals for the critical functions on the entire domain. To compute the Bayesian 95% credible intervals, we sorted the NN estimates from the smallest value to the largest value at each space-time. We then computed the 2.5 and 97.5 percentiles of the estimates as lower and upper bands of the interval, respectively. To compare our method with the non-Bayesian approach of [39], we computed the NNs-based estimates of the critical collapse functions similar to [39], treating $\omega =1.176$ with the same configuration as described above. This estimate is, henceforth, called Fixed NNs.

Figure 3 and Figure 8 show the Bayesian posterior mean and median estimates and their 95% credible intervals using Gaussian and Uniform prior distributions, respectively. Moreover, Figure 3 and Figure 8 compare the Bayesian estimates with their corresponding fixed NN-based estimates of the critical functions. One can conclude from the figures that, given all possible values for $\omega$ in the equations of motion, the true critical collapse functions, with probability 0.95, should be in the reported Bayesian credible intervals. Also, we observe that all 95% credible intervals contain the Fixed NNs estimates of [39] as one possible candidate for estimating the critical functions. We also observe that the posterior mean and median NN estimator represent the form and curvature of the critical collapse functions as the average and median of all the Bayesian NN-based candidates.

Note that the fixed NN-based method of [39] is not robust against the variability and measurement error in calculating the parameter. In other words, the Fixed NN estimates will change slightly from $w=1.176$ in estimating the parameter’s value. Unlike the fixed NNs, the Bayesian credible intervals, posterior mean and median NN estimators will remain robust against changes in measurement errors in estimating w because these Bayesian estimates have already considered all possibilities of w in the domain from the posterior distribution.

Figure 4 and Figure 9 show the posterior mean and median stacked NNs and 95% credible intervals of all the critical collapse functions together under Gaussian and Uniform prior distributions, respectively. It is apparent from Figure 4 and Figure 9 that the NN-based Bayesian credible intervals confirm that there must be a single self-similar solution, on average, in the area $z\in [1.5,4.5]$ to the elliptic equations of motion in 4 dimensions where the imaginary and real parts of the G function intersects, where the explicit form of the G function for the elliptic case in 4 dimension was given in [33]. Interestingly, this finding is compatible with the finding in the literature where [36,39] showed that there was a single self-similar solution at $z=2.56$ when the parameter of the equation is estimated to be $w=1.176$.

In order to evaluate the convergence of the proposed NN-based estimates using the posterior distribution $\pi \left(\omega \right|\mathbf{y})$, we computed the train and validation losses in the last epoch of the NN solvers for each ${\omega }_l^{*},l=1,\dots ,L$. Figure 5 and Figure 10 show the trace and box plots of the loss differences for $L=200$ realizations from the posterior distribution $\pi \left(\omega \right|\mathbf{y})$ under Gaussian and Uniform prior distributions, respectively. We clearly observe that the loss differences, on average, are less than $|3\times {10}^{-3}|$ which validates the convergence of NN-based estimates in estimating the critical collapse functions in the elliptic class of 4 dimensions.

5. Summary and Concluding Remarks

This paper proposes a Bayesian methodology to model the self-similar solutions for the spherical gravitational collapse for an elliptic case in four dimensions. In the literature on the axion-dilaton system, due to the complex and highly nonlinear equations of motion, one needs to employ various numerical techniques and constraints to make the equations of motion tractable. These constraints include, for instance, the finiteness of $f^{\prime \prime }\left(z\right)$ as $z\to z_{+}$ and the vanishing of the divergent part of $f^{\prime \prime }\left(z\right)$ producing a complex-valued constraint at $z_{+}$. Because the equations of motion are not tractable analytically, for instance, [33] used a series of numerical methods to find the solutions to the equations. Their numerical profile root-finding includes, for example, approximating the main effects of the equations by the first two orders of the Taylor expansion, applying a numerical grid search discrete optimization on the coordinates of the extended parameter space of the equations and removing the spurious roots, to name a few. This enables them to estimate the parameters of the equations and then find self-similar black hole solutions. According to the complex form of the equations and the vital role of the parameters, researchers have to ignore the measurement errors imposed in finding the parameters through numerical techniques.

To fix this problem, we have proposed a novel pipeline to take the measurement errors into the statistical models in finding the solution to the elliptic equations of motion. To do so, unlike the methods in the literature, we developed a Bayesian estimation method where the parameter is treated as a random variable. We proposed Hamiltonian Monte Carlo to derive the posterior distribution of the parameter efficiently. The posterior distribution determines all the possible outcomes in estimating the parameter of the equations of motion. The posterior distribution interestingly shows that the deterministic solution available in the literature to the elliptic equations of motion, namely $\omega =1.176$, is the true value for the parameter with a probability of almost 0.13. In addition, the posterior distribution lists all possible values for the parameter of elliptic equations of motion along with their probabilities.

Next, we utilized neural networks (NNs) to incorporate the posterior information in solving the equations of motion and estimating the critical collapse functions in the entire domain of the equations. From a probabilistic viewpoint, one can find the NN estimate of the critical function as an estimator, which would be the true trajectory of the critical collapse function whose probability distribution corresponds to the posterior distribution of the parameter. Accordingly, we applied a Bayesian model averaging strategy and found the posterior mean and median stacked NNs and the 95% Bayesian credible intervals in estimating the critical collapse functions.

Given all variability in estimating the parameter of the equations of motion, the developed Bayesian credible intervals will contain the true critical collapse functions with a probability of 0.95. Comparing our Bayesian method with the deterministic NN approach of [39], the developed Bayesian credible intervals contain the deterministic estimate as one possible candidate in estimating the critical functions. Unlike the estimate of [39], the Bayesian proposals will remain robust against measurement errors in estimating $\omega$ because these Bayesian estimates have already considered all possibilities of the parameter in the domain of the posterior distribution. It is worth mentioning that our approach will help us to reveal the value and range of the critical exponent of axion-dilaton system in a wider range in different dimensions that we would like to pursue in the near future.