# Radial Networks: Dynamic Layer Routing
for High-Performance Large Language Models

Jordan Dotzel Yash Akhauri11footnotemark: 1 Ahmed S. AbouElhamayed Carly Jiang
Mohamed Abdelfattah Zhiru Zhang
Cornell University
{dotzel, ya255, afa55, cjj43, mohamed<EMAIL_ADDRESS>Equal Contribution

###### Abstract

Large language models (LLMs) often struggle with strict memory, latency, and
power demands. To meet these demands, various forms of dynamic sparsity have
been proposed that reduce compute on an input-by-input basis. These methods
improve over static methods by exploiting the variance across individual
inputs, which has steadily grown with the exponential increase in training
data. Yet, the increasing depth within modern models, currently with hundreds
of layers, has opened opportunities for dynamic layer sparsity, which skips
the computation for entire layers. In this work, we explore the practicality
of layer sparsity by profiling residual connections and establish the
relationship between model depth and layer sparsity. For example, the residual
blocks in the OPT-66B model have a median contribution of 5% to its output. We
then take advantage of this dynamic sparsity and propose Radial Networks,
which perform token-level routing between layers guided by a trained router
module. These networks can be used in a post-training distillation from
sequential networks or trained from scratch to co-learn the router and layer
weights. They enable scaling to larger model sizes by decoupling the number of
layers from the dynamic depth of the network, and their design allows for
layer reuse. By varying the compute token by token, they reduce the overall
resources needed for generating entire sequences. Overall, this leads to
larger capacity networks with significantly lower compute and serving costs.

##  1 Introduction

Figure 1: Radial Networks – Radial Networks generalize sequential networks for
higher accuracy and performance. They take advantage of significant dynamic
layer sparsity within modern LLMs and invoke only a subset of model layers for
each token. They reduce the average network depth, lower model latency, and
provide a more scalable neural architecture.

Large language and vision models have recently achieved state-of-the-art
performance across various tasks, yet due to their large computational
requirements, they struggle with strict memory, latency, and power demands. As
these transformers grow larger, they create opportunities for dynamic layer
sparsity, which can skip individual layers on an input-by-input basis, as
shown in Figure 2. For instance, our residual block profiling in Section 4
suggests that modern state-of-the-art transformers likely have a median
contribution around 1% to the output at each block, and that these
contributions are dynamic, varying token by token. This type of sparsity was
impractical at smaller scales and with previous neural architectures. At
smaller scales, every layer contributes significantly to the computation for
each input, and with previous architectures, e.g., convolutional neural
networks (CNNs), models change their intermediate dimensions throughout their
depth and skipping causes dimensional mismatches.

This work shows that the layer contributions vary among models and tasks, and
often the earlier layers of the network contribute more than the later layers.
This indicates that early-exit methods, which dynamically prune the later
layers in the network, often focus on the wrong set of layers. This dynamic
contribution can be exploited at the token-level if it can be predicted
accurately and efficiently at runtime. This work explores the opportunities
for dynamic sparsity within the modern transformers by focusing on the OPT
family of models [9] for language and ViT models [3] for vision. It profiles
the residual blocks to quantify the importance of each intermediate layer to
its output and then highlights trends across model size and block types. Then,
it inserts oracles at every layer to calculate various accuracy proxies and
simulate greedy decisions on which layers to dynamically skip per token.

##  2 Related Work

Figure 2: Dynamic Layer Sparsity – As transformers grow larger, each layer
contributes less to the output and shows significant variation on a token by
token basis. Dynamically pruning these layers allows for models to grow
significantly without corresponding increases in model latency. Early
profiling results suggest that individual layers contribute around 1% within
modern state-of-the-art language models.

Sparsity research with deep neural networks has a long history, and broadly
can be categorized in terms of granularity, structure, and mode (static vs.
dynamic) [4]. Figure 3 shows sparsity granularity, beginning with bits that
construct parameter elements, elements that build blocks, and blocks that form
layers. As the unit becomes larger, it becomes more difficult to arbitrarily
prune without accuracy loss yet easier to accelerate with modern hardware. For
instance, unstructured element sparsity in weights leads to high compression
levels while maintaining model accuracy, yet it requires specializing sparse
accelerators to translate compression into end-to-end speedup.

In addition, the sparsity mode can either be static or dynamic. Static
sparsity leads to more regular patterns that can be optimized by compilers and
simpler architecture that do not need additional sparsity predictors, yet it
must apply to all inputs together. In contrast, dynamic sparsity can take
advantage of input-dependent characteristics to increase model accuracy at
higher levels of compression. This work focuses on dynamic layer sparsity,
which can take advantage of the recent explosion in model depth within
language models.

###  2.1 Dynamic Sparsity

Figure 3: Sparsity Granularity – Bits form the basis for elements (weight or
activations), which create blocks (rows, columns, heads), which then form
individual layers. This leads to a sparsity spectrum where the smaller units
are easier to prune without accuracy loss yet more difficult to accelerate.
Layer sparsity has the highest potential for inference speedup and largest
support within current hardware.

Multiple prior works have proposed dynamic sparsity to accelerate DNNs across
granularities. For example, Channel Gating introduced a method for dynamic
channel sparsity that reduced the compute of CNN workloads by up to $8\times$
without significant accuracy loss [5]. Precision Gating continued this line of
research by applying dynamic sparsity at the bit level to reduce the required
compute by up to $3\times$ [10]. Later, DejaVu applied a similar approach
within LLMs to induce dynamic sparsity on the channels within the FFN layer
and across the heads of the attention layer [6].

###  2.2 Early Exit

In addition to dynamic sparsity along the network width, multiple prior works
have explored sparsity in the depth dimension. For instance, early-exit DNNs
use dynamic sparsity along the depth dimension by allowing the computation to
exit prematurely at fixed points within the network [7, 1]. This process must
be trained end-to-end using a joint loss function that weights the
contributions from each early-exit layer. However, this work shows that in
many models, the earlier layers in the model often contribute more, and
therefore early-exits are significantly more difficult to apply post-training.

##  3 Layer Sparsity

Figure 4: Residual Blocks – There are two types of residual blocks within
transformers, attention (ATT) and feed-forward network (FFN). These blocks
offer natural points to profile layer strength since block inputs and outputs
are combined at a single point. To establish an upper bound on the
effectiveness of dynamic layer sparsity, oracles are inserted before each
block that know the layer contribution beforehand.

Transformer layers contain two residual blocks: attention (ATT) and feed-
forward network (FFN) [8]. These blocks each contain the main residual branch
$R(x)$, which comprises multiple individual layers, and the identity branch
$x$, which bypasses the residual branch and simply returns its input. They
combine these branch outputs together to compute $R(x)+x$, so that during
training the main branch only has to learn the function residual $R(x)-x$.

These blocks offer natural breakpoints within the model to profile and induce
layer sparsity since they already provide skip-connections that have been
trained along with the model. Figure 4 shows a lower-level view of these
blocks within two transformer layers. It shows that the main branch and skip
connections are combined at an addition node before they are passed to the
next block. This structure enables easy profiling of the blocks by measuring
the relative magnitudes into these additions.

This figure also shows the insertion of oracles that can switch on and off the
main branch using various accuracy proxies, such as the residual ratio as
defined in Section 4.1. When they are switched on, the block operates normally
combining the skip and residual branches, and when switched off, only the skip
connection is active. This work focuses on the opportunities for dynamic layer
sparsity and simulates layer skipping by allowing these oracles to have access
to future information.

##  4 Profiling

The primary proxy used by these oracles is the residual ratio, which captures
the relative importance of the main and skip branches. This section uses this
ratio to analyze the layer sparsity within OPT and ViT models with examples
taken from WikiText-2 and COCO. The WikiText-2 examples are packed together to
avoid the use of padding to simulate batch-size one inference. This batch-size
one setting is very common in practice and avoids many complications with
dynamic layer sparsity that arise when using batches of examples.

###  4.1 Residual Ratio

To profile these opportunities for dynamic sparsity, this section defines the
residual ratio $r$ as:

$\text{r}=\frac{\|R(x)\|_{2}}{\|x\|_{2}}$ (1)

This simple quantity captures the contribution of the residual branch, and
acts as an efficient post-training proxy for more expensive metrics, such as
empirical layer sensitivities. For example, a block with a 2% residual ratio
indicates the main branch provides a 2% average contribution at the output,
although there can be large element-wise variance. Therefore, skipping blocks
with ratios this small should have little overall effect on the output of the
network.

###  4.2 Model Size

Figure 5: Dynamic Depth – The deeper layers in the network contribute more
than the earlier layers, except for the very first layers. This relationship
benefits from a routed architecture as opposed to early-exit, since early-exit
skips the deeper layers. In addition, there is significant variance in the
dynamic depth of the model, allowing for token-specific sparsity. Figure 6:
Residual Ratio – As models grow larger, the residual ratio decreases for each
layer, and therefore more layers contribute less to the overall output. Each
data point represents a residual ratio for a single token for both attention
and feed-forward blocks. Figure 7: OPT-13B Routing – The majority of skipped
layers are in the earlier part of the network. A layer skip is assumed when
the residual ratio drops below 5%. ATT and FFN layers are interleaved in the
diagram.

This ratio can be used to understand the relationship between model size and
dynamic sparsity. Figure 6 explores this by plotting the residual ratio across
OPT models for the residual attention or feed-forward block (more plots shown
in Appendix A). Each data point represents a single token during the model
generation phase. It shows that as the model size grows, the ratio
distribution becomes more skewed to the left, indicating that opportunities
for layer sparsity expand with model size. For instance, while the median
residual ratio for OPT-125M is only 20%, it drops to 5.9% for OPT-66B.

In addition, the ratio seems to track the number of model parameters, not just
the number of layers. For example, OPT-2.7B and OPT-6.7B have the same number
of layers, differing only in their hidden dimensions, yet the ratio for
OPT-6.7B continues the decreasing trend. This trend likely continues for even
larger models, making dynamic layer sparsity more practical within modern
state-of-the-art models with greater than one trillion parameters.

###  4.3 Dynamic Depth

Dynamic layer sparsity leads to dynamic depth networks that adjust their depth
based on their model inputs. Figure 9 shows the residual ratio across the
layers of an OPT-13B model (all models shown in Appendix A). All values shown
are mean residual ratios taken across tokens from Wikitext-2 data using a
sequence length of 256. The ratio variance is highlighted in lighter colors
centered around the mean. It demonstrates that the earlier residual blocks
contribute more compared to the later layers, except for the first few layers.
In addition, there is significant variance across tokens across layers
suggesting the opportunity to apply dynamic layer sparsity to only the tokens
with lower ratios.

This figure additionally shows the dynamic depth caused by this layer
sparsity. It assumes oracles that threshold the residual ratio at each block
and skip the residual branch if it falls below this threshold. Since computing
the ratio requires running the residual branch, this is only used for
profiling and simulation purposes. Each data point represents an inference of
a single token using a ratio threshold of 5% is used. The figure confirms a
spread within the network depth, where most tokens only need between 40 and 70
blocks, instead of the full network at 80 blocks.

###  4.4 Routing Traces

For more detailed analysis, Figure 7 shows the routing for the OPT-13B model
across a batch of WikiText-2 examples. It reveals how the lower residual
ratios in Figure 9 lead to a significant number of skipped layers in the
beginning of the model. This again motivates the use of dynamic layer sparsity
over early-exit models, since early exit can only skip later layers, which
contribute the most to the network.

###  4.5 Vision

Figure 8: Vision Profiling – Current vision transformers are substantially
smaller than language models yet demonstrate similar trends with residual
ratio and dynamic depth. ViT-Huge (632M) shows median ratios comparable to
OPT-350M, and the dynamic depth assumes skipping layers with lower than 10%
residual ratio.

This analysis so far has focused on large language models, since they are
currently $10$ to $100\times$ larger than large vision models. Yet, recent
vision transformers have been proposed with tens of billions of parameters
[2]. These weights are not yet released, yet the trends between the smaller
language and vision models can still be aligned at smaller scales to suggest
the behavior of large vision models with billions or trillions of parameters.

Figure 8 shows a comparison for the largest released ViT model, which contains
632M parameters across 24 layers. It shows that vision transformers at this
size have comparable residual ratios to the similarly sized OPT-350M. In
addition, Appendix A lists smaller ViT versions and shows a similar trend
between model size and residual ratio, suggesting that as vision transformers
increase in size they will benefit from the same layer sparsity opportunities
as the OPT models.

##  5 Radial Networks

Given the high degree of dynamic layer sparsity, we propose a new neural
architecture that natively supports arbitrary routing between layers. As shown
in Fig 1, each token enters network at step $t_{0}$ and then is routed
dynamically to the next layer at later steps $t_{n}$. This process allows for
dynamic computation given the variable number of layers included in each token
path. As with standard transformers, the network continues auto-regressively
until an end-of-sequence token is produced.

###  5.1 Router

The router is the central component with the responsibility of directing each
embedding between layers at each time step. Building on the success of
mixture-of-expert models, we learn this router with a small multi-layer
perceptron (MLP) model. Beginning from layer $l_{t}$, it maps from
intermediate embeddings $e_{t}$ to output router logits $z_{t}$ for each
layer. These logits are then passed into a softmax function to produce
probabilities $p_{t}$. The maximum probability layer is then chosen as the
next layer in the forward pass.

$\displaystyle z_{t}$ $\displaystyle=Router(e_{t})$ $\displaystyle p_{t}$
$\displaystyle=\text{softmax}(z_{t})$ $\displaystyle l_{t+1}$
$\displaystyle=\text{argmax}(p_{t})$

The iterations stop when the model chooses the output layer, or a set maximum
number of layers are seen, which forces the output layer. This maximum number
of layers is a hyper-parameter that limits the worst-case dynamic depth in the
network.

###  5.2 Unified Cache

In standard sequential transformers, each layer activates for each token, and
the attention mechanism references the key-value cache of previous tokens.
However, given the dynamic routing of radial networks, many of the layers are
not activated for different tokens, and the key-value cache is sparse. To
solve this, we instead use shift the cache from the layers to a shared global
cache that stores all of the previous key-value pairs. Each embedding then
attends to the cached pairs of previous iterations and previous tokens.

To distinguish between the cached values for the current token and previous,
we use standard positional embeddings. These embeddings are used to encode the
give the relative position of tokens within the input sequence. A common
approach is to use sinusoidal functions; for position $p$ and dimension $i$,
the positional embedding $PE(p,i)$ can be defined as follows:

$PE(p,i)=\begin{cases}\sin(p/10000^{i/d})&\text{if }i\text{ is even}\\\
\cos(p/10000^{(i-1)/d})&\text{if }i\text{ is odd}\end{cases}$

where $d$ is the dimensionality of the embeddings.

##  6 Conclusion

In the past, dynamic layer sparsity has not been practical due to small model
sizes and incompatible neural architectures, which caused large contributions
from each layer and varying internal dimensions. For these reasons, it has
only been possible with techniques like early-exit, which requires expensive
specialized training. Yet, as language models grow in size, each layer
contributes less to output, creating opportunities for dynamic layer sparsity.
Following the trends in Figure 6, modern language models with over one
trillion parameters likely have median residual ratios less than 1%. And in
the future, as vision and multi-modal models catch up to language models,
their residual ratios should follow similar scaling trends.

To take advantage of this dynamic sparsity, we propose Radial Networks, which
perform token-level routing between layers guided by a trained router module.
These networks can be used in a post-training distillation from sequential
networks or trained from scratch to co-learn the router and layer weights.
They enable scaling to larger model sizes by decoupling the number of layers
from the dynamic depth of the network, and their design allows for layer
reuse. By varying the compute token by token, they reduce the overall
resources needed for generating entire sequences. Overall, this leads to
larger capacity networks with significantly lower compute and serving costs.

## References

  * Bolukbasi et al. [2017] Tolga Bolukbasi, Joseph Wang, Ofer Dekel, and Venkatesh Saligrama.  Adaptive neural networks for efficient inference.  _ICML_ , 2017.
  * Dehghani et al. [2023] Mostafa Dehghani, Josip Djolonga, Basil Mustafa, Piotr Padlewski, Jonathan Heek, Justin Gilmer, Andreas Steiner, Mathilde Caron, and Robert Geirhos.  Scaling vision transformers to 22 billion parameters.  _ICML_ , 2023.
  * Dosovitskiy et al. [2021] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  _ICLR_ , 2021.
  * Hao et al. [2021] Cong Hao, Jordan Dotzel, Jinjun Xiong, Luca Benini, Zhiru Zhang, and Deming Chen.  Enabling design methodologies and future trends for edge ai: Specialization and codesign.  _IEEE Design & Test_, 2021.
  * Hua et al. [2019] Weizhe Hua, Yuan Zhou, Christopher De Sa, Zhiru Zhang, and G. Edward Suh.  Channel gating neural networks.  _NeurIPS_ , 2019.
  * Liu et al. [2023] Zichang Liu, Jue Wang, Tri Dao, Tianyi Zhou, Binhang Yuan, Zhao Song, Anshumali Shrivastava, Ce Zhang, Yuandong Tian, Christopher Re, and Beidi Chen.  Deja vu: Contextual sparsity for efficient llms at inference time.  _ICML_ , 2023.
  * Teerapittayanon et al. [2016] Surat Teerapittayanon, Bradley McDanel, and H.T. Kung.  Branchynet: Fast inference via early exiting from deep neural networks.  _Intl. Conf. on Pattern Recog._ , 2016.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _NeurIPS_ , 2017.
  * Zhang et al. [2022] Susan Zhang, Stephen Roller, Naman Goyal, Mikel Artetxe, Moya Chen, Shuohui Chen, Christopher Dewan, Mona Diab, Xian Li, Xi Victoria Lin, Todor Mihaylov, Myle Ott, Sam Shleifer, Kurt Shuster, Daniel Simig, Punit Singh Koura, Anjali Sridhar, Tianlu Wang, and Luke Zettlemoyer.  Opt: Open pre-trained transformer language models, 2022.
  * Zhang et al. [2020] Yichi Zhang, Ritchie Zhao, Weizhe Hua, Nayun Xu, G. Edward Suh, and Zhiru Zhang.  Precision Gating: Improving Neural Network Efficiency with Dynamic Dual-Precision Activations.  _ICLR_ , 2020.

##  Appendix A Appendix

Figure 9: OPT Layer Residual Ratio – In general, the residual ratios decrease
with larger model sizes, implying that more layers can be more easily skipped.
Across models, there is a U-shaped distribution where the first and last
layers have the highest ratios and contribute the most to the model output.
Figure 10: OPT Residual Ratios – As models grow larger, the residual ratios
shift lower indicating more dynamic layer sparsity within the model. Each
value represents the residual ratio from a single layer and single token.
Evaluation is done on batches of sequences with length 256 sampled from the
WikiText-2 dataset. The dashed line represents the median value in the
distribution. Figure 11: ViT Family – Current vision transformers are
substantially smaller than language models yet demonstrate similar trends with
residual ratio. As model size increases, the contributions of each layer
decreases.


# Stable gravastar with large surface redshift in Einstein’s gravity with two
scalar fields

Shin’ichi Nojiri *  —<EMAIL_ADDRESS>Department of
Physics, Nagoya University, Nagoya 464-8602, Japan
&
Kobayashi-Maskawa Institute for the Origin of Particles and the Universe,
Nagoya University, Nagoya 464-8602, Japan  G. G. L. Nashed *  —
<EMAIL_ADDRESS>Centre for Theoretical Physics, The British University in
Egypt, P.O. Box 43, El Sherouk City, Cairo 11837, Egypt

(August 28, 2024)

###### Abstract

We propose a class of models, in which stable gravastar with large surface
redshift becomes a solution. In recent decades, gravastars have become a
plausible substitute for black holes. Researchers have explored stable
gravastar models in various alternative gravity theories, in addition to the
conventional framework of general relativity. In this paper, we present a
stellar model within the framework of Einstein’s gravity with two scalar
fields, in accordance with the conjecture proposed by Mazur and Mottola [Proc.
Nat. Acad. Sci. 101 (2004), 9545-9550]. In the model, the two scalar fields
become non-dynamical by imposing constraints in order to avoid ghosts. The
gravastar comprises two distinct regions, namely: (a) the interior region and
(b) the exterior region. We assume the interior region consists of the de
Sitter spacetime, and the exterior region is the Schwarzschild one. The two
regions are connected with each other by the shell region. On the shell, we
assume that the metric is given by a polynomial function of the radial
coordinate $r$. The function has six constants. These constants are fixed by
the smooth junction conditions, i.e., the interior region with the interior
layer of the shell and the exterior region with the exterior layer of the
shell. From these boundary conditions, we are able to write the coefficients
of the scalar fields in terms of the interior radius and exterior radius. To
clarify the philosophy of this study, we also give two examples of spacetimes
that asymptote as the de Sitter spacetime for small $r$ and as the
Schwarzschild spacetime for large $r$. Exploration is focused on the physical
attribute of the shell region, specifically, its proper length. The gravastar
model’s stability has frequently been examined by analyzing the relationship
between surface redshift and shell thickness, a comparison we also undertake
with previous models. Especially, we show that there exists a stable gravastar
with a large surface redshift prohibited by the instability in the previous
works. Furthermore, by checking the effective equation of state parameters, we
show that the gravastar geometry realized in this paper by using two scalar
fields could be difficult to generate with ordinary fluid.

Einstein general relativity with two scalar fields; gravastar; Junction
conditions, stability.

##  I Introduction

The final stage in the evolutionary journey of a massive celestial object
consistently captures the interest of astrophysicists. A conjecture proposes
that during the final phase of stellar evolution when the gravitational mass
of the remaining celestial body, following the creation of a planetary nebula
or the occurrence of a supernova explosion that surpasses three solar masses,
the gravitational collapse of the star leads to the creation of an
exceptionally dense entity, specifically, a black hole. Schwarzschild derived
the most straightforward black hole solution from Einstein’s vacuum field
equation in 1916, which is regarded as a significant solution for a black hole
that is static and uncharged. However, the solution of Schwarzschild exhibits
certain limitations attributed to (I) the existence of a singularity within
the confines of a black hole, along with (II) The presence of an event horizon
presents substantial, unresolved challenges. Mazur and Mottola Mitra (2002);
Mazur and Mottola (2004) were the first to introduce a novel concept for the
collapse of stellar objects, incorporating an expanded notion of the Bose-
Einstein Condensation within the gravitational system. This model proposed the
potential creation of a cold, faint, and compact entity, which they named a
gravitational vacuum condensed star or gravastar. The representation of the
gravastar must meet all the theoretical prerequisites to serve as a stable
concluding phase in stellar evolution and offer resolutions to the problems
linked to conventional black holes. Following the triumphant detection of
gravitational waves (GWs) in 2015 Abbott et al. (2016), there has been a
presumption that these GWs emanate from the collision of two substantial black
holes. Given that the observed signs do not conclusively address the
fundamental problems linked to black holes, an alternative method is necessary
to tackle these limitations. In such a framework, the gravastar might play a
pivotal role in elucidating the eventual stage of stellar evolution. Although
significant observational proof endorsing the existence of gravastars is
absent, it remains crucial to explore the notion of gravastars as a viable
alternative for addressing the conceptual dilemmas linked to understanding
black holes.

The gravastar is composed of three separate regions: the inner region, a thin
intermediate shell, and an outer region. After the suggestion made by Mazur
and Mottola Mitra (2002); Mazur and Mottola (2004), it is postulated that the
gravastar’s inner region exists in a phase resembling a de Sitter condensate
and is believed to be filled with vacuum energy. Concurrently, the outer
region is in a state of complete vacuum known as Schwarzschild vacuum, and
these two regions are divided by an ultra-relativistic matter-thin shell
characterized by an exceptionally high density. The barotropic equation of
state (EoS) is defined as $p=\rho$, with $\rho$ representing the matter-energy
density and $p$ standing for pressure. Within a gravastar, the equation of
state (EoS) parameter $\omega\equiv\frac{p}{\rho}$ varies across various
regions as depicted below:
a) In the interior region $(0\leq r<r_{1})$, the EoS parameter is $\omega=-1$,
which corresponds to $p=-\rho$, where $\rho$ is density and $p$ represents
pressure.
b) In the thin shell $(r_{1}<r<r_{2})$, the EoS parameter is $\omega=1$, which
corresponds to $p=\rho$.
c) In the exterior region $r_{2}<r$ and $\omega=0$ which corresponds to $p=0$.
In this context, $r_{1}$ and $r_{2}$ represent the radii of the inner and
outer boundaries of the gravastar, respectively. Consequently, the thickness
of the intermediate shell region can be expressed as $\delta\equiv
r_{2}-r_{1}$, where $\delta$ could be much smaller than the mass parameter of
the exterior Schwarzschild spacetime and the length parameter of the interior
de Sitter spacetime.

The concept of gravastars as an alternative to black holes has been
extensively explored in the literature Visser and Wiltshire (2004); Cattoen et
al. (2005); Carter (2005); Bilic et al. (2006); Lobo (2006); DeBenedictis et
al. (2006); Lobo and Arellano (2007); Horvat and Ilijic (2007); Chirenti and
Rezzolla (2007); Rocha et al. (2008); Horvat et al. (2009); Usmani et al.
(2008); Turimov et al. (2009); Nandi et al. (2009); Harko et al. (2009);
Usmani et al. (2011); Rahaman et al. (2012a, b); Bhar (2014); Rahaman et al.
(2015); Ghosh et al. (2017, 2018, 2019a, 2019b); Chan and da Silva (2010).
Mazur and Mottola proposed a gravastar model with a three-region structure and
examined its thermodynamic stability through entropy maximization methods.
Meanwhile, Visser Visser and Wiltshire (2004) illustrated the dynamic
stability of this system based on the three-region model, It is currently
well-established that the theoretical viability of the gravastar model,
whether in varying dimensions or with or without charge, holds significant
importance within the realm of astrophysics. Examining the Vaidya exterior
spacetime, Chan et al. Chan et al. (2011) investigated the dynamic models of
prototype radiating gravastars. Their research unveiled that the final out
could materialize as a range of entities, such as a black hole, an unsteady
gravastar, a stable gravastar, or a “bounded excursion” gravastar. The
particular result depends on variables like the evolving shell mass, the
cosmological constant, and the initial placement of the moving shell. Such
inquiries are primarily investigated within the context of Albert Einstein’s
general relativity (GR). Although it is widely acknowledged that GR remains
one of the most encouraging theories for uncovering many of nature’s enigmas,
it has become clear that an alteration of GR is necessary to overcome its
shortcomings in both the theoretical and observational realms. The
confirmation of the universe’s acceleration, along with the existence of dark
matter and dark energy, as indicated by specific observational discoveries,
presents a theoretical challenge to GR Riess et al. (1998); Perlmutter et al.
(1999); de Bernardis et al. (2000); Hanany et al. (2000); Peebles and Ratra
(2003); Padmanabhan (2003); Clifton et al. (2012); Riess et al. (2007);
Tegmark et al. (2004); Amanullah et al. (2010); Komatsu et al. (2011). In
order to account for the current acceleration of the universe, several
successive alternative gravity theories have emerged. These theories include
$f(T)$ gravity, $f(R)$ gravity, $f(R,\mathcal{T})$ gravity, and
$f(T,\mathcal{T})$ gravity, where $\mathcal{T}$ is the trace of the energy-
momentum tensor, each of which incorporates modifications either in the
geometric aspect or in the matter-energy component of the Einstein field
equations. In the majority of such gravitational theories, various functions
of $R$, the Ricci scalar, either alone or in conjunction with other scalar
quantities, are integrated into the gravitational Lagrangian within the
corresponding action. Numerous studies have been conducted on compact stars
and gravastars within the context of amended theories, as documented in Das et
al. (2015, 2016, 2017); Biswas et al. (2019); Ghosh et al. (2020); Das et al.
(2020); Sengupta et al. (2020); Banerjee et al. (2020). However, in the
majority of such gravitational theories, the energy-momentum tensor
$T_{\mu\nu}$ does not exhibit general conservation properties, in contrast to
the scenario in GR. This intriguing facet of modified gravity serves as our
motivation to conduct a comprehensive examination of gravastars within the
frame of the GR theory with two scalars, formulated in Nojiri et al. (2021).
In the model, the scalar fields often become ghosts, which have kinetic
energies unbounded below as the classical theory and generate negative norm
states in the quantum theory. In order to avoid the ghosts, we impose
constraints to make the scalar fields become non-dynamical, that is, so that
the scalar fields do not propagate. The constraints are similar to the mimetic
one proposed in Chamseddine et al. (2014), where non-dynamical dark matter is
generated. The constraints make the solution in the present model stable. Due
to this stability, we can realize the stable gravastar with large surface
redshift which was prohibited by the stability in the previous works Buchdahl
(1959); Straumann (1984); Böhmer and Harko (2006); Ivanov (2002); Barraco and
Hamity (2002); Böhmer and Harko (2006).

The structure of this current study unfolds as follows: In Section II, we
elucidate the fundamental mathematical framework of the Einstein GR theory
with two scalar fields and provide the specific expressions for the field
equations. We impose the constraints for the scalar fields so that the scalar
fields become non-dynamical. We show that any spherically symmetric solution
in the present model is always stable under any fluctuation because the two
scalar fields become non-dynamical. In Section III, we derive the model, which
generates solutions satisfying the boundary conditions between the shell
region and the interior and exterior regions of the gravastar. In Section IV,
we discuss the physics of the gravastar by driving the proper length and the
surface redshift of the shell. Here we show that the model in this paper
realizes the stable gravastar with large surface redshift prohibited due to
the instability in the previous works. In Section VI, we wrap up our study
with a summary and conclusion.

##  II Einstein gravity coupled with two scalar fields

Einstein’s GR with two scalar fields $\phi$ and $\chi$ is described by the
action as follows Nojiri et al. (2021),

$\displaystyle S_{\mathrm{GR}\phi\chi}=\int d^{4}x\sqrt{-g}$
$\displaystyle\left[\frac{R}{2\kappa^{2}}-\frac{1}{2}\,A(\phi,\chi)\partial_{\mu}\phi\partial^{\mu}\phi-B(\phi,\chi)\,\partial_{\mu}\phi\partial^{\mu}\chi\right.$
$\displaystyle\left.\quad-\frac{1}{2}\,C(\phi,\chi)\partial_{\mu}\chi\partial^{\mu}\chi-V(\phi,\chi)\right]\,.$
(1)

In this context, $g$ represents the determinant of the metric tensor
$g_{\mu\nu}$, $R$ denotes the Ricci scalar, and $V(\phi,\chi)$ represents the
potential of the scalar doublet. The functional forms of the coefficient
functions $A$, $B$, and $C$ are contingent upon the properties of the scalars.
Upon varying the action (II) with respect to the metric $g_{\mu\nu}$, we
derive the ensuing Einstein equation,

$\displaystyle\frac{1}{\kappa^{2}}\left(R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R\right)=$
$\displaystyle\,A(\phi,\chi)\partial_{\mu}\phi\partial_{\nu}\phi+B(\phi,\chi)\left(\partial_{\mu}\phi\partial_{\nu}\chi+\partial_{\nu}\phi\partial_{\mu}\chi\right)+C(\phi,\chi)\partial_{\mu}\chi\partial_{\nu}\chi$
$\displaystyle\,-g_{\mu\nu}\left[\frac{1}{2}\,A(\phi,\chi)\partial_{\rho}\phi\partial^{\rho}\phi+B(\phi,\chi)\partial_{\rho}\phi\partial^{\rho}\chi+\frac{1}{2}\,C(\phi,\chi)\partial_{\rho}\chi\partial^{\rho}\chi+V(\phi,\chi)\right]\,.$
(2)

Through the variation of action (II) concerning the scalar fields $\phi$ and
$\chi$, we acquire the following expressions,

$\displaystyle 0=$
$\displaystyle\,\frac{A_{\phi}}{2}\,\partial_{\mu}\phi\partial^{\mu}\phi+A\nabla^{\mu}\partial_{\mu}\phi+A_{\chi}\partial_{\mu}\phi\partial^{\mu}\chi+\left(B_{\chi}-\frac{1}{2}\,C_{\phi}\right)\partial_{\mu}\chi\partial^{\mu}\chi+B\nabla^{\mu}\partial_{\mu}\chi-
V_{\phi}\,,$ (3) $\displaystyle 0=$
$\displaystyle\,\left(-\frac{1}{2}\,A_{\chi}+B_{\phi}\right)\partial_{\mu}\phi\partial^{\mu}\phi+B\nabla^{\mu}\partial_{\mu}\phi+\frac{1}{2}\,C_{\chi}\partial_{\mu}\chi\partial^{\mu}\chi+C\nabla^{\mu}\partial_{\mu}\chi+C_{\phi}\partial_{\mu}\phi\partial^{\mu}\chi-
V_{\chi}\,.$ (4)

Here, let us define $A_{\phi}$ as $\partial A(\phi,\chi)/\partial\phi$, and
similarly for other derivatives. It is worth noting that Eqs. (3) and (4) can
be derived from the Bianchi identity in conjunction with Eq. (II).

In the following, we identify

$\displaystyle\phi=t\,,\quad\chi=r\,.$ (5)

As explained in the reference Nojiri et al. (2021), making the assumption (5)
does not result in any loss of generality. In the case of a spacetime with a
general spherically symmetric yet time-dependent solution, the scalar fields
$\phi$ and $\chi$ exhibit dependencies on both the time coordinate, denoted as
$t$ and the radial coordinate, denoted as $r$. In the context of a given
solution, the specific dependencies of $\phi$ and $\chi$ on both the time
variable $t$ and the radial variable $r$ are determined as functions:
$\phi=\phi(t,r)$ and $\chi=\chi(t,r)$. We may redefine the scalar fields to
replace $t$ and $r$ with new scalar fields, $\tilde{\phi}$ and $\tilde{\chi}$,
$\phi\left(\tilde{\phi},\tilde{\chi}\right)\equiv\phi\left(t=\tilde{\phi},r=\tilde{\chi}\right)$
and
$\chi\left(\tilde{\phi},\tilde{\chi}\right)\equiv\chi\left(t=\tilde{\phi},r=\tilde{\chi}\right)$.
Subsequently, we can associate the new scalar fields with the time and radial
coordinates in (5). The transformation of variables from
$\left(\phi,\chi\right)$ to $\left(\tilde{\phi},\tilde{\chi}\right)$ can be
absorbed into the redefinitions of $A$, $B$, $C$, and $V$ within the action
(II). This demonstrates that making the assumption (5) does not lead to any
loss of generality. Furthermore, as we will observe, $\phi$ is identified with
$t$ even for a static spacetime.

If we assume the static and spherically symmetric spacetime, whose metric is
given by

$\displaystyle
ds^{2}=-\mathrm{e}^{2\nu(r)}dt^{2}+\mathrm{e}^{2\mu(r)}dr^{2}+r^{2}\left(d\theta^{2}+\sin^{2}\theta\,d\phi^{2}\right)\,,$
(6)

the $\left(t,t\right)$, $\left(r,r\right)$, $\left(i,j\right)$, and
$\left(t,r\right)$ components of the Einstein equation under the
identification (5) have the following forms,

$\displaystyle\frac{\mathrm{e}^{2\left(\nu-\mu\right)}}{\kappa^{2}}\left(\frac{2\mu^{\prime}}{r}+\frac{\mathrm{e}^{2\mu}-1}{r^{2}}\right)=$
$\displaystyle\,-\mathrm{e}^{2\nu}\left(-\frac{A}{2}\mathrm{e}^{-2\nu}-\frac{C}{2}\mathrm{e}^{-2\mu}-V\right)\,,$
(7)
$\displaystyle\frac{1}{\kappa^{2}}\left(\frac{2\nu^{\prime}}{r}-\frac{\mathrm{e}^{2\mu}-1}{r^{2}}\right)=$
$\displaystyle\,\mathrm{e}^{2\mu}\left(\frac{A}{2}\,\mathrm{e}^{-2\nu}+\frac{C}{2}\,\mathrm{e}^{-2\mu}-V\right)\,,$
(8)
$\displaystyle\frac{1}{\kappa^{2}}\left[\mathrm{e}^{-2\mu}\left(r\left(\nu^{\prime}-\mu^{\prime}\right)+r^{2}\nu^{\prime\prime}+r^{2}\left(\nu^{\prime}-\mu^{\prime}\right)\nu^{\prime}\right)\right]=$
$\displaystyle\,r^{2}\left(\frac{A}{2}\mathrm{e}^{-2\nu}-\frac{C}{2}\mathrm{e}^{-2\mu}-V\right)\,,$
(9) $\displaystyle 0=$ $\displaystyle\,B\,,$ (10)

which can be solved with respect to $A$, $B$, $C$, and $V$, as follows,

$\displaystyle A(r)=$
$\displaystyle\,\frac{\mathrm{e}^{2\left(\nu-\mu\right)}}{\kappa^{2}}\left[\frac{\mathrm{e}^{2\mu}-1}{r^{2}}+\frac{\nu^{\prime}+\mu^{\prime}}{r}+\nu^{\prime\prime}+\left(\nu^{\prime}-\mu^{\prime}\right)\nu^{\prime}\right]\,,$
(11) $\displaystyle B(r)=$ $\displaystyle\,0\,,$ (12) $\displaystyle C(r)=$
$\displaystyle\,\frac{1}{\kappa^{2}}\left[-\frac{\mathrm{e}^{2\mu}-1}{r^{2}}+\frac{\nu^{\prime}+\mu^{\prime}}{r}-\nu^{\prime\prime}-\left(\nu^{\prime}-\mu^{\prime}\right)\nu^{\prime}\right]\,,$
(13) $\displaystyle V(r)=$
$\displaystyle\,\frac{\mathrm{e}^{-2\mu}}{2\kappa^{2}}\left[\frac{2\left(\mu^{\prime}-\nu^{\prime}\right)}{r}+\frac{2\left(\mathrm{e}^{2\mu}-1\right)}{r^{2}}\right]\,.$
(14)

This tells that if we consider the model where the radial coordinate $r$ in
$A$, $B$, $C$, and $V$ of Eqs. (11), (12), (13), and (14), is replaced with
$\chi$, that is, $A(r=\chi)$, $B(r=\chi)$, $C(r=\chi)$, and $V(r=\chi)$, the
spacetime given by Eq. (6) becomes a solution of the model.

We should note the product of $A$ and $C$ is always negative, which tells us
that either $\phi$ or $\chi$ is a ghost. The ghost mode has negative kinetic
energy classically and generates negative norm states as a quantum theory.
Therefore the existence of the ghost mode tells that the model is physically
inconsistent. We now eliminate ghosts by imposing constraints on $\phi$ and
$\chi$, which is similar to the mimetic constraint in Chamseddine and Mukhanov
(2013), and make the ghost modes non-dynamical. For the purpose, by using the
Lagrange multiplier fields $\lambda_{\phi}$ and $\lambda_{\chi}$, we add the
following terms to the action $S_{\mathrm{GR}\phi\chi}$ in (II) as
$S_{\mathrm{GR}\phi\chi}\to S_{\mathrm{GR}\phi\chi}+S_{\lambda}$,

$\displaystyle S_{\lambda}=\int
d^{4}x\sqrt{-g}\left[\lambda_{\phi}\left(\mathrm{e}^{-2\nu(r=\chi)}\partial_{\mu}\phi\partial^{\mu}\phi+1\right)+\lambda_{\chi}\left(\mathrm{e}^{-2\mu(r=\chi)}\partial_{\mu}\chi\partial^{\mu}\chi-1\right)\right]\,.$
(15)

The variations of $S_{\lambda}$ with respect to $\lambda_{\phi}$ and
$\lambda_{\chi}$ give the constraints111 The constraint given by Eq. (16), is
similar to the constraint of mimetic theory; see, for example, Chamseddine et
al. (2014); Nojiri and Nashed (2022); Myrzakulov et al. (2016); Nashed and
Nojiri (2022); Nashed (2023a, b); Myrzakulov et al. (2015); Vagnozzi (2017);
Nashed and Nojiri (2021); Casalino et al. (2018, 2019); Sebastiani et al.
(2017); Nashed (2021).

$\displaystyle
0=\mathrm{e}^{-2\nu(r=\chi)}\partial_{\mu}\phi\partial^{\mu}\phi+1\,,\quad
0=\mathrm{e}^{-2\mu(r=\chi)}\partial_{\mu}\chi\partial^{\mu}\chi-1\,,$ (16)

whose solutions are consistently given by (5). We should note that even in the
model modified by adding $S_{\lambda}$, $\lambda_{\phi}=\lambda_{\chi}=0$ is
always a solution. Therefore the spacetime given by Eq. (6) is a solution even
in the modified theory. The constraints in (16) make the scalar field $\phi$
and $\chi$ non-dynamical, that is, the fluctuation of $\phi$ and $\chi$ from
the background (5) do not propagate. In fact, by considering the perturbation
from (5),

$\displaystyle\phi=t+\delta\phi\,,\quad\chi=r+\delta\chi\,,$ (17)

we obtain,

$\displaystyle\partial_{t}\delta\phi=\partial_{r}\delta\chi=0\,.$ (18)

Therefore if we impose the initial condition $\delta\phi=0$, we find
$\delta\phi=0$ in the whole spacetime. On the other hand, if we impose the
boundary condition $\delta\chi\to 0$ when $r\to 0$, we find $\delta\chi=0$ in
the whole spacetime. These tell that $\phi$ and $\chi$ are non-dynamical.
Therefore the solution given by the spherically symmetric spacetime in (6) is
always stable under the perturbation of the scalar fields. Only the
propagating mode or possible fluctuation in this model is given by the
gravitational wave. Because the metric minimally couples with scalar fields,
the propagation of the gravitational wave does not change from the standard
one, and therefore the fluctuation does not cause any instability, and the
solution in (6) is always stable under any fluctuation. This situation is
similar to the original mimetic gravity theory in Chamseddine and Mukhanov
(2013), where effective dark matter appeared. The effective dark matter is not
dynamical and there does not appear any fluctuation. For example, the
effective dark matter never collapses due to gravity. Therefore the
formulation by using the constraints in (16) could be regarded as a natural
extension of the original mimetic gravity.

We often consider the spacetime where $\mu=-\nu$. For the spacetime, Eqs.
(11), (12), (13), and (14) are reduced as follows,

$\displaystyle A=$
$\displaystyle\,\frac{\mathrm{e}^{2\nu}}{\kappa^{2}}\left[-\frac{\mathrm{e}^{2\nu}-1}{r^{2}}+\frac{1}{2}\left(\mathrm{e}^{2\nu}\right)^{\prime\prime}\right]\,,\quad
B=0\,,$ $\displaystyle C=$
$\displaystyle\,\frac{\mathrm{e}^{-2\nu}}{\kappa^{2}}\left[\frac{\mathrm{e}^{2\nu}-1}{r^{2}}-\frac{1}{2}\left(\mathrm{e}^{2\nu}\right)^{\prime\prime}\right]\,,\quad
V=\frac{1}{\kappa^{2}}\left[-\frac{\left(\mathrm{e}^{2\nu}\right)^{\prime}}{r}-\frac{\mathrm{e}^{2\nu}-1}{r^{2}}\right]\,.$
(19)

We use the above expressions in the following.

##  III Realization of Gravastar geometry

We now construct the model which realizes gravastar spacetime. In the
gravastar geometry, there is an exterior region which is the Schwarzschild
spacetime where

$\displaystyle\mathrm{e}^{2\nu}=\mathrm{e}^{-2\mu}=1-\frac{2M}{r}\,,$ (20)

and the interior region is described by the de Sitter spacetime,

$\displaystyle\mathrm{e}^{2\nu}=\mathrm{e}^{-2\mu}=1-\lambda^{2}r^{2}\,.$ (21)

Here $M$ corresponds to the mass measured by the observer at infinity. On the
other hand, $\lambda^{2}$ corresponds to the cosmological constant. The
exterior region and the interior region are connected by the shell region. We
now consider the two cases. In one case,
$\mathrm{e}^{2\nu}=\mathrm{e}^{-2\mu}$ in the shell region is given by the
polynomial of $r$. In another case, we try to connect the two regions by using
smooth functions, which describe the asymptotically Schwarzschild spacetime
when $r$ is large but describe the asymptotically de Sitter spacetime when $r$
is small.

###  III.1 Shell model

We consider the shell region in $r_{1}<r<r_{2}$, with a boundary with the
interior region at $r=r_{1}$ and the boundary with the exterior region at
$r=r_{2}$. We assume that $\mathrm{e}^{2\nu}=\mathrm{e}^{-2\mu}$ in the shell
region is given by the following polynomial,

$\displaystyle\mathrm{e}^{2\nu}=\mathrm{e}^{-2\mu}=f(r)=a_{0}+a_{1}\left(r-r_{1}\right)+a_{2}\left(r-r_{1}\right)^{2}+a_{3}\left(r-r_{1}\right)^{3}+a_{4}\left(r-r_{1}\right)^{4}+a_{5}\left(r-r_{1}\right)^{5}\,.$
(22)

By imposing the continuities of $\mathrm{e}^{2\nu}$,
$\left(\mathrm{e}^{2\nu}\right)^{\prime}$, and
$\left(\mathrm{e}^{2\nu}\right)^{\prime\prime}$ at the two boundaries, we
obtain,

$\displaystyle a_{0}=$ $\displaystyle\,1-\lambda^{2}{r_{1}}^{2}>0\,,\quad
a_{1}=-2\lambda^{2}r_{1}<0\,,\quad a_{2}=-\lambda^{2}<0\,,$ $\displaystyle
a_{3}=$
$\displaystyle\,-\frac{1}{{r_{2}}^{3}\left(r_{2}-r_{1}\right)^{3}}\left\\{2\left(15{r_{2}}^{2}-6r_{2}r_{1}+{r_{1}}^{2}\right)M-{r_{2}}^{3}\left(3{r_{2}}^{2}+6r_{2}r_{1}+{r_{1}}^{2}\right)\lambda^{2}\right\\}\,,$
$\displaystyle a_{4}=$
$\displaystyle\,\frac{1}{{r_{2}}^{3}\left(r_{2}-r_{1}\right)^{4}}\left\\{2\left(24{r_{2}}^{2}-11r_{2}r_{1}+2{r_{1}}^{2}\right)M-{r_{2}}^{3}\left(3{r_{2}}^{3}+10r_{2}r_{1}+2{r_{1}}^{2}\right)\lambda^{2}\right\\}\,,$
$\displaystyle a_{5}=$
$\displaystyle\,-\frac{1}{{r_{2}}^{3}\left(r_{2}-r_{1}\right)^{5}}\left\\{2\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)M-{r_{2}}^{3}\left({r_{2}}^{2}+4r_{2}r_{1}+{r_{1}}^{2}\right)\lambda^{2}\right\\}\,.$
(23)

Because
$\frac{{r_{2}}^{2}+4r_{2}r_{1}+{r_{1}}^{2}}{10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}}>0$,
we may choose $M$ by

$\displaystyle
M=\frac{{r_{2}}^{3}\left({r_{2}}^{2}+4r_{2}r_{1}+{r_{1}}^{2}\right)\lambda^{2}}{2\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)}\,,$
(24)

so that $a_{5}=0$. In this choice, we find

$\displaystyle a_{0}=$ $\displaystyle\,1-\lambda^{2}{r_{1}}^{2}>0\,,\quad
a_{1}=-2\lambda^{2}r_{1}<0\,,\quad a_{2}=-\lambda^{2}<0\,,$ $\displaystyle
a_{3}=$
$\displaystyle\,\frac{\lambda^{2}r_{2}\left(15{r_{2}}^{2}+6r_{2}r_{1}-3{r_{1}}^{2}\right)}{\left(r_{2}-r_{1}\right)^{2}\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)}>0\,,$
$\displaystyle a_{4}=$
$\displaystyle\,-\frac{3\lambda^{2}r_{2}\left(2{r_{2}}^{2}+2r_{2}r_{1}-{r_{1}}^{2}\right)}{\left(r_{2}-r_{1}\right)^{3}\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)}<0\,.$

We also find

$\displaystyle f(r_{1})=$ $\displaystyle\,a_{0}>0\,,\quad
f^{\prime}(r_{1})=a_{1}<0\,,\quad f^{\prime\prime}(r_{1})=2a_{2}<0\,,\quad
f^{\prime\prime\prime}(r_{1})=6a_{3}>0\,.\quad
f^{\prime\prime\prime\prime}(r)=24a_{4}<0\,,$ $\displaystyle f(r_{2})=$
$\displaystyle\,1-\frac{2M}{r_{2}}>0\,,\quad
f^{\prime}(r_{2})=\frac{2M}{{r_{2}}^{2}}>0\,,\quad
f^{\prime\prime}(r_{2})=-\frac{4M}{{r_{2}}^{3}}<0\,,$ $\displaystyle
f^{\prime\prime\prime}(r_{2})=$
$\displaystyle\,6a_{3}+24a_{4}\left(r_{2}-r_{1}\right)=-\frac{54\lambda^{2}r_{2}\left({r_{2}}^{2}+2r_{2}r_{1}-{r_{1}}^{2}\right)}{\left(r_{2}-r_{1}\right)^{2}\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)}<0\,.$

This tells the function $f(r)$ has two maxima when $r<r_{1}$ and $r>r_{2}$ and
one minimum when $r_{1}<r<r_{2}$. The values of $r$ corresponding to the
extrema are given by solving the equations

$\displaystyle
0=f^{\prime}(r)=a_{1}+2a_{2}\left(r-r_{1}\right)+3a_{3}\left(r-r_{1}\right)^{2}+4a_{4}\left(r-r_{1}\right)^{3}\,.$
(27)

Because all $a_{i}$’s $\left(i=1,\cdots,4\right)$ are proportional to
$\lambda^{2}$, the solution of (27) does not depend on $\lambda$. Therefore if
we choose $\lambda$ to be sufficiently small, that is, $\lambda r_{1}$,
$\lambda r_{2}\ll 1$, only $a_{0}$ dominates $a_{0}\sim 1$ in the expression
of $f(r)$ in (22) and therefore $f(r)$ can be consistently positive. If $f(r)$
could become negative, the horizon would appear and therefore the solution
would not describe the gravastar but a black hole, which could be regular at
the origin. We should also note $f(r)<1$. On the other hand, even if $\lambda$
is not so small, that is, $\lambda r_{1}$, $\lambda r_{2}\sim 1$, when the
thickness of the shell is small enough, $\lambda\left(r_{2}-r_{1}\right)\ll
1$, $f(r)$ could be positive.

We numerically investigate typical two cases, where we choose,

  * •

Sufficiently small $\lambda$ case,

$\displaystyle r_{1}=1\,,\quad r_{2}=1.2\,,\quad\lambda=0.1\,,\quad\mbox{which
gives}\quad M=0.000665\cdots\,.$ (28)

The dimension of the mass or length is fixed by choosing $r_{1}=1$. When we
include the dimension, $\lambda r_{1,2}$ and $\lambda M$ are dimensionless.

  * •

Sufficiently thin shell case, where we choose,

$\displaystyle r_{1}=0.99\,,\quad r_{2}=1\,,\quad M=0.495\,,\quad\mbox{which
gives}\quad\lambda=1.002496757\cdots\,.$ (29)

Here the dimension of the mass or length is fixed by choosing $r_{2}=1$.

The behaviors of $f(r)$ in Eq. (22) are shown in Figure 1 0(a) for the case of
(28) and Figure 2 1(a) for the case of (29). In the case of (28) , we find
$M\ll r_{2}$ but in the case of (29), $r_{2}\sim 2M$. Note that $r_{2}$
corresponds to the radius of the surface in the gravastar. Therefore in the
case of (29), the radius of the gravastar is smaller than the radius of the
photon sphere $\sim 6M$ and the photon sphere could be observed in the case of
(29), which could be consistent with the observation of the Event Horizon
Telescope EventHorizonTelescope:2019ths . As we will show later, in the case
of (29), the surface redshift becomes very large because $r_{2}\sim 2M$ and
the redshift lies inside the region prohibited by the stability in the
previous works.

In the interior region $\left(r<r_{1}\right)$, we find

$\displaystyle A=B=C=0\,,\quad V=\frac{3\lambda^{2}}{\kappa^{2}}\,,$ (30)

which is consistent with the fact that the interior region is the de Sitter
spacetime. On the other hand, in the exterior region $\left(r>r_{2}\right)$,
we obtain,

$\displaystyle A=B=C=V=0\,,$ (31)

which corresponds to the Schwarzschild spacetime which is a vacuum solution of
the Einstein equation. In the shell region $\left(r_{1}<r<r_{2}\right)$, $A$,
$B$, $C$, and $V$ are given by substituting (22) with (III.1) or (III.1) into
the expressions in (II),

$\displaystyle A=$
$\displaystyle\,3\lambda^{2}r_{2}\left\\{1+\lambda^{2}\left(-{r_{1}}^{2}-2r_{1}\left(r-r_{1}\right)-\left(r-r_{1}\right)^{2}+\frac{3\left(5{r_{2}}^{2}+2r_{2}r_{1}-{r_{1}}^{2}\right)r_{2}\left(r-r_{1}\right)^{3}}{\left(r_{2}-r_{1}\right)^{2}\left(10{r_{2}}^{2}-5r_{2}\,r_{1}+{r_{1}}^{2}\right)}\right.\right.$
$\displaystyle\,\left.\left.\qquad-\frac{3\left(2{r_{2}}^{2}+2r_{2}r_{1}-{r_{1}}^{2}\right)r_{2}\left(r-r_{1}\right)^{4}}{\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)\left(r_{2}-r_{1}\right)^{3}}\right)\right\\}\left\\{5{r_{2}}^{3}{r_{1}}^{3}-{r_{2}}^{2}{r_{1}}^{4}-r_{2}{r_{1}}^{5}\right.$
$\displaystyle\,\qquad+\left({r_{2}}^{2}{r_{1}}^{3}-15{r_{2}}^{3}{r_{1}}^{2}+r_{2}\,{r_{1}}^{4}+{r_{1}}^{5}\right)r+\left(10{r_{2}}^{3}+10{r_{2}}^{2}r_{1}+10r_{2}{r_{1}}^{2}-6{r_{1}}^{3}\right)r^{3}$
$\displaystyle\,\left.\qquad+\left(5{r_{1}}^{2}-10{r_{2}}^{2}-10r_{2}r_{1}\right)r^{4}\right\\}\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)^{-1}\left(r_{2}-r_{1}\right)^{-3}r^{-2}\,,$
$\displaystyle C=$
$\displaystyle\,-3\lambda^{2}r_{2}\left\\{5{r_{2}}^{3}{r_{1}}^{3}-{r_{2}}^{2}{r_{1}}^{4}-r_{2}{r_{1}}^{5}+\left({r_{1}}^{5}-15{r_{2}}^{3}{r_{1}}^{2}+{r_{2}}^{2}{r_{1}}^{3}+r_{2}{r_{1}}^{4}\right)r\right.$
$\displaystyle\,\left.\qquad+\left(10{r_{2}}^{3}+10{r_{2}}^{2}r_{1}+10r_{2}{r_{1}}^{2}-6{r_{1}}^{3}\right)r^{3}+\left(-10{r_{2}}^{2}+5{r_{1}}^{2}-10r_{2}r_{1}\right)r^{4}\right\\}$
$\displaystyle\,\times\left\\{1+\lambda^{2}\left(-{r_{1}}^{2}-2r_{1}\left(r-r_{1}\right)-\left(r-r_{1}\right)^{2}+\frac{3\left(5{r_{2}}^{2}+2r_{2}r_{1}-{r_{1}}^{2}\right)r_{2}\left(r-r_{1}\right)^{3}}{\left(r_{2}-r_{1}\right)^{2}\left(10{r_{2}}^{2}-5r_{2}\,r_{1}+{r_{1}}^{2}\right)}\right.\right.$
$\displaystyle\,\left.\left.\qquad-\frac{3\left(2{r_{2}}^{2}+2r_{2}r_{1}-{r_{1}}^{2}\right)r_{2}\left(r-r_{1}\right)^{4}}{\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)\left(r_{2}-r_{1}\right)^{3}}\right)\right\\}^{-1}\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)^{-1}\left(r_{2}-r_{1}\right)^{-3}r^{-2}\,,$
$\displaystyle V=$
$\displaystyle\,3\lambda^{2}r_{2}\left\\{10{r_{2}}^{3}{r_{1}}^{3}-2{r_{2}}^{2}{r_{1}}^{4}-2r_{2}{r_{1}}^{5}+\left(-15{r_{2}}^{3}{r_{1}}^{2}+{r_{2}}^{2}{r_{1}}^{3}+r_{2}{r_{1}}^{4}+{r_{1}}^{5}\right)r\right.$
$\displaystyle\left.\qquad+\left(5{r_{2}}^{3}+5{r_{2}}^{2}r_{1}+5r_{2}{r_{1}}^{2}-3{r_{1}}^{3}\right)r^{3}+\left(-4{r_{2}}^{2}-4r_{2}r_{1}+2{r_{1}}^{2}\right)r^{4}\right\\}$
$\displaystyle\,\qquad\times\left\\{2\left(10\,{r_{2}}^{2}-5\,r_{2}\,r_{1}+{r_{1}}^{2}\right)\left(r_{2}-r_{1}\right)^{3}r^{2}\right\\}^{-1}\,.$
(32)

The above expressions are consistent with the interior and exterior solutions
in (30) and (31) because

$\displaystyle A\left(r_{1}\right)=$
$\displaystyle\,\frac{\mathrm{e}^{2\nu(r_{1})}}{\kappa^{2}}\left[-\frac{f\left(r_{1}\right)-1}{{r_{1}}^{2}}+\frac{1}{2}f^{\prime\prime}\left(r_{1}\right)\right]B\left(r_{1}\right)=0\,,$
$\displaystyle C\left(r_{1}\right)=$
$\displaystyle\,\frac{\mathrm{e}^{-2\nu\left(r_{1}\right)}}{\kappa^{2}}\left[\frac{f\left(r_{1}\right)-1}{r^{2}}-\frac{1}{2}f^{\prime\prime}\left(r_{1}\right)\right]=\frac{\mathrm{e}^{-2\nu\left(r_{1}\right)}}{\kappa^{2}}\left[\frac{1-\lambda^{2}{r_{1}}^{2}-1}{{r_{1}}^{2}}-\frac{1}{2}\left(-2\lambda^{2}\right)\right]=0\,,$
$\displaystyle V\left(r_{1}\right)=$
$\displaystyle\,\frac{1}{\kappa^{2}}\left[-\frac{f^{\prime}\left(r_{1}\right)}{r_{1}}-\frac{f\left(r_{1}\right)-1}{{r_{1}}^{2}}\right]=\frac{1}{\kappa^{2}}\left[-\frac{-2\lambda^{2}r_{1}}{r_{1}}-\frac{1-\lambda^{2}{r_{1}}^{2}-1}{{r_{1}}^{2}}\right]=\frac{3\lambda^{2}}{\kappa^{2}}\,,$
(33) $\displaystyle A\left(r_{2}\right)=$
$\displaystyle\,\frac{\mathrm{e}^{2\nu(r_{2})}}{\kappa^{2}}\left[-\frac{f\left(r_{2}\right)-1}{{r_{2}}^{2}}+\frac{1}{2}f^{\prime\prime}\left(r_{2}\right)\right]=\frac{\mathrm{e}^{2\nu(r_{2})}}{\kappa^{2}}\left[-\frac{1-\frac{2M}{r_{2}}-1}{{r_{2}}^{2}}+\frac{1}{2}\left(-\frac{4M}{{r_{2}}^{3}}\right)\right]=0\,,\quad
B\left(r_{2}\right)=0\,,$ $\displaystyle C\left(r_{2}\right)=$
$\displaystyle\,\frac{\mathrm{e}^{-2\nu\left(r_{2}\right)}}{\kappa^{2}}\left[\frac{f\left(r_{2}\right)-1}{r^{2}}-\frac{1}{2}f^{\prime\prime}\left(r_{2}\right)\right]=\frac{\mathrm{e}^{-2\nu\left(r_{2}\right)}}{\kappa^{2}}\left[\frac{1-\frac{2M}{r_{2}}-1}{{r_{2}}^{2}}-\frac{1}{2}\left(-\frac{4M}{{r_{2}}^{3}}\right)\right]=0\,,$
$\displaystyle V\left(r_{2}\right)=$
$\displaystyle\,\frac{1}{\kappa^{2}}\left[-\frac{f^{\prime}\left(r_{2}\right)}{r_{2}}-\frac{f\left(r_{2}\right)-1}{{r_{2}}^{2}}\right]=\frac{1}{\kappa^{2}}\left[-\frac{\frac{2M}{{r_{2}}^{2}}}{r_{2}}-\frac{1-\frac{2M}{r_{2}}-1}{{r_{2}}^{2}}\right]=0\,.$
(34)

Here we have used the equations (III.1) in the calculation of the above
equations. Eqs. (III.1) and (III.1) are, of course, the results of the
continuities of $\mathrm{e}^{2\nu}$,
$\left(\mathrm{e}^{2\nu}\right)^{\prime}$, and
$\left(\mathrm{e}^{2\nu}\right)^{\prime\prime}$ at the two boundaries
$r=r_{1}$ and $r=r_{2}$. The behavior of Eq. (III.1) are shown in Figures. 1
0(b), 0(c), and 0(d) for the case of (28) and in Figures. 2 1(b), 1(c), and
1(d) for the case of (29).

(a) The behavior of the polynomial function $f(r)$ given by Eq. (22)

(b) The behavior of the scalar field $A(r)$ given by Eqs. (III.1)

(c) The behavior of the scalar field $C(r)$ given by Eqs. (III.1)

(d) The behavior of the potential field $V(r)$ given by Eqs. (III.1)

Figure 1: Plots of Figures in the case of small $\lambda$ (28). 0(a) the
polynomial function given by Eq. (22) after using Eq. (III.1); Figures 0(b)
0(c), and 0(d) show the behavior of the scalars $A(r)$, $C(r)$, and the
potential $V(r)$ given by Eqs. (III.1).

(a) The behavior of the polynomial function $f(r)$ given by Eq. (22)

(b) The behavior of the scalar field $A(r)$ given by Eqs. (III.1)

(c) The behavior of the scalar field $C(r)$ given by Eqs. (III.1)

(d) The behavior of the potential field $V(r)$ given by Eqs. (III.1)

Figure 2: Plots of Figures in the case of thin shell (28). 0(a) the polynomial
function given by Eq. (22) after using Eq. (III.1); Figures 0(b) 0(c), and
0(d) show the behavior of the scalars $A(r)$, $C(r)$, and the potential $V(r)$
given by Eqs. (III.1).

###  III.2 Smooth function model

In this subsection, we try to connect the interior region and the exterior
region by using smooth functions, which describe the asymptotically
Schwarzschild spacetime for large $r$ but describes the asymptotically de
Sitter spacetime for small $r$.

####  III.2.1 Example 1

We now consider the following model

$\displaystyle\mathrm{e}^{2\nu}=\mathrm{e}^{-2\mu}=1-\frac{2Mr^{2}}{r^{3}+{r_{0}}^{3}}\,.$
(35)

The spacetime described by (35) has two horizons
$\mathrm{e}^{2\nu}=\mathrm{e}^{-2\mu}=0$ in general, which corresponds to the
Hayward black hole Hayward (2006). If we choose, however,

$\displaystyle\frac{2^{\frac{5}{3}}M}{3r_{0}}<1\,,$ (36)

$\mathrm{e}^{2\nu}$ and $\mathrm{e}^{-2\mu}$ are always positive and there
does not appear any horizon. In the following, we only consider the case that
Eq. (36) is satisfied.

When $r$ is small, $\mathrm{e}^{2\nu}=\mathrm{e}^{-2\mu}$ behaves as
$\mathrm{e}^{2\nu}=\mathrm{e}^{-2\mu}\sim 1-\frac{2M}{{r_{0}}^{3}}r^{2}$,
which is the de Sitter spacetime and we identify
$\lambda^{2}=\frac{2M}{{r_{0}}^{3}}$. On the other hand, when $r$ is large, we
find $\mathrm{e}^{2\nu}=\mathrm{e}^{-2\mu}\sim 1-\frac{2M}{r}$, which is
nothing but the Schwarzschild spacetime. Therefore the function (35) smoothly
connects the exterior region, which is the asymptotically Schwarzschild
spacetime, and the interior region, which is the asymptotically de Sitter
spacetime.

By using (II), we can find the explicit form of $A$, $B$, $C$, and $V$, as
follows,

$\displaystyle A=$
$\displaystyle\,\frac{18M{r_{0}}^{3}r^{3}}{\kappa^{2}\left(r^{3}+{r_{0}}^{3}\right)^{3}}\left(1-\frac{2Mr^{2}}{r^{3}+{r_{0}}^{3}}\right)\,,\quad
B=0\,,\quad
C=-\frac{18M{r_{0}}^{3}r^{3}}{\kappa^{2}\left(r^{3}+{r_{0}}^{3}\right)^{3}}\left(1-\frac{2Mr^{2}}{r^{3}+{r_{0}}^{3}}\right)^{-1}\,,$
$\displaystyle V=$
$\displaystyle\,\frac{6M{r_{0}}^{3}}{\kappa^{2}\left(r^{3}+{r_{0}}^{3}\right)^{2}}\,.$
(37)

We should note $A$ is always positive but $C$ is negative. Therefore $\chi$
becomes a ghost if we do not add the constraint terms (15). Due to the
constraint terms (15), both of $\phi$ and $\chi$ become non-dynamical and the
ghosts can be eliminated. Because $\phi$ is canonical although $\chi$ is not
canonical, we may impose the constraint on only $\chi$ by using

$\displaystyle S_{\lambda_{\chi}}=\int
d^{4}x\sqrt{-g}\lambda_{\chi}\left(\mathrm{e}^{-2\mu(r=\chi)}\partial_{\mu}\chi\partial^{\mu}\chi-1\right)\,.$
(38)

Instead of Eq. (15) so that only $\chi$ is not dynamical.

####  III.2.2 Example 2

Now let us use the following function

$\displaystyle\mathrm{e}^{2\nu(r)}=\mathrm{e}^{-2\mu(r)}=1-\frac{4M\arctan\left(\frac{\pi
r^{3}}{L^{3}}\right)}{\pi r}\,.$ (39)

When $r$ is large, we find that the spacetime becomes the asymptotically
Schwarzschild spacetime $\mathrm{e}^{2\nu(r)}=\mathrm{e}^{-2\mu(r)}\to
1-\frac{2M}{r}$ and when $r$ is small, we obtain the de Sitter spacetime
$\mathrm{e}^{2\nu(r)}=\mathrm{e}^{-2\mu(r)}\to 1-\frac{4\pi M}{L^{3}}r^{2}$ by
identifying $\lambda^{2}=\frac{4\pi M}{L^{3}}$. The behavior of Eq. (39) is
shown in Figure 3 2(a).

Using Eq. (39) in Eq. (II), we obtain

$\displaystyle A=$
$\displaystyle\,\frac{36Mr^{6}L^{3}\pi^{2}}{\kappa^{2}\left(L^{6}+r^{6}\pi^{2}\right)^{2}}\left(1-\frac{4M\arctan\left(\frac{r^{3}\pi}{L^{3}}\right)}{\pi
r}\right)\,,$ $\displaystyle C=$
$\displaystyle\,\frac{4M\left\\{2(L^{6}+r^{6}\pi^{2})^{2}\arctan\left(\frac{r^{3}\pi}{L^{3}}\right)-9r^{9}L^{3}\pi^{3}\right\\}}{\kappa^{2}r^{2}\left(L^{6}+r^{6}\pi^{2}\right)^{2}\left[4M\arctan\left(\frac{r^{3}\pi}{L^{3}}\right)-\pi
r\right]}\,,$ $\displaystyle V=$
$\displaystyle\,\frac{6M\left\\{\left(L^{6}+r^{6}\pi^{2}\right)\arctan\left(\frac{r^{3}\pi}{L^{3}}\right)-r^{3}L^{3}\pi\right\\}}{\kappa^{2}r^{3}\pi\left(L^{6}+r^{6}\pi^{2}\right)}\,.$
(40)

Again, we find $A$ is always positive but $C$ is negative. Therefore $\chi$
becomes a ghost and we need the constraint terms in (38). The behavior of Eq.
(III.2.2) are shown in Figures 3 2(b), 2(c), and 2(d) for $M=0.006654411765$
and $\lambda=\sqrt{\frac{4\pi M}{L^{3}}}=0.1$.

(a) The behavior of the function $\mathrm{e}^{2\nu}$ given by Eq. (39)

(b) The behavior of the scalar field $A(r)$ given by Eqs. (III.2.2)

(c) The behavior of the scalar field $C(r)$ given by Eqs. (III.2.2)

(d) The behavior of the potential field $V(r)$ given by Eqs. (III.2.2)

Figure 3: Plots of Fig. 2(a) the function given by Eq. (39); Figures 2(b)
2(c), and 2(d) show the behavior of the scalars $A(r)$, $C(r)$, and the
potential $V(r)$ given by Eq. (III.2.2).

##  IV The physical characteristic of the model

We assume the shell region that connects the interior region and the exterior
region is rather thin and therefore the matter density composed of the scalar
fiels could be high. We now investigate the proper thickness and the surface
redshift.

###  IV.1 Proper length of shell

In accordance with the proposal of Mazur and Mottola, the incompressible fluid
within the shell is located at the intersection of two separate regions. The
shell spans from radius $r_{1}$ (indicating the surface between the insdie and
the shell) to radius $r_{2}$ (indicating the surface between the shell and the
external spacetime), with the assumption that
$\lambda\left(r_{2}-r_{1}\right)\ll 1$. This allows for the calculation of the
proper spacing separating these two interfaces using the following equation:

$\displaystyle l=\int_{r_{1}}^{r_{2}}\mathrm{e}^{\mu(r)}dr\,.$ (41)

In this subsection, we only consider the case that the shell is given by the
polynomial in (22) because the boundaries of the shell are well defined to be
$r=r_{1}$ and $r=r_{2}$.

By using Eq. (22) in Eq. (41), we obtain,

$\displaystyle l\approx$
$\displaystyle\,\frac{r_{2}-r_{1}}{\mathcal{N}}-\frac{3\lambda^{2}\left({r_{2}}^{2}-{r_{1}}^{2}\right)r_{2}{r_{1}}^{2}\left(-15{r_{2}}^{3}+{r_{2}}^{2}r_{1}+r_{2}{r_{1}}^{2}+{r_{1}}^{3}\right)}{4\mathcal{N}\mathcal{M}}$
$\displaystyle\,+\frac{\left({r_{2}}^{3}-{r_{1}}^{3}\right)}{3\mathcal{N}}\left(\frac{\lambda^{2}\left(-10{r_{2}}^{5}-10{r_{2}}^{4}r_{1}-55{r_{2}}^{3}{r_{1}}^{2}+19{r_{2}}^{2}{r_{1}}^{3}+r_{2}{r_{1}}^{4}+{r_{1}}^{5}\right)}{2\mathcal{M}}\right.$
$\displaystyle\,\left.+\frac{27\lambda^{4}{r_{2}}^{2}{r_{1}}^{4}\left(-15{r_{2}}^{3}+{r_{2}}^{2}r_{1}+r_{2}{r_{1}}^{2}+{r_{1}}^{3}\right)^{2}}{8\mathcal{M}^{2}}\right)+\frac{\left({r_{2}}^{4}-{r_{1}}^{4}\right)}{4\mathcal{N}}\left(-\frac{3\lambda^{2}r_{2}\left(-5{r_{2}}^{3}-5{r_{2}}^{2}r_{1}-5r_{2}{r_{1}}^{2}+3{r_{1}}^{3}\right)}{2\mathcal{M}}\right.$
$\displaystyle\,-\frac{9\lambda^{4}r_{2}{r_{1}}^{2}\left(-15{r_{2}}^{3}+{r_{2}}^{2}r_{1}+r_{2}{r_{1}}^{2}+{r_{1}}^{3}\right)\left(-10{r_{2}}^{5}-10{r_{2}}^{4}r_{1}-55{r_{2}}^{3}{r_{1}}^{2}+19{r_{2}}^{2}{r_{1}}^{3}+r_{2}{r_{1}}^{4}+{r_{1}}^{5}\right)}{4\mathcal{M}^{2}}$
$\displaystyle\,\left.-\frac{135\lambda^{6}{r_{2}}^{3}{r_{1}}^{6}\left(-15{r_{2}}^{3}+{r_{2}}^{2}r_{1}+r_{2}{r_{1}}^{2}+{r_{1}}^{3}\right)^{3}}{16\mathcal{M}^{3}}\right)+\frac{\left({r_{2}}^{5}-{r_{1}}^{5}\right)}{5\mathcal{N}}\left(\frac{3\lambda^{2}r_{2}\,\left(-2{r_{2}}^{2}-2r_{2}r_{1}+{r_{1}}^{2}\right)}{2\mathcal{M}}\right.$
$\displaystyle\,+\frac{27\lambda^{4}{r_{2}}^{2}{r_{1}}^{2}\left(-15{r_{2}}^{3}+{r_{2}}^{2}r_{1}+r_{2}{r_{1}}^{2}+{r_{1}}^{3}\right)\left(-5{r_{2}}^{3}-5{r_{2}}^{2}r_{1}-5r_{2}{r_{1}}^{2}+3{r_{1}}^{3}\right)}{4\mathcal{M}^{2}}$
$\displaystyle\,+\frac{3\lambda^{4}\left(-10{r_{2}}^{5}-10{r_{2}}^{4}r_{1}-55{r_{2}}^{3}{r_{1}}^{2}+19{r_{2}}^{2}{r_{1}}^{3}+r_{2}{r_{1}}^{4}+{r_{1}}^{5}\right)^{2}}{8\mathcal{M}^{2}}$
$\displaystyle\,+\frac{135\lambda^{6}{r_{2}}^{2}{r_{1}}^{4}\left(-15{r_{2}}^{3}+{r_{2}}^{2}r_{1}+r_{2}{r_{1}}^{2}+{r_{1}}^{3}\right)^{2}\left(-10{r_{2}}^{5}-10{r_{2}}^{4}r_{1}-55{r_{2}}^{3}{r_{1}}^{2}+19{r_{2}}^{2}{r_{1}}^{3}+r_{2}{r_{1}}^{4}+{r_{1}}^{5}\right)}{16\mathcal{M}^{3}}$
$\displaystyle\,\left.+\frac{2835\lambda^{8}{r_{2}}^{4}{r_{1}}^{8}\left(-15{r_{2}}^{3}+{r_{2}}^{2}r_{1}+r_{2}{r_{1}}^{2}+{r_{1}}^{3}\right)^{4}}{128\mathcal{M}^{4}}\right)\,.$
(42)

Here

$\displaystyle\mathcal{N}\equiv$
$\displaystyle\,\sqrt{-\frac{\mathcal{M}}{\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)\left(r_{2}-r_{1}\right)^{3}}}\,,$
$\displaystyle\mathcal{M}\equiv$
$\displaystyle\,\lambda^{2}{r_{2}}^{2}{r_{1}}^{3}\left(15{r_{2}}^{2}-3r_{2}r_{1}-3{r_{1}}^{2}\right)-10{r_{2}}^{5}+35{r_{2}}^{4}r_{1}-46{r_{2}}^{3}{r_{1}}^{2}+28{r_{2}}^{2}{r_{1}}^{3}-8r_{2}{r_{1}}^{4}-{r_{1}}^{5}\,.$
(43)

For example, in the case of (28), we find $\lambda l=0.0201\cdots$, which is
not so changed from $\lambda\left(r_{2}-r_{1}\right)=0.02$.

###  IV.2 Surface redshift

In this particular subsection, we will examine the surface redshift. The
examination of gravastar surface redshift stands as a crucial means for
assessing its detection and the surface redshift has been often used to
investigate the stability of the gravastar configuration. Surface
gravitational redshift, denoted as $Z_{s}$, is defined by using the wavelength
$\lambda_{e}$ of the emitted signal at the surface of the gravastar and the
wavelength $\lambda_{0}$ of the observed signal, that is,
$Z_{s}=\frac{\lambda_{0}-\lambda_{e}}{\lambda_{e}}$.

Buchdahl’s claim, founded on the conditions of isotropy, staticity, and a
perfect fluid distribution, suggests that the surface redshift should remain
below 2, denoted as $Z_{s}<2$ Buchdahl (1959); Straumann (1984); Böhmer and
Harko (2006). Conversely, Ivanov Ivanov (2002) argued that in the case of an
anisotropic fluid distribution, it could increase to as high as 3.84.
Furthermore, Barraco and Hamity Barraco and Hamity (2002) showed that in the
absence of a cosmological constant, for an isotropic fluid distribution,
$Z_{s}$ stays at a value less than or equal to 2. Nevertheless, Bohmer and
Harko Böhmer and Harko (2006) demonstrated that in an anisotropic star with
the influence of a cosmological constant, $Z_{s}$ can reach as high as 5. In
our model, it is not always possible to represent the scalar fields using the
fluid approximation. The constraints outlined in (16) were employed to emulate
the behavior of dark matter, as described in Chamseddine et al. (2014). Once
more, it’s essential to emphasize that effective dark matter is not actual
dark matter, and it does not undergo gravitational collapse. In the model
presented in this paper, the arrangement remains stable when subjected to
perturbations. To compare with earlier studies employing conventional fluids,
we assess the surface redshift within the current model. To calculate the
surface redshift, we utilize the subsequent equation:

$\displaystyle Z_{s}=-1+\frac{1}{\sqrt{\left|g_{tt}\right|}}\,,$ (44)

which is evaluated at the surface of the gravastar. Therefore in the
polynomial model (22), we find

$\displaystyle Z_{s}=-1+\frac{1}{\sqrt{1-\frac{2M}{r_{2}}}}\,,$ (45)

because the surface is given by $r=r_{2}$. In the case of small $\lambda$ in
(28), we obtain

$\displaystyle Z_{s}\sim\frac{M}{r_{2}}\fallingdotseq 0.0054\cdots\,,$ (46)

which is very small and therefore the value is in the stable region. On the
other hand, in the case of thin shell in (29), we find,

$\displaystyle Z_{s}=9\,,$ (47)

which is beyond the stable region in the previous works Buchdahl (1959);
Straumann (1984); Böhmer and Harko (2006); Ivanov (2002); Barraco and Hamity
(2002); Böhmer and Harko (2006) and in the region prohibited by the
instability. We should note, again, that the solution in the present model is
always stable under fluctuations. Therefore we can realize the stable
gravastar with large surface redshift.

In the case of smooth function models in (35) and (39), the position of the
surface is not well-defined. We now evaluate the redshift in some region of
$r$. In the case of (35), by using Eq. (44) for general $r$, we obtain

$\displaystyle Z_{s}=-1+\frac{1}{\sqrt{\frac{2Mr^{2}}{r^{3}+{r_{0}}^{3}}}}\,,$
(48)

and for Eq. (39), we obtain the following redshift as,

$\displaystyle Z_{s}=-1+\frac{1}{\sqrt{1-\frac{4M\arctan\left(\frac{\pi
r^{3}}{L^{3}}\right)}{\pi r}}}\,.$ (49)

The behavior of Eq. (48) is shown in Figure 4 for $M=0.006654411765$ and
$\lambda=\sqrt{\frac{4\pi M}{L^{3}}}=0.1$, again, which ensures that the
redshift in our model is in the stable region given in Buchdahl (1959);
Straumann (1984); Böhmer and Harko (2006). The behavior of redshift given by
Eq. (49) has the same pattern as Eq. (48) which tells us that the redshift is
in the stable region.

Figure 4: Plot of the redshift given by Eq. (48).

##  V Effective equation of state parameter

The matter that composes the shell of the gravastar obeys a rather excentric
equation of state (EoS). In our model, although the shell is composed of
scalar fields, we investigate what kind of matter corresponds to the shell in
the polynomial model (22).

From Eqs. (7), (8), and (9), one can find the energy density $\rho$, radial
and tangential pressures $p_{r}$ and $p_{t}$, as follows,

$\displaystyle\rho=$
$\displaystyle\,\frac{\kappa^{2}}{2}\left[2\mathrm{e}^{2\nu\left(r\right)}V\left(r\right)+A\left(r\right)+C\left(r\right)\mathrm{e}^{2\nu\left(r\right)-2\mu\left(r\right)}\right]\,,$
$\displaystyle p_{r}=$
$\displaystyle\,\frac{\kappa^{2}}{2}\left[C\left(r\right)+A\left(r\right)\mathrm{e}^{2\mu\left(r\right)-2\nu\left(r\right)}-2\mathrm{e}^{2\mu\left(r\right)}V\left(r\right)\right]\,,$
$\displaystyle p_{t}=$
$\displaystyle\,\frac{\kappa^{2}}{2}\left[A\left(r\right)\mathrm{e}^{2\mu\left(r\right)-2\nu\left(r\right)}-2\mathrm{e}^{2\mu\left(r\right)}V\left(r\right)-C\left(r\right)\right]\,.$
(50)

We use the above expressions for the models proposed in this paper.

In the case that the shell is described by the polynomial function (22), by
substituting (22) with (III.1) and Eq. (III.1), we obtain the radial and
tangential EoS parameters $\omega_{r}$ and $\omega_{t}$ as follows,

$\displaystyle\omega_{r}\equiv$
$\displaystyle\,\frac{p_{r}}{\rho}=-\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)^{2}\left(r_{2}-r_{1}\right)^{6}\left[10{r_{2}}^{5}-35{r_{2}}^{4}r_{1}+46{r_{2}}^{3}{r_{1}}^{2}-28{r_{2}}^{2}{r_{1}}^{3}+8r_{2}{r_{1}}^{4}-{r_{1}}^{5}\right.$
$\displaystyle\,+\lambda^{2}\left\\{\left(-6{r_{2}}^{3}-6{r_{2}}^{2}r_{1}+3r_{2}{r_{1}}^{2}\right)r^{4}+\left(15{r_{2}}^{4}+15{r_{2}}^{3}r_{1}+15{r_{2}}^{2}{r_{1}}^{2}-9r_{2}{r_{1}}^{3}\right)r^{3}\right.$
$\displaystyle\,+\left(-10{r_{2}}^{5}-10{r_{2}}^{4}r_{1}-55{r_{2}}^{3}{r_{1}}^{2}+19{r_{2}}^{2}{r_{1}}^{3}+r_{2}{r_{1}}^{4}+{r_{1}}^{5}\right)r^{2}+\left(+45{r_{2}}^{4}{r_{1}}^{2}-3{r_{2}}^{3}{r_{1}}^{3}-3{r_{2}}^{2}{r_{1}}^{4}-3r_{2}{r_{1}}^{5}\right)r$
$\displaystyle\,\left.\left.-15{r_{2}}^{4}{r_{1}}^{3}+3{r_{2}}^{3}{r_{1}}^{4}+3{r_{2}}^{2}{r_{1}}^{5}\right\\}\right]^{-2}\,,$
$\displaystyle\omega_{t}\equiv$
$\displaystyle\,\frac{p_{t}}{\rho}=r\left(r_{2}-r_{1}\right)^{6}\left(10{r_{2}}^{2}-5r_{2}r_{1}+{r_{1}}^{2}\right)^{2}\left\\{\left(-16{r_{2}}^{2}-16r_{2}r_{1}+8{r_{1}}^{2}\right)r^{2}+\left(15{r_{2}}^{3}-{r_{2}}^{2}r_{1}-r_{2}{r_{1}}^{2}-{r_{1}}^{3}\right)r\right.$
$\displaystyle\left.+15{r_{2}}^{3}r_{1}-{r_{2}}^{2}{r_{1}}^{2}-r_{2}{r_{1}}^{3}-{r_{1}}^{4}\right\\}\left(r-r_{1}\right)^{-1}\left\\{\left(-4{r_{2}}^{2}-4r_{2}r_{1}+2{r_{1}}^{2}\right)r^{2}+\left(5{r_{2}}^{3}-3{r_{2}}^{2}r_{1}-3r_{2}{r_{1}}^{2}+r{r_{1}}^{3}\right)r\right.$
$\displaystyle\,\left.+10{r_{2}}^{3}r_{1}-2{r_{2}}^{2}{r_{1}}^{2}-2r_{2}{r_{1}}^{3}\right\\}^{-1}\left[10{r_{2}}^{5}-35{r_{2}}^{4}r_{1}+46{r_{2}}^{3}{r_{1}}^{2}-28{r_{2}}^{2}{r_{1}}^{3}+8r_{2}{r_{1}}^{4}-{r_{1}}^{5}\right.$
$\displaystyle\,+\lambda^{2}\left\\{\left(-6{r_{2}}^{3}-6{r_{2}}^{2}r_{1}+3r_{2}{r_{1}}^{2}\right)r^{4}+\left(15{r_{2}}^{4}+15{r_{2}}^{3}r_{1}+15{r_{2}}^{2}{r_{1}}^{2}-9r_{2}{r_{1}}^{3}\right)r^{3}\right.$
$\displaystyle\,+\left(-10{r_{2}}^{5}-10{r_{2}}^{4}r_{1}-55{r_{2}}^{3}{r_{1}}^{2}+19{r_{2}}^{2}{r_{1}}^{3}+r_{2}{r_{1}}^{4}+{r_{1}}^{5}\right)r^{2}$
$\displaystyle\,\left.\left.+\left(45{r_{2}}^{4}{r_{1}}^{2}-3{r_{2}}^{3}{r_{1}}^{3}-3{r_{2}}^{2}{r_{1}}^{4}-3r_{2}{r_{1}}^{5}\right)r-15{r_{2}}^{4}{r_{1}}^{3}+3{r_{2}}^{3}{r_{1}}^{4}+3{r_{2}}^{2}{r_{1}}^{5}\right\\}\right]^{-2}\,.$
(51)

The behaviors of Eqs. (V) are shown in Figures 5 4(a) and 4(b) for
sufficiently small $\lambda$ case in (28) and in Figures 6 5(a) and 5(b) for
sufficiently thin shell case in (29). As clear from these Figures, the
signature of the radial EoS parameter $\omega_{r}$ is different from that of
the tangential one $\omega_{t}$ and especially the value of the tangential EoS
parameters are very huge. This could tell that the gravastar geometry realized
in this paper by using two scalar fields could be difficult to be generated by
using ordinary fluid.

(a) The behavior of the radial EoS given by Eq. (V)

(b) The behavior of the tangential EoS given by Eq. (V)

Figure 5: Plots of Figures for sufficiently small $\lambda$ case in (28), 4(a)
the radial EoS given by Eq. (V); 4(b) the tangential EoS given by Eq. (V).

(a) The behavior of the radial EoS given by Eq. (V)

(b) The behavior of the tangential EoS given by Eq. (V)

Figure 6: Plots of Figures for sufficiently thin shell case in (29), 4(a) the
radial EoS given by Eq. (V); 4(b) the tangential EoS given by Eq. (V).

From Eqs. (35) and (III.2.1) for Example 1 in the smooth function model, one
can find the energy density $\rho$, radial and tangential pressures $p_{r}$
and $p_{t}$, as follows,

$\displaystyle\rho=$
$\displaystyle\,-\frac{3\left(\left(-\lambda^{2}+r^{2}\right)M-\frac{r^{3}}{2}\right)M^{2}\lambda^{2}}{\left(M\lambda^{2}+\frac{r^{3}}{2}\right)^{3}}\,,$
$\displaystyle p_{r}=$
$\displaystyle\,-\frac{12\lambda^{2}M^{2}}{\left(r^{3}+2M\lambda^{2}\right)\left(r^{3}+2M\lambda^{2}-2Mr^{2}\right)}\,,$
$\displaystyle p_{t}=$
$\displaystyle\,\,{\frac{3M^{2}\lambda^{2}\left(M\lambda^{2}-r^{3}\right)}{\left(M\lambda^{2}+\frac{r^{3}}{2}\right)^{2}\left(\left(-\lambda^{2}+r^{2}\right)M-\frac{r^{3}}{2}\right)}}\,.$
(52)

Here $\lambda=\sqrt{\frac{2M}{{r_{0}}^{3}}}$. Using Eq. (V), we obtain the
radial and tangential EoS parameters as follows,

$\displaystyle\omega_{r}=$
$\displaystyle\,\frac{\left(4M\lambda^{2}+\frac{r^{3}}{2}\right)^{3}}{\left(r^{3}+2M\lambda^{2}\right)\left(r^{3}+2M\lambda^{2}-2Mr^{2}\right)\left(\left(-\lambda^{2}+r^{2}\right)M-\frac{r^{3}}{2}\right)}\,,$
$\displaystyle\omega_{t}=$
$\displaystyle\,-\frac{\left(M\lambda^{2}-r^{3}\right)\left(M\lambda^{2}+\frac{r^{3}}{2}\right)}{\left(\left(-\lambda^{2}+r^{2}\right)M-\frac{r^{3}}{2}\right)^{2}}\,.$
(53)

The behavior of Eqs. (V) are shown in Figures 7 6(a) and 6(b) for
$M=0.006654411765$ and $\lambda=\sqrt{\frac{2M}{{r_{0}}^{3}}}=0.1$, which
satisfies the condition (36). Although the signature of the radial EoS
parameter is different from that of the tangential one, the values are not
surprisingly large compared with the polynomial cases in Figures 5 and 6.

(a) The behavior of the radial EoS given by Eq. (V)

(b) The behavior of the tangential EoS given by Eq. (V)

Figure 7: Plots of Fig. 6(a) the radial EoS given by Eq. (V); 6(b) the
tangential EoS given by Eq. (V).

From Eqs. (39), and (III.2.2) for Example 2 in the smooth function model, one
can find the energy density $\rho$, radial and tangential pressures $p_{r}$,
and $p_{t}$, as follows,

$\displaystyle\rho=$
$\displaystyle\,\frac{8M}{r^{4}\left(L^{6}+r^{6}\pi^{2}\right)^{2}\pi^{2}}\left(\left(L^{6}+r^{6}\pi^{2}\right)^{2}\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)+15r^{9}L^{3}\pi^{3}-3r^{3}L^{9}\pi\right)\left(\frac{\pi
r}{4}-M\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)\right)\kappa^{2}\,,$
$\displaystyle p_{r}=$
$\displaystyle\,-\frac{10M\kappa^{2}}{r^{2}\left(\pi\,r-4\,M\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)\right)\left(L^{6}+r^{6}\pi^{2}\right)^{2}}\left(\left(L^{6}+r^{6}\pi^{2}\right)^{2}\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)-\frac{3L^{3}r^{3}\pi}{5}\left(7r^{6}\pi^{2}+L^{6}\right)\right)\,,$
$\displaystyle p_{t}=$ $\displaystyle\,\frac{2M\kappa^{2}}{\left(-\pi
r+4M\arctan\left(\frac{r^{3}\pi}{L^{3}}\right)\right)r^{2}\left(L^{6}+r^{6}\pi^{2}\right)}\left(\left(L^{6}+r^{6}\pi^{2}\right)\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)-3r^{3}L^{3}\pi\right)\,.$
(54)

Using Eq. (V), we obtain the radial and tangential EoS parameters as follows,

$\displaystyle\omega_{r}=$
$\displaystyle\,-\frac{\left(5\left(L^{6}+r^{6}\pi^{2}\right)^{2}\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)-3L^{3}r^{3}\pi\left(7r^{6}\pi^{2}+L^{6}\right)\right)r^{2}\pi^{2}}{\left(\pi
r-4M\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)\right)^{2}\left(\left(L^{6}+r^{6}\pi^{2}\right)^{2}\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)+15r^{9}L^{3}\pi^{3}-3r^{3}L^{9}\pi\right)}\,,$
$\displaystyle\omega_{t}\equiv$
$\displaystyle\,-\frac{\left(\left(L^{6}+r^{6}\pi^{2}\right)\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)-3r^{3}L^{3}\pi\right)r^{2}\left(L^{6}+r^{6}\pi^{2}\right)\pi^{2}}{\left(\pi
r-4M\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)\right)^{2}\left(\left(L^{6}+r^{6}\pi^{2}\right)^{2}\arctan\left({\frac{r^{3}\pi}{L^{3}}}\right)+15r^{9}L^{3}\pi^{3}-3r^{3}L^{9}\pi\right)}\,.$
(55)

The behavior of Eqs. (V) are shown in Figures 8 7(a) and 7(b) for
$M=0.006654411765$ and $\lambda=\sqrt{\frac{4\pi M}{L^{3}}}=0.1$. In this
case, both EoS parameters are negative.

(a) The behavior of the radial EoS given by Eq. (V)

(b) The behavior of the tangential EoS given by Eq. (V)

Figure 8: Plots of Fig. 7(a) the radial EoS given by Eq. (V); 7(b) the
tangential EoS given by Eq. (V).

##  VI Summary and discussion

The intention behind the gravastar solution, incorporating a modified matter
source within the framework of Einstein’s GR, is to investigate alternative
gravastar models and examine the consequences they may entail. Gravastars are
theoretical entities suggested as potential substitutes for black holes. These
objects are believed to consist of exotic material capable of obstructing the
creation of an event horizon, a characteristic typically associated with black
holes. In contrast, gravastars feature a surface referred to as a
“gravitational vacuum star”, which can replicate certain attributes of a black
hole but lacks the singularity at its core. In this particular scenario, the
background is assumed to be Einstein’s GR with two scalar fields.

Applying the spherically symmetric spacetime to the field equation of
Einstein’s GR with two scalar fields, we derive the coefficients $A(r)$ and
$C(r)$ that characterize the scalar functions. We show that the product of
these two functions is always negative, which means that either one of the
scalars is a ghost. We use constraints similar to those of mimetic theory to
remove such ghosts. The constraints also make the solution in the model stable
under the perturbation.

We assume that the interior of the gravastar is characterized by the de Sitter
geometry, which is characterized by the cosmological constant $\lambda$, and
the exterior by the Schwarzschild one, which is characterized by the
gravitational mass $M$. On the shell, we assume a polynomial of fifth order
that has six unknown constants (22). These constants are fixed through the
smooth junction condition in the interior geometry and the interior of the
shell and by the smooth junction condition between the exterior geometry and
the exterior shell. Through these fixations, we show that the polynomial
function has two maxima when $r<r_{1}$ and $r>r_{2}$, respectively, and one
minimum when $r_{1}<r<r_{2}$, where $r_{1}$ and $r_{2}$ are the radii of the
boundaries between the shell region and the interior and exterior regions,
respectively. Through these constraints, we succeeded in writing the two
functions $A$ and $C$ and the potential $V$ in terms of $r_{1}$, $r_{2}$, and
$\lambda$. We show that on the boundaries of the shell $r=r_{1}$ and
$r=r_{2}$, the junction conditions are satisfied also for these functions,
$A$, $C$, and $V$.

To clarify our philosophy, we give two additional examples: For small $r$
asymptotes of de Sitter spacetime and for large $r$ asymptotes of
Schwarzschild geometry. The first example was given by the geometry similar to
the Hayward black hole solution but without horizons. The second one was the
one given by equation (39). These two examples give the de Sitter for small
$r$ and Schwarzschild geometry for large $r$. We derive the explicative forms
of $A$, $C$, and $V$ of these two models and show their behaviors in Figure 1,
Figure 2, and Figure 3, which ensures our analytical derivation. We calculated
the proper length and (surface) redshift of these models as given by Eqs.
(22), (35), and (39) and we have shown that there exists a stable gravastar
with a large surface redshift prohibited by the instability in the previous
works. Furthermore, we check the behavior of the effective EoS parameters and
show that the gravastar geometry could be difficult to be realized by the
ordinary fluid.

In conclusion, our study of gravastars within the framework of Einstein’s GR,
enriched by the inclusion of two scalar fields, offers a promising avenue for
furthering our understanding of exotic astrophysical objects. We have shown
that these scalar fields have a substantial impact on the properties and
behavior of gravastars, making them more than just theoretical constructs.

Our research not only reaffirms the potential of gravastars as alternatives to
black holes but also highlights the versatility and complexity introduced by
scalar fields. The interplay between these fields and the gravastar’s
structure opens up a wide array of possibilities for future investigations in
the fields of astrophysics and theoretical physics.

As we continue to explore the intricacies of gravastars with scalar fields, we
look forward to unveiling new layers of understanding in the ever-evolving
narrative of our universe. This study represents a vital contribution to the
ongoing quest to unravel the mysteries of compact astrophysical objects and
their role in the cosmos.

Is the procedure carried out in this study applicable in the frame of $f(R)$
theory with two scalar fields or within the frame of Gauss-Bonnet theory with
two scalar fields? The answer to these questions will be found elsewhere.

## References

  * Mitra (2002) A. Mitra, Found. Phys. Lett. 15, 439 (2002), eprint astro-ph/0207056.
  * Mazur and Mottola (2004) P. O. Mazur and E. Mottola, Proc. Nat. Acad. Sci. 101, 9545 (2004), eprint gr-qc/0407075.
  * Abbott et al. (2016) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), eprint 1602.03837.
  * Visser and Wiltshire (2004) M. Visser and D. L. Wiltshire, Class. Quant. Grav. 21, 1135 (2004), eprint gr-qc/0310107.
  * Cattoen et al. (2005) C. Cattoen, T. Faber, and M. Visser, Class. Quant. Grav. 22, 4189 (2005), eprint gr-qc/0505137.
  * Carter (2005) B. M. N. Carter, Class. Quant. Grav. 22, 4551 (2005), eprint gr-qc/0509087.
  * Bilic et al. (2006) N. Bilic, G. B. Tupper, and R. D. Viollier, JCAP 02, 013 (2006), eprint astro-ph/0503427.
  * Lobo (2006) F. S. N. Lobo, Class. Quant. Grav. 23, 1525 (2006), eprint gr-qc/0508115.
  * DeBenedictis et al. (2006) A. DeBenedictis, D. Horvat, S. Ilijic, S. Kloster, and K. S. Viswanathan, Class. Quant. Grav. 23, 2303 (2006), eprint gr-qc/0511097.
  * Lobo and Arellano (2007) F. S. N. Lobo and A. V. B. Arellano, Class. Quant. Grav. 24, 1069 (2007), eprint gr-qc/0611083.
  * Horvat and Ilijic (2007) D. Horvat and S. Ilijic, Class. Quant. Grav. 24, 5637 (2007), eprint 0707.1636.
  * Chirenti and Rezzolla (2007) C. B. M. H. Chirenti and L. Rezzolla, Class. Quant. Grav. 24, 4191 (2007), eprint 0706.1513.
  * Rocha et al. (2008) P. Rocha, R. Chan, M. F. A. da Silva, and A. Wang, JCAP 11, 010 (2008), eprint 0809.4879.
  * Horvat et al. (2009) D. Horvat, S. Ilijic, and A. Marunovic, Class. Quant. Grav. 26, 025003 (2009), eprint 0807.2051.
  * Usmani et al. (2008) A. A. Usmani, P. P. Ghosh, U. Mukhopadhyay, P. C. Ray, and S. Ray, Mon. Not. Roy. Astron. Soc. 386, L92 (2008), eprint 0801.4529.
  * Turimov et al. (2009) B. V. Turimov, B. J. Ahmedov, and A. A. Abdujabbarov, Mod. Phys. Lett. A 24, 733 (2009), eprint 0902.0217.
  * Nandi et al. (2009) K. K. Nandi, Y. Z. Zhang, R. G. Cai, and A. Panchenko, Phys. Rev. D 79, 024011 (2009), eprint 0809.4143.
  * Harko et al. (2009) T. Harko, Z. Kovacs, and F. S. N. Lobo, Class. Quant. Grav. 26, 215006 (2009), eprint 0905.1355.
  * Usmani et al. (2011) A. A. Usmani, F. Rahaman, S. A. Rakib, S. Ray, K. K. Nandi, P. K. F. Kuhfittig, and Z. Hasan, Phys. Lett. B 701, 388 (2011), eprint 1012.5605.
  * Rahaman et al. (2012a) F. Rahaman, S. Ray, A. A. Usmani, and S. Islam, Phys. Lett. B 707, 319 (2012a), eprint 1108.4824.
  * Rahaman et al. (2012b) F. Rahaman, A. A. Usmani, S. Ray, and S. Islam, Phys. Lett. B 717, 1 (2012b), eprint 1205.6796.
  * Bhar (2014) P. Bhar, Astrophys. Space Sci. 354, 2109 (2014).
  * Rahaman et al. (2015) F. Rahaman, S. Chakraborty, S. Ray, A. A. Usmani, and S. Islam, Int. J. Theor. Phys. 54, 50 (2015), eprint 1209.6291.
  * Ghosh et al. (2017) S. Ghosh, F. Rahaman, B. K. Guha, and S. Ray, Phys. Lett. B 767, 380 (2017), eprint 1511.05417.
  * Ghosh et al. (2018) S. Ghosh, S. Ray, F. Rahaman, and B. K. Guha, Annals Phys. 394, 230 (2018), eprint 1701.01046.
  * Ghosh et al. (2019a) S. Ghosh, S. Biswas, F. Rahaman, B. K. Guha, and S. Ray, Annals Phys. 411, 167968 (2019a).
  * Ghosh et al. (2019b) S. Ghosh, D. Shee, S. Ray, F. Rahaman, and B. K. Guha, Results Phys. 14, 102473 (2019b).
  * Chan and da Silva (2010) R. Chan and M. F. A. da Silva, JCAP 07, 029 (2010), eprint 1005.3703.
  * Chan et al. (2011) R. Chan, M. F. A. da Silva, J. F. V. da Rocha, and A. Wang, JCAP 10, 013 (2011), eprint 1109.2062.
  * Riess et al. (1998) A. G. Riess et al. (Supernova Search Team), Astron. J. 116, 1009 (1998), eprint astro-ph/9805201.
  * Perlmutter et al. (1999) S. Perlmutter et al. (Supernova Cosmology Project), Astrophys. J. 517, 565 (1999), eprint astro-ph/9812133.
  * de Bernardis et al. (2000) P. de Bernardis et al. (Boomerang), Nature 404, 955 (2000), eprint astro-ph/0004404.
  * Hanany et al. (2000) S. Hanany et al., Astrophys. J. Lett. 545, L5 (2000), eprint astro-ph/0005123.
  * Peebles and Ratra (2003) P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559 (2003), eprint astro-ph/0207347.
  * Padmanabhan (2003) T. Padmanabhan, Phys. Rept. 380, 235 (2003), eprint hep-th/0212290.
  * Clifton et al. (2012) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Physics reports 513, 1 (2012).
  * Riess et al. (2007) A. G. Riess et al., Astrophys. J. 659, 98 (2007), eprint astro-ph/0611572.
  * Tegmark et al. (2004) M. Tegmark et al. (SDSS), Phys. Rev. D 69, 103501 (2004), eprint astro-ph/0310723.
  * Amanullah et al. (2010) R. Amanullah et al., Astrophys. J. 716, 712 (2010), eprint 1004.1711.
  * Komatsu et al. (2011) E. Komatsu et al. (WMAP), Astrophys. J. Suppl. 192, 18 (2011), eprint 1001.4538.
  * Das et al. (2015) A. Das, F. Rahaman, B. K. Guha, and S. Ray, Astrophys. Space Sci. 358, 36 (2015), eprint 1507.04959.
  * Das et al. (2016) A. Das, F. Rahaman, B. Guha, and S. Ray, The European Physical Journal C 76, 1 (2016).
  * Das et al. (2017) A. Das, S. Ghosh, B. K. Guha, S. Das, F. Rahaman, and S. Ray, Phys. Rev. D 95, 124011 (2017), eprint 1702.08873.
  * Biswas et al. (2019) S. Biswas, S. Ghosh, S. Ray, F. Rahaman, and B. K. Guha, Annals Phys. 401, 1 (2019), eprint 1803.00442.
  * Ghosh et al. (2020) S. Ghosh, A. D. Kanfon, A. Das, M. J. S. Houndjo, I. G. Salako, and S. Ray, Int. J. Mod. Phys. A 35, 2050017 (2020), eprint 2003.03194.
  * Das et al. (2020) A. Das, S. Ghosh, D. Deb, F. Rahaman, and S. Ray, Nucl. Phys. B 954, 114986 (2020), eprint 2003.10842.
  * Sengupta et al. (2020) R. Sengupta, S. Ghosh, S. Ray, B. Mishra, and S. K. Tripathy, Phys. Rev. D 102, 024037 (2020), eprint 2004.01519.
  * Banerjee et al. (2020) S. Banerjee, S. Ghosh, N. Paul, and F. Rahaman, Eur. Phys. J. Plus 135, 185 (2020).
  * Nojiri et al. (2021) S. Nojiri, S. D. Odintsov, and V. Faraoni, Phys. Rev. D 103, 044055 (2021), eprint 2010.11790.
  * Chamseddine and Mukhanov (2013) A. H. Chamseddine and V. Mukhanov, JHEP 11, 135 (2013), eprint 1308.5410.
  * Chamseddine et al. (2014) A. H. Chamseddine, V. Mukhanov, and A. Vikman, JCAP 06, 017 (2014), eprint 1403.3961.
  * Nojiri and Nashed (2022) S. Nojiri and G. G. L. Nashed, Phys. Lett. B 830, 137140 (2022), eprint 2202.03693.
  * Myrzakulov et al. (2016) R. Myrzakulov, L. Sebastiani, S. Vagnozzi, and S. Zerbini, Class. Quant. Grav. 33, 125005 (2016), eprint 1510.02284.
  * Nashed and Nojiri (2022) G. G. L. Nashed and S. Nojiri, JCAP 05, 011 (2022), eprint 2110.08560.
  * Nashed (2023a) G. G. L. Nashed, Nucl. Phys. B 993, 116264 (2023a), eprint 2307.03199.
  * Nashed (2023b) G. G. L. Nashed, Gen. Rel. Grav. 55, 63 (2023b), eprint 2305.09694.
  * Myrzakulov et al. (2015) R. Myrzakulov, L. Sebastiani, and S. Vagnozzi, Eur. Phys. J. C 75, 444 (2015), eprint 1504.07984.
  * Vagnozzi (2017) S. Vagnozzi, Class. Quant. Grav. 34, 185006 (2017), eprint 1708.00603.
  * Nashed and Nojiri (2021) G. G. L. Nashed and S. Nojiri, Phys. Rev. D 104, 044043 (2021), eprint 2107.13550.
  * Casalino et al. (2018) A. Casalino, M. Rinaldi, L. Sebastiani, and S. Vagnozzi, Phys. Dark Univ. 22, 108 (2018), eprint 1803.02620.
  * Casalino et al. (2019) A. Casalino, M. Rinaldi, L. Sebastiani, and S. Vagnozzi, Class. Quant. Grav. 36, 017001 (2019), eprint 1811.06830.
  * Sebastiani et al. (2017) L. Sebastiani, S. Vagnozzi, and R. Myrzakulov, Adv. High Energy Phys. 2017, 3156915 (2017), eprint 1612.08661.
  * Nashed (2021) G. G. L. Nashed, Astrophys. J. 919, 113 (2021), eprint 2108.04060.
  * (64) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 875 (2019) no.1, L4 doi:10.3847/2041-8213/ab0e85 [arXiv:1906.11241 [astro-ph.GA]].
  * Hayward (2006) S. A. Hayward, Phys. Rev. Lett. 96, 031103 (2006), eprint gr-qc/0506126.
  * Buchdahl (1959) H. A. Buchdahl, Physical Review 116, 1027 (1959).
  * Straumann (1984) N. Straumann, _GENERAL RELATIVITY AND RELATIVISTIC ASTROPHYSICS_ (1984).
  * Böhmer and Harko (2006) C. Böhmer and T. Harko, Classical and Quantum Gravity 23, 6479 (2006).
  * Ivanov (2002) B. V. Ivanov, Physical Review D 65, 104011 (2002).
  * Barraco and Hamity (2002) D. Barraco and V. H. Hamity, Phys. Rev. D 65, 124028 (2002).