# Unified Interacting Quark Stars in Regularized 4D Einstein Gauss-Bonnet
Gravity

Michael Gammon<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Waterloo, Waterloo, Ontario, Canada, N2L 3G1  Sarah Rourke
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of
Waterloo, Waterloo, Ontario, Canada, N2L 3G1  Department of Physics, McGill
University, Montreal, Quebec, Canada  Robert B. Mann<EMAIL_ADDRESS>Department of Physics and Astronomy, University of Waterloo, Waterloo,
Ontario, Canada, N2L 3G1

###### Abstract

Since the derivation of a well-defined $D\rightarrow 4$ limit for 4D Einstein
Gauss-Bonnet (4DEGB) gravity hennigar_2020_on  coupled to a scalar field,
there has been interest in testing it as an alternative to Einstein’s general
theory of relativity. Using the Tolman-Oppenheimer-Volkoff (TOV) equations
modified for 4DEGB gravity, we model the stellar structure of quark stars
using a novel interacting quark matter equation of state zhang_2021_unified .
We find that increasing the Gauss-Bonnet coupling constant $\alpha$ or the
interaction parameter $\lambda$ both tend to increase the mass-radius profiles
of quark stars described by this theory, allowing a given central pressure to
support larger quark stars in general. These results logically extend to cases
where $\lambda<0$, in which increasing the magnitude of the interaction
effects instead diminishes masses and radii. We also analytically identify a
critical central pressure in both regimes, below which no quark star solutions
exist due to the pressure function having no roots. Most interestingly, we
find that quark stars can exist below the general relativistic Buchdahl bound
and Schwarzschild radius $R=2M$, due to the lack of a mass gap between black
holes and compact stars in 4DEGB. Even for small $\alpha$ well within current
observational constraints, we find that quark star solutions in this theory
can describe Extreme Compact Objects (ECOs), objects whose radii are smaller
than what is allowed by general relativity.

##  I Introduction

Modified theories of gravity continue to attract attention despite the
empirical success of general relativity (GR). Such theories are motivated by a
variety of problems, including addressing issues in in modern cosmology,
quantizing gravity, eliminating singularities, and, perhaps most importantly,
finding viable phenomenological competitors against which GR can be tested in
the most stringent manner possible.

Higher curvature theories (or HCTs) are amongst the most popular
modifications. An HCT modifies the assumed linear relationship in GR between
the curvature and the stress-energy, replacing the former with an arbitrary
sum of powers of the curvature tensor (appropriately contracted to two
indices). Such modifications could conceivably improve the empirical success
of GR, while also making new testable predictions.

Many quantum gravity proposals ahmed2017  and cosmological puzzles (such as
dark energy, dark matter, and early-time inflation bueno2016 ) suggest that
HCTs could play an important role in physics. Lovelock theories lovelock1971
have long been at the forefront of this search, since they posses the
distinctive feature of having 2nd order differential equations of motion. The
physical significance of such theories has been unclear, however, since their
higher order terms yield non-trivial contributions to the equations of motion
only in more than four spacetime dimensions ($D>4$).

Recently this restriction was circumvented in the quadratic case, or what is
better known as “Einstein-Gauss-Bonnet” gravity. The Gauss-Bonnet (GB)
contribution to the gravitational action is

$S_{D}^{GB}=\alpha\int
d^{D}x\sqrt{-g}\left[R^{\mu\nu\rho\tau}R_{\mu\nu\rho\tau}-4R^{\mu\nu}R_{\mu\nu}+R^{2}\right]\equiv\alpha\int
d^{D}x\sqrt{-g}\mathcal{G}$ (1)

(where $R_{\mu\nu\rho\tau}$ is the Riemann curvature tensor), which becomes
the integral of a total derivative in $D=4$, and thus cannot contribute to a
system’s gravitational dynamics in less than five dimensions. For this reason
it is often referred to as a “topological term” having no relevance to
physical problems. Indeed, the Lovelock theorem lovelock1971  ensures that a
$D=4$ dimensional metric theory of gravity must incorporate additional fields
in order to have second order equations of motion and diffeomorphism
invariance.

Recently it was noted glavan2020  that several exact solutions to
$D$-dimensional Einstein-Gauss Bonnet gravity have a sensible limit under the
rescaling

$\lim_{D\to 4}(D-4)\alpha\rightarrow\alpha,$ (2)

of the Gauss-Bonnet coupling constant. Using this approach a variety of
4-dimensional metrics can be obtained, including cosmological glavan2020 ;
li2020 ; kobayashi2020 , spherical black hole glavan2020 ; kumar2022 ;
fernandes2020 ; kumar2020 ; kumar2022_2 , collapsing malafarina2020 , star-
like doneva2021 ; charmousis2022 , and radiating ghosh2020  metrics, each
carrying imprints of the quadratic curvature effects of their $D>4$
counterparts. However a number of objections to this approach were
subsequently raised gurses2020 ; ai2020 ; shu2020 , based on the fact that the
existence of a limiting solution does not imply the existence of a well-
defined 4D theory whose field equations have that solution. This shortcoming
was quickly addressed when it was shown that the $D\to 4$ limit in (2) can be
taken in the gravitational action hennigar_2020_on ; Fernandes:2020nbq ,
generalizing an earlier procedure employed in obtaining the $D\to 2$ limit of
GR Mann:1992ar . One can also compactify $D$-dimensional Gauss-Bonnet gravity
on a $(D-4)$-dimensional maximally symmetric space and then use (2) to obtain
a $D=4$ HCT Lu:2020iav . This approach yields the same result (up to trivial
field redefinitions), in addition to terms depending on the curvature of the
maximally symmetric $(D-4)$-dimensional space. Taking this to vanish yields

$\displaystyle S_{4}^{GB}=$ $\displaystyle\alpha\int
d^{4}x\sqrt{-g}\left[\phi\mathcal{G}+4G_{\mu\nu}\nabla^{\mu}\phi\nabla^{\nu}\phi-4(\nabla\phi)^{2}\square\phi+2(\nabla\phi)^{4}\right]$
(3)

where we see that an additional scalar field $\phi$ appears. Surprisingly, the
spherically symmetric black hole solutions to the field equations match those
from the naïve $D\rightarrow 4$ limit of solutions glavan2020 . The resultant
4D scalar-tensor theory is a particular type of Horndeski theory horndeski1974
, and solutions to its equations of motion can be obtained without ever
referencing a higher dimensional spacetime hennigar_2020_on .

We are interested here in what is called 4D Einstein-Gauss-Bonnet gravity
(4DEGB), whose action is given by (4) plus the Einstein-Hilbert term:

$\displaystyle S=S^{GR}+S_{4}^{GB}=$ $\displaystyle\int
d^{4}x\sqrt{-g}\left[R+\alpha\left\\{\phi\mathcal{G}+4G_{\mu\nu}\nabla^{\mu}\phi\nabla^{\nu}\phi-4(\nabla\phi)^{2}\square\phi+2(\nabla\phi)^{4}\right\\}\right]$
(4)

which has been shown to be an interesting phenomenological competitor to GR
Clifton:2020xhc . Despite much exploration Fernandes:2022zrq , it is unclear
whether these higher curvature terms play an important role in real
gravitational dynamics. One important arena for testing such theories against
standard general relativity is via observations of compact astrophysical
objects like neutron stars. The correct theory should be able to accurately
describe recent gravitational wave observations of astrophysical objects
existing in the mass gap between the heaviest compact stars and the lightest
black holes.

Modern observational astrophysics is rich in findings of compact objects and
as such our understanding of highly dense gravitational objects is rapidly
advancing. However there is as of yet no strong consensus on their underlying
physics. A number of such objects have been recently observed that are
inconsistent with standard GR and a simple neutron star equation of state. It
was recently shown zhang_2021_unified ; zhang_2021_stellar  that in standard
general relativity, the secondary component of the merger GW190814 could
feasibly be a quark star with an interacting equation of state governed by a
single parameter $\lambda$. This parameterization of strong interaction
effects was inspired by another recent theoretical study showing that non-
strange quark matter could feasibly be the ground state of baryonic matter at
sufficient density and temperature holdom_2018_quark . Similar analyses with a
different equation of state and/or QM phase miao_2021_bayesian ;
lopes_2022_onthe ; oikonomou_2023_colourflavour  found similarly promising
results. This same object was subsequently shown to be well described as a
slowly-rotating neutron star in the 4DEGB theory without resorting to exotic
quark matter EOSs, while also demonstrating that the equilibrium sequence of
neutron stars asymptotically matches the black hole limit, thus closing the
mass gap between NS/black holes of the same radius charmousis2022 . More
recently, some groups have also been interested in modelling ECOs
zhang_2023_rescaling  as well as unusually light compact stars
horvath_2023_alight ; oikonomou_2023_colourflavour  (like that in the gamma
ray remnant J1731-347, which is inconsistent with minimum mass calculations of
neutron stars generated by iron cores) as quark stars to explain their unusual
properties.

To further illuminate the range of possibilities, we consider in this paper
quark star solutions to the 4DEGB theory. Inspired by the unified interacting
quark matter equation of state derived in zhang_2021_unified , we present a
simple model for interacting quark stars in the regularized 4DEGB theory in
which all of the strong interaction effects (including corrections from
perturbative QCD, colour superconductivity, and the strange quark mass) are
characterized by a single parameter ${\lambda}$ in the equation of state. In
doing so, we solve the modified TOV equations for a number of combinations of
Gauss-Bonnet coupling and QM interaction strengths (constrained by
observational limits). Although quark star solutions have previously been
considered in the context of 4DEGB banerjee_2021_strange ; banerjee_2021_quark
; pretel_2022 , a different (less general) equation of state was used, and the
upper limit of the coupling $\alpha$ was taken to be much smaller than that
allowed by current observational constraints Fernandes:2022zrq ;
Clifton:2020xhc  and the 4DEGB Buchdahl bound chakraborty2020_limits . In
considering a unified, interacting equation of state zhang_2021_unified  and
values of $\alpha$ up to presently allowed bounds, we obtain a number of
interesting novel results.

Our most intriguing result is that quark stars in 4DEGB can be Extreme Compact
Objects (ECOs), objects whose radii are smaller than that allowed by the
Buchdahl bound in GR. Indeed, there exist quark stars in 4DEGB whose radii are
smaller than that of a corresponding black hole of the same mass in GR111This
phenomenon was also shown to be present for neutron stars charmousis2022 ,
though the relationship with the Buchdahl bound was not noted.. Observations
of these latter objects, apart from indicating a new class of astrophysical
phenomena Mann:2021mnc , would provide strong evidence for 4DEGB as a physical
theory. These ECOs respect a generalization of the Buchdahl bound, whose
small-radius limit is that of the horizon radius of the corresponding black
hole.

We find in general that for a given central pressure, quark stars in 4DEGB
have larger mass and radius than their GR counterparts with the same pressure.
This can be attributed to the ‘less attractive’ nature of gravity in 4DEGB:
for $\alpha>0$ the gradient of the effective potential yields a weaker force
than is the case in GR. We consequently find a larger maximal mass for a given
value of $\lambda$ in 4DEGB than in GR.

Surprisingly we also find a critical central pressure for certain
$\alpha/\lambda$ combinations, below which quark stars cannot form, and above
which the mass and radius are greater than or equal to those of a black hole
in the 4DEGB theory. With this, it is clear that for such parameter
combinations no stable quark star solutions exist. We analytically derive an
expression for this critical quantity in terms of the couplings $\alpha$ and
$\lambda$. Moreover, we find that both increasing the 4DEGB coupling constant
$\alpha$ and increasing the interaction parameter $\lambda$ make the mass-
radius profiles of the stars larger (and consequently lead to a larger maximum
mass point for a given parameter set). The converse is also true for smaller
values of these parameters, with the trend continuing as $\lambda<0$.

The outline of our paper is as follows: In section II we introduce the basic
theory underlying 4DEGB gravity, as well as the unified, interacting quark
matter equation of state that we make use of. Following this, a perfect fluid
stress-energy tensor is employed to derive the 4DEGB TOV equations, and
current observational constraints on the coupling constant are briefly
discussed. Section III outlines the results of our numerical calculations,
with the first part of the section covering $\lambda>0$ solutions (positive
interaction strength), and the second part seeing how these results change
when $\lambda<0$. We conclude our results with a brief analysis of the
stability of unified, interacting 4DEGB quark stars. Section V summarizes our
key findings and suggests topics for future study.

##  II Theory

###  II.1 4D Einstein-Gauss-Bonnet Gravity

The field equations of 4DEGB are obtained from a straightforward variational
principle applied to the action (4). Variation with respect to the scalar
$\phi$ yields

$\displaystyle\mathcal{E}_{\phi}=$
$\displaystyle-\mathcal{G}+8G^{\mu\nu}\nabla_{\nu}\nabla_{\mu}\phi+8R^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi-8(\square\phi)^{2}+8(\nabla\phi)^{2}\square\phi+16\nabla^{a}\phi\nabla^{\nu}\phi\nabla_{\nu}\nabla_{\mu}\phi$
(5)
$\displaystyle\qquad+8\nabla_{\nu}\nabla_{\mu}\phi\nabla^{\nu}\nabla^{\mu}\phi$
$\displaystyle=$ $\displaystyle\;0$

and the variation with respect to the metric gives

$\displaystyle\mathcal{E}_{\mu\nu}$ $\displaystyle=\Lambda
g_{\mu\nu}+G_{\mu\nu}+\alpha\left[\phi
H_{\mu\nu}-2R\left[\left(\nabla_{\mu}\phi\right)\left(\nabla_{\nu}\phi\right)+\nabla_{\nu}\nabla_{\mu}\phi\right]+8R_{(\mu}^{\sigma}\nabla_{\nu)}\nabla_{\sigma}\phi+8R_{(\mu}^{\sigma}\left(\nabla_{\nu)}\phi\right)\left(\nabla_{\sigma}\phi\right)\right.$
(6)
$\displaystyle-2G_{\mu\nu}\left[(\nabla\phi)^{2}+2\square\phi\right]-4\left[\left(\nabla_{\mu}\phi\right)\left(\nabla_{\nu}\phi\right)+\nabla_{\nu}\nabla_{\mu}\phi\right]\square\phi-\left[g_{\mu\nu}(\nabla\phi)^{2}-4\left(\nabla_{\mu}\phi\right)\left(\nabla_{\nu}\phi\right)\right](\nabla\phi)^{2}$
$\displaystyle+8\left(\nabla_{(\mu}\phi\right)\left(\nabla_{\nu)}\nabla_{\sigma}\phi\right)\nabla^{\sigma}\phi-4g_{\mu\nu}R^{\sigma\rho}\left[\nabla_{\sigma}\nabla_{\rho}\phi+\left(\nabla_{\sigma}\phi\right)\left(\nabla_{\rho}\phi\right)\right]+2g_{\mu\nu}(\square\phi)^{2}$
$\displaystyle-4g_{\mu\nu}\left(\nabla^{\sigma}\phi\right)\left(\nabla^{\rho}\phi\right)\left(\nabla_{\sigma}\nabla_{\rho}\phi\right)+4\left(\nabla_{\sigma}\nabla_{\nu}\phi\right)\left(\nabla^{\sigma}\nabla_{\mu}\phi\right)$
$\displaystyle\left.-2g_{\mu\nu}\left(\nabla_{\sigma}\nabla_{\rho}\phi\right)\left(\nabla^{\sigma}\nabla^{\rho}\phi\right)+4R_{\mu\nu\sigma\rho}\left[\left(\nabla^{\sigma}\phi\right)\left(\nabla^{\rho}\phi\right)+\nabla^{\rho}\nabla^{\sigma}\phi\right]\right]$
$\displaystyle=\;T_{\mu\nu}$

where

$\displaystyle
H_{\mu\nu}=2\Big{[}RR_{\mu\nu}-2R_{\mu\alpha\nu\beta}R^{\alpha\beta}+R_{\mu\alpha\beta\sigma}R_{\nu}^{\alpha\beta\sigma}-2R_{\mu\alpha}R_{\nu}^{\alpha}-\frac{1}{4}g_{\mu\nu}\mathcal{G}\Big{]}$
(7)

is the Gauss-Bonnet tensor. These field equations satisfy the following
relationship

$g^{\mu\nu}T_{\mu\nu}=g^{\mu\nu}\mathcal{E}_{\mu\nu}+\frac{\alpha}{2}\mathcal{E}_{\phi}=4\Lambda-R-\frac{\alpha}{2}\mathcal{G}$
(8)

which can act as a useful consistency check to see whether prior solutions
generated via the Glavin/Lin method are even possible solutions to the theory.
For example, using (8) it is easy to verify that the rotating metrics
generated from a Newman-Janis algorithm kumar2020 ; wei2020  are not solutions
to the field equations of the scalar-tensor 4DEGB theory.

###  II.2 Unified Interacting Quark Matter Equation of State

Compact stars described in terms of deconfined quark degrees of freedom are
often modelled using the simple, non-interacting quark matter equation of
state glendenningbook

$p(r)=\frac{1}{3}(\rho(r)-4B_{\mathrm{eff}})$ (9)

where $p$ and $\rho$ are pressure and mass density (respectively), and
$B_{\mathrm{eff}}$ is the effective bag constant from the MIT bag model for
quark confinement. The bag constant values associated with conventional non-
interacting strange quark matter (SQM) or up-down quark matter (udQM) are not
small enough to explain the large compact star masses found in recent binary
merger events using GR alone zhang2020 ; ren2020 ; zhou_2018 .

Accordingly, Zhang and Mann zhang_2021_unified  investigated whether a
strongly interacting equation of state would allow a quark star to fit the
observational constraints in the context of GR. This EOS was inspired by a
recent theoretical development holdom_2018_quark  showing that udQM could be
the grount state of baryonic matter at sufficient density and temperature. In
employing this, they also provided a unified framework for all strongly
interacting phases of quark matter, condensing the entirety of strong
interaction effects into a single parameter $\lambda$. This greatly simplifies
the problem as it is no longer necessary to solve the same equations for
different phases of quark matter (namely the two-flavour superconducting
without/with strange quarks (2SC, 2SC+s), and colour-flavour locking (CFL)
phases), different values of the gap parameter from colour superconductivity,
different perturbative contributions etc. as all these effects are unified in
the model. They found that the interaction strength parameter can be written
explicitly as zhang_2021_unified

$\lambda=\frac{\xi_{2a}\Delta^{2}-\xi_{2b}m_{s}^{2}}{\sqrt{\xi_{4}a_{4}}}$
(10)

where $\Delta$ is the gap parameter, $m_{s}$ accounts for corrections from the
finite strange quark mass (if applicable), $a_{4}$ represents the perturbative
Quantum Chromodynamics (pQCD) contribution from one-gluon exchange, and the
constant coefficients

$\left(\xi_{4},\xi_{2a},\xi_{2b}\right)=\begin{cases}\left(\left(\left(\frac{1}{3}\right)^{\frac{4}{3}}+\left(\frac{2}{3}\right)^{\frac{4}{3}}\right)^{-3},1,0\right)&\text{
2SC phase }\\\ (3,1,3/4)&\text{ 2SC+s phase }\\\ (3,3,3/4)&\text{ CFL phase
}\end{cases}$ (11)

account for the various possible QM phases. It is easy to see from here how
different combinations of parameters may result in the same unified
interaction strength. The parameter $\lambda$ takes into account quark star
physics through its inclusion in the following unified, interacting QM EOS
zhang_2021_unified :

$p(r)=\frac{1}{3}\left(\rho(r)-4B_{\text{eff
}}\right)+\frac{4\lambda^{2}}{9\pi^{2}}\left(-1+\operatorname{sgn}(\lambda)\sqrt{1+3\pi^{2}\frac{\left(\rho(r)-B_{\text{eff
}}\right)}{\lambda^{2}}}\right)$ (12)

where $\rm{sgn}(\lambda)$ represents the sign of $\lambda$. This can be
further generalized by dividing out the $B_{\rm{eff}}$ dependence, leaving us
with fully dimensionless equations. Performing the rescalings

$\bar{\rho}=\frac{\rho}{4B_{\mathrm{eff}}},\bar{p}=\frac{p}{4B_{\mathrm{eff}}},$
(13)

and

$\bar{\lambda}=\frac{\lambda^{2}}{4B_{\text{eff
}}}=\frac{\left(\xi_{2a}\Delta^{2}-\xi_{2b}m_{s}^{2}\right)^{2}}{4B_{\text{eff
}}\xi_{4}a_{4}},$ (14)

allows us to rewrite the EOS (12) in terms of these dimensionless parameters
as

$\bar{p}(\bar{r})=\frac{1}{3}(\bar{\rho}(\bar{r})-1)+\frac{4}{9\pi^{2}}\bar{\lambda}\left(-1+\operatorname{sgn}(\lambda)\sqrt{1+\frac{3\pi^{2}}{\bar{\lambda}}\left(\bar{\rho}(\bar{r})-\frac{1}{4}\right)}\right)$
(15)

where

$\bar{m}=m\sqrt{4B_{\mathrm{eff}}}\qquad\bar{r}=r\sqrt{4B_{\mathrm{eff}}}\qquad\bar{\alpha}=\alpha\cdot
4B_{\text{eff}}$ (16)

and in the sequel we shall assume a characteristic value of
$B_{\mathrm{eff}}=60\;\rm{MeV/fm^{3}}$.

In the limit $\bar{\lambda}\rightarrow 0$, (15) reduces back to the expected
non-interacting EOS $\bar{p}(\bar{r})=\frac{1}{3}(\bar{\rho}(\bar{r})-1)$,
whereas in the limit of extreme positive interaction strength
($\lambda\to+\infty$), (15) approaches the form

$\left.\bar{p}\right|_{\bar{\lambda}\rightarrow\infty}=\bar{\rho}-\frac{1}{2}$
(17)

or, equivalently, $p(r)=\rho(r)-2B_{\mathrm{eff}}$. In effect, this means that
the strong interaction can reduce the surface mass density of the quark star
from $\rho_{0}=4B_{\mathrm{eff}}$ to $\rho_{0}=2B_{\mathrm{eff}}$, and
increase the speed of sound in QM from $\frac{1}{3}c$ to $c$ maximally. On the
other hand, when $\lambda$ is negative no well-defined large magnitude limit
exists.

###  II.3 4DEGB TOV Equations

The standard Tolman-Oppenheimer-Volkoff (TOV) equations for stellar structure
are well-known in GR. To model the structure of a quark star in 4DEGB gravity,
these relations need to be re-derived. In the following we do this, starting
with a static, spherically symmetric metric ansatz in natural units ($G=c=1$):

$ds^{2}=-e^{2\Phi(r)}c^{2}dt^{2}+e^{2\Lambda(r)}dr^{2}+r^{2}d\Omega^{2}.$ (18)

As usual hennigar_2020_on ; gammon_2022 , so long as $e^{2\Phi}=e^{2\Lambda}$
outside the star, the combination $\mathcal{E}_{0}^{0}-\mathcal{E}_{1}^{1}$ of
the field equations can be used to derive the following equation for the
scalar field:

$\left(\phi^{\prime
2}+\phi^{\prime\prime}\right)\left(1-\left(r\phi^{\prime}-1\right)^{2}e^{2\Lambda}\right)=0.$
(19)

which, apart from the irrelevant $\phi=\ln\left(\frac{r-r_{0}}{l}\right)$
(with $r_{0}$ and $l$ being constants of integration), has the solution

$\phi_{\pm}=\int\frac{e^{\Lambda}\pm 1}{r}dr$ (20)

where $\phi_{-}$ falls off as as $\frac{1}{r}$ in asymptotically flat
spacetimes. Choosing $\phi=\phi_{-}$ ensures that (5) is automatically
satisfied.

Modelling the quark matter by a perfect fluid matter source, the stress-energy
tensor is

$T_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu},$ (21)

from which we obtain the equations

$\displaystyle\frac{2}{r}\frac{d\Lambda}{dr}=e^{2\Lambda}\left[\frac{8\pi
G}{c^{4}}\rho-\frac{1-e^{-2\Lambda}}{r^{2}}\left(1-\frac{\alpha\left(1-e^{-2\Lambda}\right)}{r^{2}}\right)\right]\left[1+\frac{2\alpha\left(1-e^{-2\Lambda}\right)}{r^{2}}\right]^{-1},$
(22) $\displaystyle\frac{2}{r}\frac{d\Phi}{dr}=e^{2\Lambda}\left[\frac{8\pi
G}{c^{4}}P+\frac{1-e^{-2\Lambda}}{r^{2}}\left(1-\frac{\alpha\left(1-e^{-2\Lambda}\right)}{r^{2}}\right)\right]\left[1+\frac{2\alpha\left(1-e^{-2\Lambda}\right)}{r^{2}}\right]^{-1},$
(23) $\displaystyle\frac{dp}{dr}=-(\rho+p)\frac{d\Phi}{dr},$ (24)

in 4DEGB, matching those found previously banerjee_2021_quark .

Asymptotic flatness as a condition means that
$\Phi(\infty)=\Lambda(\infty)=0$, and regularity at the center of the star
implies $\Lambda(0)=0$. If, in the usual way, we define the gravitational mass
$m$ within a sphere of radius $r$ through the relation hennigar_2020_on

$e^{-2\Lambda}=1+\frac{r^{2}}{2\alpha}\left(1-\sqrt{1+\frac{8\alpha
m(r)}{r^{3}}}\right)$ (25)

we arrive at the 4DEGB modified TOV equations, namely

$\displaystyle\frac{dp}{dr}=\frac{(p+\rho)\left[r^{3}(\Gamma+8\pi\alpha
p-1)-2\alpha m\right]}{r^{2}\Gamma\left[r^{2}(\Gamma-1)-2\alpha\right]}$ (26)
$\displaystyle\frac{dm}{dr}=4\pi r^{2}\rho$ (27)

along with the EOS (15), where $\Gamma=\sqrt{1+\frac{8\alpha m}{r^{3}}}$. Note
that this expression differs from the equivalent equation in
banerjee_2021_quark ; banerjee_2021_strange  due to our inclusion of the 4DEGB
coupling in the definition of gravitational mass, however we do find agreement
with pretel_2022 .

The vacuum solution is given by $m(r)=M$, where $M$ is constant, implying that
$\Phi=-\Lambda$. Writing $e^{-2\Lambda}=1+2\varphi(r)$, we can compute the
gravitational force in 4DEGB due to a spherical body

$\vec{F}=-\frac{d\varphi}{dr}\hat{r}=-\frac{r}{2\alpha}\left(1-\frac{r^{3}+2\alpha
M}{r^{3}+8\alpha M}\sqrt{1+\frac{8\alpha M}{r^{3}}}\right)\hat{r}$ (28)

which is smaller in magnitude than its Newtonian $\alpha=0$ counterpart
($\vec{F}_{N}=-\frac{M}{r^{2}}\hat{r}$) for $\alpha>0$. The force in (28)
vanishes at $r=(\alpha M)^{1/3}$, but this is always at a smaller value of $r$
than the outer horizon $R_{h}=M-\sqrt{M^{2}-\alpha}$ of the corresponding
black hole. Hence the gravitational force outside of any spherical body, while
weaker than that in GR, is always attractive provided $\alpha>0$. If
$\alpha<0$ then the corresponding gravitational force is more attractive than
in GR.

Rescaling the various quantities using (13), (14), and (16) we obtain the
unitless equations

$\displaystyle\frac{d\bar{p}}{d\bar{r}}=\frac{(\bar{p}+\bar{\rho})\left[\bar{r}^{3}(\Gamma+8\pi\bar{\alpha}\bar{p}-1)-2\bar{\alpha}\bar{m}\right]}{\bar{r}^{2}\Gamma\left[\bar{r}^{2}(\Gamma-1)-2\bar{\alpha}\right]}$
(29) $\displaystyle\frac{d\bar{m}}{d\bar{r}}=4\pi\bar{r}^{2}\bar{\rho}$ (30)

which may be solved numerically. In the limit $\alpha\to 0$, the above
equations reduce back to the well-known TOV equations for a static,
spherically symmetric gravitating body in GR.

To solve (29) and (30) numerically we impose the boundary conditions

$m(0)=0,\quad\rho(0)=\rho_{\mathrm{c}},$ (31)

where the star’s surface radius $R$ is defined via $\bar{p}(\bar{R})=0$,
namely the radius at which the pressure goes to 0 (i.e. $p(R)=0$). We
similarly define the total mass of the star to be $M=m(R)$. Numerical
solutions can thus be obtained by scanning through a range of values of
$\rho_{c}$ and solving for the star’s total mass and radius.

Before proceeding to solve the TOV equations, we consider the behaviour of the
scalar field $\phi$ in the interior. Inserting the interior solution (25) into
(19), we find

$\displaystyle\lim_{r\to 0}\phi^{\prime}\approx\sqrt{\frac{m(0)}{2\alpha
r}}+\frac{3m(0)}{4\alpha}+\frac{m(0)\sqrt{r}\left(2\alpha
m^{\prime}(0)+5m(0)^{2}\right)}{4\sqrt{2}(\alpha m(0))^{3/2}}$ (32)
$\displaystyle+\frac{r\left(8\alpha\left(3m^{\prime}(0)-1\right)+35m(0)^{2}\right)}{32\alpha^{2}}+\mathcal{O}\left(r^{3/2}\right).$

Furthermore, provided $m(r)$ vanishes at least quadratically in $r$ for small
$r$ (which is ensured from (27) for the boundary conditions (31)) we find that
near the origin

$\displaystyle\lim_{r\to 0}\phi^{\prime}\sim-\frac{r}{4\alpha}\approx
0,\;\lim_{r\to 0}\phi\approx K-\frac{r^{2}}{8\alpha}\approx K,$ (33)

(where $K$ is a constant) and thus regularity of the scalar at the origin is
ensured.

Finally we note that the effective bag constant $B_{\mathrm{eff}}=60$
MeV/$\mathrm{fm}^{3}$ can be converted to units of inverse length squared
(‘gravitational units’) with the factor $G/c^{4}$, yielding

$B_{\mathrm{eff}}=7.84\times 10^{-5}\;\mathrm{km}^{-2}.$ (34)

###  II.4 Observational Constraints on the 4DEGB Coupling Constant

A recent study of the observational constraints on the coupling $\alpha$
yielded Clifton:2020xhc ; Fernandes:2022zrq

$-10^{-30}\textrm{m}^{2}<\alpha<10^{10}\;\textrm{m}^{2}$ (35)

where the lower bound comes from “early universe cosmology and atomic nuclei”
charmousis2022 , and the upper bound follows from LAGEOS satellite
observations. Regarding the lower bound as negligibly close to zero, the
dimensionless version of (35) reads

$0<\bar{\alpha}\lesssim 3.2.$ (36)

We note that inclusion of preliminary calculations on recent GW data suggest
these constraints could potentially tighten to $0<\alpha\lesssim
10^{7}\;\textrm{m}^{2}$, or alternatively $0<\bar{\alpha}\lesssim 0.0032$.
This would mean that deviations from GR due to 4DEGB would only be detectable
in extreme environments such as in the very early universe or near the surface
of extremely massive objects like black holes. Even tighter bounds were
assumed in previous studies of quark stars, where only solutions with $\alpha$
below 6 $\mathrm{km}^{2}$ ($\bar{\alpha}\leq 0.0019$) were considered
banerjee_2021_strange ; banerjee_2021_quark ; pretel_2022 . Adopting such a
tight bound would make compact stars near the upper end of the mass gap an
ideal candidate for investigation the effects of 4DEGB theory.

At this point in time such tighter bounds are not warranted. A proper study of
the effects of gravitational radiation in 4DEGB has yet to be carried out. In
view of this we shall assume the bound (36), which has strong observational
support Clifton:2020xhc ; Fernandes:2022zrq .

##  III Results

In this section we numerically obtain the total quark star mass as a function
of both total star radius and central density, with the former plots being
supplemented by the GR/4DEGB black hole horizon radii hennigar_2020_on ;
glavan_2020_einsteingaussbonnet  and the GR/4DEGB Buchdahl limits
chakraborty2020_limits .

Recall that the vacuum solution to the field equations is given by (25) with
$m(r)=M$, a constant, and $\Phi=-\Lambda$. Assuming no star, this solution
describes a black hole with two horizons provided
$M>\sqrt{\alpha}=M^{\mathrm{BH}}_{\mathrm{min}}$. The outer horizon radius is
hennigar_2020_on ; glavan_2020_einsteingaussbonnet ; gammon_2022

$R_{h}=M+\sqrt{M^{2}-\alpha}\;.$ (37)

The Buchdahl bound has been derived in 4DEGB chakraborty2020_limits  using
arguments similar to those in GR, namely that the pressure of the star remains
positive throughout and for a given mass, the radius of the star must be
larger than that of the (outer) horizon of a black hole of the same mass. This
yields

$\sqrt{1-\mu
R^{2}}\left(1+\alpha\mu\right)>\frac{1}{3}\left(1-\alpha\mu\right)$ (38)

where

$\mu\equiv\frac{1}{2\alpha}(\sqrt{1+\frac{8M\alpha}{R^{3}}}-1)$

for a star of radius $R$ and mass $M$. For small $\alpha$, (38) becomes

$\frac{M}{R}\leq\frac{4}{9}+\frac{16}{27}\left(\frac{\alpha}{R^{2}}\right)$
(39)

and we see that the GR Buchdahl bound is recovered in the $\alpha\to 0$ limit.
This latter relation holds provided spherical symmetry and isotropy are valid,
regardless of whether the internal density distribution is constant or a
function of the radial coordinate chakraborty2020_limits .

Unlike GR, the 4DEGB theory lacks a mass gap between a compact star and black
hole of the same radius. This feature was first observed in for neutron stars,
where the $M$ vs. $R$ curves were found to asymptote to the outer black hole
horizon for sufficiently large $\alpha$ charmousis2022 . To find the point at
which the Buchdahl bound intersects the black hole horizon, we substitute (37)
into (38) and solve for the mass $M$. Doing so, we find

$M_{\mathrm{int}}=\sqrt{\alpha}=M^{\mathrm{BH}}_{\mathrm{min}},$ (40)

or that the Buchdahl bound asymptotically approaches the smallest black hole
mass allowed by the theory.

These results imply that it is possible to have stable compact objects in
4DEGB whose radii are smaller than that of the GR Buchdahl bound $R\geq 9M/4$
or even that of the Schwarzschild radius $R=2M$. We shall demonstrate below
that this situation is indeed realized, and that quark stars can have radii
that are arbitrarily close to that of a black hole of the same mass.

###  III.1 Results for $\lambda\geq 0$

Here we present the solutions to the 4DEGB TOV equations for a positive
coupling constant $\lambda\geq 0$. This is equivalent to the condition

$\xi_{2a}\Delta^{2}>\xi_{2b}m_{s}^{2},$ (41)

and clearly from (14) as $\lambda$ increases in magnitude, so does
$\bar{\lambda}$. Our results are illustrated for a range of values of
$\bar{\alpha}$ and $\bar{\lambda}$ in Figures 1 and 2.

We find that in general, both increasing $\bar{\alpha}$ and increasing
$\bar{\lambda}$ expand the curves to larger mass and radii for the same choice
of central pressure, yielding larger maximum mass points. This is due to
gravitational attraction in 4DEGB being increasingly weaker as $\alpha>0$
increases, as well as a larger $\lambda$ mapping to a stiffer equation of
state as a consequence of the strong interaction effects. We see that for any
given $\alpha$, the hook-shaped $M$ vs. $R$ curves move upward and to the
right as $\lambda$ increases. Likewise, for any given $\lambda$, the curves
again move upward and rightward, with the upper part more rapidly asymptoting
to the Buchdahl/horizon limit as $\alpha$ increases.

It is quite interesting to note that even for small $\alpha$, with a high
enough central pressure there are stable solutions that are smaller in radius
than both the GR Buchdahl bound and the Schwarzschild black hole. We
illustrate this in Fig. 3, which shows that, even for $\bar{\alpha}=0.001$
there is a narrow range for small mass, small radius stars that are within not
only the GR Buchdahl bound, but also the $2M$ Schwarzschild radius. Similar
results can be seen in the 2nd panel of Figure 6c for $\bar{\alpha}=0.0001$
when $\lambda<0$. The stability of these objects remains an interesting
subject for investigation.

Our results are commensurate with those of Banerjee et al. banerjee_2021_quark
in the limited overlapping regions of parameter space, in that increasing the
coupling strength to the higher curvature terms tends to increase the mass-
radius profiles, in turn yielding a larger local maximal mass point. In our
case not all curves have this local extrema, but in the small $\alpha$ and
$\rho_{c}$ regime we see the same behaviour. In banerjee_2021_quark , rather
than varying coupling strength directly, different values of the MIT bag
constant are considered, with a smaller $B_{\mathrm{eff}}$ corresponding to
larger M-R profiles (which in our case is correlated with a larger positive
interaction strength).

(a)  (b)  (c)  (d)

(e)  (f)

Figure 1: Mass vs. radius/central density curves for unified interacting quark
stars when $\lambda>0$. Each plot shows results for a unique value of the
4DEGB coupling. In order from blue to purple curves, the $\bar{\lambda}$
values considered are 0, 0.5, 10, 100, $\infty$ respectively. The solid and
dashed black lines correspond to the 4DEGB equivalent of the Schwarzschild and
Buchdahl limits, respectively (with their red counterparts marking the
equivalent bounds in GR). In general a larger $\alpha$ and/or $\lambda$ tends
to increase the mass and radius of a solution.

(a)  (b)  (c)  (d)

(e)  (f)

Figure 2: Additional mass vs. radius/central density curves for unified
interacting quark stars when $\lambda>0$. In order from blue to purple curves,
the $\bar{\lambda}$ values considered are 0, 0.5, 10, 100, $\infty$
respectively. The solid and dashed black lines correspond to the 4DEGB
equivalent of the Schwarzschild and Buchdahl limits, respectively (with their
red counterparts marking the equivalent bounds in GR). In general a larger
$\alpha$ and/or $\lambda$ tends to increase the mass and radius of a solution.
At $\alpha/\lambda$ combinations large enough for $p_{\mathrm{crit}}>0$ we
observe solutions which start at the black hole horizon (as shown in the lower
two panels represented by coloured dots). These solutions are unstable for all
choices of central pressure.

(a)  (b)

(c)

Figure 3: Mass vs. radius and central density for $\bar{\alpha}=.001$ when
$\lambda>0$. We see that even for very small alpha there exists a range of
solutions yielding quark stars smaller than the GR Buchdahl bound (dashed red)
and the $2M$ Schwarzschild radius (solid red); the corresponding bounds for
4DEGB are in black.

Furthermore, for sufficiently large $\alpha/\lambda$, no quark star solutions
exist outside of the black hole horizon, since the pressure profiles smaller
than criticality diverge to $\infty$ with no real roots. This in turn implies
that there should be a well-defined central pressure/density at which
$p(r)=p(0)=\mathrm{const}$, and below which quark star solutions do not exist.
This value can be derived analytically (see Appendix A for a detailed
derivation), and we find that below the critical central pressure:

$\displaystyle p_{\mathrm{crit}}=\frac{1}{2\pi(9\pi-32\alpha\lambda)^{2}}$
$\displaystyle\Bigg{[}256\alpha(2\pi\alpha+3)\lambda^{2}+8\pi(2\pi\alpha(8\pi\alpha+21)-9)\lambda+9\pi^{3}(4\pi\alpha-3)$
(42)
$\displaystyle-\mathrm{sgn}(\lambda)6\sqrt{\lambda(32\alpha\lambda+\pi(8\pi\alpha-3))^{2}\left(16(2\pi\alpha+1)\lambda+\pi^{2}(8\pi\alpha+3)\right)}\Bigg{]},$

there are no real roots of $p(r)$ with which to define the star’s surface.

This central pressure is never realized physically. We find that all $M$ vs.
$R$ curves with $p_{\mathrm{crit}}>0$ begin at the 4DEGB black hole horizon
and extend to values within the horizon, or in other words, in a disallowed
region of parameter space. We illustrate this in Figure 2c, where we notice
that the $M$ vs. $R$ curves with $p_{\mathrm{crit}}>0$ begin at the 4DEGB
black hole horizon (as coloured points). The first threshold of interest is to
consider is the smallest value of $\bar{\alpha}$ at $\bar{\lambda}=0$ for
which none of the solutions lie below the Buchdahl threshold – for such
parameter combinations no physical quark star solutions exist.

To investigate this further, consider the pressure profile of a star with a
central pressure just past criticality (ie. $p(0)=p_{\mathrm{crit}}+\delta$
where $\delta<<1$). As shown in Appendix A, the pressure profile is well-
approximated by a step function times a constant, allowing us to approximately
solve (30) analytically:

$\bar{m}(\bar{R})^{\mathrm{min}}=\bar{M}^{\mathrm{min}}\approx
4\pi\int_{0}^{\bar{R}}\Theta(\bar{r}-\bar{R})\bar{\rho}(0)\bar{r}^{2}d\bar{r}=\frac{4\pi}{3}\bar{R}^{3}\bar{\rho}(0).$
(43)

Then for the special case of vanishing interaction strength,

$\bar{M}_{\lambda\to
0}^{\mathrm{min}}\approx\frac{4\pi}{3}\bar{R}^{3}(3\bar{p}(0)+1).$ (44)

Since $p(0)\approx p_{\mathrm{crit}}$ and $p_{\mathrm{crit}}^{\lambda\to
0}=\frac{2\pi\bar{\alpha}}{9}-\frac{1}{6}$, this can be further simplified
into the following form:

$\bar{M}_{\lambda\to
0}^{\mathrm{min}}\approx\frac{2}{3}\pi(\frac{4\pi}{3}\bar{\alpha}+1)\bar{R}^{3}.$
(45)

Inserting this minimal mass solution into (29) and finding where
$\bar{p}^{\prime}(\bar{r})$ diverges, we obtain the corresponding radius

$\bar{R}_{\lambda\to 0}^{\mathrm{min}}\approx\sqrt{\frac{3}{4\pi}}$ (46)

of this minimal mass point, which is (interestingly) independent of
$\bar{\alpha}$. Consequently

${\bar{M}_{\lambda\to
0}^{\mathrm{min}}}\approx\frac{3+4\pi\bar{\alpha}}{4\sqrt{3\pi}}.$ (47)

Inserting this into the full 4DEGB Buchdahl bound and solving for
$\bar{\alpha}$, we find that (38) is saturated when
$\bar{\alpha}=\frac{3}{4\pi}$, which consequently is also the solution to
$\lim_{\lambda\to 0}p_{\mathrm{crit}}(\bar{\alpha})=0$ (the smallest
$\bar{\alpha}$ with a non-negative critical pressure in the limit of vanishing
QM interaction), implying that $\bar{M}_{\lambda\to
0}^{\mathrm{min}}|_{\bar{\alpha}=\frac{3}{4\pi}}=\bar{M}_{\mathrm{int}}$.

For a more general choice of $\bar{\alpha}$ it is easy to see that by
inserting (46) into the 4DEGB black hole mass
$M_{BH}=\frac{\alpha+R_{H}^{2}}{2R_{H}}$ that $M_{BH}(\bar{R}_{\lambda\to
0}^{\mathrm{min}})=\bar{M}_{\lambda\to 0}^{\mathrm{min}}$. In other words, all
$\bar{\lambda}=0$ solution curves with non-zero critical pressure must start
at the black hole horizon (or the Buchdahl bound/horizon intersection, in the
case where $\bar{\alpha}=\frac{3}{4\pi}$), and thus none of these solutions
are stable. This can also be checked numerically by considering
$\bar{\alpha}=\frac{3}{4\pi}+\delta$ (where $\delta<<1$) and plotting the
lowest pressure solution possible, which indeed intersects the minimal mass
point of the Buchdahl curve.

A similar analysis can be repeated in the limit $\lambda\to\infty$. Doing so
we find an equivalent expression for the mass of the threshold point:

$\bar{M}_{\lambda\to\infty}^{\mathrm{min}}\approx\frac{\left(1+2\pi\bar{\alpha}-\sqrt{2\pi\bar{\alpha}+1}\right)}{2\bar{\alpha}}\bar{R}^{3}$
(48)

and radius:

$\bar{R}_{\lambda\to\infty}^{\mathrm{min}}\approx\sqrt{\frac{2\bar{\alpha}}{\sqrt{8\pi\bar{\alpha}-4\sqrt{2\pi\bar{\alpha}+1}+5}-1}},$
(49)

the latter of which now has an explicit $\bar{\alpha}$ dependence. Combining
the two we find

$\bar{M}_{\lambda\to\infty}^{\mathrm{min}}\approx\sqrt{2\bar{\alpha}}\left(2\pi\bar{\alpha}+1-\sqrt{2\pi\bar{\alpha}+1}\right)\left(\frac{1}{\sqrt{8\pi\bar{\alpha}-4\sqrt{2\pi\bar{\alpha}+1}+5}-1}\right)^{3/2}.$
(50)

As before, we insert this mass/radius pair into the full 4DEGB Buchdahl bound
(38), and find that the expression is saturated when
$\bar{\alpha}=\frac{3}{2\pi}$, once again the solution to
$\lim_{\lambda\to\infty}p_{\mathrm{crit}}(\bar{\alpha})=0$. Inserting
$\bar{\alpha}=\frac{3}{2\pi}$ in to (50) we confirm once again that
$\bar{M}_{\lambda\to\infty}^{\mathrm{min}}|_{\bar{\alpha}=\frac{3}{2\pi}}=M_{\mathrm{int}}$,
implying that solutions with critical pressure will start at the intersection
point of the 4DEGB Buchdahl bound and black hole horizon, and are thus always
unphysical. For an arbitrary $\bar{\alpha}>\frac{3}{2\pi}$ it is easy to show
that again
$M_{BH}(\bar{R}_{\lambda\to\infty}^{\mathrm{min}})=\bar{M}_{\lambda\to\infty}^{\mathrm{min}}$,
indicating that the lowest pressure solution intersects the black hole
horizon. As with the non-interacting case, the $\bar{\lambda}\to\infty$
solutions are also unstable for any central pressure above the critical value,
and thus for any $\bar{\alpha}/\bar{\lambda}$ combination with a positive
critical pressure. These results are manifest in Figures 1 and 2 as points
lying along the solid black line. The above analysis also breaks the parameter
space up into three different regions, which can be seen in Figure 4.

Figure 4: The solid black line (bordering the gray region) shows
$p_{\mathrm{crit}}(\alpha)$ for a test value of $\bar{\lambda}=0.1$ (and
$\lambda>0$). For a fixed central pressure of $\bar{p}_{0}=$ 0.15 the
parameter space has three distinct regions. From $\bar{\alpha}=0$ to the value
that satisfies $p_{\mathrm{crit}}(\alpha)=0$ (green) there exist valid,
physical solutions at fixed $p_{0}$. As $\bar{\alpha}$ increases past this the
red region becomes manifest, having parameter combinations with
$p_{0}>p_{\mathrm{crit}}$, and thus with unphysical radii inside the black
hole horizon. the gray region has parameter combinations with
$p_{0}<p_{\mathrm{crit}}$ (and thus radial pressure profiles that diverge to
$\infty$). The boundary between the gray and red regions corresponds with
solutions lying directly at the black hole horizon - these are included in our
M-R plots as coloured points.

(a)  (b)

Figure 5: Mass vs. radius and central density where the interaction strength
is fixed at 0, and $\bar{\alpha}=$ 0.006, 0.0065, 0.007, 0.0072, 0.0074,
0.0075, 0.0076, 0.0078, 0.0079, 0.00795, 0.008, 0.0085 (lower $\bar{\alpha}$
corresponding to lower maximal mass points). The purpose of this plot is to
show the progression of the maximum mass point as $\alpha$ changes slowly.

###  III.2 Results for $\lambda\leq 0$

(a)  (b)  (c)  (d)

(e)  (f)

Figure 6: Mass-radius and mass-central density relations in the case where
$\lambda<0$. Each $\bar{\alpha}$ corresponds to a different set of
$\bar{\lambda}$ for which the solutions are plotted, due to the extreme
sensitivity of the critical pressure with respect to negative $\lambda$. In
all cases the ratio $\frac{\bar{\lambda}}{\bar{\lambda}_{\mathrm{crit}}}$
takes on the same set of values, namely 0, 1/25, 1/10, 1/2, 1 respectively
from the blue to purple curves. The effect of increasing $\bar{\lambda}$
(decreasing $\lambda$) is to suppress the mass and radius, a logical extension
of the $\lambda>0$ results.

(a)  (b)  (c)  (d)

(e)  (f)

Figure 7: Mass-radius and mass-central density relations in the case where
$\lambda<0$. Each $\bar{\alpha}$ corresponds to a different set of
$\bar{\lambda}$ for which the solutions are plotted, due to the extreme
sensitivity of the critical pressure with respect to negative $\lambda$. In
all cases the ratio $\frac{\bar{\lambda}}{\bar{\lambda}_{\mathrm{crit}}}$
takes on the same set of values, namely 0, 1/25, 1/10, 1/2, 1 respectively
from the blue to purple curves. The effect of increasing $\bar{\lambda}$
(decreasing $\lambda$) is to suppress the mass and radius, a logical extension
of the $\lambda>0$ results.

Here we summarize the results for $\lambda<0$, which is equivalent to the
condition

$\xi_{2a}\Delta^{2}<\xi_{2b}m_{s}^{2}.$ (51)

In this case (and so throughout this section) a larger $\bar{\lambda}$
corresponds to a negative $\lambda$ that is larger in magnitude (more
negative). As is shown in detail in appendix A, the critical pressure now has
a divergence at $\bar{\alpha}=\frac{9\pi}{32\bar{\lambda}}$ in this regime.
This defines a critical relationship between the 4D Gauss-Bonnet coupling and
the strong interaction parameter at which the central pressure must approach
infinity for a valid quark star solution. This also means that the allowed
range of interaction strengths inside the star vary based on the choice of
coupling to the higher curvature term.

We illustrate our results in figures 6 and 7. In order to keep the results
manageable and easy to interpret, we choose to never use
$\bar{\alpha}-\bar{\lambda}$ combinations that result in a critical centre
pressure larger than $\bar{p}_{\mathrm{crit}}=10$ (and define
$\bar{\lambda}_{\mathrm{crit}}$ to be that which corresponds with this
solution set, since $\bar{\lambda}_{\mathrm{max}}$ has a divergent critical
pressure - see Figure 13). By doing this, we can more easily compare a range
of $\bar{\alpha}$ solutions with the same choices of
$\frac{\bar{\lambda}}{\bar{\lambda}_{\mathrm{crit}}}$, and view them on axes
of comparable magnitudes.

As expected we see that decreasing $\lambda$ (and hence increasing
$\bar{\lambda}$) causes the mass and radius of a given solution to trend
downward, whereas increasing $\bar{\alpha}$ yields the same upward trend as
before. In other words, in the negative lambda regime strong interaction
effects tend to somewhat counteract the effects of non-zero coupling to the
Gauss-Bonnet theory. We also notice that the effects of the critical pressure
are first manifest on the curves with a larger negative magnitude, and thus
the end of the solution space occurs when $\alpha$ is large enough to trigger
criticality for the case of vanishing interaction effects ($\bar{\lambda}=0$).
As in the previous section this will occur when $\bar{\alpha}=\frac{3}{4\pi}$,
and thus a stronger coupling to 4DEGB theory will result in empty plots as
seen in Figures 7e and 7f.

Finally, we note from Figure 6 that there are again regions where solutions
exist that violate the GR Buchdahl bounds and Schwarzschild radius, even for
very small values of $\bar{\alpha}$. The range of solutions gets large as
$\bar{\alpha}$ increases, as shown in Figure 7.

##  IV Stability Analysis

We now consider the stability of the 4DEGB unified interacting quark stars. A
necessary but insufficient condition for an uncharged compact star’s stability
in Einstein gravity is $dM/d\rho_{c}<0$ glendenning1998_nonidentical ;
arbanil2016_radial ; zhang_2021_stellar , corresponding to the parts of the
solution curves before maximum mass points are reached. In Einstein-Maxwell
theory a net charge can offset stability from the maximum mass point in either
direction zhang_2021_unified . In a similar vein, when the 4DEGB theory
coupling is non-zero, it is not obvious whether the coincidence of stability
and maximum mass point will hold. Leaving a thorough analysis of the
fundamental radial oscillation modes for future work, we note that stable
compact objects must obey the causality condition that speed of sound $c_{s}$
never exceeds the speed of light $c$. Inserting the equation of state (12)
directly into the definition for $c_{s}$, we find

$c_{s}=\left(3\;-\;\frac{\textrm{sgn}({\lambda})4}{\sqrt{\frac{4\lambda+4\pi^{2}\bar{p}(\bar{r})+\pi^{2}}{\lambda}}}\right)^{-1/2}\;.$
(52)

To attain the limit $c_{s}\to 1$ from below implies that $\bar{p}(r)=-1/4$
(and $\lambda>0$), a value that is never part of the quark star solution
space. In Figure 8 we plot the sound speed as a function of the interaction
strength $\lambda$ for various fixed values of the pressure. We see that for
positive pressures the $c_{s}\to 1$ limit is asymptotically approached from
below as $\lambda\to\infty$; the sound speed is never superluminal. If
$\lambda<0$ instead, the analogous solutions curve downward and asymptotically
approach $c_{s}=1/\sqrt{5}$ as expected from zhang_2021_unified ;
zhang_2021_stellar .

(a) Zoomed in. (b) Zoomed out.

Figure 8: Sound speed inside the unified interacting quark star as a function
of interaction strength $\bar{\lambda}$ when pressure is fixed (note that this
quantity is independent of the coupling strength $\alpha$). The solid lines
correspond to $\lambda>0$ and the dashed lines $\lambda<0$.

Another indicator of stability is the effective adiabatic index $\gamma$ of
perturbations. For adiabatic oscillations, this can be defined via the speed
of sound through the quark-gluon plasma:

$\gamma_{\mathrm{eff}}\equiv\left(1+\frac{\rho}{P}\right)\left(\frac{dp}{d\rho}\right)_{S}$
(53)

where the subscript $S$ indicates that we consider the sound speed at constant
specific entropy (recalling that $c_{s}=\sqrt{(\frac{\partial
p}{\partial\rho})_{S}}$). One can consider this quantity to be a bridge
“between the relativistic structure of a spherical static object and the
equation of state of the interior fluid” moustakidis2017_stability . In
principle a critical value for $\langle\gamma_{\mathrm{eff}}\rangle$ exists,
below which configurations are unstable against radial perturbations. In
standard general relativity this critical value can be written
chandrasekhar1964 ; moustakidis2017_stability
$\gamma_{cr}=\frac{4}{3}+\frac{19}{42}\beta$, where $\beta=2M/R=R_{S}/R$ is
the compactness parameter. If $\beta\to 0$ the well-known classical Newtonian
limit is recovered as expected
($\langle\gamma_{\mathrm{eff}}\rangle\geq\frac{4}{3}$).

An equivalent bound has not yet been derived for the 4DEGB theory. Despite
this, it is common practice to plot the adiabatic index of the star relative
to the Newtonian critical value banerjee_2021_quark ; banerjee_2021_strange ;
hansraj_2020_isotropic ; singh_2022_anisotropic . In general $\gamma$ will
depend on $\lambda,\alpha,\rho_{0}$ and $r$. Given the size of our parameter
space, plotting this would be cumbersome and not particularly illuminating
(although an example case is shown in Figure 9). We already know due to the
form of the equations that $\lim_{r\to R}\gamma\to\infty$, and due to the
monotonically decreasing nature of our pressure profiles,
$\frac{\partial\gamma}{\partial r}>0$. With this, we are mostly interested in
the lower limit of $\gamma_{r\to 0}$ (ie. where the curves in Figure 9 start
from). This information is laid out in table 1 where, noting that
$\lim_{p\to\infty}\gamma_{\mathrm{eff}}$ is the lower limit on $\gamma_{r=0}$,
we can see that in all cases $\gamma_{r=0}\geq\frac{4}{3}$, and also that
$\gamma$ always diverges at the surface of the star. Since there is no well
defined upper limit on interaction when $\lambda<0$, the first column of the
third row cannot be further simplified. The results of this table imply that
$\langle\gamma\rangle>\gamma_{\mathrm{crit}}^{\mathrm{GR}}$ in all cases.

(a) $\bar{\lambda}=0$ (b) $\bar{\lambda}=\infty$

Figure 9: Adiabatic index $\gamma_{\mathrm{eff}}$ plotted over the surface of the star when $\bar{\alpha}=0.0001$ and $\lambda>0$. Different coloured curves represent different central pressures (with the smallest central pressures having the largest $\gamma_{\mathrm{eff}}$ at the origin), and the black curve is the Newtonian limit 4/3.  Table 1: Limits on the adiabatic index $\gamma$. | $\gamma_{\mathrm{eff}}$ | $\lim_{p\to\infty}\gamma_{\mathrm{eff}}$ |  $\lim_{p\to 0}\gamma_{\mathrm{eff}}$ ($\gamma_{r\to R}$)
---|---|---|---
$\lambda=0$ | $\frac{1}{3}(4+\frac{1}{p(r)})$ | $\frac{4}{3}$ | $\infty$
$\lambda\to\infty$ | $2(1+\frac{1}{4p(r)})$ | $2$ | $\infty$
$\lambda<0$ | $\frac{\sqrt{\frac{4\lambda+4\pi^{2}p(r)+\pi^{2}}{\lambda}}\left(4\lambda+2\lambda\sqrt{\frac{4\lambda+4\pi^{2}p(r)+\pi^{2}}{\lambda}}+4\pi^{2}p(r)+\pi^{2}\right)}{\pi^{2}p(r)\left(3\sqrt{\frac{4\lambda+4\pi^{2}p(r)+\pi^{2}}{\lambda}}+4\right)}$ | $\frac{4}{3}$ | $\infty$

##  V Summary

In this paper we have investigated the stellar structure of strongly
interacting quark stars in the regularized 4D Einstein Gauss-Bonnet theory of
gravity for different combinations of the 4DEGB coupling constant $\alpha$ and
the unified strong interaction parameter $\lambda$ in our interacting quark
matter equation of state. In accord with the lack of a mass gap in the 4DEGB
theory charmousis2022 , we find that – for both signs of the coupling
$\lambda$ – even for small $\alpha$ the quark star solutions asymptotically
approach the 4DEGB black hole horizon radius, and thus have solutions with
smaller radii than the GR Buchdahl/Schwarzschild limits. It is worth noting
that in the $\lambda<0$ case a much larger central pressure is required to
approach this limit as compared to the analogous $\lambda>0$ case. In general,
larger $\alpha$ and $\lambda$ tend to increase the mass-radius profile of
quark stars, and a large negative $\lambda$ has the opposite effect, pushing
back against the contributions of non-zero coupling to the higher curvature
gravity terms. These findings are generally consistent with what was found in
the regime of weak coupling to the 4DEGB theory banerjee_2021_quark .

We have found many additional features in the unexplored regions of parameter
space, the most striking of which is that 4DEGB quark stars can exist not only
below the GR Buchdahl bound, but also smaller than the $2M$ Schwarzschild
radius. These objects obey the criterion that their sound speed is subluminal
and have an average effective adiabatic index above what is required in GR. An
interesting avenue for future work would be to study the stability properties
of these objects, and of compact 4DEGB stars in general.

We have also analytically identified a critical centre pressure (which is a
function of $\alpha$/$\lambda$ only), below which quark star solutions do not
exist. This critical pressure cannot be physically realized; rather it
corresponds to quark stars whose radii are either equal to or smaller than the
corresponding black hole of the same mass. In principle, the $M-R$ pairs
associated with this critical pressure can be used to constrain the constants
of the theory (for instance, if a particular $\alpha$/$\lambda$ pair does not
allow for $M-R$ pairs found observationally).

In addition to a detailed study of the stability of compact stars in 4DEGB
gravity, our work should be extended to consider the effects of charge on
stellar structure, since the confining nature of the strong interaction makes
quark stars one of the few remaining candidates for stellar objects which are
not electrically neutral zhang_2021_stellar ; pretel_2022 ; gammon_2022_beyond
. Moreover, real astrophysical objects generally also have a net angular
momentum, which can alter the stellar structure. Our calculations could be
extended to a slowly rotating metric ansatz adair2020 ; gammon_2022  to see
how the introduction of angular momentum changes the mass-radius relations of
our interacting 4DEGB quark star.

## Acknowledgements

This work was supported in part by the Natural Sciences and Engineering
Research Council of Canada. We are grateful to Sharon Morsink for helpful
correspondence.

## References

  * [1] Robie A. Hennigar, David Kubizňák, Robert B. Mann, and Christopher Pollack.  On taking the d → 4 limit of gauss-bonnet gravity: theory and solutions.  Journal of High Energy Physics, 2020, 07 2020.
  * [2] Chen Zhang and Robert B. Mann.  Unified interacting quark matter and its astrophysical implications.  Physical Review D, 103, 03 2021.
  * [3] Jamil Ahmed, Robie A. Hennigar, Robert B. Mann, and Mozhgan Mir.  Quintessential Quartic Quasi-topological Quartet.  JHEP, 05:134, 2017.
  * [4] P. Bueno & P. Cano.  Einsteinian cubic gravity.  Phys. Rev. D, 94, 2016.
  * [5] D. Lovelock.  The einstein tensor and its generalizations.  Journal of Mathematical Physics, 12, 1971.
  * [6] Drazen Glavan and Chunshan Lin.  Einstein-gauss-bonnet gravity in four-dimensional spacetime.  Physical Review Letters, 124(8), Feb 2020.
  * [7] Shou-Long Li, Puxun Wu, and Hongwei Yu.  Stability of the Einstein Static Universe in $4D$ Gauss-Bonnet Gravity, 4 2020.
  * [8] Tsutomu Kobayashi.  Effective scalar-tensor description of regularized lovelock gravity in four dimensions.  Journal of Cosmology and Astroparticle Physics, 2020(07):013–013, jul 2020.
  * [9] Arun Kumar, Rahul Kumar Walia, and Sushant G. Ghosh.  Bardeen black holes in the regularized 4d einstein gauss bonnet gravity.  Universe, 8(4):232, apr 2022.
  * [10] Pedro G.S. Fernandes.  Charged black holes in AdS spaces in 4d einstein gauss-bonnet gravity.  Physics Letters B, 805:135468, jun 2020.
  * [11] Rahul Kumar and Sushant G. Ghosh.  Rotating black holes in 4d einstein-gauss-bonnet gravity and its shadow.  Journal of Cosmology and Astroparticle Physics, 2020(07):053–053, jul 2020.
  * [12] Arun Kumar, Dharmanand Baboolal, and Sushant G. Ghosh.  Nonsingular black holes in 4d einstein–gauss–bonnet gravity.  Universe, 8(4):244, apr 2022.
  * [13] Daniele Malafarina, Bobir Toshmatov, and Naresh Dadhich.  Dust collapse in 4d einstein–gauss–bonnet gravity.  Physics of the Dark Universe, 30:100598, dec 2020.
  * [14] Daniela D. Doneva and Stoytcho S. Yazadjiev.  Relativistic stars in 4d einstein-gauss-bonnet gravity.  Journal of Cosmology and Astroparticle Physics, 2021(05):024, may 2021.
  * [15] C. Charmousis, A. Lehébel, E. Smyrniotis, and N. Stergioulas.  Astrophysical constraints on compact objects in 4d einstein-gauss-bonnet gravity.  Journal of Cosmology and Astroparticle Physics, 2022(02):033, Feb 2022.
  * [16] Sushant G. Ghosh and Sunil D. Maharaj.  Radiating black holes in the novel 4D Einstein–Gauss–Bonnet gravity.  Phys. Dark Univ., 30:100687, 2020.
  * [17] Metin Gürses, Tahsin Çağrı Şişman, and Bayram Tekin.  Is there a novel einstein–gauss–bonnet theory in four dimensions?  The European Physical Journal C, 80(7), jul 2020.
  * [18] Wen-Yuan Ai.  A note on the novel 4d einstein–gauss–bonnet gravity.  Communications in Theoretical Physics, 72(9):095402, aug 2020.
  * [19] Fu-Wen Shu.  Vacua in novel 4d einstein-gauss-bonnet gravity: Pathology and instability?  Physics Letters B, 811:135907, dec 2020.
  * [20] Pedro G. S. Fernandes, Pedro Carrilho, Timothy Clifton, and David J. Mulryne.  Derivation of Regularized Field Equations for the Einstein-Gauss-Bonnet Theory in Four Dimensions.  Phys. Rev. D, 102(2):024025, 2020.
  * [21] Robert B. Mann and S. F. Ross.  The D —$>$ 2 limit of general relativity.  Class. Quant. Grav., 10:1405–1408, 1993.
  * [22] H. Lu and Yi Pang.  Horndeski gravity as $D\rightarrow 4$ limit of Gauss-Bonnet.  Phys. Lett. B, 809:135717, 2020.
  * [23] Gregory Walter Horndeski.  Second-order scalar-tensor field equations in a four-dimensional space.  Int. J. Theor. Phys., 10:363–384, 1974.
  * [24] Timothy Clifton, Pedro Carrilho, Pedro G. S. Fernandes, and David J. Mulryne.  Observational Constraints on the Regularized 4D Einstein-Gauss-Bonnet Theory of Gravity.  Phys. Rev. D, 102(8):084005, 2020.
  * [25] Pedro G. S. Fernandes, Pedro Carrilho, Timothy Clifton, and David J. Mulryne.  The 4D Einstein–Gauss–Bonnet theory of gravity: a review.  Class. Quant. Grav., 39(6):063001, 2022.
  * [26] Chen Zhang, Michael Gammon, and Robert B. Mann.  Stellar structure and stability of charged interacting quark stars and their scaling behaviour.  Physical Review D, 104, 12 2021.
  * [27] Bob Holdom, Jing Ren, and Chen Zhang.  Quark matter may not be strange.  Physical Review Letters, 120, 05 2018.
  * [28] Zhiqiang Miao, Jin-Liang Jiang, Ang Li, and Lie-Wen Chen.  Bayesian inference of strange star equation of state using the GW170817 and GW190425 data.  The Astrophysical Journal Letters, 917(2):L22, aug 2021.
  * [29] Luiz L. Lopes and Debora P. Menezes.  On the nature of the mass-gap object in the GW190814 event.  The Astrophysical Journal, 936(1):41, aug 2022.
  * [30] P. T. Oikonomou and Ch. C. Moustakidis.  Colour-flavour locked quark stars in light of the compact object in hess j1731-347 and the gw190814 event, 2023.
  * [31] Chen Zhang, Yong Gao, Cheng-Jun Xia, and Renxin Xu.  Rescaling strangeon stars and its implications on gravitational-wave echoes, 2023.
  * [32] J. E. Horvath, L. S. Rocha, L. M. de Sá , P. H. R. S. Moraes, L. G. Barão, M. G. B. de Avellar, A. Bernardo, and R. R. A. Bachega.  A light strange star in the remnant HESS j1731-347: Minimal consistency checks.  Astronomy & Astrophysics, 672:L11, apr 2023.
  * [33] Ayan Banerjee, Takol Tangphati, and Phongpichit Channuie.  Strange quark stars in 4d einstein–gauss–bonnet gravity.  The Astrophysical Journal, 909(1):14, mar 2021.
  * [34] Ayan Banerjee, Takol Tangphati, Daris Samart, and Phongpichit Channuie.  Quark stars in 4d einstein–gauss–bonnet gravity with an interacting quark equation of state.  The Astrophysical Journal, 906:114, 01 2021.
  * [35] Juan M. Z. Pretel, Ayan Banerjee, and Anirudh Pradhan.  Electrically charged quark stars in 4d einstein–gauss–bonnet gravity.  The European Physical Journal C, 82(2), feb 2022.
  * [36] Sumanta Chakraborty and Naresh Dadhich.  Limits on stellar structures in lovelock theories of gravity.  Physics of the Dark Universe, 30:100658, dec 2020.
  * [37] Robert B. Mann, Sebastian Murk, and Daniel R. Terno.  Black holes and their horizons in semiclassical and modified theories of gravity.  Int. J. Mod. Phys. D, 31(09):2230015, 2022.
  * [38] Shao-Wen Wei and Yu-Xiao Liu.  Testing the nature of Gauss-Bonnet gravity by four-dimensional rotating black hole shadow.  Eur. Phys. J. Plus, 136(4):436, 2021.
  * [39] Norman K. Glendenning.  Compact Stars.  1996\.
  * [40] Yu-Peng Zhang, Shao-Wen Wei, and Yu-Xiao Liu.  Spinning test particle in four-dimensional einstein gauss bonnet black holes.  Universe, 6(8):103, jul 2020.
  * [41] Jing Ren and Chen Zhang.  Quantum nucleation of up-down quark matter and astrophysical implications.  Physical Review D, 102(8), oct 2020.
  * [42] En-Ping Zhou, Xia Zhou, and Ang Li.  Constraints on interquark interaction parameters with GW170817 in a binary strange star scenario.  Physical Review D, 97(8), apr 2018.
  * [43] Michael Gammon and Robert Mann.  Slowly rotating black holes in 4d gauss-bonnet gravity, 2022.
  * [44] Dražen Glavan and Chunshan Lin.  Einstein-gauss-bonnet gravity in four-dimensional spacetime.  Physical Review Letters, 124, 02 2020.
  * [45] Norman K. Glendenning and Christiane Kettner.  Nonidentical neutron star twins.  Astron. Astrophys., 353:L9, 2000.
  * [46] José D.V. Arbañil and M. Malheiro.  Radial stability of anisotropic strange quark stars.  Journal of Cosmology and Astroparticle Physics, 2016(11):012, nov 2016.
  * [47] Ch. C. Moustakidis.  The stability of relativistic stars and the role of the adiabatic index.  Gen. Rel. Grav., 49(5):68, 2017.
  * [48] S. Chandrasekhar.  The Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity.  Astrophys. J. , 140:417, August 1964.
  * [49] Sudan Hansraj, Ayan Banerjee, Lushen Moodly, and M K Jasim.  Isotropic compact stars in 4d einstein–gauss–bonnet gravity.  Classical and Quantum Gravity, 38(3):035002, dec 2020.
  * [50] Ksh. Newton Singh, S. K. Maurya, Piyali Bhar, and Riju Nag.  Anisotropic solution for polytropic stars in 4d einstein–gauss–bonnet gravity.  The European Physical Journal C, 82(9), sep 2022.
  * [51] Michael Gammon.  Beyond general relativity: New structure and stability predictions from modified gravity and high-density interaction models.  Master’s thesis, University of Waterloo, 2022.
  * [52] C. Adair et al.  Slowly rotating black holes in einsteinian cubic gravity.  Phys. Rev. D, 102, 2020.

##  Appendix A Critical Centre Pressure

The quark star solution space is filled out numerically by scanning through
test values for the star’s centre pressure $p_{\mathrm{c}}$ and finding the
radius at which pressure vanishes to define the stellar surface. Typical
compact star solutions have pressure profiles reminiscent of Figure 10a where
pressure decreases monotonically starting from $r=0$. For particular
$\alpha-\lambda$ combinations, no quark star solutions exist below some
critical center pressure $p_{\mathrm{crit}}$. For $p_{c}$ below this critical
value we find that the pressure blows up toward infinity and thus has no
roots. In theory there should exist some $p_{c}=p_{\mathrm{crit}}$ whose
pressure plot is constant, straight line.

To derive this critical pressure, we replace $\bar{p}(r)$ with the constant
$\bar{p}_{\mathrm{crit}}$ in the EOS (15). The differential equation (30) for
$m(r)$ becomes

$\bar{m}^{\prime}(\bar{r})=4\pi\bar{r}^{2}\left(-\frac{4\bar{\lambda}\left(\sqrt{\frac{\pi^{2}\left(\bar{p}_{\mathrm{crit}}+\frac{1}{4}\right)}{\lambda}+1}-1\right)}{\pi^{2}}+3\bar{p}_{\mathrm{crit}}+1\right)$
(54)

which is straightforwardly integrated, yielding

$\bar{m}(\bar{r})=\frac{4\bar{r}^{3}\left(\pi^{2}(3\bar{p}_{\mathrm{crit}}+1)-2\bar{\lambda}\left(\sqrt{\frac{4\bar{\lambda}+4\pi^{2}\bar{p}_{\mathrm{crit}}+\pi^{2}}{\bar{\lambda}}}-2\right)\right)}{3\pi}+c_{1}.$
(55)

where the boundary condition $\bar{m}(0)=0$ implies $c_{1}=0$. Upon inserting
this into the differential equation (29) for $\bar{p}(\bar{r})$ (and setting
$\bar{p}^{\prime}(\bar{r})=0$), we arrive at a condition corresponding to
$\bar{p}(\bar{r})=\bar{p}_{c}=\bar{p}_{\mathrm{crit}}$:

$\displaystyle\bar{p}_{\mathrm{crit}}=\frac{1}{2\pi(9\pi-32\bar{\alpha}\bar{\lambda})^{2}}$
$\displaystyle\Bigg{[}256\bar{\alpha}(2\pi\bar{\alpha}+3)\bar{\lambda}^{2}+8\pi(2\pi\bar{\alpha}(8\pi\bar{\alpha}+21)-9)\bar{\lambda}+9\pi^{3}(4\pi\bar{\alpha}-3)$
(56)
$\displaystyle-\textrm{sgn}(\lambda)6\sqrt{\bar{\lambda}(32\bar{\alpha}\bar{\lambda}+\pi(8\pi\bar{\alpha}-3))^{2}\left(16(2\pi\bar{\alpha}+1)\bar{\lambda}+\pi^{2}(8\pi\bar{\alpha}+3)\right)}\Bigg{]}.$

For fixed $\lambda>0$, the critical pressure takes the form shown in Figure
11. If $\lambda=0$ then $\bar{p}_{crit}=\frac{4\pi\bar{\alpha}-3}{18}$, so the
limiting intersection point is $\bar{\alpha}=\frac{3}{4\pi}=0.238$.

For a given $\alpha-\lambda$ pair only central pressures above $p_{crit}$ give
valid quark star solutions. Since pressure must always be non-negative, a
negative critical pressure implies that all choices of positive central
pressure yield valid quark star solutions, and so to find solutions we can
start scanning from an arbitrarily small positive pressure. However once
$\alpha$ is sufficiently large, there is a threshold value of the central
pressure below which there are no solutions. The solutions in these cases when
$p>p_{\mathrm{crit}}$ are found to start at the black hole horizon and thus
are never physically relevant. More simply, $\bar{\alpha}/\bar{\lambda}$
combinations for which $p_{\mathrm{crit}}>0$ have no physical quark star
solutions, and thus in a sense this condition marks the physically relevant
boundary of the parameter space. The requirement of a larger $p_{c}$ to avoid
pressure divergence is somewhat counter-intuitive, although one can see from
equation (54) that a larger critical pressure corresponds to a larger mass
which is more gravitationally attractive.

To demonstrate this, for test values $\alpha=1,\lambda=1$ we find that
$p_{0}=0.11249510580691144$. We plot $p(r)$ for this trio of parameters in
Figure 10b, where as expected we see a flat line. If we instead we plot $p(r)$
for $p_{0}=0.112$ and $p_{0}=0.113$, (figures 12b and 12a), we can see what
happens on either side of criticality.

(a) A “typical” monotonically decreasing pressure profile for a quark star
with a well-defined surface at $\bar{p}(\bar{r})=0$.  (b) $p_{0}=p_{c}$ for
our test case when $\lambda>0$. As expected we find a constant, flat line.

Figure 10: Critical center pressure analysis for $\lambda>0$. Figure 11:
Critical pressure as a function of $\bar{\alpha}$ when $\lambda>0$, plotted
for $\bar{\lambda}=0.1,0.5,1,10,100$ (with the darkest lines corresponding to
the smallest $\bar{\lambda}$).

(a) $p_{0}=p_{\mathrm{crit}}-\delta$. (b) $p_{0}=p_{\mathrm{crit}}+\delta$.

Figure 12: Either side of criticality for our test case when $\lambda>0$.

On the other hand, when $\lambda<0$ the critical pressure solution has a
divergence at $\bar{\alpha}=\frac{9\pi}{32\bar{\lambda}}$, which can be seen
in Figure 13. For values of $\bar{\alpha}$ smaller than the divergence point
the critical pressure is defined analogously to the $\lambda>0$ case (namely
only values above the curve are allowed). Past the divergence point there are
no quark star solutions for any central pressures, and thus in the $\lambda<0$
regime we can say that
$\bar{\alpha}_{\mathrm{max}}=\frac{9\pi}{32\bar{\lambda}}$.

Figure 13: Critical pressure as a function of $\alpha$ when $\lambda<0$,
plotted for $\bar{\lambda}=0.1,0.5,1,2,10$ (with the darkest lines
corresponding to the smallest $\bar{\lambda}$). Note that mathematically these
functions past their discontinuity, though we find that this region does not
produce quark star solutions and thus the corresponding curves have been
removed (ie. the value of $\alpha$ at which $p_{\mathrm{crit}}$ diverges is
the end of the solution space).


# Bifurcation of dividing surfaces constructed from a pitchfork bifurcation of
periodic orbits in a symmetric potential energy surface with a post-
transition-state bifurcation

Matthaios Katsanikas<EMAIL_ADDRESS>School of
Mathematics, University of Bristol,
Fry Building, Woodland Road, Bristol, BS8 1UG, United Kingdom.  Makrina
Agaoglou<EMAIL_ADDRESS>Instituto de Ciencias Matemáticas, CSIC,
C/Nicolás Cabrera 15, Campus Cantoblanco, 28049 Madrid, Spain  Stephen Wiggins
<EMAIL_ADDRESS>School of Mathematics, University of Bristol,
Fry Building, Woodland Road, Bristol, BS8 1UG, United Kingdom.

###### Abstract

In this work we analyze the bifurcation of dividing surfaces that occurs as a
result of a pitchfork bifurcation of periodic orbits in a two degrees of
freedom Hamiltonian System. The potential energy surface of the system that we
consider has four critical points:two minima, a high energy saddle and a lower
energy saddle separating two wells (minima). In this paper we study the
structure, the range, and the minimum and maximum extent of the periodic orbit
dividing surfaces of the family of periodic orbits of the lower saddle as a
function of the total energy.

Keywords: Bifurcation, periodic orbit dividing surfaces, Valley Ridge
Inflection Point Potential, Phase space structure, Chemical reaction dynamics.

##  I Introduction

The aim of this paper is to study a bifurcation of periodic orbit dividing
surfaces that occurs as a result of a pitchfork bifurcation of periodic
orbits. For this reason we choose a pitchfork bifurcation of periodic orbits
that takes place in a Hamiltonian system with two degrees of freedom. The
potential energy surface (PES) that we use is symmetric about the $x$-axis and
it has four critical points: two index-1 saddles and two minima. One index-1
saddle (the upper index-1 saddle) has higher energy than the other index-1
saddle (the lower index-1 saddle) and the lower index-1 saddle separates the
two potential wells (the upper well and the lower well). Between the two
index-1 saddles is a valley ridge inflection (VRI) point, where the shape of
the potential near the VRI point changes from a valley to a ridge. This
potential has been studied recently ([1, 2, 3, 4, 5, 6, 7]).

Dividing surfaces are surfaces of one less dimension than that of the
potential energy surface. These surfaces are important for the Transition
State Theory (TST) (see [8, 9]) in chemical reaction dynamics. In Hamiltonian
systems with two degrees of freedom, the periodic orbit dividing surfaces are
2-dimensional surfaces in the energy surface. The construction of the periodic
orbit dividing surfaces is constructed through the classical method of [10,
11, 12, 13, 14] in Hamiltonian systems with two degrees of freedom. The
construction of the periodic orbit dividing surfaces in Hamiltonian with three
or more degrees of freedom can be done through the method of [15, 16].

An invariant, closed, oriented, co-dimension-2 submanifold of an energy level
that can be spanned by two compact codimension-1 submanifolds of
unidirectional flux whose union, a dividing surface, separates locally the
energy level into two components and it has no local recrossing is called a
transition state (see [17]). The bifurcations of transition states (unstable
periodic orbits and NHIMs) in chemical reaction dynamics have been used in
order to understand: A) The breakdown of transition state theory in collinear
$HgI_{2}\rightarrow HgI+I$ reaction [18] and the hydrogen exchange reaction
(see [19]). B) The influence of no-return transition states in the reaction
rate of hydrogen exchange reaction [20]. C) The loss and regain of normal
hyperbolicity of the transition states in a three-dimensional Hamiltonian
model of the hydrogen exchange reaction [21]. D) The bridge between reactant
and product channels in the collinear $FH_{2}$ reactive system [22]. E) The
study of systems in which as the energy varies the transition state becomes
singular and then regains the normal hyperbolicity with a change in
differentiability as described through Morse bifurcations (see [17]). F) The
separation of reactants and products and the flux through the dividing surface
(see [23]). G) The transport of trajectories from the region of the higher
index-1 saddle to the region of the top and bottom well in a symmetric
potential energy surface with a post-transition-state bifurcation through a
pitchfork bifurcation (see [2, 3]).

The basic and initial motivation of this work is to continue our previous work
([2, 3]) by investigating the effects of a pitchfork bifurcation on the
structure of the dividing surfaces of the top and bottom unstable periodic
orbits (that bifurcate from the family of the lower saddle through a pitchfork
bifurcation) that are responsible for the transport of the trajectories from
the region of the higher index-1 saddle to the region of the top and bottom
well (see [2, 3]). In all the papers that were mentioned in the previous
paragraph except ([17],[23]), the studies were focused on the applications of
the bifurcations of periodic orbits (see [24] for a review of the many types
of bifurcations in Hamiltonian systems with two and three degrees of freedom
in chemical reaction dynamics) and NHIMs in different settings of chemical
reactions, but they did not focus on the effects of these bifurcations on the
structure and geometry of dividing surfaces. In this paper, we deal with this
problem and in particular with the evolution of the geometry of dividing
surfaces as the energy increases and the family of periodic orbits of the
lower saddle undergoes a pitchfork bifurcation in a symmetric potential energy
surface with a post-transition-state bifurcation ([2, 3]). These surfaces are
spheres or ellipsoids before and after a bifurcation, as is also shown for
integrable cases of Hamiltonian systems [23]). They are presented as disks in
the 3D projection $(x,y,p_{x})$ of the phase space, before and after a
pitchfork bifurcation in a non-integrable Hamiltonian system with two degrees
of freedom ([25]). In our case, we will study in a more systematic way the
change of the geometry of the dividing surfaces, before and after a pitchfork
bifurcation, investigating the structure of the dividing surfaces in all 3D
projections of the phase space of a non-integrable Hamiltonian system with two
degrees of freedom. We will also describe in detail the structure of the
dividing surfaces in a system in which the structure of dividing surfaces is
more complicated (the dividing surfaces are initially not disks in the 3D
projection $(x,y,p_{x})$ but filamentary structures) than that of the dividing
surfaces that are found in the Hamiltonian model of [25]. Furthermore, we will
follow the complex distortion of these surfaces that are initially ellipsoids
or filamentary structures in the 3D projections of the phase space as we
change the energy. We emphasize the fact that we investigated the distortion
of the dividing surfaces (that are initially ellipsoids or filamentary
structures) through a pitchfork bifurcation of periodic orbits, and not
situations where the dividing surfaces change topology and become tori or
other topological structures though a Morse bifurcation (see [17]). The study
of topological changes of this kind (through a Morse bifurcation) is an
interesting topic but out of the scope of this paper. We study, also the
evolution of other geometrical characteristics, such as the range and the
maximum and minimum value of the dividing surfaces.

We describe the PES of our system in section II. We describe the pitchfork
bifurcation of periodic orbits of the lower saddle in section III. Then in
section IV, we describe our results about the structure of the diving
surfaces. Finally, we discuss our results and we present our conclusions in
the last section V.

##  II Model

The analytical form of the PES that was inspired by [1] is the following:

$V(x,y)=\frac{8x^{3}}{3}-4x^{2}+\frac{y^{2}}{2}+xy^{2}(y^{2}-2).$ (1)

The PES has an upper index-1 saddle (exit/entrance channel) and a lower
index-1 saddle, which is an energy barrier separating two potential wells.
Moreover, the PES is symmetric about the x axis. A three dimensional plot of
the PES is given in Fig. 1. The two dimensional (2D) Hamiltonian is given by

$H(x,y,p_{x},p_{y})=\frac{p_{x}^{2}}{2}+\frac{p_{y}^{2}}{2}+V(x,y),$ (2)

and the corresponding Hamiltonian vector field (i.e. Hamilton’s equations of
motion) is:

$\begin{cases}\dot{x}=p_{x},\\\ \dot{y}=p_{y},\\\
\dot{p}_{x}=8x(1-x)+y^{2}(2-y^{2}),\\\
\dot{p}_{y}=y\left(4x(1-y^{2})-1\right).\end{cases}$ (3)

Hamilton’s equations conserve the total energy, which we will denote as $E$
throughout the text.

Figure 1: Plot of the PES given in Eq. (1).

##  III Periodic Orbits

In this section, we present a pitchfork bifurcation of the family of the lower
saddle. As we explained in the previous section, we have two index-1 saddles,
one for higher values of energy and one for lower values of energy. According
to the Lyapunov subcenter theorem we have a family of periodic orbits
associated with every index-1 saddle (see [26, 27, 28]). In this section we
follow the family of periodic orbits associated with the lower index-1 saddle,
increasing the energy above the energy of the lower index-1 saddle. The
periodic orbits of this family lie on the $x$-axis (see Fig. 8).

In the characteristic diagram (a diagram that shows the evolution of a
coordinate of periodic orbits of a family in a Poincaré section versus the
energy of the system - see Fig. 2) we see that the family of periodic orbits
associated with the lower saddle undergoes a pitchfork bifurcation at
$E=-0.00056$. At this value the family of periodic orbits associated with the
lower index-1 saddles, that was initially unstable, becomes stable. At this
point we have two new unstable families of periodic orbits which we refer to
as the top and bottom families of periodic orbits (see more details for the
role of the these periodic orbits in the dynamics of the system in [2, 3, 6,
5]). We observe, in Figs. 9, 10, that by increasing the energy the periodic
orbits of these families become more elongated in the direction of the
$y$-axis. In these figures we see that the POs of the top and bottom families
are converging to the line $y=0$ from above and below, respectively.

Figure 2: The evolution of the coordinate $y$ of periodic orbits on the
Poincaré section $x=0.05$ versus energy (characteristic diagram). The family
of periodic orbits of the lower saddle is depicted by black color and its
bifurcating families, the families of the top and bottom unstable periodic
orbits, are depicted by red and green color, respectively.

##  IV Results

Initially (for values of energy above the energy of the lower saddle, $E=-4/3$
and before the pitchfork bifurcation for $E=-0.00056$) we have the family of
periodic orbits of the lower saddle that gives rise to a family of dividing
surfaces. For the construction of the periodic orbit dividing surfaces we use
he classical algorithm of [10, 11, 14, 29, 13], that is presented in our
previous papers (see ([6, 30]). After the pitchfork bifurcation of the
periodic orbits of the family of the lower saddles (see the previous section)
we have two new bifurcating families, the top and bottom unstable periodic
orbits. These families give rise to two new families of dividing surfaces. In
this section, we study the evolution of the structure and the range and of the
minimum and maximum of the dividing surfaces that are constructed from the
periodic orbits of the lower saddle and its bifurcations as energy is varied:

  1. 1.

The structure of the dividing surfaces: The periodic orbit dividing surfaces
that are constructed from the periodic orbits of the lower saddles are
presented as ellipsoids in the $(x,p_{x},p_{y})$ projection (see Fig. 3). This
ellipsoid becomes more stretched and sharp as we increase the energy (see the
2D projections $(x,p_{x})$ of the dividing surfaces in Fig. 11). Furthermore,
these dividing surfaces are presented as filamentary structures, that lie on
the plane $y=0$, in the $(x,y,p_{x})$ and $(x,y,p_{y})$ subspaces of the phase
space (see panel B of Fig. 3 and panel A of Fig. 5) and in the
$(y,p_{x},p_{y})$ subspace of the phase space (see panel B of Fig. 5). This
filamentary structure of these surfaces is because of the morphology of the
periodic orbits of the lower saddles in the configuration space (it lies on
the line $y=0$ \- see the previous section). This is the reason that the 2D
projections $(y,p_{y})$ of the dividing surfaces of the family of lower saddle
(Fig. 14) is the line $y=0$ for all energies.

As the energy is increased above that of the pitchfork bifurcation, in
addition to the dividing surfaces of the family of the lower saddles, we have
two other families of dividing surfaces. One family corresponds to the
dividing surfaces of the top unstable periodic orbits and the other to the
dividing surfaces of the bottom unstable periodic orbits. The dividing
surfaces of these families are on both sides of the dividing surfaces of the
family corresponding to the lower saddles in the 3D projections $(x,y,p_{x})$,
$(x,y,p_{y})$ and $(y,p_{x},p_{y})$ (see panels C and D of Fig. 5). The
dividing surfaces of the bifurcating families in the 3D projections,
$(x,y,p_{x})$ and $(x,y,p_{y})$, are presented as filamentary structures that
are curved in the direction of positive $y$-values (the dividing surfaces of
the top unstable periodic orbits) or on the direction of negative $y$-values
(the dividing surfaces of the bottom unstable periodic orbits). The structure
of these dividing surfaces in the 3D subspace $(y,p_{x},p_{y})$ is similar to
that of an ellipsoid that is a little sharper in one direction (see for
example Fig. 4). In addition, the dividing surfaces of the family of the lower
saddle coincide with the dividing surfaces of the two bifurcating families in
the $(x,p_{x},p_{y})$ subspace of the phase space and they have the same
structure (an ellipsoid) as that of panel A of Fig. 3.

For values of energy larger than 0 (the energy of the upper saddle) the family
of periodic orbits of the lower saddle vanishes. Consequently, the family of
dividing surface of the lower saddle vanishes too. This means that the
dividing surfaces of the family of the lower saddle do not exist between the
dividing surfaces of the bifurcating families in the 3D projections
$(x,y,p_{x})$, $(x,y,p_{y})$ and $(y,p_{x},p_{y})$ (see the panels E and F of
Fig. 5). We notice also that the distance between the dividing surfaces of the
top and bottom unstable periodic orbits, in the 3D projections $(x,y,p_{x})$,
$(x,y,p_{y})$ and $(y,p_{x},p_{y})$, increases as the energy increases
(compare the panels C and D of Fig. 5 with the panels E and F of Fig. 5).

As we increase the energy, the dividing surfaces of the bifurcating families
in the 3D subspace $(x,p_{x},p_{y})$ become more round/ellipsoidal as the
energy increases (see the 2D projections $(x,p_{x})$ of the dividing surfaces
of the top and bottom unstable periodic orbits in Figs. 12, 13). Furthermore,
we observe a decrease of the sharpness of the ellipsoidal dividing surfaces of
the bifurcating families in the 3D subspace $(y,p_{x},p_{y})$, as we increase
the energy (see the 2D projections $(y,p_{y})$ of the dividing surfaces of the
top and bottom unstable periodic orbits in Figs. 15, 16).

  2. 2.

Maximum and Minimum of the dividing surfaces: In Fig. 6 we present in the left
column the minimum of the dividing surfaces with respect to the $x$
coordinate, $y$ coordinate and $p_{x}$ coordinate (panels A, C and E
respectively) of the family of the lower saddle (black line), top family (red
line) and bottom family (green line), respectively. Similarly, in the right
column of the figure we show the maximum of the dividing surfaces with respect
to the $x$ coordinate, $y$ coordinate and $p_{x}$ coordinate (panels B, D and
F respectively) of the family of the lower saddle (black line), top family
(red line) and bottom family (green line), respectively. Note that we did not
show the minimum and maximum of $p_{y}$ coordinate because it is determined
from the other three coordinates through the Hamiltonian of the system. We
observe that in the panel A the minimum of the family of the lower saddle is
decreasing as the energy increases, as does does the minimum of the top and
bottom families so that they coincide. In the panels B, E and F we notice that
there is a continuation of the minima and maxima of the family of the lower
saddle to the other two bifurcating families (they coincide with each other).
Finally in the panel C we show that the minimum of the family of the lower
saddle is always zero whilst the minimum of the top family is increasing and
the minimum of the bottom family is decreasing. The same happens in the panel
D for the maximum.

  3. 3.

The range of the dividing surfaces: In Figure 7 we illustrate the range of the
$x$ coordinate, $y$ coordinate and $p_{x}$ coordinate (panels A, B and C,
respectively) of the dividing surfaces of the family of the lower saddle
(black line), top family (red line) and bottom family (green line),
respectively. Note that we did not present the range of $p_{y}$ coordinate
because it is determined from the other three coordinates through the
Hamiltonian of the system. In the panels A, B and C the range of the
coordinates of the top family coincide with the range of the same coordinates
of the bottom family. Moreover the range of the $y$ coordinate of the family
of the lower saddle is zero for all energies (see panel B). Finally, the
continuation of the range of the $p_{x}$ coordinate of the families is evident
as the energy varies (see panel C).

A) B)

Figure 3: 3D projections $(x,p_{x},p_{y})$ (panel A) and $(x,y,p_{x})$ (panel
B) of the periodic orbit dividing surface of the family of the lower saddle
for $E=-0.2$.  Figure 4: 3D projection $(x,p_{x},p_{y})$ of the dividing
surface that is constructed by the family of the top unstable periodic orbits
for $E=-0.00055$.

A) B)
C) D)
E) F)

Figure 5: 3D projections $(x,y,p_{x})$ (left column of the panels) and
$(y,p_{x},p_{y})$ (right column of the panels) of the periodic orbit dividing
surface of the family of the lower saddle (with black color) and of the
dividing surfaces that are constructed from the periodic orbits of the top
(with red color) and bottom (with green color) unstable periodic orbits for
$E=-0.2$ (first row), $E=-0.00055$ (second row) and $E=0.2$ (third row). We
encountered the same structures, as these in the 3D subspace $(x,y,p_{x})$, in
the 3D subspace $(x,y,p_{y})$.

A) B)
C) D)
E) F)

Figure 6: Diagrams for the A) mininum and B) maximum of the X coordinate of
the central (black), top (red) and bottom (green) families versus the energy,
C) mininum and D) maximum of the y coordinate of the central, top and bottom
families versus the energy and E) mininum and F) maximum of the PX coordinate
of the central, top and bottom families versus the energy.

A)
B)
C)

Figure 7: A) Range of the X coordinate of the central family (black line), top
family (red line) and bottom family (green line) versus the energy, B) Range
of the Y coordinate of the central family (black line), top family (red line)
and bottom family (green line) versus the energy and C) Range of the PX
coordinate of the central family (black line), top family (red line) and
bottom family (green line) versus the energy

A) B)
C) D)

Figure 8: 2D projections of the periodic orbits of the central family in the
configuration space for energies A) $-1.2$, B) $-0.6$, C) $-0.2$ and D)
$-0.00055$

A) B)
C) D)

Figure 9: 2D projections of the periodic orbits of the top family in the
configuration space for energies A) $-0.00055$, B) $0.2$, C) $0.4$ and D)
$0.65$

A) B)
C) D)

Figure 10: 2D projections of the periodic orbits of the bottom family in the
configuration space for energies A) $-0.00055$, B) $0.2$, C) $0.4$ and D)
$0.65$

A) B)
C) D)

Figure 11: 2D projections of the dividing surfaces of the central family in
the X-PX space for energies A) $-1.2$, B) $-0.6$, C) $-0.2$ and D) $-0.00055$

A) B)
C) D)

Figure 12: 2D projections of the dividing surfaces of the top family in the
X-PX space for energies A) $-0.00055$, B) $0.2$, C) $0.4$ and D) $0.65$

A) B)
C) D)

Figure 13: 2D projections of the dividing surfaces of the bottom family in the
X-PX space for energies A) $-0.00055$, B) $0.2$, C) $0.4$ and D) $0.65$

A) B)
C) D)

Figure 14: 2D projections of the dividing surfaces of the central family in
the Y-PY space for energies A) $-1.2$, B) $-0.6$, C) $-0.2$ and D) $-0.00055$

A) B)
C) D)

Figure 15: 2D projections of the dividing surfaces of the top family in the
Y-PY space for energies A) $-0.00055$, B) $0.2$, C) $0.4$ and D) $0.65$

A) B)
C) D)

Figure 16: 2D projections of the dividing surfaces of the bottom family in the
Y-PY space for energies A) $-0.00055$, B) $0.2$, C) $0.4$ and D) $0.65$

##  V Conclusions

We studied a bifurcation of periodic orbit dividing surfaces that occurs
because of a pitchfork bifurcation of periodic orbits in a nonlinear
Hamiltonian system with two degrees of freedom. We investigated the structure,
the range and the minimum and maximum of the periodic orbit dividing surfaces
of the family of periodic orbits of the lower saddle of our system and of its
bifurcating families of periodic orbits. Our main conclusions are:

  1. 1.

The dividing surfaces of the family of periodic orbits of the lower saddle are
presented as ellipsoids and filamentary structures in the 3D subspaces of the
phase space. The ellipsoids become sharper as we increase the energy.

  2. 2.

The dividing surfaces of the bifurcating families are presented as ellipsoids
with a sharp edge or as curved filamentary structures. The sharpness of the
ellipsoids decreases as we increase the energy.

  3. 3.

The dividing surfaces of the bifurcating families coincide with the dividing
surfaces of the family of the lower saddle in the $(x,p_{x},p_{y})$ subspace
of the phase space.

  4. 4.

The dividing surfaces of the bifurcating families are on both sides of the
dividing surfaces of the family of the lower saddle in the 3D subspaces of the
phase space $(x,y,p_{x})$, $(x,y,p_{y})$ and $(y,p_{x},p_{y})$. This happens
until the disappearance of the dividing surfaces of the family of lower
saddles for values of energy from 0 (for which we have the appearance of the
high energy saddle) and above. As we increase the energy, the dividing
surfaces of the bifurcating families move away from each other.

  5. 5.

The ranges of the dividing surfaces of the bifurcating families coincide with
each other and they follow the evolution of the range of the dividing surfaces
of the lower saddle (except the range in the $y$ coordinate where the range of
the dividing surfaces of the lower saddle is 0).

  6. 6.

The minimum and the maximum of the dividing surfaces coincide with each other
with respect to the $x$ and $p_{x}$ coordinate (the $p_{y}$ coordinate is
determined through the Hamiltonian of the system) and they follow, in these
cases, the evolution of the minimum and maximum of the dividing surfaces of
the lower saddles.

## Acknowledgments

The authors acknowledge the financial support provided by the EPSRC Grant No.
EP/P021123/1 and MA also acknowledges support from the grant “Juan de la
Cierva Incorporación” (IJC2019-040168-I).

## References

  * Collins  _et al._ [2013] P. Collins, B. K. Carpenter, G. S. Ezra, and S. Wiggins, Nonstatistical dynamics on potentials exhibiting reaction path bifurcations and valley-ridge inflection points, The Journal of Chemical Physics 139, 154108 (2013).
  * Katsanikas  _et al._ [2020] M. Katsanikas, V. J. García-Garrido, M. Agaoglou, and S. Wiggins, Phase space analysis of the dynamics on a potential energy surface with an entrance channel and two potential wells, Physical Review E 102, 012215 (2020).
  * Agaoglou  _et al._ [2020] M. Agaoglou, V. J. García-Garrido, M. Katsanikas, and S. Wiggins, The phase space mechanism for selectivity in a symmetric potential energy surface with a post-transition-state bifurcation, Chemical Physics Letters 754, 137610 (2020).
  * García-Garrido  _et al._ [2020] V. García-Garrido, M. Katsanikas, M. Agaoglou, and S. Wiggins, Tuning the branching ratio in a symmetric potential energy surface with a post-transition state bifurcation using external time dependence, Chemical Physics Letters 754, 137714 (2020).
  * Agaoglou  _et al._ [2021] M. Agaoglou, M. Katsanikas, and S. Wiggins, The influence of a parameter that controls the asymmetry of a potential energy surface with an entrance channel and two potential wells, Regular and Chaotic Dynamics , (was submitted) (2021).
  * Douglas  _et al._ [2021] H. Douglas, M. Katsanikas, M. Agaoglou, and S. Wiggins, The time evolution of the trajectories after the selectivity in a symmetric potential energy surface with a post-transition-state bifurcation, Regular and Chaotic Dynamics , (was submitted) (2021).
  * Crossley  _et al._ [2021] R. Crossley, M. Agaoglou, M. Katsanikas, and S. Wiggins, From Poincaré maps to Lagrangian descriptors: The case of the valley ridge inflection point potential, Regular and Chaotic Dynamics 26, 147 (2021).
  * Wigner [1938] E. Wigner, The transition state method, Transactions of the Faraday Society 34, 29 (1938).
  * Waalkens  _et al._ [2007] H. Waalkens, R. Schubert, and S. Wiggins, Wigner’s dynamical transition state theory in phase space: classical and quantum, Nonlinearity 21, R1 (2007).
  * Pechukas and McLafferty [1973] P. Pechukas and F. J. McLafferty, On transition-state theory and the classical mechanics of collinear collisions, The Journal of Chemical Physics 58, 1622 (1973).
  * Pechukas and Pollak [1977] P. Pechukas and E. Pollak, Trapped trajectories at the boundary of reactivity bands in molecular collisions, The Journal of Chemical Physics 67, 5976 (1977).
  * Pechukas and Pollak [1979] P. Pechukas and E. Pollak, Classical transition state theory is exact if the transition state is unique, The Journal of Chemical Physics 71, 2062 (1979).
  * Pollak [1985] E. Pollak, Periodic orbits and the theory of reactive scattering, Theory of Chemical Reaction Dynamics 3, 123 (1985).
  * Pechukas [1981] P. Pechukas, Transition state theory, Annual Review of Physical Chemistry 32, 159 (1981).
  * Katsanikas and Wiggins [2021a] M. Katsanikas and S. Wiggins, The generalization of the periodic orbit dividing surface for Hamiltonian systems with three or more degrees of freedom in chemical reaction dynamics - I, International Journal of Bifurcation and Chaos , in press (2021a).
  * Katsanikas and Wiggins [2021b] M. Katsanikas and S. Wiggins, The generalization of the periodic orbit dividing surface for Hamiltonian systems with three or more degrees of freedom in chemical reaction dynamics - II, International Journal of Bifurcation and Chaos , in press (2021b).
  * MacKay and Strub [2014] R. S. MacKay and D. C. Strub, Bifurcations of transition states: Morse bifurcations, Nonlinearity 27, 859 (2014).
  * Burghardt and Gaspard [1995] I. Burghardt and P. Gaspard, The molecular transition state: From regular to chaotic dynamics, The Journal of Physical Chemistry 99, 2732 (1995).
  * Inarrea  _et al._ [2011] M. Inarrea, J. F. Palacian, A. I. Pascual, and J. P. Salas, Bifurcations of dividing surfaces in chemical reactions, J. Chem. Phys. 135, 014110 (2011).
  * Li  _et al._ [2009] C. B. Li, M. Toda, and T. Komatsuzaki, Bifurcation of no-return transition states in many-body chemical reactions, J. Chem. Phys. 130, 124116 (2009).
  * Allahem and Bartsch [2012] A. Allahem and T. Bartsch, Chaotic dynamics in multidimensional transition states, J. Chem. Phys. 137, 214310 (2012).
  * Founargiotakis  _et al._ [1997] M. Founargiotakis, S. Farantos, C. Skokos, and G. Contopoulos, Bifurcation diagrams of periodic orbits for unbound molecular systems: Fh2, Chemical Physics Letters 277, 456 (1997).
  * Mauguiere  _et al._ [2013] F. A. Mauguiere, P. Collins, G. S. Ezra, and S. Wiggins, Bifurcations of normally hyperbolic invariant manifolds in analytically tractable models and consequences for reaction dynamics, International Journal of Bifurcation and Chaos 23, 1330043 (2013).
  * Farantos  _et al._ [2009] S. C. Farantos, R. Schinke, H. Guo, and M. Joyeux, Energy localization in molecules, bifurcation phenomena, and their spectroscopic signatures: The global view, Chemical Reviews 109, 4248 (2009).
  * Lyu  _et al._ [2021] W. Lyu, S. Naik, and S. Wiggins, Hamiltonian pitchfork bifurcation in transition across index-1 saddles, arXiv preprint arXiv:2102.10933 (2021).
  * Rabinowitz [1982] P. H. Rabinowitz, Periodic solutions of Hamiltonian systems: a survey, SIAM Journal on Mathematical Analysis 13, 343 (1982).
  * Weinstein [1973] A. Weinstein, Normal modes for nonlinear Hamiltonian systems, Inventiones Mathematicae 20, 47 (1973).
  * Moser [1976] J. Moser, Periodic orbits near an equilibrium and a theorem by alan weinstein, CPAM 29, 727 (1976).
  * Pollak and Pechukas [1978] E. Pollak and P. Pechukas, Transition states, trapped trajectories, and classical bound states embedded in the continuum, The Journal of Chemical Physics 69, 1218 (1978).
  * Ezra and Wiggins [2018] G. S. Ezra and S. Wiggins, Sampling phase space dividing surfaces constructed from normally hyperbolic invariant manifolds (NHIMs), The Journal of Physical Chemistry A 122, 8354 (2018).