Charged black holes from a family of Born-Infeld-type electrodynamics models

Sebastián<EMAIL_ADDRESS>and Leonardo
<EMAIL_ADDRESS>
Departamento de Ciencias Físicas,
Facultad de Ingeniería y Ciencias
Universidad de La Frontera, Casilla 54-D
Temuco, Chile.

###### Abstract

We present a family of nonlinear electrodynamics models that are free of the
infinite self-energy of the point charge. Each model is dependent on a
dimensional nonlinearity parameter and is determined by the integer value of a
dimensionless parameter $n$. The Born-Infeld model is recovered when $n=1$.
Some of the characteristics of this family are studied. In addition, we study
the solutions of electrically charged AdS black holes, which result from
coupling these electrodynamic models to General Relativity. A Smarr formula
consistent with the first law of thermodynamics is also obtained in an
extended phase space.

Key words: Static black hole; Born-Infeld type electrodynamics; Smarr formula.

##  1 Introduction

The model of Born-Infeld electrodynamics [1] is free of singularities. The
electric field derived from this model decreases from its maximum value at the
origin (equivalent to the nonlinearity parameter $b$ of the model) until it
behaves like the electric field of model Maxwell at large distances. One can
also see that the self-energy of a point charge is finite. The parameter $b$
has been related to the string tension [2, 3]. Possible constraints for the
parameter $b$ of this model have been analyzed in Refs. [4, 5, 6, 7, 8, 9,
10].

Unlike the Euler-Heisenberg electrodynamics [11], the Born-Infeld model does
not exhibit the phenomenon of vacuum birefringence in the presence of external
electric field [12, 13]. The causality and unitarity principles (tachyons and
ghosts are absent, respectively) are held in the Born-Infeld theory [14].

Born-Infeld electrodynamics has inspired other models free of singularities,
which, in general, have the same characteristics: the Soleng model [15],
various models presented by S. Kruglov [16, 17, 18, 19, 20, 21], others
proposed by Gaete et al. [22, 23, 24], the exponential model [25], the models
presented by Mazharimousavi and Halilsoy [26, 27] and the double-logarithmic
model [28].

Some models that generalize the BI model have been recently introduced, that
is, at some limit of the model parameters, the Born-Infeld electrodynamics is
recovered: the Born-Infeld like model with two parameter[29], the Born-Infeld
type model with three parameters [30], the ModMax model [31], the ModMax
generalized model [32].

The first extension of the Reissner-Nordström solution was obtained using
Born-Infeld electrodynamics [33, 34]. In contrast to Bardeen-type black holes
(for example Refs. [35, 36, 37, 38, 39, 40]) whose metric functions and
curvature invariants are regular in everywhere, a charged black hole obtained
with electrodynamics of Born-Infeld, presents singularities when evaluating
these quantities. Additionally, depending on the values of the mass, the
charge and the parameter $b$ of the Born-Infeld black hole, the solutions
obtained can be categorized as Reissner-Nordström type or Schwarzschild type
[41]. The first ones can have two horizons or one degenerate horizon; the
second ones have a single horizon.

The energy-momentum tensor of the Born-Infeld model has a non-zero trace,
because of this a Smarr formula consistent with the first law of black hole
thermodynamics cannot be found without considering a new thermodynamic
quantity [42] (see also Ref. [43]). In the study of the thermodynamics of
Born-Infeld Anti-de Sitter black holes in the extended phase space [44], where
the cosmological constant represents the thermodynamic pressure [45, 46, 47],
the parameter $b$ has been considered as a thermodynamic variable, which,
together with its corresponding conjugate quantity [48], allows obtaining a
Smarr formula consistent with the first law of thermodynamics. In this
framework, the quantity conjugated to the nonlinearity parameter has been
interpreted as vacuum polarization per unit volume [44]. Other works [49, 50,
51, 52, 53, 54, 55] have also considered thermodynamics in the extended phase
space, also profiting the nonlinearity parameter of the model as an extra
thermodynamic quantity, to study some features of Born-Infeld or Born-Infeld-
type black holes.

In this paper we propose a family of nonlinear electrodynamics models that are
also dependent on a nonlinearity parameter $b$, but where the different
solutions of this family are obtained by specifying the value of a
dimensionless parameter $n=1,2,...$, which determines the order of the field
strength that follow the term obtained from Maxwell electrodynamics when
$r\rightarrow\infty$. In this family when $n=1$ the Born-Infeld model is
obtained and in all cases when $b\rightarrow\infty$ we recover the Maxwell
model. We study some of the characteristics that these models present,
including the phenomenon of birefringence in the presence of an external
magnetic field and later in the presence of an external electric field. In
addition, we obtain static electrically charged AdS black hole solutions by
coupling Einstein equations to these electrodynamics models and from here we
obtain algebraic expressions for both the electric field and the metric
function. We also show that for each value of $n$, a black hole solution is
obtained that has a different behavior depending on the mass, the electric
charge and the parameter $b$, which is manifested in the number of horizons it
has, such as is the case for the Born-Infeld black hole. Finally, considering
$b$ as a new thermodynamic quantity in the extended phase space allows us to
obtain a Smarr formula consistent with the corresponding first law of
thermodynamics for any value of $n$.

The paper is organized as follows. In section 2, the family of nonlinear
electrodynamics models considered is presented and some of its main features
are studied. In section 3, we present the black hole solutions that are
obtained by coupling General Relativity with these electrodynamic models.
Section 4 is devoted to describe the different types of behavior of black hole
solutions for each value of $n$. Section 5 presents the obtaining of a Smarr
formula consistent with the first law of black hole thermodynamics.
Conclusions are given in section 6.

##  2 Family of Born-Infeld-type electrodynamics models

In this section we propose a family of new nonlinear electrodynamics models.
The motivation for introducing this family is that it has two features of
interest. First, each of the models in this family is free of singularities in
self-energy of the point charge. Second, the experimental interest in
measuring the effect of birefringence (see for example Ref. [56]), which is
present in these models, except for the case $n=1$.

Let us summarize the family of new nonlinear electrodynamics models in the
following Lagrangian

$\mathcal{L}(\mathcal{F},\mathcal{G})=-\left(\mathcal{F}-\frac{1}{2b^{2}}\mathcal{G}^{2}\right){}_{2}F_{1}\left(\frac{1}{2n},\frac{1}{n};1+\frac{1}{n};\frac{(-2)^{n}}{b^{2n}}\left(\mathcal{F}-\frac{1}{2b^{2}}\mathcal{G}^{2}\right)^{n}\right)\,,$
(1)

where $n>0$ is an integer,
$\mathcal{F}=F_{\mu\nu}F^{\mu\nu}/4=(B^{2}-E^{2})/2$,
$\mathcal{G}=F_{\mu\nu}\tilde{F}^{\mu\nu}/4=\mathbf{E}\cdot\mathbf{B}$, $b$ is
a nonlinearity parameter of the theory and ${}_{2}F_{1}\left(a,b;c;z\right)$
is the Gauss hypergeometric function. In the limit $b\rightarrow\infty$ (in
the weak field limit) the Lagrangian behaves as $\mathcal{L}=-\mathcal{F}$.
Note that when $n=1$ in the Lagrangian (1), we recover the Lagrangian of Born-
Infeld electrodynamics [1]

$\mathcal{L}_{BI}=b^{2}\left(1-\sqrt{1+\frac{2\mathcal{F}}{b^{2}}-\frac{\mathcal{G}^{2}}{b^{4}}}\right)\,,$
(2)

The form of Lagrangian (1) is inspired by that of the Born-Infeld model when
written using the Gauss hypergeometric function.

If we assume $\mathbf{B}=0$ and $\mathbf{E}=E(r)\hat{r}$, that is,
$\mathcal{F}=-E^{2}(r)/2$ and $\mathcal{G}=0$, then from Eq. (1) we obtain the
following field equations

$\nabla_{\mu}(F^{\mu\nu}\mathcal{L}_{\mathcal{F}})=0\,,$ (3)

where $\mathcal{L}_{\mathcal{F}}$ denotes differentiation of $\mathcal{L}$
with respect to $\mathcal{F}$, which yields

$\frac{E}{\left[1-(-1)^{-n}b^{-2n}\left(-E^{2}\right)^{n}\right]^{\frac{1}{2n}}}=\frac{q}{r^{2}}\,\,\,.$
(4)

The solution of this equation provides an expression for the electric field

$E(r)=\frac{b\,q}{(q^{2n}+b^{2n}r^{4n})^{\frac{1}{2n}}}\,\,\,.$ (5)

Here, it can be noted that the electric field asymptotically behaves as the
Maxwell field, i.e.,

$E(r)=\frac{q}{r^{2}}-\frac{q^{2n+1}}{2n\,b^{2n}\,r^{2(2n+1)}}+\frac{(2n+1)q^{4n+1}}{8n^{2}b^{4n}r^{2(4n+1)}}+O\left(\frac{1}{r^{2(6n+1)}}\right)\,\,\,.$
(6)

For its part, expanding the electric field about $r=0$ yields

$E(r)=b-\frac{b^{2n+1}}{2nq^{2n}}r^{4n}+\frac{(2n+1)b^{4n+1}}{8n^{2}q^{4n}}r^{8n}+O\left(r^{12n}\right)\,\,\,.$
(7)

In addition one can see that the parameter $b$ represents the maximum value of
the electromagnetic field strength for any value of $n$, just as is known to
occur for the case $n=1$. The shape of the electric field (5) for $n=1,2,3$
and $4$ is shown in Fig. 1 where the starting point is at the value of $b$.
For comparison purposes, the figure also shows the Maxwell case.

Figure 1: The electric field is plotted for some values of $n$. Maxwell
electric field is also depicted. Here, in all cases, the parameters are chosen
as $q=0.7$ and $b=2$.

###  2.1 Self-energy

In Ref. [57] the authors conjecture that for theories belonging to the Born-
Infeld family in 3+1-dimensions, the value of the self-energy for charged
particles is finite. In our case, we calculate these values and obtain a
result that is finite for each $n$, as it is shown in the following.

Self-energy in 3+1 dimensions is defined as

$\varepsilon=4\pi\int^{\infty}_{0}T_{0}^{0}\,r^{2}dr=-\int^{\infty}_{0}\left(E^{2}\mathcal{L}_{\mathcal{F}}+\mathcal{L}\right)\,r^{2}dr\,.$
(8)

Looking at Eq. (1), one readily obtains

$\mathcal{L}_{\mathcal{F}}=-\left(1-\frac{(-2)^{n}\mathcal{F}^{n}}{b^{2n}}\right)^{-\frac{1}{2n}}$
(9)

which together with the expression for the electric field allows us to
calculate the self-energy

$\varepsilon(n)=\frac{2\Gamma\left(\frac{1}{4n}\right)\Gamma\left(\frac{4n+1}{4n}\right)}{3\Gamma\left(\frac{1}{2n}\right)}\sqrt{bq^{3}}$
(10)

where $\varepsilon(1)=1.23605\sqrt{bq^{3}}$ and it can reach a maximum value
of $\varepsilon(n\rightarrow\infty)=4\sqrt{bq^{3}}/3$.

###  2.2 Causality and unitarity

Here we analyze the causality and unitarity of the model when $n$ is odd and
when it is even.

The principles of causality and unitarity are satisfied if
$\mathcal{L}_{\mathcal{F}}\leq 0$, $\mathcal{L}_{\mathcal{FF}}\geq 0$,
$\mathcal{L}_{\mathcal{GG}}\geq 0$,
$\mathcal{L}_{\mathcal{F}}+2\mathcal{F}\mathcal{L}_{\mathcal{FF}}\leq 0$ and
$2\mathcal{F}\mathcal{L}_{\mathcal{GG}}-\mathcal{L}_{\mathcal{F}}\geq 0$ [14].

If we consider Lagrangian (1) and the abbreviation
$\mathcal{K}=-2b^{2}\mathcal{F}+\mathcal{G}^{2}$, we get

$\mathcal{L}_{\mathcal{F}}=-\frac{1}{\left(1-b^{-4n}\mathcal{K}^{n}\right)^{\frac{1}{2n}}}$
(11)

$\mathcal{L}_{\mathcal{F}\mathcal{F}}=\frac{b^{4}\mathcal{K}^{n-1}}{\left(b^{4n}-\mathcal{K}^{n}\right)^{1+\frac{1}{2n}}}$
(12)

$\mathcal{L}_{\mathcal{G}\mathcal{G}}=\frac{b^{4n}+2b^{2}\mathcal{F}\mathcal{K}^{n-1}}{\left(b^{4n}-\mathcal{K}^{n}\right)^{1+\frac{1}{2n}}}$
(13)

$\mathcal{L}_{\mathcal{F}}+2\mathcal{F}\mathcal{L}_{\mathcal{F}\mathcal{F}}=-\frac{b^{2}\left(b^{4n}-\mathcal{G}^{2}\mathcal{K}^{n-1}\right)}{\left(b^{4n}-\mathcal{K}^{n}\right)^{1+\frac{1}{2n}}}$
(14)

$2\mathcal{F}\mathcal{L}_{\mathcal{G}\mathcal{G}}-\mathcal{L}_{\mathcal{F}}=\frac{\left[b^{4n}\left(b^{2}+2\mathcal{F}\right)-b^{2}\mathcal{K}^{n}+4b^{2}\mathcal{F}^{2}\mathcal{K}^{n-1}\right]}{\left(b^{4n}-\mathcal{K}^{n}\right)^{1+\frac{1}{2n}}}$
(15)

According to these results, we obtain the following table showing the cases in
which the inequalities are satisfied (✓) or not (✕).

Condition | $n$ | $B=0$ | $E>B$ | $E=0$ | $B>E$
---|---|---|---|---|---
$\mathcal{L}_{\mathcal{F}}\leq 0$ | odd | ✓ | ✓ | ✓ | ✓
even | ✓ | ✓ | ✓* | ✓*
$\mathcal{L}_{\mathcal{FF}}\geq 0$ | odd | ✓ | ✓ | ✓ | ✓
even | ✓ | ✓ | ✕ | ✕
$\mathcal{L}_{\mathcal{GG}}\geq 0$ | odd | ✓ | ✓ | ✓ | ✓
even | ✓ | ✓ | ✓* | ✓*
$\mathcal{L}_{\mathcal{F}}+2\mathcal{F}\mathcal{L}_{\mathcal{FF}}\leq 0$ | odd | ✓ | ✓ | ✓ | ✓
even | ✓ | ✓ | ✓* | ✓*
$2\mathcal{F}\mathcal{L}_{\mathcal{GG}}-\mathcal{L}_{\mathcal{F}}\geq 0$ | odd | ✓ | ✓ | ✓ | ✓
even | ✓ | ✓ | ✓* | ✓*

Table 1: Causality and unitarity. Here $\mathcal{G}=0$ and the asterik
indicate that the condition is satisfied only if $\mathcal{F}\leq b^{2}$.

We can notice in the table above that, just like for Born-Infeld
electrodynamics [14], for all other values of $n$ odd with $\mathcal{G}=0$ the
principles of causality and unitarity are always satisfied. For its part, for
$n$ even, only when $\mathcal{F}<0$ and considering $\mathcal{G}=0$, the
conditions are fully satisfied.

###  2.3 Birefringence

Another feature that can be considered for our model is the birefringence
[58]. The Born-Infeld model does not present the phenomenon of birefringence
[59]. Some authors have considered the study of the phenomenon of
birefringence in different models of non-linear electrodynamics: Euler-
Heisenberg theory [60], Soleng electrodynamics [22], Born-Infeld-type
electrodynamics [30], Gaete-Helayël-Neto model [23], inverse electrodynamics
[61], logarithmic electrodynamics [24], Euler-Heisenberg-like electrodynamics
[62] and electrodynamics with three parameters [64] and ModMax model [31, 63].

To analyze this phenomenon in the case of our model, and following Ref. [12],
we consider the propagation of the electromagnetic plane wave
$(\mathbf{E}_{p},\mathbf{B}_{p})$ in a constant external magnetic field which
is along $x$ axis $\overline{\mathbf{B}}=(\overline{B},0,0)$ where

$\mathbf{E}_{p}(\textbf{r},t)=\mathbf{E}_{0}\,e^{-i(\omega
t-\textbf{k}\cdot\textbf{r})}\,\,\,,\,\,\,\,\,\mathbf{B}_{p}(\textbf{r},t)=\mathbf{B}_{0}\,e^{-i(\omega
t-\textbf{k}\cdot\textbf{r})}\,.$ (16)

We assume the wave propagates in the $z$-direction and the background magnetic
field is much stronger than the amplitudes of electric and magnetic fields.
Then the resultants electromagnetic fields become
$\mathbf{B}=\mathbf{B}_{p}+\overline{\mathbf{B}}$ and
$\mathbf{E}=\mathbf{E}_{p}$. Thus substituting into the Lagrangian (1) we get

$\mathcal{L}(\mathbf{E}_{p},\mathbf{B}_{p}+\overline{\mathbf{B}})=-\left(F_{\mathbf{B}}\right)\,_{2}F_{1}\left(\frac{1}{2n},\frac{1}{n};1+\frac{1}{n};-2^{n}\left(F_{\mathbf{B}}\right)^{n}b^{-2n}\right)\,,$
(17)

where $F_{\mathbf{B}}$ is defined as

$F_{\mathbf{B}}=\frac{1}{2}\left(\frac{\left(\mathbf{B}_{p}+\overline{\mathbf{B}}\right)^{2}}{2}-\frac{\mathbf{E}_{p}^{2}}{2}\right)-\frac{\left(\left(\mathbf{B}_{p}+\overline{\mathbf{B}}\right)\cdot\mathbf{E}_{p}\right)^{2}}{2b^{2}}\,.$
(18)

One can obtain the electric displacement and magnetic field by using
$\mathbf{D}=\partial\mathcal{L}/\partial\mathbf{E}$ and
$\mathbf{H}=-\partial\mathcal{L}/\partial\mathbf{B}$, respectively and taking
into account only the linear terms in $\mathbf{E}_{p}$ and $\mathbf{B}_{p}$

$\mathbf{D}_{p}=\frac{1}{\mathcal{R}}\left(\mathbf{E}_{p}+\frac{\overline{\mathbf{B}}\,(\overline{\mathbf{B}}\cdot\mathbf{E}_{p})}{b^{2}}\right)$
(19)
$\mathbf{H}_{p}=\frac{1}{\mathcal{R}}\left(\mathbf{B}_{p}+\left(-1\right)^{n}\frac{\overline{\mathbf{B}}^{2(n-1)}\overline{\mathbf{B}}\,(\overline{\mathbf{B}}\cdot\mathbf{B}_{p})}{b^{2n}\mathcal{R}^{2n}}\right)\,,$
(20)

where we have introduced the abbreviation

$\mathcal{R}=\left(1+\left(-1\right)^{n}\frac{\overline{\mathbf{B}}^{2n}}{b^{2n}}\right)^{\frac{1}{2n}}\,.$
(21)

We may calculate the electrical permittivity and the inverse of the magnetic
permeability tensors, using the relations $D_{i}=\varepsilon_{ij}(E_{p})_{j}$
y $H_{i}=\left(\mu^{-1}\right)_{ij}(B_{p})_{j}$

$\varepsilon_{ij}=\frac{1}{\mathcal{R}}\left(\delta_{ij}+\frac{\overline{B}_{i}\overline{B}_{j}}{b^{2}}\right)$
(22)
$\left(\mu^{-1}\right)_{ij}=\frac{1}{\mathcal{R}}\left(\delta_{ij}+\left(-1\right)^{n}\frac{\overline{\mathbf{B}}^{2(n-1)}\overline{B}_{i}\overline{B}_{j}}{(b\,\mathcal{R})^{2n}}\right)\,.$
(23)

From Maxwell equations written in terms $\mathbf{k}$ and $\omega$, we obtain
the wave equation for the electric field $\mathbf{E}_{p}$ [12]

$\left[\varepsilon_{ijl}\,\varepsilon_{abc}\,k_{j}\left(\mu^{-1}\right)_{la}k_{b}+\omega^{2}\varepsilon_{ic}\right](E_{p})_{c}=0\,,$
(24)

where $\varepsilon_{ijl}$ is the Levi-Civita symbol.

Substituting $\varepsilon_{ij}$ and $\left(\mu^{-1}\right)_{ij}$ in this last
equation, we get the following equation in matrix notation,

$\frac{\omega^{2}}{\mathcal{R}}\left[AI+C\overline{\textbf{B}}\cdot\overline{\textbf{B}}\right]\textbf{E}_{p}=0\,$
(25)

where $I$ is the identity matrix,

$A=1-n^{2}-(-1)^{n}\frac{n^{2}\overline{\mathbf{B}}^{2n}}{b^{2n}\mathcal{R}^{2n}}\,\,\,\,,\,\,\,\,\,\,\,\,C=\frac{1}{b^{2}}+(-1)^{n}\frac{\overline{\mathbf{B}}^{2(n-1)}}{b^{2n}\mathcal{R}^{2n}}$
(26)

and

$\overline{\textbf{B}}\cdot\overline{\textbf{B}}=\left(\begin{array}[]{ccc}\overline{B}_{x}^{2}&\overline{B}_{x}\overline{B}_{y}&\overline{B}_{x}\overline{B}_{z}\\\
\overline{B}_{x}\overline{B}_{y}&\overline{B}_{y}^{2}&\overline{B}_{y}\overline{B}_{z}\\\
\overline{B}_{x}\overline{B}_{z}&\overline{B}_{y}\overline{B}_{z}&\overline{B}_{z}^{2}\\\
\end{array}\right)\,.$ (27)

Then we can find the following eigenvalues

$\lambda_{1}=1-n^{2}\left(1+(-1)^{n}\frac{\overline{\mathbf{B}}^{2n}}{(b\mathcal{R})^{2n}}\right),\,\,\,\,\lambda_{2}=1-n^{2}+\frac{\overline{\mathbf{B}}^{2}}{b^{2}}\,.$
(28)

By setting $\lambda_{1}=0$ y $\lambda_{2}=0$, we arrive at the indices of
refraction when the polarization of the electromagnetic wave is perpendicular
and parallel to the external magnetic field $\overline{\mathbf{B}}$,
respectively

$n_{\perp}=\sqrt{1-\left(-1\right)^{n}\frac{\overline{\mathbf{B}}^{2n}}{b^{2n}}}\,,\,\,\,\,\,n_{\parallel}=\sqrt{1+\frac{\overline{\mathbf{B}}^{2}}{b^{2}}}\,,$
(29)

or if we consider $b$ large we can write them as

$n_{\perp}\approx
1-(-1)^{n}\frac{\overline{\mathbf{B}}^{2n}}{2b^{2n}}\,,\,\,\,\,n_{\parallel}\approx
1+\frac{\overline{\mathbf{B}}^{2}}{2b^{2}}\,.$ (30)

Note that when $n\neq 1$ the phenomenon of birefringence is present. In other
words, electromagnetic waves polarized in perpendicular directions travel at
different velocities. Since the velocity of the electromagnetic wave can be
found using the relation $v_{\parallel(\perp)}=c/n_{\parallel(\perp)}$ (here
we explicitly write the constant $c$), it can also be noted, even when the
higher order contributions are very small, that if $n$ is odd then
$n_{\perp}>1$, whereas if $n$ is even we would have non-physics values, that
is $n_{\perp}<1$.

Analogously, we can now calculate the refractive indices in the presence of
only a purely external electric field. If we consider that the direction of
the background electric field is along the $x$ axis
$\overline{\mathbf{E}}=(\overline{E},0,0)$ and the propagation direction of
the wave in the $z$ direction as in the previous case, then in this case we
arrive at the following refractive indices

$n_{\perp}=\left(1-\frac{\overline{\mathbf{E}}^{2}}{b^{2}}\right)^{-\frac{1}{2}}\,,\,\,\,\,n_{\parallel}=\left(1-\frac{\overline{\mathbf{E}}^{2n}}{b^{2n}}\right)^{-\frac{1}{2}}\,.$
(31)

If we consider $b$ large, then

$n_{\perp}\approx
1+\frac{\overline{\mathbf{E}}^{2}}{2b^{2}}\,,\,\,\,\,n_{\parallel}\approx
1+\frac{\overline{\mathbf{E}}^{2n}}{2b^{2n}}\,.$ (32)

We can notice that, unlike the previous case, the refractive indices are
greater than 1 for $n$ even or odd. Here also the phenomenon of birefringence
is present, except for the case where $n=1$.

##  3 AdS Black holes with Born-Infeld-type electrodynamics

In this and the next two sections we present solutions of electrically charged
AdS black holes that are obtained from the nonlinear electrodynamics model
discussed above. We study the three categories of solutions that are obtained
for each value of $n$. We also consider thermodynamics in extended phase space
and obtain a Smarr formula consistent with the first law of thermodynamics.

The action of the 3+1-dimensional Einstein theory coupled with a nonlinear
electrodynamics model is described by the action

$S=\int d^{4}x\sqrt{-g}\left[\frac{(R-2\Lambda)}{16\pi
G}+\mathcal{L(F)}\right]\,.$ (33)

Here, $g$ is the determinant of the metric tensor, $R$ is the Ricci scalar,
$G$ is the gravitational constant, $\Lambda=-\frac{3}{l^{2}}$ is the
cosmological constant and $\mathcal{L(F)}$ is the Lagrangian given in Eq. (1).

The Einstein fields equations resulting from action (33) are given by

$G_{\mu\nu}+\Lambda g_{\mu\nu}=8\pi T_{\mu\nu}$ (34) $4\pi
T_{\mu\nu}=g_{\mu\nu}\mathcal{L(F)}-F_{\mu\rho}F_{\nu}^{\,\,\rho}\mathcal{L}_{\mathcal{F}}\,$
(35)

and

$\nabla_{\mu}(F^{\mu\nu}\mathcal{L}_{\mathcal{F}})=0$ (36)

Let us consider static spherically symmetric black holes where the line
elements read

$ds^{2}=-f(r)dt^{2}+f(r)^{-1}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})\,.$ (37)

Using the Einstein equations (34) and the Lagrangian (1) allows us to obtain
the metric function given by

$\displaystyle f(r)$ $\displaystyle=$ $\displaystyle
1-\frac{2M}{r}+\frac{r^{2}}{l^{2}}+\frac{4\sqrt{bq^{3}}}{3r}\frac{\Gamma\left(\frac{1}{4n}\right)\Gamma\left(\frac{4n+1}{4n}\right)}{\Gamma\left(\frac{1}{2n}\right)}+\frac{2b^{2}r^{2}}{3}\frac{\Gamma\left(\frac{1}{n}\right)\Gamma\left(\frac{2n-1}{2n}\right)}{\Gamma\left(\frac{1}{2n}\right)}$
(38)
$\displaystyle-\frac{2bq}{3}\,_{2}F_{1}\left(-\frac{1}{2n},\frac{1}{2n};1-\frac{1}{2n};-\frac{b^{2n}r^{4n}}{q^{2n}}\right)$
$\displaystyle-\frac{4bq}{3}\,_{2}F_{1}\left(\frac{1}{4n},\frac{1}{2n};1+\frac{1}{4n};-\frac{b^{2n}r^{4n}}{q^{2n}}\right)\,\,\,,$

where $M$ represents the ADM mass, $q$ the electric charge of the black hole.
This function can also be written as

$\displaystyle f(r)$ $\displaystyle=$ $\displaystyle
1-\frac{2M}{r}+\frac{r^{2}}{l^{2}}+\frac{2b^{2}r^{2}}{3}\frac{\Gamma\left(\frac{1}{n}\right)\Gamma\left(\frac{2n-1}{2n}\right)}{\Gamma\left(\frac{1}{2n}\right)}$
(39)
$\displaystyle-\frac{2bq}{3}\,_{2}F_{1}\left(-\frac{1}{2n},\frac{1}{2n};1-\frac{1}{2n};-\frac{b^{2n}r^{4n}}{q^{2n}}\right)$
$\displaystyle+\frac{4q^{2}}{3r^{2}}\,_{2}F_{1}\left(\frac{1}{4n},\frac{1}{2n};1+\frac{1}{4n};-\frac{q^{2n}}{b^{2n}r^{4n}}\right)\,\,\,\,.$

It can be seen that the asymptotic behavior of this metric function is

$f(r)=1+\frac{r^{2}}{l^{2}}-\frac{2M}{r}+\frac{q^{2}}{r^{2}}+O\left(\frac{1}{r^{4n+2}}\right)\,,$
(40)

Figure 2: Various metric functions corresponding to different values of $n$
are plotted. As in the Reissner-Nordström and Born-Infeld cases ($n=1$), the
solutions for other values of $n$ can also have two horizons, for values of
$q$ less than those listed in table 2. However, as in the Born-Infeld case,
the solutions can also have only one horizon, depending on the value of $q$,
as illustrated below. Here, in all cases, the parameters are chosen as $M=1$,
$q=0.7$, $\Lambda=0$ and $b=2$.

In Fig. 2, the metric functions for some values of $n$ together with the
Reissner-Nordström metric function are displayed.

Figure 3: Metric functions of extremal cases for $n=1,2,3$ and $4$. The circle
amplifies the region where the outer and inner horizons coincide for the
values of $n$ that we are displaying. In all cases $M=1$, $\Lambda=0$ and
$b=2$, for the electric charge we use the values indicated in table 2.

It is straightforward to find the mass function $m(r)$ from the metric
function, this allows us to obtain

$\lim_{r\rightarrow
0}m(r)=M-\frac{2\sqrt{bq^{3}}\,\Gamma\left(1+\frac{1}{4n}\right)\Gamma\left(\frac{1}{4n}\right)}{3\Gamma\left(\frac{1}{2n}\right)}\,.$
(41)

To ensure that this last result is positive, it is necessary that the electric
charge $q$ satisfies the following relation

$q\leq\frac{\left(\frac{3}{2}\right)^{2/3}M^{2/3}\Gamma\left(\frac{1}{2n}\right)^{2/3}}{\sqrt[3]{b}\,\Gamma\left(1+\frac{1}{4n}\right)^{2/3}\Gamma\left(\frac{1}{4n}\right)^{2/3}}=q_{ex}\,,$
(42)

where $q_{ex}$ is the electric charge of the extremal black hole. In Fig. (3)
extremal solutions are shown for different values of $n$.

Extremal electric charge
---
$n$ |  $q_{ext}$
1 | 1.00626550322001
2 | 1.00028890439369
3 | 1.00002502213484
4 | 1.00000287157123

Table 2: Electric charge of the extremal black holes
$q_{ext}$ for $n=1,2,3,4$, with $\Lambda=0$, $M=1$ and $b=2$.

##  4 Categories of solutions for each value of $n$.

As we shall see, for each value of $n$ the metric function can have zero, one
or two horizons, depending on the values of $M$, $q$ and $b$. Here we follow
the analysis as given in Ref. [41] or Refs. [65, 66].

To classify the black hole solutions that can be obtained, it is useful to
define the extreme mass. For this purpose, let us consider the metric function
(38) with fixed values of $M$, $l$ and $n$. Then one can numerically find the
extremal values $q_{ext}$ y $r_{ext}$ for which $f(r)=0$ and $f^{\prime}(r)=0$
at the horizon. Hence, the extremal mass of the black hole is defined as

$\displaystyle M_{ext}(q)$ $\displaystyle=$
$\displaystyle\frac{b^{2}r^{3}_{ext}\Gamma\left(1-\frac{1}{2n}\right)\Gamma\left(\frac{1}{n}\right)}{3\Gamma\left(\frac{1}{2n}\right)}+\frac{2\sqrt{bq^{3}}\,\Gamma\left(1+\frac{1}{4n}\right)\Gamma\left(\frac{1}{4n}\right)}{3\Gamma\left(\frac{1}{2n}\right)}$
(43)
$\displaystyle+\frac{1}{2}\left(\frac{r^{3}_{ext}}{l^{2}}+r_{ext}\right)-\frac{bqr_{ext}}{3}\bigg{[}\,_{2}F_{1}\left(-\frac{1}{2n},\frac{1}{2n};1-\frac{1}{2n};-\frac{b^{2n}r^{4n}_{ext}}{q^{2n}}\right)$
$\displaystyle+2\,_{2}F_{1}\left(\frac{1}{4n},\frac{1}{2n};1+\frac{1}{4n};-\frac{b^{2n}r^{4n}_{ext}}{q^{2n}}\right)\bigg{]}\,.$

Note that $M_{ext}(q_{ext})=M$.

At the turn, in the limit $r\rightarrow 0$, the metric function (38) behaves
as

$f(r)=1-\frac{2}{r}\left(M-\frac{2\sqrt{bq^{3}}}{3}\frac{\Gamma\left(\frac{1}{4n}\right)\Gamma\left(\frac{4n+1}{4n}\right)}{\Gamma\left(\frac{1}{2n}\right)}\right)-2bq+\frac{r^{2}}{l^{2}}+\frac{2b^{2}r^{2}}{3}\frac{\Gamma\left(\frac{1}{n}\right)\Gamma\left(\frac{2n-1}{2n}\right)}{\Gamma\left(\frac{1}{2n}\right)}+O(r^{4n})\,.$
(44)

and using the definition

$\alpha\equiv\frac{2\sqrt{bq^{3}}}{3}h(n)\,,$ (45)

where

$h(n)\equiv\frac{\Gamma\left(\frac{1}{4n}\right)\Gamma\left(\frac{4n+1}{4n}\right)}{\Gamma\left(\frac{1}{2n}\right)}\,,$
(46)

with $1.854<h(1)<h(2)<h(3)....<2$, one finds the following possible cases:

If $M>\alpha$ the solution behaves as a Schwarzschild black hole (S-type) and
has only one horizon.

If $M<\alpha$ the solution behaves as a Reissner-Nordström black hole (RN-
type) and this can have:

\- a naked singularity, when also $M<M_{ext}$.

\- two horizon, when also $M>M_{ext}$.

\- a degenerate horizon when also $M=M_{ext}$.

If $M=\alpha$, then $f(0)=1-2bq$ and one can distinguish two cases:

\- the solution behaves as a Schwarzschild black hole and has only one
horizon, when $bq>1/2$.

\- the solution has a naked singularity, when $bq<1/2$.

Figure 4: Behaviors of the different categories of black hole solutions for
$n=4$ with $M=1$, $b=2$ and $\Lambda=0$. The black curves represent the three
possible cases of RN-type solutions. The upper curve represents a solution
with a naked singularity ($M<M_{ext}$). The following curves represent the
extremal case ($M=M_{ext}$) and a solution with two horizons
($M~{}>~{}M_{ext}$), respectively. The upper red curve shows the marginal
case, which is the one that does not present singularity in the origin. The
other red curve represents a solution with a single horizon. The values of the
$q_{i}$ that allow to obtain each curve are displayed in table 3.

Note that since $h(n)$ can have the indicated range of values, all these
categories of solutions can be found for every $n$. In Fig. 4 the first five
solutions mentioned above are shown. For this we consider the case $n=4$ and
choose $M=1$, $b=2$ and $\Lambda=0$, where the values of electric charges are
given according to the correspondence indicated in table 3.

Electric charge $q$ | $\times M$
---|---
$q_{1}$ | 1.2
$q_{2}$ (extremal) | $\approx 1.0000029$
$q_{3}$ | 0.8
$q_{4}$ | $\approx 0.6578$
$q_{5}$ | 0.5

Table 3: Values of the electric charge considered in Fig. 4.

In Ref. [44] the parameter $\alpha$ is called marginal mass, hence the case
where $M=\alpha$ is called the marginal case. In Fig. 5 we illustrate
solutions that are inside the marginal case, for different values of the
parameter $b$, comparing the cases when $n=1$ and $n=4$. The solutions that
intersect the vertical axis above the horizontal axis correspond to solutions
with naked singularity, that is, solutions in which $bq<1/2$. Solutions such
that $bq=1/2$ are those in which the metric function at $r=0$ is zero. The
solutions under the horizontal axis are those that behave as Schwarzschild
black holes, i.e., the solutions have a single horizon. Note that if we fix
the value of $bq$, then the solutions for any $n$ have the same behavior when
$r\rightarrow 0$. Meanwhile, in Fig. 6 some metric functions are represented
for different values of $n$, illustrating the marginal case for $bq>0$, that
is, solutions that have a single horizon. In this figure, in each case the
values of the mass $M=\alpha$ and $l$ are left constant, and the parameter $b$
varies according to the correspondence given in table 4.

$n$ | parameter $b$
---|---
$1$ | $1.9082$
$2$ | $1.7129$
$3$ | $1.6737$
$4$ | $1.6593$

Table 4: Values of the parameter
$b$ used in the graph of Fig. 6. Figure 5: Comparison between the metric
functions for the marginal case when $n=1$ (Born-Infeld) and $n=4$ for
different values of $b$. The first two upper curves (which are closely
overlapping) represent solutions with naked singularity. Meanwhile, solutions
below the horizontal axis represent S-type solutions, which only have a single
horizon. We use $q=1$, $l=1$ and the values of $M=\alpha$ are obtained from
Eq. (45). Figure 6: Again the marginal case is illustrated. Here we compare
the metric functions of the Schwarzschild-type for different values of $n$,
where the values of $q$ and the marginal mass remain fixed. We use $q=0.7$ and
$M=\alpha=1$ and accordingly the values of the parameter $b$ vary for each $n$
as shown in table 4.

Carga eléctrica marginal
---
n | carga eléctrica q
1 | 0.68912714805
2 | 0.6647656184347
3 | 0.659652555844
4 | 0.65776354642

Table 5: Values of the electric charge that produce
marginal cases, considering $\Lambda=0$, $M=1$ y $b=2$.

Finally, in table 5, we include the values of the electric charge that produce
a marginal black hole for different values of $n$. This allows us to compare
the behavior of the marginal case for black holes with the same mass $M=1$ and
$\Lambda=0$, as shown in Fig. 7.

Figure 7: Marginal cases of the Schwarzschild-type for $n=1,2,3$ and $4$. Here
we consider fixed $\Lambda=0$ and $b=2$, thus $q$ varies according to what is
displayed in table 5. With these quantities, the values of $M=\alpha$ for each
$n$ are obtained from Eq. (45). The oval amplifies the region where the curves
intersect before reaching the vertical axis.

##  5 Smarr formula and first law of thermodynamics

In this section, we will derive thermodynamic quantities, first law and the
Smarr formula for the Born-Infeld-type black hole. We will be studying
thermodynamics in the extended phase space where the thermodynamic pressure is
given by

$P=-\frac{\Lambda}{8\pi}=\frac{3}{8\pi l^{2}}\,.$ (47)

The Hawking temperature on the event horizon $r_{+}$ can be obtained by using
the surface gravity $\kappa$ as

$\displaystyle T=\frac{\kappa}{2\pi}=\frac{1}{4\pi}\frac{df(r_{+})}{dr}$
$\displaystyle=$ $\displaystyle\frac{1}{4\pi r_{+}}+\frac{3r_{+}}{4\pi
l^{2}}+\frac{b^{2}r_{+}}{2\pi}\frac{\Gamma\left(1-\frac{1}{2n}\right)\Gamma\left(\frac{1}{n}\right)}{\Gamma\left(\frac{1}{2n}\right)}$
(48) $\displaystyle-\frac{bq}{2\pi
r_{+}}\,_{2}F_{1}\left(-\frac{1}{2n},\frac{1}{2n};1-\frac{1}{2n};-\frac{b^{2n}r_{+}^{4n}}{q^{2n}}\right)\,.$

The electric potential $\Phi$ on the event horizon can be calculated by using
the expression for the electric field given in Eq. (5)

$\Phi(r_{+})=\int_{r_{+}}^{\infty}E(r)dr=\frac{\sqrt{bq}\,\Gamma\left(1+\frac{1}{4n}\right)\Gamma\left(\frac{1}{4n}\right)}{\Gamma\left(\frac{1}{2n}\right)}-br_{+}\,{}_{2}F_{1}\left(\frac{1}{4n},\frac{1}{2n};1+\frac{1}{4n};-\frac{b^{2n}r_{+}^{4n}}{q^{2n}}\right)\,$
(49)

The black hole entropy is related to the event horizon as

$S=\pi r_{+}^{2}\,$ (50)

On the horizon $r_{+}$ the metric function satisfies $f(r_{+})=0$. From this
one obtains the mass as a function of entropy $S$, pressure $P$, electric
charge $q$ and nonlinearity parameter $b$

$\displaystyle M(S,P,q,b)$ $\displaystyle=$
$\displaystyle\frac{b^{2}S^{3/2}\Gamma\left(1-\frac{1}{2n}\right)\Gamma\left(\frac{1}{n}\right)}{3\pi^{3/2}\Gamma\left(\frac{1}{2n}\right)}+\frac{2\sqrt{bq^{3}}\,\Gamma\left(1+\frac{1}{4n}\right)\Gamma\left(\frac{1}{4n}\right)}{3\Gamma\left(\frac{1}{2n}\right)}+\frac{4PS^{3/2}}{3\sqrt{\pi}}$
(51)
$\displaystyle+\frac{1}{2}\sqrt{\frac{S}{\pi}}-\frac{bq}{3}\sqrt{\frac{S}{\pi}}\,_{2}F_{1}\left(-\frac{1}{2n},\frac{1}{2n};1-\frac{1}{2n};-\left(\frac{bS}{\pi
q}\right)^{2n}\right)$
$\displaystyle-\frac{2bq}{3}\sqrt{\frac{S}{\pi}}\,_{2}F_{1}\left(\frac{1}{4n},\frac{1}{2n};1+\frac{1}{4n};-\left(\frac{bS}{\pi
q}\right)^{2n}\right)\,$

The black hole entropy is conjugate to the temperature, then

$\displaystyle T=\left(\frac{\partial M}{\partial S}\right)_{q,P,b}$
$\displaystyle=$
$\displaystyle\frac{1}{4\sqrt{\pi}\sqrt{S}}+\frac{2P\sqrt{S}}{\sqrt{\pi}}+\frac{b^{2}}{2}\sqrt{\frac{S}{\pi^{3}}}\frac{\Gamma\left(1-\frac{1}{2n}\right)\Gamma\left(\frac{1}{n}\right)}{\Gamma\left(\frac{1}{2n}\right)}$
(52)
$\displaystyle-\frac{bq}{2\sqrt{\pi}\sqrt{S}}\,_{2}F_{1}\left(-\frac{1}{2n},\frac{1}{2n};1-\frac{1}{2n};-\left(\frac{bS}{\pi
q}\right)^{2n}\right)$

One can easily show, using Eqs. (47) and (50), that the two expressions (48)
and (52) for the black hole temperature are equivalent.

The electric potential on the event horizon is conjugate to the electric
charge $q$, so it can also be calculated as

$\displaystyle\Phi(r_{+})=\left(\frac{\partial M}{\partial q}\right)_{S,P,b}$
$\displaystyle=$
$\displaystyle\frac{\sqrt{bq}\,\Gamma\left(1+\frac{1}{4n}\right)\Gamma\left(\frac{1}{4n}\right)}{\Gamma\left(\frac{1}{2n}\right)}$
(53)
$\displaystyle-\sqrt{\frac{S}{\pi}}b\,_{2}F_{1}\left(\frac{1}{4n},\frac{1}{2n};1+\frac{1}{4n};-\left(\frac{bS}{\pi
q}\right)^{2n}\right)\,.$

Using Eq. (50), it is also easy to show that the result obtained for the
electric potential with Eq. (49) is equivalent to the one arrived at using Eq.
(53).

The differentiation of the mass $M(S,P,q,b)$, leads to the first law

$dM=\left(\frac{\partial M}{\partial S}\right)_{q,P,b}dS+\left(\frac{\partial
M}{\partial P}\right)_{S,q,b}dP+\left(\frac{\partial M}{\partial
q}\right)_{S,P,b}dq+\left(\frac{\partial M}{\partial
b}\right)_{S,q,P}db\,\,\,,$ (54)

where the black hole volume $V$ is the conjugate quantity to pressure $P$

$V=\left(\frac{\partial M}{\partial P}\right)_{S,q,b}\,,$ (55)

Following Refs. [48, 44], for the Born-Infeld black hole, we have introduced
the quantity $B$ as the conjugate quantity to $b$. Thus, we can calculate the
variable $B$ as

$\displaystyle B=\left(\frac{\partial M}{\partial b}\right)_{S,q,P}$
$\displaystyle=$
$\displaystyle\sqrt{\frac{q^{3}}{b}}\frac{\Gamma\left(1+\frac{1}{4n}\right)\Gamma\left(\frac{1}{4n}\right)}{3\Gamma\left(\frac{1}{2n}\right)}+\frac{2br_{+}^{3}\Gamma\left(1-\frac{1}{2n}\right)\Gamma\left(\frac{1}{n}\right)}{3\Gamma\left(\frac{1}{2n}\right)}$
(56)
$\displaystyle-\frac{2}{3}qr_{+}\,{}_{2}F_{1}\left(-\frac{1}{2n},\frac{1}{2n};1-\frac{1}{2n};-\frac{b^{2n}r_{+}^{4n}}{q^{2n}}\right)$
$\displaystyle-\frac{1}{3}qr_{+}\,{}_{2}F_{1}\left(\frac{1}{4n},\frac{1}{2n};1+\frac{1}{4n};-\frac{b^{2n}r_{+}^{4n}}{q^{2n}}\right)\,,$

where we have replaced $S$ using Eq. (50).

Consequently, these quantities satisfy the first law of thermodynamics

$dM=TdS+VdP+\Phi dq+Bdb\,\,\,.$ (57)

Considering the parameter $b$ as a thermodynamic quantity helps us to obtain a
Smarr formula [67] consistent with the first law of thermodynamics obtained
above. Following the scaling argument given in Ref. [45], we perform the
dimensional analysis, by using the following scaling relations $S\propto
L^{2}$, $P\propto L^{-2}$, $q\propto L^{1}$ and $b\propto L^{-1}$, where
$L=[\mbox{lenght}]$. Hence, as the mass is a homogeneous function of degree 1,
the Euler theorem implies

$M=(2)\left(\frac{\partial M}{\partial S}\right)S+(-2)\left(\frac{\partial
M}{\partial P}\right)P+(1)\left(\frac{\partial M}{\partial
q}\right)q+(-1)\left(\frac{\partial M}{\partial b}\right)b\,\,\,,$ (58)

Thus, we can write the Smarr formula as

$M=2TS-2VP+\Phi q-Bb\,\,\,.$ (59)

Note that using the Komar integral [68], we can find the same expression for
Smarr formula. In this case, following, for example, the procedure given in
Ref. [43], we arrive at the following expression

$M=\frac{\kappa A}{4\pi}-2VP+\Phi q-\int_{v}dv\,\omega\,\,\,,$ (60)

where $A$ is the surface area of the black hole, $\kappa$ is the surface
gravity, $\omega=-1/2(\mbox{trace }T)$ is the work density, $T$ represents the
energy-momentum tensor of the considered electrodynamics model and the
integral is evaluated outside surface of the black hole. Here it can be
computed the following result

$\int_{v}dv\,\omega=Bb\,\,\,.$ (61)

Using a similar analysis as in Ref. [44], and since here too the parameter $b$
has units of electric field for all $n$, then $B$ can also be interpreted as
polarization per unit volume, or as the work due to the nonlinearity of the
electrodynamic model.

##  6 Conclusions

In this study, we have presented a family of nonlinear electrodynamics models
that depend on a parameter $b$, which in all cases represents the maximum
value reached by the electric field, and on a dimensionless parameter $n$
(positive integers), that determines the different models of electrodynamics.
Each member of this family of solutions includes the invariants
$\mathcal{F}=F_{\mu\nu}F^{\mu\nu}$ and
$\mathcal{G}=F_{\mu\nu}\tilde{F}^{\mu\nu}$. Similar to the Born-Infeld
electrodynamics, when $b\rightarrow\infty$ each model becomes Maxwell theory
and the solution with $n=1$ corresponds to the Born-Infeld model. The value of
the self-energy for charged particles is finite in all cases. We have analyzed
the causality and unitarity for our family of electrodynamic models
considering $\mathcal{F}>0$ and $\mathcal{F}<0$ separately with
$\mathcal{G}=0$. Only in the cases where $n$ is even and $\mathcal{F}>0$ the
conditions are not fully satisfied.

When considering the phenomenon of birefringence neglecting the contributions
of the higher orders in Eqs. (30) and (32), we basically obtain that all the
models quantitatively present the same birefringence. That is, except for
$n=1$, in all cases the Cotton-Mouton effect is exhibited, that is, an
external magnetic field induces birefringence, which is represented in the
following way

$\Delta
n_{CM}=n_{\parallel}-n_{\perp}\approx\frac{\overline{\mathbf{B}}^{2}}{2b^{2}}\,\,\,.$
(62)

In a similar way, we obtain that for all cases, excluding $n=1$, the Kerr
effect is present, that is, in the presence of an external electric field
birefringence is induced, which is characterized as

$\Delta
n_{K}=n_{\parallel}-n_{\perp}\approx-\frac{\overline{\mathbf{E}}^{2}}{2b^{2}}\,\,\,.$
(63)

According to Ref. [58] for the Cotton-Moutton and Kerr birefringences, $\Delta
n_{CM}=k_{CM}\mathbf{B}^{2}$ and $\Delta n_{K}=k_{K}\mathbf{E}^{2}$,
respectively, where $k_{CM}\approx 4.03\times 10^{-24}T^{-2}$ and
$k_{K}\approx-4.4\times 10^{-41}m^{2}V^{-2}$.

We have found a family of solutions for charged black holes coupling General
Relativity with the electrodynamic models that we have presented. As in the
case of a Born-Infeld black hole, we have also classified, according to their
behaviour, the different solutions obtained for each $n$. All cases exhibit
the same categories, which are defined according the horizons that each
solution presents. In the transition from a Schwarzschild-type black hole to a
Reissner-Nordström-type black hole, a solution (marginal case) is found whose
metric function is regular at the origin, as illustrated in Figs. 5, 6 and 7.
All the cases illustrated here have a horizon, except those corresponding to
$b=0.01$ in Fig. 5.

Even though it is possible to estimate its values, the idea of considering the
cosmological constant [69] or the nonlinearity parameter $b$ [48] as dynamic
variables and from this as thermodynamic quantities has been a subject of
great interest [45, 44, 70]. In this approach, the cosmological constant
contributes to the first law as a thermodynamic pressure in the term $VdP$,
and the nonlinearity parameter contributes as a term $Bdb$ that allows us to
solve the problem of finding a Smarr formula consistent with the first law.
The first of these quantities has been considered as a conserved charge in the
form $C=-\sqrt{|\Lambda|}/4\pi G$ (where the minus sign corresponds to the AdS
case) [71] with the chemical potential $\Theta$ as its thermodynamic
conjugate, contributing to the first law with the term $\Theta dC$. Perhaps
also the parameter $b$ deserves equal attention from this point of view.

## Acknowledgments

L.B. was supported by DIUFRO through the project DI22-0026.

## References

  * [1] M. Born and L. Infeld, “Foundations of the new field theory,” Proc. Roy. Soc. Lond. A 144, no. 852, 425 (1934).
  * [2] G. W. Gibbons, “Aspects of Born-Infeld theory and string / M theory,” AIP Conf. Proc. 589, no.1, 324-350 (2001).
  * [3] E. S. Fradkin and A. A. Tseytlin, “Nonlinear Electrodynamics from Quantized Strings,” Phys. Lett. B 163, 123-130 (1985).
  * [4] J. M. Davila, C. Schubert and M. A. Trejo, “Photonic processes in Born-Infeld theory,” Int. J. Mod. Phys. A 29, 1450174 (2014).
  * [5] J. Ellis, N. E. Mavromatos and T. You, “Light-by-Light Scattering Constraint on Born-Infeld Theory,” Phys. Rev. Lett. 118, no.26, 261802 (2017).
  * [6] P. Niau Akmansoy and L. G. Medeiros, “Constraining Born–Infeld-like nonlinear electrodynamics using hydrogen’s ionization energy,” Eur. Phys. J. C 78, no.2, 143 (2018).
  * [7] P. Niau Akmansoy and L. G. Medeiros, “Constraining nonlinear corrections to Maxwell electrodynamics using $\gamma\gamma$ scattering,” Phys. Rev. D 99, no.11, 115005 (2019).
  * [8] M. J. Neves, J. B. de Oliveira, L. P. R. Ospedal and J. A. Helayël-Neto, “Dispersion relations in nonlinear electrodynamics and the kinematics of the Compton effect in a magnetic background,” Phys. Rev. D 104, no.1, 015006 (2021).
  * [9] M. J. Neves, L. P. R. Ospedal, J. A. Helayël-Neto and P. Gaete, “Considerations on anomalous photon and Z-boson self-couplings from the Born–Infeld weak hypercharge action,” Eur. Phys. J. C 82, no.4, 327 (2022).
  * [10] P. De Fabritiis, P. C. Malta and J. A. Helayël-Neto, “Phenomenology of a Born-Infeld extension of the $U(1)_{Y}$ sector at lepton colliders,” Phys. Rev. D 105, no.1, 016007 (2022).
  * [11] W. Heisenberg and H. Euler, “Consequences of Dirac’s theory of positrons,” Z. Phys. 98, no.11-12, 714-732 (1936).
  * [12] S. I. Kruglov, “Vacuum birefringence from the effective Lagrangian of the electromagnetic field,” Phys. Rev. D 75, 117301 (2007).
  * [13] J. Y. Kim, “Deflection of light by a Coulomb charge in Born–Infeld electrodynamics,” Eur. Phys. J. C 81, no.6, 508 (2021).
  * [14] A. E. Shabad and V. V. Usov, “Effective Lagrangian in nonlinear electrodynamics and its properties of causality and unitarity,” Phys. Rev. D 83, 105006 (2011).
  * [15] H. H. Soleng, “Charged black points in general relativity coupled to the logarithmic U(1) gauge theory,” Phys. Rev. D 52, 6178 (1995).
  * [16] S. I. Kruglov, “A model of nonlinear electrodynamics,” Annals Phys. 353, 299-306 (2014).
  * [17] S. I. Kruglov, “Nonlinear arcsin-electrodynamics,” Annalen Phys. 527, 397-401 (2015).
  * [18] S. I. Kruglov, “Modified nonlinear model of arcsin-electrodynamics,” Commun. Theor. Phys. 66, no.1, 59-65 (2016).
  * [19] S. I. Kruglov, “Acceleration of Universe by Nonlinear Electromagnetic Fields,” Int. J. Mod. Phys. D 25, no.11, 1640002 (2016).
  * [20] S. I. Kruglov, “Magnetized black holes and nonlinear electrodynamics,” Int. J. Mod. Phys. A 32, no.23n24, 1750147 (2017).
  * [21] S. I. Kruglov, “Nonlinear Electrodynamics and Magnetic Black Holes,” Annalen Phys. 529, no.8, 1700073 (2017).
  * [22] P. Gaete and J. Helayël-Neto, “Finite Field-Energy and Interparticle Potential in Logarithmic Electrodynamics,” Eur. Phys. J. C 74, no.3, 2816 (2014).
  * [23] P. Gaete and J. A. Helayël-Neto, “A note on nonlinear electrodynamics,” EPL 119, no.5, 51001 (2017).
  * [24] P. Gaete, J. A. Helayël-Neto and L. P. R. Ospedal, “Coulomb’s law modification driven by a logarithmic electrodynamics,” EPL 125, no.5, 51001 (2019).
  * [25] S. H. Hendi, “Asymptotic charged BTZ black hole solutions,” JHEP 03, 065 (2012).
  * [26] S. H. Mazharimousavi and M. Halilsoy, “Electric Black Holes in a Model of Nonlinear Electrodynamics,” Annalen Phys. 531, no.12, 1900236 (2019).
  * [27] S. H. Mazharimousavi and M. Halilsoy, “Electric and magnetic black holes in a new nonlinear electrodynamics model,” Annals Phys. 433, 168579 (2021).
  * [28] I. Gullu and S. H. Mazharimousavi, “Double-logarithmic nonlinear electrodynamics,” Phys. Scripta 96, no.4, 045217 (2021).
  * [29] P. Gaete and J. Helayël-Neto, “Remarks on nonlinear Electrodynamics,” Eur. Phys. J. C 74, no.11, 3182 (2014).
  * [30] S. I. Kruglov, “Notes on Born–Infeld-type electrodynamics,” Mod. Phys. Lett. A 32, no.36, 1750201 (2017).
  * [31] I. Bandos, K. Lechner, D. Sorokin and P. K. Townsend, “A non-linear duality-invariant conformal extension of Maxwell’s equations,” Phys. Rev. D 102, 121703 (2020).
  * [32] S. I. Kruglov, “On generalized ModMax model of nonlinear electrodynamics,” Phys. Lett. B 822, 136633 (2021).
  * [33] B. Hoffmann, “Gravitational and electromagnetic mass in the Born-Infeld electrodynamics,” Phys. Rev. 47, no. 11, 877 (1935).
  * [34] B. Hoffmann and L. Infeld, “On the choice of the action function in the new field theory,” Phys. Rev. 51, no.9, 765-773 (1937).
  * [35] E. Ayon-Beato and A. Garcia, Regular black hole in general relativity coupled to nonlinear electrodynamics, Phys. Rev. Lett. 80, 5056 (1998).
  * [36] E. Ayon-Beato and A. Garcia, New regular black hole solution from nonlinear electrodynamics, Phys. Lett. B 464, 25 (1999).
  * [37] K. A. Bronnikov, Regular magnetic black holes and monopoles from nonlinear electrodynamics, Phys. Rev. D 63, 044005 (2001).
  * [38] I. Dymnikova, Regular electrically charged structures in nonlinear electrodynamics coupled to general relativity, Class. Quant. Grav. 21, 4417 (2004).
  * [39] L. Balart and E. C. Vagenas, Regular black hole metrics and the weak energy condition, Phys. Lett. B 730, 14 (2014).
  * [40] L. Balart and E. C. Vagenas, Regular black holes with a nonlinear electrodynamics source, Phys. Rev. D 90, no. 12, 124045 (2014).
  * [41] S. Fernando and D. Krug, “Charged black hole solutions in Einstein-Born-Infeld gravity with a cosmological constant,” Gen. Rel. Grav. 35, 129-137 (2003).
  * [42] D. A. Rasheed, “Nonlinear electrodynamics: Zeroth and first laws of black hole mechanics,” [arXiv:hep-th/9702087 [hep-th]].
  * [43] L. Balart and S. Fernando, “A Smarr formula for charged black holes in nonlinear electrodynamics,” Mod. Phys. Lett. A 32, no.39, 1750219 (2017).
  * [44] S. Gunasekaran, R. B. Mann and D. Kubiznak, “Extended phase space thermodynamics for charged and rotating black holes and Born-Infeld vacuum polarization,” JHEP 11, 110 (2012).
  * [45] D. Kastor, S. Ray and J. Traschen, Enthalpy and the Mechanics of AdS Black Holes, Class. Quant. Grav. 26, 195011 (2009).
  * [46] B. P. Dolan, “The cosmological constant and the black hole equation of state,” Class. Quant. Grav. 28, 125020 (2011).
  * [47] B. P. Dolan, “Pressure and volume in the first law of black hole thermodynamics,” Class. Quant. Grav. 28, 235017 (2011).
  * [48] W. Yi-Huan, “Energy and first law of thermodynamics for Born-Infeld-anti-de-Sitter black hole,” Chin. Phys. B 19, 090404 (2010).
  * [49] D. C. Zou, S. J. Zhang and B. Wang, “Critical behavior of Born-Infeld AdS black holes in the extended phase space thermodynamics,” Phys. Rev. D 89, no.4, 044002 (2014).
  * [50] S. H. Hendi, S. Panahiyan and B. Eslam Panah, “P–V criticality and geometrical thermodynamics of black holes with Born–Infeld type nonlinear electrodynamics,” Int. J. Mod. Phys. D 25, no.01, 1650010 (2015).
  * [51] N. Kumar, S. Bhattacharyya and S. Gangopadhyay, “Phase transitions in Born-Infeld AdS black holes in D-dimensions,” Gen. Rel. Grav. 52, no.2, 20 (2020).
  * [52] X. X. Zeng and H. Q. Zhang, “Thermodynamics and weak cosmic censorship conjecture in Born-Infeld-anti-de Sitter black holes,” Chin. Phys. C 45, no.2, 025112 (2021).
  * [53] I. Gullu and S. H. Mazharimousavi, “Black holes in double-Logarithmic nonlinear electrodynamics,” Phys. Scripta 96, no.9, 095213 (2021).
  * [54] L. Balart and S. Fernando, “Thermodynamics and heat engines of black holes with Born–Infeld-type electrodynamics,” Mod. Phys. Lett. A 36, no.15, 15 (2021).
  * [55] N. Kumar, S. Sen and S. Gangopadhyay, “Phase transition structure and breaking of universal nature of central charge criticality in a Born-Infeld AdS black hole,” Phys. Rev. D 106, no.2, 026005 (2022).
  * [56] J. Béard, J. Agil, R. Battesti and C. Rizzo, “A novel pulsed magnet for magnetic linear birefringence measurements,” Rev. Sci. Instrum. 92, no.10, 104710 (2021).
  * [57] A. Dehghani, M. R. Setare and S. Zarepour, “Self-energy problem, vacuum polarization, and dual symmetry in Born–Infeld-type U(1) gauge theories,” Eur. Phys. J. Plus 137, no.7, 859 (2022) [erratum: Eur. Phys. J. Plus 137, no.9, 1090 (2022)].
  * [58] R. Battesti and C. Rizzo, “Magnetic and electric properties of quantum vacuum,” Rept. Prog. Phys. 76, no.1, 016401 (2013).
  * [59] G. Boillat, “Nonlinear electrodynamics - Lagrangians and equations of motion,” J. Math. Phys. 11, no.3, 941-951 (1970).
  * [60] Z. Bialynicka-Birula and I. Bialynicki-Birula, “Nonlinear effects in Quantum Electrodynamics. Photon propagation and photon splitting in an external field,” Phys. Rev. D 2, 2341-2345 (1970).
  * [61] P. Gaete and J. A. Helayël-Neto, “Remarks on inverse electrodynamics,” Eur. Phys. J. C 81, no.10, 899 (2021).
  * [62] P. Gaete, “Some Remarks on Nonlinear Electrodynamics,” Adv. High Energy Phys. 2016, 2463203 (2016).
  * [63] D. P. Sorokin, “Introductory Notes on Non-linear Electrodynamics and its Applications,” Fortsch. Phys. 70, no.7-8, 2200092 (2022).
  * [64] S. I. Kruglov, “Nonlinear electrodynamics with birefringence, Phys. Lett. B 379, 623-625 (2015).
  * [65] S. Fernando, “Thermodynamics of Born-Infeld-anti-de Sitter black holes in the grand canonical ensemble,” Phys. Rev. D 74, 104032 (2006).
  * [66] N. Bretón, T. Clark and S. Fernando, “Quasinormal modes and absorption cross-sections of Born–Infeld–de Sitter black holes,” Int. J. Mod. Phys. D 26, no.10, 1750112 (2017).
  * [67] L. Smarr, “Mass formula for Kerr black holes,” Phys. Rev. Lett. 30, 71 (1973). Erratum: [Phys. Rev. Lett. 30, 521 (1973)].
  * [68] A. Komar, “Covariant conservation laws in general relativity,” Phys. Rev. 113, 934-936 (1959).
  * [69] C. Teitelboim, “The cosmological constant as a thermodynamic black hole parameter,” Phys. Lett. B 158, 293-297 (1985).
  * [70] D. Kubiznak, R. B. Mann and M. Teo, “Black hole chemistry: thermodynamics with Lambda,” Class. Quant. Grav. 34, no.6, 063001 (2017).
  * [71] D. Chernyavsky and K. Hajian, “Cosmological constant is a conserved charge,” Class. Quant. Grav. 35, no.12, 125012 (2018).


# Orbital Stability and Secular Dynamics of the Proxima Centauri Planetary
System

Joseph R. Livesey Department of Astronomy, University of Washington, Seattle,
WA 98195  Department of Astronomy, University of Wisconsin–Madison, Madison,
WI 53706  Rory Barnes Department of Astronomy, University of Washington,
Seattle, WA 98195  NASA Virtual Planetary Laboratory Lead Team, USA  Russell
Deitrick School of Earth and Ocean Sciences, University of Victoria, Victoria,
BC V8P 3E6

###### Abstract

The two innermost planets of the Proxima Centauri system are separated by just
0.02 au, inducing strong gravitational interactions between them. We assess
this interaction by leveraging fast orbital stability indicators and find that
orbital stability is very likely if the initial eccentricities of planets b
and d are less than $\sim 0.2$, but cannot confirm stability at larger values.
We find that stability is not strongly affected by the true masses of the
planets or by the distant planet c. However, mutual inclinations between
$95^{\circ}$ and $142^{\circ}$ often result in unstable motion. We further
explore the long-term evolution of the orbits in these stable regions of
parameter space and find that circularization can take over 5 Gyr. This tidal
evolution could support surface energy fluxes in excess of 1 W m-2 for over 1
Gyr, possibly affecting planet b’s habitability.

Astrobiology (74), Exoplanet dynamics (490), $N$-body simulations (1083)

††software: Rebound (Rein & Liu, 2012), VPLanet (Barnes et al., 2020), Astropy
(Astropy Collaboration et al., 2022), Matplotlib (Hunter, 2007), Numpy (Harris
et al., 2020), Pandas (Pandas Development Team, 2020; McKinney, 2010), Seaborn
(Waskom, 2021).

##  1 Introduction

Proxima Centauri is the nearest stellar neighbor to the Sun and is known to
host at least 3 planets, including one $\sim$Earth-mass planet in the
habitable zone (Anglada-Escudé et al., 2016; Damasso et al., 2020; Suárez
Mascareño et al., 2020; Faria et al., 2022). Given the system’s proximity and
the possible habitability of the middle planet, b, this system is uniquely
valuable in the study of exoplanet formation, evolution, and habitability.
Herein we describe a series of experiments to test the orbital stability of
the Proxima Centauri system over a range of plausible orbital architectures.

The orbital evolution of Proxima’s planetary system could affect the
habitability of Proxima b in many ways. Orbital oscillations could affect the
seasons and trigger snowball states (Spiegel et al., 2009; Armstrong et al.,
2014; Deitrick et al., 2018b; Wilhelm et al., 2022), while tidal damping can
drive long term decay of orbital semi-major axis and eccentricity (Barnes et
al., 2008; Barnes, 2017). Additionally, tidal friction in Proxima b’s interior
could create large amounts of energy that affect the generation of magnetic
fields (Driscoll & Barnes, 2015), vigorous volcanic activity (Jackson et al.,
2008; Barnes et al., 2009), and possibly even drive a runaway greenhouse
(Barnes et al., 2013). Thus, defining the physically plausible range of
orbital architectures in this system could significantly improve our
understanding of Proxima b’s habitability.

Not surprisingly, then, this study is not the first to consider the stability
of this system. In particular, Meng et al. (2019) estimated the parameter
space that permits orbital stability, but they used an older architecture that
included only planets b and c. In this study, we update their results by
considering the currently observed planetary system.

In particular, the recently detected planet d (Faria et al., 2022) demands a
reanalysis of this system’s stability. The orbits of Proxima d and b are
sufficiently near so that gravitational perturbations between the two could
drive the two into instability, i.e., the ejection of a planet from the
system. Such cases should therefore be deemed implausible starting conditions
for the system, i.e., evolutionary models can safely avoid those unstable
parameter combinations. Using semi-analytic and $N$-body simulations, we
systematically searched for unstable configurations of this system. We also
identify architectures that give rise to chaotic dynamics, which do not
necessarily imply instability. Identifying these chaotic configurations is
useful for investigating long-term behavior, as the dynamics of systems
exhibiting orbital chaos are ill-described by a disturbing function containing
only secular terms (e.g., that employed in DistOrb; see §2.5 and Barnes et
al., 2020).

In §2, we outline the theoretical models we use to assess stability and the
parameter space we will explore. In §3, we present stability maps of the
system obtained through thousands of simulations varying the initial
conditions. In §4, we discuss the major outcomes of this work and implications
for future studies of this planetary system. Note that code for reproducing
the results presented in this paper is available
online.111https://github.com/jrlivesey/ProximaStability

##  2 Methods

In this section we describe theories of orbital stability, the numerical
methods we employed to identify unstable cases, and the range of initial
conditions we explored. Throughout this study we consider only the three-body
system consisting of the host star and planets d and b. We tested the role of
planet c, with an orbital period of $\sim$ 1,900 days, and found that it has
negligible impact on the stability of the system, so we will not present
detailed results of its role.

Proxima Centauri is a member of the $\alpha$ Centauri triple star system,
orbiting the binary $\alpha$ Centauri AB at a distance of about 8700 au
(Kervella et al., 2017). Asteroseismological analysis has revealed that the
binary is $5.3\pm 0.3$ Gyr old (Joyce & Chaboyer, 2018). We assume that
Proxima Centauri formed in the same cluster as $\alpha$ Centauri AB and
therefore has the same age. While Proxima’s planetary system could be modified
by the binary (Barnes et al., 2016), we will ignore that possibility in our
studies.

###  2.1 Stability Criteria

In this study we evaluate the Lagrange stability of the Proxima system. In a
Lagrange stable planetary system, all planets remain gravitationally bound to
the star and the ordering of the planets is preserved. There is no analytic
formula to determine whether a system is Lagrange stable, and so we resort to
$N$-body simulations.

We first consider, however, the possibility that the three-body system is Hill
stable (Szebehely & McKenzie, 1977; Marchal & Bozis, 1982; Gladman, 1993;
Barnes & Greenberg, 2006), which is determined analytically, but can only
approximate orbital stability. The definition of Hill stability is that the
ordering of the bodies cannot change. In other words, if the outer planet
escapes to infinity, the system would still be Hill stable. Although this
definition is very limited, it does reduce to a single equation that can be
quickly evaluated. In the present study, we assess the agreement between the
“Hill stability boundary” of the system in orbital element space and the
stability boundary obtained through simulations. Such an analysis has been
performed previously for systems of two $\sim$Jupiter-mass planets (Barnes &
Greenberg, 2006, 2007), but not, to our knowledge, for systems of terrestrial
planets.

A system is guaranteed to be Hill stable if it satisfies the inequality
$\beta>\beta_{\text{crit}}$ for

$\beta=-\frac{2M}{G^{2}\mathcal{M}^{3}}CL^{2}$ (1)

and

$\beta_{\text{crit}}\simeq
1+\frac{3^{4/3}m_{1}m_{2}}{m_{\star}^{2/3}(m_{1}+m_{2})^{4/3}}-\frac{m_{1}m_{2}(11m_{1}+7m_{2})}{3m_{\star}(m_{1}+m_{2})^{2}},$
(2)

where $m_{\star}$, $m_{1}$, and $m_{2}$ are the masses of the star, inner
planet, and outer planet, respectively; $M=\sum_{i}m_{i}$ is the total mass of
the system; $\mathcal{M}=(1/2)\sum_{i\neq j}m_{i}m_{j}$; $C$ is the total
orbital energy of the system; and $L$ is its total angular momentum (Marchal &
Bozis, 1982). A three-body system for which $\beta/\beta_{\text{crit}}<1$ is
likely but not certain to be Hill unstable.

To evaluate Hill stability, we generated $10^{6}$ possible orbital
configurations and calculated $\beta/\beta_{\text{crit}}$. We randomly sampled
the longitudes of ascending node, arguments of periastron, and mean anomalies
of both planets from $\mathcal{U}[0,360)^{\circ}$. The inclinations of the
planets were sampled from $\mathcal{U}[0,15]^{\circ}$. The inclination of the
planets relative to the sky plane, used to calculate their true masses, was
sampled from $\mathcal{U}[10,90]^{\circ}$. For both of the orbital
eccentricities, we used $e\sim\mathcal{U}[0,0.9]$. We assumed the orbits were
coplanar for these calculations.

###  2.2 MEGNO

A hybrid approach to estimating orbital stability is to compute a quantity
called the MEGNO (Mean Exponential Growth of Nearby Orbits; Cincotta & Simó,
2000), which measures the divergence of two trajectories in phase space. If
$\boldsymbol{\delta}(t)$ denotes the separation between two initially
infinitesimally close particles in phase space, then the MEGNO $Y$ is defined
as

$Y(t)=\frac{2}{t-t_{0}}\int_{t_{0}}^{t}\frac{\dot{{\bf\it\delta}}\dotproduct{\bf\it\delta}}{{\bf\it\delta}\dotproduct{\bf\it\delta}}t^{\prime}\>\differential
t^{\prime}$ (3)

and its time-averaged value as

$\langle
Y\rangle(t)=\frac{1}{t-t_{0}}\int_{t_{0}}^{t}Y(t^{\prime})\>\differential
t^{\prime}.$ (4)

For periodic (stable) evolution, $\langle Y\rangle\to 2$ as $t\to\infty$. For
a chaotic system, $\langle Y\rangle\to(\gamma/2)t$, where $\gamma$ is the
maximum Lyapunov characteristic number (Cincotta & Simó, 2000; Goździewski et
al., 2001; Morbidelli, 2002; Cincotta et al., 2003; Hinse et al., 2010). The
MEGNO is defined similarly to the Lyapunov characteristic number, but is
weighted with time, amplifying stochastic behavior in the evolution and
allowing for an earlier detection of chaos (Cincotta & Simó, 2000; Goździewski
et al., 2001; Morbidelli, 2002). We follow standard practice, and label the
motion as chaotic if $|\langle Y\rangle-2|>10^{-2}$, indicating the system may
be unstable on Gyr timescales.

###  2.3 General Relativistic Effects

Relativistic (GR) effects can modify the precessional frequency of the orbit
by an amount

$\dot{\varpi}_{\text{GR}}\simeq\frac{3Gm_{\star}n}{c^{2}a(1-e^{2})}=\frac{3(Gm_{\star})^{3/2}}{c^{2}a^{5/2}(1-e^{2})},$
(5)

where $n$ is the mean motion of the planet in question (Sterne, 1939). If we
apply this relation to the most extreme case studied here, i.e., that with
$e=0.9$ for both planets, we obtain GR precession periods of about 22,000
years and 82,000 years for planets d planet b, respectively. Since previous
work has shown that $N$-body integrations up to $10^{6}$ orbits of the
outermost body are sufficient to identify most unstable configurations of a
planetary system (Barnes & Quinn, 2004). We therefore run our simulations for
about $10^{6}$ orbits of planet b, which is $\sim$ 31,000 yr. Therefore, GR
will only induce one extra circulation over the length of our simulation and
should therefore not be a significant effect. To decrease computational cost,
we therefore do not include a post-Newtonian correction to our model.

By the same token, we may safely assume that the rate of precession due to the
Lense–Thirring effect is also negligible (Will & Nordtvedt, 1972; Iorio,
2008).

###  2.4 N-Body Simulations

Table 1 lists the current best fits and uncertainties for the physical and
orbital elements of the system. We assume the host star has a mass of
$0.12{M_{\odot}}$ (Kervella et al., 2017). However, we are also interested in
the range of plausible architectures at the end of planet formation as tidal
damping can cause long-term changes to the planetary orbits (e.g., Barnes et
al., 2008; Barnes, 2017; Meng et al., 2019). We will thus consider a wider
parameter space than allowed by observational uncertainties, focusing on
larger eccentricities and semi-major axes.

Table 1: Best-fit parameters for the Proxima Centauri system. Body | $a$ | $e$ | $m\sin i_{\text{obs}}$ | $i_{\text{obs}}$
---|---|---|---|---
| (au) |  |  | (deg)
staraaKervella et al. (2017) |  |  | $0.1221\pm 0.0022{M_{\odot}}$ |
Planet dbbFaria et al. (2022, TV RM analysis) | $0.02885^{+0.00019}_{-0.00022}$ | $0.04^{+0.15}_{-0.04}$ | $0.26\pm 0.05{M_{\oplus}}$ |  $\cdots$
Planet bbbFaria et al. (2022, TV RM analysis) | $0.04856^{+0.00030}_{-0.00030}$ | $0.02^{+0.04}_{-0.02}$ | $1.07\pm 0.06{M_{\oplus}}$ |  $\cdots$
Planet cccDamasso et al. (2020) | $1.48^{+0.08}_{-0.08}$ | $0.04^{+0.01}_{-0.01}$ | $5.8\pm 1.9{M_{\oplus}}$ |  $133\pm 1$ddBenedict & McArthur (2020)

To gauge the dependence of the system’s stability on the system architecture,
we performed three sets of 10,000 $N$-body simulations with the publicly
available REBOUND code222Available at https://github.com/hannorein/rebound.
Our simulations use v3.28.3. (Rein & Liu, 2012; Rein & Tamayo, 2015; Rein &
Spiegel, 2015; Rein & Tamayo, 2016). In Set I we varied the initial
eccentricities of both planets; in Set II we varied the planets’ masses; and
in Set III we varied the initial inclinations with respect to the fundamental
plane.

In addition to evaluating the MEGNO, in our Set I simulations we identify the
parameter space in which the orbits cross, and in which we therefore
anticipate chaos. As there are two planets included in our $N$-body
simulations, planet d’s apocenter coincides with planet b’s pericenter, and
thus the orbits cross where $a_{b}(1-e_{b})<a_{d}(1+e_{d})$. We also
incorporate the two-planet chaos criterion given by Hadden & Lithwick (2018)
in terms of the relative eccentricity, which in our case simplifies somewhat
to

$\displaystyle e_{d}+e_{b}$
$\displaystyle>\left(\frac{a_{b}-a_{d}}{a_{d}}\right)$
$\displaystyle\quad\times\exp\left[-\frac{11}{5}\left(\frac{m_{d}+m_{b}}{m_{\star}}\right)^{1/3}\left(\frac{a_{b}}{a_{b}-a_{d}}\right)^{4/3}\right].$
(6)

In our Set II simulations, we test the stability of the system with
combinations of masses up to $10{M_{\oplus}}$ for both planets. Previous
studies have specified a critical orbital separation in a two-planet system
with low $e$ and $i$, at which the system becomes unstable. The critical
separation is $\Delta a=2\sqrt{3}R_{\text{H}}$, where

$R_{\text{H}}=\frac{a+a^{\prime}}{2}\left(\frac{m+m^{\prime}}{3m_{\star}}\right)^{1/3}$
(7)

is the Hill radius (Hasegawa & Nakazawa, 1990; Gladman, 1993; Chambers et al.,
1996). Primed quantities refer to the outer body. In the most extreme case
considered here, both masses are $10{M_{\oplus}}$, with which $\Delta
a=36.3R_{\text{H}}$. We therefore expect every experiment in Set II to result
in stable motion; we perform them and include the results here for
completeness.

The mutual inclination (or relative inclination) $\Psi$ between the two
planets’ orbits is given by

$\cos(\Psi)=\sin(i)\sin(i^{\prime})\cos(\Omega-\Omega^{\prime})+\cos(i)\cos(i^{\prime}).$
(8)

Note that because there are only two planets included in our $N$-body
simulations, $\cos(\Omega-\Omega^{\prime})\simeq-1$ for our purposes. The
resulting distribution of simulated cases is weighted toward higher values of
$\Psi$, since nearly coplanar, prograde orbits are already expected to be
stable.

For each set of simulations, we identify the Hill stability boundary in the
relevant parameter space as an additional analytic stability criterion. We
compute this boundary for each set because the Hill stability parameter
depends upon the eccentricities, masses, and inclinations of both planets
according to Eq. (1).

The simulations are integrated for $\sim 10^{6}$ orbits of Proxima b. This
timescale is longer than the conventional integration times used for
identifying chaotic motion with the MEGNO, which are typically $\sim 10^{3-4}$
times the longest dynamical timescale in the system (Cincotta & Simó, 2000;
Hinse et al., 2010; Meng et al., 2019). We perform our integrations over this
longer timescale because we find some cases that are surprisingly stable (see
Section 3.1).

The initial conditions for all of our REBOUND simulations are provided in
Table 2. The semi-major axes used are slightly greater than the measured
values to account for tidal effects over time. The approximate mass values
used in these simulations for both planets were obtained by setting
$i_{\text{obs}}=133^{\circ}$, the mean observational inclination of Proxima c
obtained astrometrically by Benedict & McArthur (2020).

Table 2: Initial conditions for the REBOUND simulations. Set | Planet | $a$ | $e$ | $m\sin i_{\text{obs}}$ | $i$ | $\Omega$ | $\varpi$
---|---|---|---|---|---|---|---
|  | (au) |  | (${M_{\oplus}}$) | (deg) | (deg) | (deg)
I | d | 0.029 | 0–0.9 | 0.26 | 0 | 0 | 0
| b | 0.049 | 0–0.9 | 1.07 | 0 | 180 | 180
II | d | 0.029 | 0 | 0.21–10 | 0 | 0 | 0
| b | 0.049 | 0 | 1.01–10 | 0 | 180 | 180
III | d | 0.029 | 0 | 0.26 | 0–180 | 0 | 0
| b | 0.049 | 0 | 1.07 | 0–180 | 180 | 180

We performed two varieties of $N$-body simulations with REBOUND. We first used
WHFAST, a symplectic Wisdom–Holman integrator (see Rein & Tamayo, 2015; Wisdom
& Holman, 1991) to integrate the simulations. WHFAST is fast and efficient,
but does not fare well when the planets’ orbits are strongly perturbed (i.e.,
when they become significantly non-Keplerian). In cases where the relative
energy error exceeds $10^{-6}$, we re-run the simulations using IAS15, a 15th
order Gauss–Radau integrator (Rein & Spiegel, 2015). IAS15 employs adaptive
timestepping to resolve close encounters. At the end of each simulation, we
find $\langle Y\rangle$ using built-in functions in REBOUND that employ
variational equations to calculate the trajectories of nearby orbits (Rein &
Tamayo, 2016).

###  2.5 Secular Orbital Evolution

The $N$-body experiments discussed in the previous section yield a set of
points in initial parameter space that give rise to quasi-periodic motion,
where $\langle Y\rangle\to 2$, and are thus suitable starting points for long-
term simulations of the secular evolution.

We performed six such simulations using the modular planetary evolution code
VPLanet333Available at https://github.com/VirtualPlanetaryLaboratory/ vplanet.
Our simulations use v2.3.28. (Barnes et al., 2020). In particular, we used the
EqTide and DistOrb modules to simultaneously simulate evolution of the orbits
due to equilibrium stellar tides (the “constant phase lag” model of Goldreich
& Soter, 1966) and due to the disturbing function parameterizing gravitational
interactions between the planets (Murray & Dermott, 1999). For detailed
descriptions and validation of these models, the reader is directed to Barnes
et al. (2020) and Deitrick et al. (2018a).

Each body in the system is given a value of $k_{2}$ and $Q$. $k_{2}$ is the
Love number of degree 2 and $Q$ is the “tidal quality factor” that
parameterizes the rate of tidal dissipation in a body’s interior. Earth has
$Q=12$ (Williams et al., 1978), but a terrestrial planet without large oceans
may have a quality factor an order of magnitude larger. In each of our
simulations, $k_{2}=0.5$ for the star, $k_{2}=0.3$ for the planets. The star’s
$Q$ is set to $10^{6}$, while the planet’s $Q$ is sampled in the range
$\mathcal{U}[10,10^{3}]$. Planetary radii are necessary for calculating the
tidal torques, and are determined according to the power law $R\propto
m^{0.274}$ (Sotin et al., 2007). We can only estimate these radii, as no
planet in this system is transiting (Jenkins et al., 2019). The values we use
fall within the range established by the examination of Proxima b’s internal
structure given in Brugger et al. (2016).

In addition to evaluating the time dependence of the orbital elements, we
calculate the tidal heating of planet b throughout each evolution according to
$\dot{E}_{\text{tide}}=-(\dot{E}_{\text{orb}}+\dot{E}_{\text{rot}})$, where

$\dot{E}_{\text{orb}}=\frac{Z}{8}\left[4\zeta-\left(20\zeta+76\right)e^{2}\right]$
(9)

and

$\dot{E}_{\text{rot}}=\frac{Z}{8}\frac{\omega}{n}\left[4\zeta-\left(20\zeta+48\right)e^{2}\right]$
(10)

are the power dissipated from the orbit and rotation of the planet,
respectively.444Our Equations 9 and 10 are simplified from the CPL tidal
heating equations in Barnes et al. (2020, see their Appendix E). We set
$\omega/n=1$ at the start of each simulation, and assume that $\omega/n\simeq
1$ throughout. Thus, we take $\operatorname{sgn}(2\omega-3n)=-1$,
$\operatorname{sgn}(2\omega-n)=+1$ and $\operatorname{sgn}(n)=+1$. Proxima b
has zero obliquity in each of our simulations, and so we discard terms related
to obliquity tides. Here, $n$ and $\omega$ are the orbital and rotational
frequencies, $\zeta=\operatorname{sgn}(2\omega-2n)$, and

$Z=3G^{2}k_{2}m_{\star}^{2}(m_{\star}+m)\frac{R^{5}}{a^{9}}\frac{1}{nQ}.$ (11)

Initial conditions were sampled randomly from the permitted region of
parameter space using VSPACE, a tool for initializing VPLanet’s parameter
sweeps. They are presented in Table 4. In these simulations, we include planet
c and assume for the sake of simplicity that Proxima c’s orbit begins with
zero inclination and does not evolve tidally. Therefore, values of $R$,
$\Omega$, and $Q$ are not provided for this planet.

Table 3: Distributions from which the initial conditions of our VPLanet simulations are sampled. Planet | Quantity | Distribution
---|---|---
d | $P_{\text{orb}}$ |  $\mathcal{U}(5.122,5.4)$ days
| $Q$ | $\mathcal{N}(100,250)$
| $\Omega$ |  $\mathcal{U}(0,360)$ deg
| $\omega$ |  $\mathcal{U}(0,360)$ deg
b | $P_{\text{orb}}$ |  $\mathcal{U}(11.187,11.5)$ days
| $Q$ | $\mathcal{N}(100,250)$
| $\Omega$ |  $\mathcal{U}(0,360)$ deg
| $\omega$ |  $\mathcal{U}(0,360)$ deg
c | $P_{\text{orb}}$ |  $\mathcal{N}(1902,24)$ days
| $\Omega$ |  $\mathcal{U}(0,360)$ deg
| $\omega$ |  $\mathcal{U}(0,360)$ deg
Table 4: Initial conditions for the VPLanet simulations, numbered 1–6. Simulation | Planet | $P$ | $e$ | $m$ | $R$ | $i$ | $\Omega$ | $\varpi$ | $Q$
---|---|---|---|---|---|---|---|---|---
|  | (days) |  | (${M_{\oplus}}$) | (${R_{\oplus}}$) | (deg) | (deg) | (deg) |
1 | d | 5.23736 | 0.12679 | 0.51653 | 0.83443 | 5.35445 | 214.29919 | 301.19731 | 389.31695
| b | 11.49553 | 0.13508 | 2.18208 | 1.23837 | 8.52253 | 278.74258 | 287.98729 | 85.28869
| c | 1,881.12305 | 0.01 | 6 |  | 0 |  | 375.5932 |
2 | d | 5.1691 | 0.09071 | 0.39862 | 0.777237 | 5.35445 | 273.62635 | 211.36366 | 322.54608
| b | 11.3132 | 0.0925 | 1.71255 | 1.158826 | 8.52253 | 308.83944 | 221.8431 | 466.03801
| c | 1,881.91953 | 0.01 | 6 |  | 0 |  | 71.65173 |
3 | d | 5.24638 | 0.12754 | 0.35092 | 0.75056 | 2.0301 | 30.27576 | 62.72857 | 250.40924
| b | 11.27027 | 0.0711 | 1.23215 | 1.05887 | 8.76222 | 53.66543 | 195.10852 | 220.5822
| c | 1,885.90255 | 0.01 | 6 |  | 0 |  | 136.16203 |
4 | d | 5.34536 | 0.17036 | 0.442 | 0.80 | 2.0301 | 7.32011 | 72.85377 | 277.98179
| b | 11.32216 | 0.09628 | 1.13222 | 1.03461 | 8.76222 | 130.82045 | 197.2157 | 646.02112
| c | 1,914.12162 | 0.01 | 6 |  | 0 |  | 360.19573 |
5 | d | 5.37027 | 0.18371 | 0.33433 | 0.74067 | 2.0301 | 283.00882 | 211.43098 | 243.22927
| b | 11.4511 | 0.12222 | 1.23139 | 1.05869 | 8.76222 | 264.70184 | 242.43391 | 15.94921
| c | 1,892.59296 | 0.01 | 6 |  | 0 |  | 329.46207 |
6 | d | 5.35688 | 0.17726 | 0.32236 | 0.73331 | 2.0301 | 62.43061 | 24.80976 | 452.43231
| b | 11.33074 | 0.08986 | 1.37201 | 1.09053 | 8.76222 | 174.8711 | 299.8717 | 124.37175
| c | 1,866.94917 | 0.01 | 6 |  | 0 |  | 185.42062 |

##  3 Results

###  3.1 Orbital Stability

We find that the orbits of the two inner planets are Hill stable in just 7.2%
of the parameter space we explored. The Hill stable configurations are those
with the lowest eccentricities for both planets and the boundary is shown with
the black curve in Figure 1.

Figure 1: Panel (a) shows a stability map of Proxima d and b, with both
eccentricities varying between 0 and 0.9. Blue cells indicate periodic
evolution according to the MEGNO criterion, red cells indicate one planet was
ejected, and white cells indicate chaotic (but not necessarily unstable)
evolution. The solid black curve marks the Hill stability boundary in this
plane, with Hill stable cases lying below the curve. The dashed-dotted black
line indicates the limit of crossing orbits, i.e., above this line the
apocenter distance of planet d is initially greater than the pericenter
distance of planet b. The dotted line indicates the stability limit of Hadden
& Lithwick (2018). Yellow stars indicate the simulations shown in detail in
Figure 2. Green dots indicate the locations of the secular dynamics
simulations detailed in Section 2.5. Grey triangles indicate the locations of
the longer (1 Myr) integrations described in Section 3.1. Panel (b) shows a
MEGNO map of the same parameter space. Most points in this map are similarly
colored, as these simulations concluded with $\langle Y\rangle\geq 10$ or with
an ejection, in which case an arbitrarily high MEGNO value was assigned.

The results of Set I are also shown, with colors denoting values of $\langle
Y\rangle$. Blue points indicate periodic (stable) configurations ($|\langle
Y\rangle-2|\leq 10^{-2}$), white points are chaotic ($|\langle
Y\rangle-2|>10^{-2}$, but both planets remained bound), and red points are
cases in which one of the planets escaped to infinity. A clear trend is
present in which larger eccentricities lead to decreasing stability. There are
many cases in which $\langle Y\rangle<2$ — or even becomes negative — at the
end of a simulation. Such trajectories lie in the neighborhood of stable
periodic orbits (e.g., Breiter et al., 2005). There are cases in which the
variational equations fail to converge at the end of a simulation, in which
case REBOUND returns no value for the MEGNO. This outcome means that the
separation in phase space has grown very large, and thus we label these cases
as chaotic without having an exact MEGNO value.

We find that most periodic configurations lie within the Hill stability
boundary, with $e_{d}\lesssim 0.25$ and $e_{b}\lesssim 0.2$. Most
configurations with $e_{d}\lesssim 0.55$ and $e_{b}\lesssim 0.45$ give rise to
motion that is chaotic, but does not become unstable within the length of
these simulations. Beyond these boundaries, the system is very likely to be
unstable, with one of the planets escaping to infinity within the length of
these simulations.

We also find “islands” of configurations at higher initial eccentricities that
are mostly classified as chaotic, but with some periodic results, e.g., near
$(0.55,0.4)$ and $(0.85,0.55)$. The observational uncertainties suggest the
system cannot currently exist in these Hill unstable islands (see Table 1). We
return to this possibility in §4.

The outlying stable regions shown in Figure 1 could correspond to mean motion
resonances, however, we find no evidence that MMRs are responsible for
stabilizing these configurations. For a selection of simulated high-$e$ stable
cases, we calculated the time evolution of resonant arguments, which are

$\phi=j_{1}\lambda^{\prime}+j_{2}\lambda+j_{3}\varpi^{\prime}+j_{4}\varpi$
(12)

for the $j_{1}:j_{2}$ MMR, where the coefficients satisfy the d’Alembert
relation $\sum_{i}j_{i}=0$. We found no libration in $\phi$ up to $j_{1}=30$.
We also re-ran five cases for 1 Myr, and found that all cases initially
designated as periodic were reclassified as chaotic (i.e., $\langle
Y\rangle\not\to 2$). We therefore conclude that these outlying islands are in
fact unstable and that the MEGNO method was not able to identify instabilities
over the first $10^{6}$ orbits. This conjecture is supported by the fact that
these anomalous cases lie in the crossing orbits region of parameter space,
and thus chaotic evolution is expected. Note that simulating the system for
its age would require at least 5,000 CPU hours per simulation with WHFAST on a
modern workstation, and so is currently intractable.

Figure 2 shows the eccentricity evolution for three example cases: a periodic
case, a chaotic case, and an unstable case. The locations of these cases in
Figure 1 are marked by stars. Recall that Set I assumes coplanar orbits, so
there is no evolution with respect to the inclinations. The stable case shows
classic periodic motion in which both eccentricities oscillate with just one
frequency. The chaotic case appears regular, but contains at least two
frequencies, probably resulting from the larger eccentricities. The unstable
case ejects planet d in less than 250 years.

Figure 2: Example evolutions of Proxima d and b. Periodic evolution
originating from the initial conditions $e_{d}=e_{b}=0.05$ is shown in panel
(a). An evolution originating from $e_{d}=0.3$ and $e_{b}=0.2$ is shown in
panel (b). The time-averaged MEGNO for this case diverges significantly from
2, signifying that the motion is chaotic. At least 3 distinct frequencies are
seen to be present in both planets’ eccentricities, which is greater than the
number predicted by secular theory. An unstable case originating from
$e_{d}=0.45$ and $e_{b}=0.3$ is shown in panel (c). The ejection of planet d
from the system is evident at $t=222$ yr, where $e_{d}>1$.

Next we turn to Set II, which focused on the role of planetary masses on
stability. In each of these simulations, the initial orbits are circular and
coplanar. We found that every configuration we tested resulted in periodic
motion. We conclude that stability is not strongly influenced by the planetary
masses, which is consistent with the $m^{1/3}$ dependence in the Hill
stability criterion.

Finally, we present results for Set III, which focused on the inclinations. As
shown in Figure 3, we find that the motion is periodic for $\Psi<50^{\circ}$
and $\Psi>150^{\circ}$. With few exceptions, systems in between evolve
chaotically. 80% of systems with $95^{\circ}<\Psi<142^{\circ}$ resulted in
ejections. This finding is consistent with previous results that show
stability increases as the orbits approach coplanarity (e.g., Meng et al.,
2019). Our obtained distribution of chaotic configurations in mutual
inclination space, with a wide chaotic region in the middle skewed toward
greater $\Psi$, agrees with the stability analyses of two-planet systems
performed in Veras & Ford (2010) in terms of the maximum eccentricity
variation. Our Set III data are presented again in two dimensions — $i_{d}$
and $i_{b}$ — in Figure 4. Figure 5 shows the eccentricity evolution for three
example cases of various initial mutual inclination.

Figure 3: MEGNO values of the system over all possible initial mutual
inclinations. The dashed black line indicates $\langle Y\rangle=2$. Blue dots
are cases with periodic or quasi-periodic motion, grey dots are cases with
chaotic motion, and the red ticks indicate values of $\Psi$ for which the
system is unstable (and therefore any MEGNO value is meaningless). All cases
but one between 50 and $140^{\circ}$ resulted in chaotic motion. 165
simulations resulted in ejections, all with mutual inclinations between 95 and
$142^{\circ}$. Figure 4: Similar to Figure 1, but with the initial
inclinations varied instead of eccentricities. The black lines in panel (a)
indicate the Hill stability boundary. Yellow stars indicate the simulations
presented in detail in Figure 5, and the green dots point out two combinations
of inclinations used in the secular dynamics simulations. Figure 5: Sample
evolutions of Proxima d and b, beginning with circular orbits. Panel (a) shows
periodic evolution, where initially $\Psi=3^{\circ}$. Panel (b) shows chaotic
evolution, where initially $\Psi=60^{\circ}$. Panel (c) shows an unstable case
originating from $\Psi=110^{\circ}$. The ejection of planet d from the system
is evident at $\sim 63$ kyr.

###  3.2 Long-term Orbital Evolution

Using initial orbital configurations that give rise to regular motion, we
simulated the long-term orbital evolution of the system. One result is shown
in Figure 6. Notably, in this scenario tidal evolution persists longer than
the 7 Gyr for which the simulation was run, and much longer than the expected
value for the age of the system.

The evolution of Proxima d and b in terms of $e$ and $a$ in each of the
VPLanet simulations is displayed in Figure 7. Five of the six secular
evolutions result in tidal heating that remains above the value for modern
Earth for up to 7 Gyr. In the remaining simulation, Proxima b has $Q\simeq
16$. This result suggests that if Proxima b has $Q\sim 10$ (on the order of
Earth), then the surface flux of tidal heating is not significant after the
first few Gyr because the orbit has circularized. If, however, $Q\sim 10^{2}$,
then tidal interactions render a first-order effect on the heat delivered to
the world’s surface over its entire history.

Figure 6: Secular orbital evolution of all three planets in Simulation 1. In
this case, we use $e_{d}=0.127$, $e_{b}=0.135$, $m_{d}=0.517{M_{\oplus}}$,
$m_{b}=2.182{M_{\oplus}}$, and $\Psi=7.859^{\circ}$. Figure 7: Secular orbital
evolution of Proxima d in each simulated case is displayed in panels (a) and
(b), and Proxima b in panels (c) and (d), with the tidal heat flux through the
planet b’s surface shown in panel (e). Three simulations end with a tidal heat
flux of $\sim 0.1$ W m-2, two end with a near-Earth value of $\sim 0.01$ W
m-2, and one with a far lower value. The present tidal heat flux at Io’s and
Earth’s surfaces are indicated for reference (Moore et al., 2007; Munk &
Wunsch, 1998).

##  4 Discussion

Our stability analysis has revealed key constraints in the orbital properties
of Proxima d and b. In particular, eccentricities in excess of 0.5 are
strongly disfavored, as are mutual inclinations between 95 and 142∘. We found
cases in which the system is not Hill stable, and yet exhibits quasi-periodic
motion. However, most quasi-periodic configurations of the system reside
within the Hill stability boundary, as we should expect. Our analysis also
identified regions where our simulations did not produce ejections, despite
having crossing orbits initially. Although these configurations appear stable,
a subset of simulations re-run for a longer integration time results in
instabilities, suggesting that most of the region where
$\beta<\beta_{\text{crit}}$ is in fact unstable.

The separated islands of stable motion in Figure 1, however, are very likely
inaccessible. The large region of instability separating them from the low
eccentricity “mainland” is probably too large for the system to successfully
cross. As shown in Figure 2, instabilities can arise within a few centuries,
which is short compared to the tidal damping timescale. We thus conclude that
Proxima b and d could only have formed in the low eccentricity region of
stability.

Our simulations of secular orbital evolution in the system reveal that,
beginning from the “stability mainland” of parameter space determined with our
$N$-body integrations, the tidal orbital evolution of Proxima b can easily
continue beyond 5 Gyr. In each case simulated with VPLanet, we see the
eccentricity behaves as a damped harmonic oscillator, with secular amplitude
decreasing rapidly with time. As such, we find no new constraint on the
present-day eccentricity of this world, which has a value at $+1\sigma$ of
0.06 according to Faria et al. (2022). This result suggests that tidal heating
in Proxima b is likely ongoing, providing an additional source of energy to
drive surface and atmospheric processes, perhaps including active biology. On
the other hand, the increased volcanism could outgas water that is ultimately
destroyed by high energy photons, resulting in a desiccated, and therefore
uninhabitable world (Luger & Barnes, 2015). Regardless, our research shows
that tidal heating should not be ignored for planet b unless the upper limit
of its eccentricity can be significantly reduced.

Models of the evolution and habitability of Proxima b can use our results to
target their studies in the physically allowed regions we have identified in
Figures 1 and 4. Climate studies (e.g. Turbet et al., 2016; Del Genio et al.,
2017; Colose et al., 2021) could consider eccentricities for Proxima b up to
0.45, but should probably focus on the blue triangle in the lower left corner
of Figure 1. Studies of the planetary interior could also consider the role of
tidal heating over the same range (see e.g., Henning et al., 2009; Barnes et
al., 2013; Driscoll & Barnes, 2015; Barnes et al., 2016). These results may
also prove useful in dynamical modeling of the Proxima system from radial
velocity data, providing limits on the maximum achievable orbital
eccentricities and inclinations of planets d and b.

## Acknowledgments

J.R.L. gratefully acknowledges partial support from a Mary Gates Research
Scholarship. R.B. acknowledges support from NASA grant number 80NSSC20K0229
and the NASA Virtual Planetary Laboratory under grant number 80NSSC18K0829.
The authors thank Gongjie Li, Hanno Rein, and Juliette Becker for valuable
insights and discussion, and thank the anonymous referee whose thoughtful
recommendations substantially improved this manuscript. This work was
facilitated through the use of advanced computational, storage, and networking
infrastructure provided by the Hyak supercomputer system and funded by the
Student Technology Fund at the University of Washington.

## References

  * Anglada-Escudé et al. (2016) Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437, doi: 10.1038/nature19106
  * Armstrong et al. (2014) Armstrong, J. C., Barnes, R., Domagal-Goldman, S., et al. 2014, Astrobiology, 14, 277, doi: 10.1089/ast.2013.1129
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Barnes (2017) Barnes, R. 2017, Celestial Mechanics and Dynamical Astronomy, 129, 509, doi: 10.1007/s10569-017-9783-7
  * Barnes & Greenberg (2006) Barnes, R., & Greenberg, R. 2006, ApJ, 647, L163, doi: 10.1086/507521
  * Barnes & Greenberg (2007) —. 2007, ApJ, 665, L67, doi: 10.1086/521144
  * Barnes et al. (2009) Barnes, R., Jackson, B., Greenberg, R., & Raymond, S. N. 2009, ApJ, 700, L30, doi: 10.1088/0004-637X/700/1/L30
  * Barnes et al. (2013) Barnes, R., Mullins, K., Goldblatt, C., et al. 2013, Astrobiology, 13, 225, doi: 10.1089/ast.2012.0851
  * Barnes & Quinn (2004) Barnes, R., & Quinn, T. 2004, ApJ, 611, 494, doi: 10.1086/421321
  * Barnes et al. (2008) Barnes, R., Raymond, S. N., Jackson, B., & Greenberg, R. 2008, Astrobiology, 8, 557, doi: 10.1089/ast.2007.0204
  * Barnes et al. (2016) Barnes, R., Deitrick, R., Luger, R., et al. 2016.  https://arxiv.org/abs/1608.06919
  * Barnes et al. (2020) Barnes, R., Luger, R., Deitrick, R., et al. 2020, Publications of the Astronomical Society of the Pacific, 132, 1, doi: 10.1088/1538-3873/ab3ce8
  * Benedict & McArthur (2020) Benedict, G. F., & McArthur, B. E. 2020, Research Notes of the American Astronomical Society, 4, 86, doi: 10.3847/2515-5172/ab9ca9
  * Breiter et al. (2005) Breiter, S., Wytrzyszczak, I., & Melendo, B. 2005, Advances in Space Research, 35, 1313, doi: 10.1016/j.asr.2005.02.033
  * Brugger et al. (2016) Brugger, B., Mousis, O., Deleuil, M., & Lunine, J. I. 2016, ApJ, 831, L16, doi: 10.3847/2041-8205/831/2/L16
  * Chambers et al. (1996) Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996, Icarus, 119, 261, doi: 10.1006/icar.1996.0019
  * Cincotta et al. (2003) Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica D Nonlinear Phenomena, 182, 151, doi: 10.1016/S0167-2789(03)00103-9
  * Cincotta & Simó (2000) Cincotta, P. M., & Simó, C. 2000, A&AS, 147, 205, doi: 10.1051/aas:2000108
  * Colose et al. (2021) Colose, C. M., Haqq-Misra, J., Wolf, E. T., et al. 2021, ApJ, 921, 25, doi: 10.3847/1538-4357/ac135c
  * Damasso et al. (2020) Damasso, M., Sordo, F. D., Anglada-Escudé, G., et al. 2020, Science Advances, 6, 1, doi: 10.1126/sciadv.aax7467
  * Deitrick et al. (2018a) Deitrick, R., Barnes, R., Quinn, T. R., et al. 2018a, AJ, 155, 60, doi: 10.3847/1538-3881/aaa301
  * Deitrick et al. (2018b) Deitrick, R., Barnes, R., Bitz, C., et al. 2018b, AJ, 155, 266, doi: 10.3847/1538-3881/aac214
  * Del Genio et al. (2017) Del Genio, A. D., Way, M. J., Amundsen, D. S., et al. 2017, arXiv e-prints, arXiv:1709.02051.  https://arxiv.org/abs/1709.02051
  * Driscoll & Barnes (2015) Driscoll, P. E., & Barnes, R. 2015, Astrobiology, 15, 739, doi: 10.1089/ast.2015.1325
  * Faria et al. (2022) Faria, J. P., Suárez Mascareño, A., Figueira, P., et al. 2022, A&A, 658, A115, doi: 10.1051/0004-6361/202142337
  * Gladman (1993) Gladman, B. 1993, Icarus, 106, 247, doi: 10.1006/icar.1993.1169
  * Goldreich & Soter (1966) Goldreich, P., & Soter, S. 1966, Icarus, 5, 375, doi: 10.1016/0019-1035(66)90051-0
  * Goździewski et al. (2001) Goździewski, K., Bois, E., Maciejewski, A. J., & Kiseleva-Eggleton, L. 2001, A&A, 378, 569, doi: 10.1051/0004-6361:20011189
  * Hadden & Lithwick (2018) Hadden, S., & Lithwick, Y. 2018, AJ, 156, 95, doi: 10.3847/1538-3881/aad32c
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hasegawa & Nakazawa (1990) Hasegawa, M., & Nakazawa, K. 1990, A&A, 227, 619
  * Henning et al. (2009) Henning, W. G., O’Connell, R. J., & Sasselov, D. D. 2009, ApJ, 707, 1000, doi: 10.1088/0004-637X/707/2/1000
  * Hinse et al. (2010) Hinse, T. C., Christou, A. A., Alvarellos, J. L. A., & Goździewski, K. 2010, MNRAS, 404, 837, doi: 10.1111/j.1365-2966.2010.16307.x
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Iorio (2008) Iorio, L. 2008, Scholarly Research Exchange, 2008, 5235, doi: 10.3814/2008/105235
  * Jackson et al. (2008) Jackson, B., Barnes, R., & Greenberg, R. 2008, MNRAS, 391, 237, doi: 10.1111/j.1365-2966.2008.13868.x
  * Jenkins et al. (2019) Jenkins, J. S., Harrington, J., Challener, R. C., et al. 2019, MNRAS, 487, 268, doi: 10.1093/mnras/stz1268
  * Joyce & Chaboyer (2018) Joyce, M., & Chaboyer, B. 2018, ApJ, 864, 99, doi: 10.3847/1538-4357/aad464
  * Kervella et al. (2017) Kervella, P., Thévenin, F., & Lovis, C. 2017, A&A, 598, L7, doi: 10.1051/0004-6361/201629930
  * Luger & Barnes (2015) Luger, R., & Barnes, R. 2015, Astrobiology, 15, 119, doi: 10.1089/ast.2014.1231
  * Marchal & Bozis (1982) Marchal, C., & Bozis, G. 1982, Celestial Mechanics, 26, 311, doi: 10.1007/BF01230725
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a
  * Meng et al. (2019) Meng, T., Ji, J., & Dong, Y. 2019, MNRAS, 482, 372, doi: 10.1093/mnras/sty2682
  * Moore et al. (2007) Moore, W. B., Schubert, G., Anderson, J. D., & Spencer, J. R. 2007, in Io After Galileo: A New View of Jupiter’s Volcanic Moon, ed. R. M. C. Lopes & J. R. Spencer, 89, doi: 10.1007/978-3-540-48841-5_5
  * Morbidelli (2002) Morbidelli, A. 2002, Modern celestial mechanics: Aspects of solar system dynamics (CRC Press)
  * Munk & Wunsch (1998) Munk, W., & Wunsch, C. 1998, Deep Sea Research Part I: Oceanographic Research, 45, 1977, doi: 10.1016/S0967-0637(98)00070-3
  * Murray & Dermott (1999) Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge University Press)
  * Pandas Development Team (2020) Pandas Development Team. 2020, pandas-dev/pandas: Pandas, latest, Zenodo, doi: 10.5281/zenodo.3509134
  * Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Rein & Spiegel (2015) Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424, doi: 10.1093/mnras/stu2164
  * Rein & Tamayo (2015) Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376, doi: 10.1093/mnras/stv1257
  * Rein & Tamayo (2016) —. 2016, MNRAS, 459, 2275, doi: 10.1093/mnras/stw644
  * Sotin et al. (2007) Sotin, C., Grasset, O., & Mocquet, A. 2007, Icarus, 191, 337, doi: 10.1016/j.icarus.2007.04.006
  * Spiegel et al. (2009) Spiegel, D. S., Menou, K., & Scharf, C. A. 2009, ApJ, 691, 596, doi: 10.1088/0004-637X/691/1/596
  * Sterne (1939) Sterne, T. E. 1939, MNRAS, 99, 451, doi: 10.1093/mnras/99.5.451
  * Suárez Mascareño et al. (2020) Suárez Mascareño, A., Faria, J. P., Figueira, P., et al. 2020, Astronomy & Astrophysics, 639, doi: 10.1051/0004-6361/202037745
  * Szebehely & McKenzie (1977) Szebehely, V., & McKenzie, R. 1977, AJ, 82, 79, doi: 10.1086/112011
  * Turbet et al. (2016) Turbet, M., Leconte, J., Selsis, F., et al. 2016, A&A, 596, A112, doi: 10.1051/0004-6361/201629577
  * Veras & Ford (2010) Veras, D., & Ford, E. B. 2010, ApJ, 715, 803, doi: 10.1088/0004-637X/715/2/803
  * Waskom (2021) Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021, doi: 10.21105/joss.03021
  * Wilhelm et al. (2022) Wilhelm, C., Barnes, R., Deitrick, R., & Mellman, R. 2022, , 3, 13, doi: 10.3847/PSJ/ac3b61
  * Will & Nordtvedt (1972) Will, C. M., & Nordtvedt, Kenneth, J. 1972, ApJ, 177, 757, doi: 10.1086/151754
  * Williams et al. (1978) Williams, J. G., Sinclair, W. S., & Yoder, C. F. 1978, Geophys. Res. Lett., 5, 943, doi: 10.1029/GL005i011p00943
  * Wisdom & Holman (1991) Wisdom, J., & Holman, M. 1991, AJ, 102, 1528, doi: 10.1086/115978