#  Defrosting frozen stars:
spectrum of internal fluid modes

Ram Brustein(1), A.J.M. Medved(2,3), Tom Shindelman(1)
${}^{\textrm{\normalsize(1)\ Department of Physics, Ben-Gurion University,
Beer-Sheva 84105, Israel}}$ ${}^{\textrm{\normalsize(2)\ Department of Physics
\& Electronics, Rhodes University, Grahamstown 6140, South Africa}}$
${}^{\textrm{\normalsize(3)\ National Institute for Theoretical Physics
(NITheP), Western Cape 7602, South Africa}}$<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The frozen star model provides a classical description of a regularized black
hole and is based upon the idea that regularizing the singularity requires
deviations from the Schwarzschild geometry which extend over horizon-sized
scales, as well as maximally negative radial pressure as an equation of state.
The frozen star has also been shown to be ultra-stable against perturbations;
a feature that can be attributed to the equation of state and corresponds to
this model mimicking a black hole in the limit $\hbar\to 0$ or, equivalently,
the limit of infinite Newton’s constant. Here, we “defrost” the frozen star by
allowing its radial pressure to be perturbatively less negative than maximal.
This modification to the equation of state is implemented by appropriately
deforming the background metric so as to allow the frozen star to mimic a
quantum black hole at finite $\hbar$ and Newton’s constant. As a consequence,
the defrosted star acquires a non-trivial spectrum of oscillatory
perturbations. To show this, we first use the Cowling approximation to obtain
generic equations for the energy density and pressure perturbations of a
static, spherically symmetric background with an anisotropic fluid. The
particular setting of a deformed frozen star is then considered, for which the
dispersion relation is obtained to leading order in terms of the deviation
from maximal pressure. The current results compare favorably with those
obtained earlier for the collapsed polymer model, whose strongly non-classical
interior is argued to provide a microscopic description of the frozen and
defrosted star geometries.

##  1 Introduction

The name “frozen star”, which was first used by Ruffini [1] to describe
objects that are now commonly known as black holes (BHs), has recently been
resurrected [2, 3] to describe a particular model for a regularized BH
geometry [4, 5]. It is a highly compact object, whose compactness is
parametrically close to that of a BH, and yet it manages to evade singularity
theorems [6, 7] and “Buchdahl-like” matter bounds [8, 9, 10, 11, 12] by
invoking a maximally negative (radial) pressure, $\;p=-\rho\;$, 111Here, we
use $p$ to denote radial pressure, $q$ for the transverse pressure components
and $\rho$ is the energy density. while not possessing a trapped surface. It
also avoids certain inconsistencies by deviating from the standard general-
relativistic solution over horizon-sized scales [13, 14]. As far as we know,
the frozen star is the first completely regular and stable solution of
Einstein’s equations whose compactness is arbitrarily close to that of a
Schwarzschild BH. Furthermore, the solution obeys the only universally agreed-
upon energy condition, the null energy condition. It is free of any pathology,
unlike most, if not all, of the other proposals for regular compact objects.

Having such a solution at hand is important as it provides a model for a
compact object — differing from a general-relativistic BH — that can be used
to estimate the validity of general relativity by looking at the resulting
dynamics and gravitational-waves emissions during a binary-BH merger event.
Additionally, it provides a final state in which to address the process of
gravitational collapse (a paper that includes such a discussion is already in
progress).

An unusual feature of the frozen star that distinguishes it from the gravastar
[15] and others of the negative-pressure class (e.g.) [16, 17] is that two of
its metric component are identically vanishing, $\;g_{tt}=g^{rr}=0\;$.
However, we do not mean to imply that astrophysical BHs are made of an
anisotropic classical fluid with maximally negative radial pressure. The
anisotropy and maximal negative pressure turned out to be needed so as to
connect this solution to the collapsed polymer BH [18], which is another model
for a regularized compact object that is based, in large part, on the notion
of having a maximally entropic interior [19]. The “fluid” that supports this
maximal entropy was found to consist of long, closed, fundamental, highly
excited, interacting strings. From a geometric point of view, maximal entropy
translates into a mass profile that saturates the Schwarzschild limit at all
interior radii, $\;2m(r)=r\;$, and then $\;g_{tt}=g^{rr}=0\;$ follows from
Einstein’s equations. The catch is that the polymer model lacks a
semiclassical geometry, as maximally entropic is synonymous with strongly non-
classical. Hence, the desire for a formal classical description, and so a
frozen star model was born. The frozen star is, therefore, merely a classical
proxy of a quantum model, the polymer BH, which has strictly non-negative
pressure. The main purpose of the frozen star model is to facilitate some
calculations that would present a challenge to someone using the quantum
model.

The frozen star and polymer models are then meant to provide complementary
pictures of the same model for a regular, ultra-compact object. The former
allows for precise calculations, whereas the latter provides a microscopic
description with connections to string theory. With that in mind, $g^{rr}$ and
$|g_{tt}|$ should be small enough to maintain the frozen stars’ character, but
not that much changes if these components are alternatively assigned a small
but finite constant, which we denote as $\;\varepsilon\ll 1\;$. What is
maintained is the frozen star’s ultra-stability against small perturbations of
the metric and fluid densities [4, 2]. This stability is due to the equation
of state $\;p+\rho=0\;$. To see this, one can inspect the stability
calculation in [2], where the vanishing of the metric components (but not the
equation-of-state relation) was relaxed for a portion of the interior.
Additionally, the limit $\;p+\rho\to 0\;$ can be taken in the current analysis
to verify this conclusion. This ultra-stability was also found for the polymer
model in the absence of quantum effects, which corresponds in that case to the
absence of closed-string coupling or $g_{s}^{2}$ corrections [20, 21]. The
reason being that, in the limit $\;g_{s}^{2}\to 0\;$ or, equivalently,
$\;\hbar\to 0\;$, the polymer model behaves just like a hairless classical BH.
It was also shown in [21] that, once stringy quantum effects are taken into
account, the polymer model does exhibit a spectrum of excitations and some
quantum hair.

The main objective of this paper is to determine the dispersion relation for
the oscillatory modes of a appropriately deformed version of the frozen star;
the “defrosted star”. To have such modes at all, it is clear that the
background must be suitably modified to support fluctuations. Technically
speaking, this means that the original equation of state needs to be changed
such that the radial pressure is no longer exactly equal to the negative of
the energy density and the transverse pressure components — which are
identically zero when the star is truly frozen — no longer vanish. The change
that we consider is parametrically small and is supposed to model some
microscopic corrections to the equation of state of the constituents. As just
discussed, in the polymer model, these corrections can be attributed to
quantum perturbations that are proportional to the square of the string
coupling strength [21].

The spectrum of oscillatory modes for the defrosted star could be relevant to
BH merger events. If astrophysical BHs are indeed described by such objects,
the internal modes of the star could be excited during a merger event and then
act as sources for gravitational-wave emissions. After all, it is the success
of the LIGO and Virgo collaboration at detecting and then analyzing the waves
from such mergers — starting with the now world-famous event GW150914 [22] and
the data analysis thereof [23] — that has led to an intensive theoretical
effort to understand how departures from general relativity would imprint on
the gravitational-wave spectrum of a BH. In this regard, each candidate for a
regular BH mimicker can be expected to have its own unique and experimentally
verifiable or falsifiable signature (see [24] for further discussion and a
catalogue of relevant models).

We also have in mind a previous study that determined the mode spectrum for
the collapsed polymer model [21] (also see [25]). We would like to verify that
the two models produce compatible results, even though the polymer analysis
was, at times, heuristic in nature because of the absence of an interior
geometry. It is already notable that, in both cases, non-trivial mode
solutions require some out-of-equilibrium physics [26].

The spectral analysis for the collapsed polymer model utilized the Cowling
approximation [27], which assumes, for a calculation of the fluid-mode
spectrum, that the spacetime modes have fully decoupled. For the sake of
consistency and simplicity, we will use the same approximation here. Although
it is, in any event, a justifiable assumption insofar as the Cowling
approximation only gets more accurate as the star gets more compact [28].

The rest of the paper is arranged as follows: We start by deriving the
required perturbative formalism for an anisotropic fluid, with our only
assumptions being that the background geometry is static and spherically
symmetric. 222We also insist on $3+1$ spacetime dimensions for concreteness.
This part of the analysis culminates in a coupled pair of equations that
describe the dynamics of two types of modes: one is associated with
fluctuations of the radial velocity and the other with those of the transverse
velocity. Some brief background material on the frozen star geometry is then
provided; after which the deformations are introduced and explained in some
detail. We then return to the dynamical equations for the modes but now with
the focus on the solutions for the defrosted or deformed frozen star. The
scaling behavior of the modes and their dispersion relation are identified and
discussed. The paper ends with a brief overview followed by an appendix that
explains in a thorough way how to obtain one of the key equations in the
generic portion of the perturbative analysis.

##  2 Perturbation equations: Generic formalism

In the following, the perturbation equations for the energy density and
pressure of a star containing anisotropic fluid are reproduced in the Cowling
approximation, meaning that the fluctuations of the metric are neglected. We
will be closely following and reviewing the analysis of [29], while providing
additional details as needed to ensure compatibility of the generic formalism
with that of the defrosted star model.

###  2.1 Variation of the stress tensor

The background under consideration is, for the meantime, described by a
static, spherically symmetric but otherwise generic geometry,

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

The corresponding stress (energy–momentum) tensor can be expressed as

$T_{\mu\nu}\;=\;\rho
u_{\mu}u_{\nu}+pk_{\mu}k_{\nu}+q\left(g_{\mu\nu}+u_{\mu}u_{\nu}-k_{\mu}k_{\nu}\right)\;,$
(2)

where the transverse pressure $q$ is generally different from the radial
pressure $p$, so that $\;\sigma=p-q\neq 0\;$. Also, $u^{\mu}$ is the fluid
4-velocity with standard normalization, $\;u^{\mu}u_{\mu}=-1\;$, and $k^{\mu}$
is a radial unit vector that is defined by $\;k^{\mu}k_{\mu}=+1\;$ such that
$\;u^{\mu}k_{\mu}=0\;$.

The oscillatory modes of the fluid will be obtained from the perturbation
equations that arise from the variation of the conservation equation
$\;\nabla_{\nu}T_{\mu}^{\leavevmode\nobreak\ \nu}=0\;$. In the Cowling
approximation, the variation of interest is

$\nabla_{\nu}\delta T_{\mu}^{\leavevmode\nobreak\ \nu}\;=\;0\;,$ (3)

and the variation of Eq. (2) under this approximation yields

$\displaystyle\delta T_{\mu}^{\nu}\;=\;$ $\displaystyle\left(\delta\rho+\delta
q\right)u_{\mu}u^{\nu}+\left(\rho+q\right)\left(\delta
u_{\mu}u^{\nu}+u_{\mu}\delta u^{\nu}\right)$ (4) $\displaystyle+$
$\displaystyle\delta q\delta_{\mu}^{\nu}+\delta\sigma
k_{\mu}k^{\nu}+\sigma\delta k_{\mu}k^{\nu}+\sigma k_{\mu}\delta k^{\nu}\;.$
(5)

The velocity $u^{\mu}$ satisfies some useful identities,

$\nabla_{\nu}u^{\nu}\;=\;u^{\mu}\nabla_{\nu}u_{\mu}\;=\;u^{\mu}\delta
u_{\mu}\;=\;0\;,$ (6)

and similarly for the radial unit vector $k^{\mu}$. By choosing, without loss
of generality, to work in the comoving frame, we also have
$\;u^{\mu}=u^{t}\delta_{t}^{\mu}\;$, $\;k^{\mu}=k^{r}\delta_{r}^{\mu}\;$ and
$\;u^{t}u_{t}=-k^{r}k_{r}=-1\;$.

Equation (3) leads to an independent pair of perturbation equations when it is
projected, respectively, parallel and perpendicular to $u^{\mu}$; the latter
by way of the projection operator
$\;\mathcal{P}_{\alpha}^{\mu}=\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\;$.
Taking the covariant derivative of Eq. (5) and projecting onto $u^{\mu}$, one
obtains a scalar equation of the form

$\displaystyle 0\;=\;u^{\mu}\nabla_{\nu}\delta T_{\mu}^{\nu}\;=\;$
$\displaystyle-\nabla_{\nu}\delta\rho
u^{\nu}-\nabla_{\nu}\left[\left(\left(\rho+q\right)\delta_{\mu}^{\nu}+\sigma
k^{\nu}k_{\mu}\right)\delta u^{\mu}\right]$
$\displaystyle-\left(\rho+q\right)a_{\mu}\delta
u^{\mu}-\nabla_{\nu}u^{\mu}\delta\left(\sigma k^{\nu}k_{\mu}\right)\;,$ (7)

where $\;a_{\mu}=u^{\nu}\nabla_{\nu}u_{\mu}\;$ is the 4-acceleration. Whereas,
in the direction perpendicular to $u^{\mu}$, one obtains the following vector
equation (see Appendix A for details):

$\displaystyle 0\;=\;\mathcal{P}_{\alpha}^{\mu}\nabla_{\nu}\delta
T_{\mu}^{\nu}\;=\;$
$\displaystyle\delta\left(\rho+q\right)a_{\alpha}+\left(\rho+q\right)u^{\nu}\left(\nabla_{\nu}\delta
u_{\alpha}-\nabla_{\alpha}\delta u_{\nu}\right)$
$\displaystyle+\nabla_{\nu}\delta
q\delta_{\alpha}^{\nu}+u^{\nu}u_{\alpha}\nabla_{\nu}\delta
q+\mathcal{P}_{\alpha}^{\mu}\nabla_{\nu}\delta\left(\sigma
k_{\mu}k^{\nu}\right)\;.$ (8)

###  2.2 Density and pressure perturbations

It proves to be useful if the perturbations of the spacelike components of the
velocity vector are expressed in terms of a displacement vector $\xi^{i}$ such
that

$\frac{\partial\xi^{i}}{\partial t}\;=\;\frac{\delta u^{i}}{u^{t}}={v}_{i}\;,$
(9)

where $\;i=\\{r,\theta,\phi\\}\;$ and $v_{i}$ is a component of the fluid’s
3-velocity.

As the velocity is a function of $t$ and $r$ only,
$\;a_{\theta}=a_{\phi}=0\;$. Then the angular components of Eq. (7) reduce to

$\left(\rho+q\right)\left(u^{t}\right)^{2}\partial_{t}^{2}\xi_{\theta}+\partial_{\theta}\delta
q\;=\;0\;,$ (10)
$\left(\rho+q\right)\left(u^{t}\right)^{2}\partial_{t}^{2}\xi_{\phi}+\partial_{\phi}\delta
q\;=\;0\;.$ (11)

Differentiating the former equation by $\phi$ and the latter by $\theta$, one
is led to a consistency condition,

$\partial_{\phi}\xi_{\theta}\;=\;\partial_{\theta}\xi_{\phi}\;,$ (12)

which implies that $\xi_{\theta}$ and $\xi_{\phi}$ can be expressed as angular
derivatives of a multipole expansion [30],

$\displaystyle\xi_{\theta}$ $\displaystyle\;=\;-\sum_{\ell,m}V_{\ell
m}\left(r,t\right)\partial_{\theta}Y_{\ell m}\left(\theta,\phi\right)\;,$ (13)
$\displaystyle\xi_{\phi}$ $\displaystyle\;=\;-\sum_{\ell,m}V_{\ell
m}\left(r,t\right)\partial_{\phi}Y_{\ell m}\left(\theta,\phi\right)\;.$ (14)

The definition of $\xi^{i}$ enables one to obtain the variation in the energy
density by integrating Eq. (7) over time, which eventually leads to

$\displaystyle\delta\rho\;=\;$
$\displaystyle-\frac{1}{\sqrt{-g}}\partial_{\nu}\left[\sqrt{-g}\left(\left(\rho+q\right)\delta_{i}^{\nu}+\sigma
k^{\nu}k_{i}\right)\xi^{i}\right]$
$\displaystyle-\left[\left(\left(\rho+q\right)\delta_{i}^{\nu}+\sigma
k^{\nu}k_{i}\right)\xi^{i}\right]\partial_{\nu}\left(\ln
u^{t}\right)-\left(\rho+q\right)a_{i}\xi^{i}\;.$ (15)

A convenient form for the radial component of the displacement vector is the
following:

$\xi^{r}\;=\;\sum_{\ell,m}e^{-\Lambda}\frac{W_{\ell
m}\left(r,t\right)}{r^{2}}Y_{\ell m}\;.$ (16)

To avoid clutter, we will henceforth be suppressing the subscripts $\ell$ and
$m$ on $W$ and $V$, as well as the accompanying factors of $Y_{\ell m}$ and
the summation over $\ell,m$ in any equation involving $W$ and/or $V$. These
indices, factors and summations can be easily restored at the end of the
calculation.

We next employ $\;u^{t}=e^{-\Phi\left(r\right)}\;$,
$\;a_{\mu}=a_{r}\delta_{\leavevmode\nobreak\ \mu}^{r}=\partial_{r}\Phi\;$, the
previous forms for $\xi^{i}$ and the equation for energy conservation,

$\frac{dp}{dr}\;=\;-\partial_{r}\Phi\left(\rho+p\right)-\frac{2\sigma}{r}\;,$
(17)

to recast Eq. (15) as

$\displaystyle\delta\rho\;=\;$
$\displaystyle-\left(\rho+p\right)\left[e^{-\Lambda}\frac{W^{\prime}}{r^{2}}+\frac{\ell\left(\ell+1\right)}{r^{2}}V\right]-\frac{d\rho}{dr}e^{-\Lambda}\frac{W}{r^{2}}+\frac{2\sigma}{r^{3}}e^{-\Lambda}W$
$\displaystyle+\sigma\frac{\ell\left(\ell+1\right)}{r^{2}}V\;,$ (18)

where a prime denotes a differentiation by $r$.

The perturbation of the radial pressure can be obtained by first recalling the
relation between the Eulerian ($\delta\rho$) and Lagrangian ($\Delta\rho$)
variations of the energy density,

$\delta\rho\;=\;\Delta\rho-\xi^{r}\partial_{r}\rho\;.$ (19)

Given an equation of state of the form $\;p=p\left(\rho\right)\;$, the
Lagrangian perturbations for the radial pressure and density can be directly
related,

$\Delta
p\;=\;\frac{dp}{d\rho}\Delta\rho\;=\;\frac{dp}{d\rho}\left(\delta\rho+\xi^{r}\partial_{r}\rho\right)\;,$
(20)

and the Eulerian variation of the pressure goes similarly as

$\displaystyle\delta p$ $\displaystyle=$ $\displaystyle\Delta
p-\xi^{r}\partial_{r}p$ (21) $\displaystyle=$
$\displaystyle\frac{dp}{d\rho}\left(\delta\rho+\xi^{r}\partial_{r}\rho\right)-\xi^{r}\partial_{r}p$
$\displaystyle=$ $\displaystyle\frac{dp}{d\rho}\;\delta\rho\;,$

where the last line follows from
$\;\left(dp/d\rho\right)\partial_{r}\rho=\partial_{r}p\;$.

More explicitly,

$\displaystyle\delta p\;=\;$
$\displaystyle\frac{dp}{d\rho}\left(-\left(\rho+p\right)\left[e^{-\Lambda}\frac{W^{\prime}}{r^{2}}+\frac{\ell\left(\ell+1\right)}{r^{2}}V\right]+\frac{2\sigma}{r^{3}}e^{-\Lambda}W+\sigma\frac{\ell\left(\ell+1\right)}{r^{2}}V\right)$
$\displaystyle-\frac{dp}{dr}e^{-\Lambda}\frac{W}{r^{2}}\;.$ (22)

Although the Eulerian variation $\delta p$ is the quantity that determines the
form of the oscillatory modes, the Lagrangian variation $\Delta p$ is the
quantity that is directly relevant to the boundary conditions. One can see
from Eq. (16) and the first line of Eq. (21) that the difference between the
two variations is the same as the sole term in the second line of Eq. (22).
Hence, the Lagrangian variation $\Delta p$ can be obtained directly from Eq.
(22) by simply dropping the second line.

###  2.3 Dynamical equations for oscillatory modes

We are interested in oscillating modes of the form
$\;W\left(r,t\right)=W\left(r\right)e^{i\omega t}\;$ and likewise for $V$. The
dynamical equations for $W$ and $V$ follow from Eq. (8) for $\;\mu=r\;$ and
$\;\mu=\theta\;$, respectively, and take on the corresponding forms,

$\displaystyle 0\;=\;$
$\displaystyle-\omega^{2}\left(\rho+p\right)e^{\Lambda-2\Phi}\frac{W}{r^{2}}+\partial_{r}{\delta
p}+\left({\delta\rho}+{\delta
p}\right)\partial_{r}\Phi+\frac{2}{r}{\delta\sigma}\;$ (23)

and

$0\;=\;\left(\rho+q\right)e^{-2\Phi}\omega^{2}V+{\delta q}\;,$ (24)

where, just like for the $W$ and $V$ modes, the $\ell$, $m$ indices on
$\delta\rho$, the associated summations and accompanying spherical harmonics
are all implied in the above and subsequent equations. So $\delta\rho$ here is
meant to represent ${\delta\rho}_{\ell m}Y_{\ell m}$ such that the total
Eulerian density perturbation is given by
$\;\sum\limits_{\ell,m}{\delta\rho}_{\ell m}Y_{\ell m}\;$.

The next step in [29] and similar discussions relies on a formal equation-of-
state relation between $\rho$ and $p$ and also between $q$ (or $\sigma$) and
$p$, as these enable one to calculate $\frac{\partial\rho}{\partial p}$ and
then $\;\delta\rho=\frac{\partial\rho}{\partial p}\delta p\;$ (and similarly
for $\delta q$). For the frozen star model, one has $\;\rho+p=0\;$ but no such
relation exists between $q$ and $p$ because the transverse pressure is
identically vanishing. Nevertheless, the situation changes for the defrosted
star, as the $\rho$–$p$ relation picks up a perturbative correction and we are
now able to relate $q$ and $p$ by varying (in an Euler–Lagrange sense) the
conservation equation (17), which leads to $\;\frac{\partial q}{\partial p}\;$
and thus $\;\delta q=\frac{\partial q}{\partial p}\delta p\;$.

Expressing ${\delta\rho}$ and ${\delta q}$ in terms of ${\delta p}$ and
background quantities allows us to rewrite both of the dynamical equations
(23,24) in terms of $\delta p$, which itself is determined by Eq. (22). The
final result is then a pair of coupled equations in terms of a pair of unknown
functions, $W(r)$ and $V(r)$, for any given values of $\ell$ and $m$. Thanks
to the appearances of $\partial_{r}{\delta p}$ in Eq. (23) and $W^{\prime}$ in
Eq. (22), one can, after diagonalizing, expect to obtain a second-order
differential equation for $W(r)$ and first-order equation for $V(r)$.

The dynamical equations are to be supplemented by a pair of boundary
conditions involving the Lagrangian variation $\Delta p$ of the radial
pressure. In particular, the Lagrangian variation should vanish both at the
star’s outer surface $\;r=R\;$ and at the star’s center $\;r=0\;$.
Importantly, this variation should also remain finite as it approaches the
center of the star; that is, for $\;r\ll R\;$.

##  3 Frozen star: A refresher

To formulate the frozen star model, one starts with a metric and stress tensor
just like those in Eqs. (1) and (2). The first key ingredient is that the
interior is endowed with a maximally negative radial pressure, $\;p=-\rho\;$,
throughout. Consistency between the Einstein field equations and the energy-
conservation equation (17) is what dictates the form of the transverse
pressure $q$.

Next, let us define a mass function $m(r)$ in the standard way,

$m(r)\;=\;4\pi\int\limits_{0}^{r}dx\,x^{2}\rho(x)\;\;\;{\rm for}\;\;\;r\leq
R\;,$ (25)

for which Einstein’s equations indicate that

$e^{-2\Lambda}\;=\;1-\frac{2m(r)}{r}\;.$ (26)

It is the choice of $m(r)$ that specifies the exact nature of the model within
this negative-pressure class. The most popular option is to choose $m$ such
that $\rho$ is constant, with the result that $\;q=p=-\rho\;$; what is known
as the gravastar model [15]. In its originally prescribed form, the frozen
star was specified by the profile $\;m(r)=r/2\;$, as this choice ensured that
each spherical shell within the star saturated the Schwarzschild bound. This
would be the second key ingredient and, just like the first, it was chosen to
make the frozen star a proxy to the (strongly non-classical) collapsed polymer
model. However, as this same mass profile also leads to a pair of apparently
singular metric components, $\;e^{-2\Lambda}=e^{2\Phi}=0\;$, we have more
recently opted to follow a different path: Namely, we now set
$\;m(r)=\frac{r}{2}(1-\varepsilon)\;$ for a small dimensionless and constant
$\varepsilon$, and so $\;e^{-2\Lambda}=e^{2\Phi}=\varepsilon\ll 1\;$. It is
implied that $\varepsilon$ is the smallest dimensionless scale in the model.

In this relaxed version of the frozen star, the energy density and pressure
components go as

$\displaystyle 8\pi G\rho$ $\displaystyle=$
$\displaystyle\frac{1-(rf)^{\prime}}{r^{2}}\;=\;\frac{1-\varepsilon}{r^{2}}\;,$
(27) $\displaystyle 8\pi Gp$ $\displaystyle=$
$\displaystyle-\frac{1-(rf)^{\prime}}{r^{2}}\;=\;-\frac{1-\varepsilon}{r^{2}}\;,$
(28) $\displaystyle 8\pi Gq$ $\displaystyle=$
$\displaystyle\frac{(rf)^{\prime\prime}}{2r}\;=\;0\;.$ (29)

Not much has changed from the $\;\varepsilon=0\;$ limit. In fact, a stability
study in [2] reveals that the ultra-stability of the model against
$r,t$-dependent perturbations persists irrespective of the inclusion of non-
zero metric components. The critical feature for ensuring stability is rather
the condition $\;p=-\rho\;$. Whether $\varepsilon$ is zero or small but finite
is besides the point.

Some final notes about the frozen star: It is necessary to modify the metric
over a thin translational layer near (and including) the outer boundary so
that the star’s metric and its first two derivatives match smoothly to the
Schwarzschild metric in the exterior [2]. Additionally, the metric near the
center of the star needs to be regularized so as to ensure that the energy
density remains finite at $\;r=0\;$ [3]. This regularization process is
mandatory for the frozen star because, unlike a Schwarzschild BH, any would-be
singularity is not protected by cosmic censorship.

The width of the transitional surface layer is taken to be very small, of
order of the string or Planck length scale, so that its inclusion in our
analysis would lead to highly suppressed corrections. We can then safely
ignore the outer layer except that its presence is relevant to the outer-
surface boundary condition. In our case, the value of $\Delta p$ should be set
equal to some non-zero value at $\;r=R\;$, thus allowing it to decay to zero
only after passing through the translational layer where it would match on to
the external geometry. The mathematical details of this process are
inconsequential to the current analysis but would follow along the lines of
[2].

Similarly, the regularized core has a similarly sized width and could also be
ignored except for the following observation (again pertaining to boundary
conditions): A particularly relevant feature of the core regularization
process in [3] is that $\rho$, $-p$ and $-q$ all tend to the same constant
value at $\;r=0\;$, yet the rate in which they approach that value is somewhat
arbitrary. We then have the freedom to send the combinations $\;\sigma=p-q\;$
and $\;\rho+p\;$ of the deformed metric (see the next section) to zero at an
arbitrarily fast rate. This is pertinent to the boundary condition at the
center because, as one can see from the first line of Eq. (22), every term in
the Lagrangian variation $\Delta p$ is preceded by just such a combination.
Meaning that the $\;r=0\;$ regularity condition can always be easily satisfied
for the defrosted star model. 333In [3], the regularization method was applied
to the undeformed version of the frozen star model, but it is a
straightforward exercise to include deformations and obtain similar results.

###  3.1 Briefly on stability

Before discussing the deformation of the frozen star, let us consider how the
dynamical equations (23) and (24) can be used to learn about the star’s
stability in the undeformed case. Starting with Eq. (23), one immediately sees
that the $W$ term vanishes because of the factor of $\;\rho+p=0\;$ and the
$a_{r}$ term because the equation of state dictates that
$\;{\delta\rho}+{\delta p}\;$=0 or, more simply, because $a_{r}$ itself is
vanishing for constant $\varepsilon$.

As we are assuming an object of fixed mass, it follows that
$\;{\delta\rho}=0\;$ and then likewise for the equal and opposite variation
$\;{\delta p}\;$. To see that $\;{\delta q}\;$ also vanishes, one can use
$\;\rho+p=0\;$ to rewrite the conservation equation (17) as

$0\;=\;r^{2}\partial_{r}p+2r\sigma\;=\;\partial_{r}(r^{2}p)-2rq\;,$ (30)

which confirms that both $q$ and $\;{\delta q}\;$ vanish because $\;r^{2}p\;$
must be constant by virtue of Eq. (28) and $\delta p=0\;$.

Since $\;\delta q=0\;$, Eq. (24) reduces to

$0\;=\;\left(\rho+q\right)e^{-2\Phi}\omega^{2}V\;,$ (31)

meaning that either $\omega$ or $V$ is vanishing. No matter which, the angular
components of the velocity perturbations are vanishing, as these are given by
the time derivative of $V$.

This leaves the radial component of the velocity perturbations, which are
determined by the time derivative of $W$. However, the dynamical equations
tell us nothing about $W$ for the undeformed star because it carries a
vanishing factor of $\rho+p$. The way around this is to actually deform the
star (as we do next) and then consider the limit as the deformation parameter
goes to zero. What one finds is that $\omega$ does indeed go to zero in this
limit, meaning that the time derivative of $W$ and thus the radial component
of the velocity perturbations must also tend to zero. Note, though, that it
can never be established that $W$ is itself vanishing. Fortunately, there is
no reason that it has to. Only the time derivative of $W$ has physical
meaning.

##  4 Defrosting a frozen star

The frozen star geometry is unusually stable against perturbations, even when
$\;\varepsilon=|g_{tt}|=g^{rr}\;$ differs from zero. For the star to sustain
dynamical modes, we have found that at least two modifications are required.
The first requirement is that the equation of state deviates from
$\;p=-\rho\;$, which necessitates that $\;|g_{tt}|\neq g^{rr}\;$. The second
is that at least one of $g_{tt}$ and $g^{rr}$ have some radial dependence. To
implement these requirements, we introduce a perturbative, dimensionless and
positive parameter $\;\gamma\ll 1\;$ which controls the difference between the
frozen and defrosted stars.

###  4.1 Deformed background geometry

Let us first re-express the relevant metric components as

$\displaystyle-g_{tt}$ $\displaystyle=$
$\displaystyle\varepsilon+\gamma\left(\frac{r}{R}\right)^{a}\;,$ (32)
$\displaystyle g^{rr}$ $\displaystyle=$
$\displaystyle\varepsilon+\gamma\left(\frac{r}{R}\right)^{b}\;,$ (33)

where the relation between the star’s outer surface $R$ and mass $M$ will be
determined shortly. It will be shown later that the constants $a$ and $b$ are
fixed, for consistency, such that $\;a=2\;$ and $\;b=0\;$.

We assume that $\;\varepsilon\ll\gamma\;$ and so neglect $\varepsilon$ in what
follows, leaving us with

$\displaystyle-g_{tt}$ $\displaystyle=$
$\displaystyle\gamma\left(\frac{r}{R}\right)^{a}\;,$ (34) $\displaystyle
g^{rr}$ $\displaystyle=$ $\displaystyle\gamma\left(\frac{r}{R}\right)^{b}\;.$
(35)

Furthermore, only the leading-order terms in $\gamma$ will ever be considered.
We assume that regularization procedures similar to those used in [2] and [3]
have been implemented but, as already discussed, these are expected to have no
tangible affect on our analysis except to trivialize the enforcement of the
central boundary condition and modify the outer boundary condition.

Let us now consider two components of Einstein’s equations,

$\;r^{2}\rho\;=\;\left[r\left(1-g^{rr}\right)\right]^{\prime}\;,$ (36)
$\;r^{2}p\;=\;g^{rr}-1+rg^{rr}\left[\ln{|g_{tt}|}\right]^{\prime}\;,$ (37)

where the convention $\;8\pi G=1\;$ has now been adopted. Then, for the
deformed frozen star,

$r^{2}\rho\;=\;1-(1+b)\leavevmode\nobreak\
\gamma\left(\frac{r}{R}\right)^{b}\;,$ (38)
$r^{2}p\;=\;-1+(1+a)\leavevmode\nobreak\
\gamma\left(\frac{r}{R}\right)^{b}\;,$ (39)

which can be combined into

$p\;=\;\rho\left(-1+(a-b)\leavevmode\nobreak\
\gamma\left(\frac{r}{R}\right)^{b}\right)\;$ (40)

and

$r^{2}(\rho+p)\;=\;(a-b)\leavevmode\nobreak\
\gamma\left(\frac{r}{R}\right)^{b}\;,$ (41)

thus making the deviation from $\;\rho+p=0\;$ quite clear. Satisfying the null
energy condition requires that $\;a\geq b\;$, which we will assume.

Via the energy conservation equation (17), it can be shown that $q$ is no
longer vanishing,

$r^{2}q\;=\;\frac{1}{4}\left(a^{2}+ab+2b\right)\leavevmode\nobreak\
\gamma\left(\frac{r}{R}\right)^{b}\;,$ (42)

and consequently that

$q\;=\;\frac{1}{4}\left(\frac{a^{2}+ab+2b}{a-b}\right)(\rho+p)\;.\;$ (43)

We can use the above equation to determine $\frac{\partial p}{\partial q}$,
which will be needed later on for the calculation of the perturbation $\delta
q$,

$\displaystyle\frac{\partial q}{\partial p}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(\frac{a^{2}+ab+2b}{a-b}\right)\left(\frac{\partial\rho}{\partial
p}+1\right)$ (44) $\displaystyle=$
$\displaystyle-\frac{1}{4}\left(a^{2}+ab+2b\right)\leavevmode\nobreak\
\gamma\left(\frac{r}{R}\right)^{b}\;,$ (45)

where the variation $\frac{\partial\rho}{\partial p}$ can be obtained by
inverting Eq. (40).

The following relation will also be required:

$\partial_{r}\Phi\;=\;\frac{a}{r}\;.$ (46)

The radius of a defrosted star of mass $M$ is larger by a factor of order
$\frac{\gamma}{R}$ from its Schwarzschild size. The way to see this is to
match $\;g^{rr}=\gamma\left(\frac{r}{R}\right)^{b}\;$ to its Schwarzschild
form $\;g^{rr}=1-\frac{2GM}{r}\;$ at $\;r=R\;$, which yields
$\;R=\frac{2GM}{1-\gamma}\simeq{2GM}(1+\gamma)\;$. One can use this to indeed
verify that the total mass of the star is the same mass $M$ as its undeformed
counterpart. With $E$ denoting the deformed star’s mass,

$\displaystyle E$ $\displaystyle=$ $\displaystyle\int\limits^{R}_{0}dr\;4\pi
r^{2}\;\rho(r)\;=\;\frac{1}{2G}\int\limits^{R}_{0}dr\left(1-(b+1)\gamma\left(\frac{r}{R}\right)^{b}\right)$
(47) $\displaystyle=$ $\displaystyle\frac{1-\gamma}{2G}R\;=\;M\;,$

where Eq. (38) (with a restored factor of $8\pi G$ on its left side) has been
used in the top line.

###  4.2 Perturbations

We start this discussion by recalling the pressure perturbation in Eq. (22),

$\displaystyle r^{2}\delta
p\;=\;(\rho+p)e^{-\Lambda}{W^{\prime}}-\left[\partial_{r}(r^{2}p)-2rq\right]e^{-\Lambda}\frac{W}{r^{2}}-\leavevmode\nobreak\
p\leavevmode\nobreak\ \ell\left(\ell+1\right){V}\;.$ (48)

To obtain this form from the previous one, we have used that the factors in
front of $e^{-\Lambda}W^{\prime}$ and $e^{-\Lambda}W$ are linear in $\gamma$
(the latter by using the conservation equation) and that, as shown later,
444See the scaling relations that immediately follow Eq. (64). the leading-
order contribution to $V$ is also linear in $\gamma$. These observations
allowed us, in particular, to set $\frac{\partial p}{\partial\rho}=-1\;$,
recalling that our calculations are performed to leading order in $\gamma$.

Substituting the background quantities for the deformed geometry from the
previous section, we find that, to leading order in $\gamma$,

$\displaystyle r^{2}\delta
p\;=\;(a-b)\gamma\left(\frac{r}{R}\right)^{b}\leavevmode\nobreak\
\left[e^{-\Lambda}\frac{W^{\prime}}{r^{2}}+\frac{a}{2}e^{-\Lambda}\frac{W}{r^{3}}\right]+\leavevmode\nobreak\
\ell\left(\ell+1\right)\frac{V}{r^{2}}\;.$ (49)

The two other perturbations $\delta\rho$ and $\delta q$ can be related to
$\delta p$ by using Eqs. (40), (41) and (45). From Eq. (40),

$\delta\rho\;=\;\frac{\partial\rho}{\partial p}\delta p\;\simeq\;-\delta p\;,$
(50)

and then from Eq. (41),

$\delta\rho+\delta p\;\simeq\;-(a-b)\gamma\left(\frac{r}{R}\right)^{b}\delta
p\;,$ (51)

and finally, from Eq. (45),

$\delta q\;=\;\frac{\partial q}{\partial p}\delta
p\;\simeq\;-\frac{1}{4}\left(a^{2}+ab+2b\right)\gamma\left(\frac{r}{R}\right)^{b}\delta
p\;.$ (52)

As expected, the perturbation in the transverse pressure is subleading to the
radial perturbation, $\;\delta q\sim\gamma\leavevmode\nobreak\ \delta p\;$,
same as for the deformed background components,
$\;q\sim\gamma\leavevmode\nobreak\ p\;$.

##  5 Frozen star: Dynamical equations and dispersion relations

The initial step here is to rewrite the dynamical equations (23,24) for $W$
and $V$ in terms of the deformed frozen star geometry. Recall the suppression
of the angular indices, accompanying spherical harmonics and summations.

First, multiplying Eq. (23) by $r^{2}$ and then rewriting, we find that

$\displaystyle-$
$\displaystyle\frac{\omega^{2}}{\gamma^{2}\left(\frac{r}{R}\right)^{a+b}}\left((a-b)\leavevmode\nobreak\
\gamma\left(\frac{r}{R}\right)^{b}\right)e^{-\Lambda}\frac{W}{r^{2}}$ (53)
$\displaystyle+$ $\displaystyle\partial_{r}(r^{2}\delta
p)+\frac{1}{2}ar\left({\delta\rho}+{\delta p}\right)-2r\delta q\;=\;0\;.\;$
(54)

Meanwhile, for Eq. (24), the resulting expression is

$\frac{\omega^{2}}{\gamma^{2}\left(\frac{r}{R}\right)^{a}}\;\gamma\;\frac{V}{r^{2}}+{\delta
q}\;=\;0\;.$ (55)

We may now substitute the expressions for $\delta q$ from Eq. (52) and for
${\delta\rho}+{\delta p}$ from Eq. (51) into the two previous equations,
giving

$\displaystyle-$
$\displaystyle\frac{\omega^{2}}{\gamma^{2}\left(\frac{r}{R}\right)^{a+b}}\left((a-b)\leavevmode\nobreak\
\gamma\left(\frac{r}{R}\right)^{b}\right)e^{-\Lambda}\frac{W}{r^{2}}$
$\displaystyle+$ $\displaystyle\partial_{r}(r^{2}\delta
p)-\frac{1}{2}a(a-b)\gamma\left(\frac{r}{R}\right)^{b}r\delta
p+\frac{1}{2}\left(a^{2}+ab+2b\right)\gamma\left(\frac{r}{R}\right)^{b}r\delta
p\;=\;0\;$

and

$\frac{\omega^{2}}{\gamma^{2}\left(\frac{r}{R}\right)^{a}}\;\gamma\;\frac{V}{r^{2}}-\frac{1}{4}\left(a^{2}+ab+2b\right)\gamma\left(\frac{r}{R}\right)^{b}\delta
p\;=\;0\;.$ (57)

An inspection of Eq. (49) reveals the scaling $\;r^{2}\delta p\sim V/r^{2}\;$,
whereas Eq. (57) rather implies $\;r^{a+b}\delta p\sim V/r^{2}\;$. The
conclusion is that

$a+b\;=\;2\;.$ (58)

Similarly, by comparing the radial dependence of the term
$\partial_{r}(r^{2}\delta p)$ in Eq. (5) with any of the other terms involving
$\delta p$, one observes that $\;r\sim r^{b+1}\;$, so that

$\displaystyle a$ $\displaystyle=$ $\displaystyle 2\;,$ (59) $\displaystyle b$
$\displaystyle=$ $\displaystyle 0\;.$ (60)

These values for $a$ and $b$ simplify the dynamical equations and the
expression for $\delta p$ by a considerable amount. Respectively,

$\displaystyle-$ $\displaystyle
2\frac{\omega^{2}R^{2}}{\gamma^{2}}{\gamma}^{3/2}\frac{W}{r^{4}}+\partial_{r}(r^{2}\delta
p)\;=\;0\;,\;$ (61)
$\frac{\omega^{2}R^{2}}{\gamma^{2}}\frac{V}{r^{2}}-r^{2}\delta p\;=\;0\;,$
(62) $\displaystyle r^{2}\delta
p\;=\;2\gamma^{3/2}\left[\frac{W^{\prime}}{r^{2}}+\frac{W}{r^{3}}\right]+\leavevmode\nobreak\
\ell\left(\ell+1\right)\frac{V}{r^{2}}\;.$ (63)

Defining a dimensionless frequency $\widetilde{\omega}$,

$\omega^{2}\;=\;\gamma^{2}\tfrac{1}{R^{2}}\widetilde{\omega}^{2}$ (64)

and rescaling the perturbations as $\;\widetilde{W}=\sqrt{\gamma}W\;$,
$\gamma\widetilde{V}=V\;$ and $\;{\gamma}\widetilde{\delta p}=r^{2}{\delta
p}\;$, we then have

$\displaystyle\partial_{r}\widetilde{\delta
p}\;=\;2\widetilde{\omega}^{2}\frac{\widetilde{W}}{r^{4}}\;,\;$ (65)
$\widetilde{\delta
p}\;=\;\widetilde{\omega}^{2}\frac{\widetilde{V}}{r^{2}}\;,$ (66)
$\displaystyle\widetilde{\delta
p}\;=\;2\frac{\partial_{r}\left(r\widetilde{W}\right)}{r^{3}}+\leavevmode\nobreak\
\ell\left(\ell+1\right)\frac{\widetilde{V}}{r^{2}}\;.$ (67)

We conclude that the the oscillation frequencies are non-relativistic,
governed by a scale of $\gamma/R$ which is parametrically smaller than the
relativistic frequency scale; that is,
$\;\omega\sim\frac{\gamma}{R}\ll\frac{1}{R}\;$.

We can use Eqs. (65-67) to obtain a single “wave equation” for
$\widetilde{\delta p}$ (or, equivalently, for $\widetilde{W}$ or
$\widetilde{V}$),

$-\widetilde{\omega}^{2}\widetilde{\delta
p}+\frac{1}{r^{3}}\partial_{r}\left(r^{5}\partial_{r}\widetilde{\delta
p}\right)+\ell(\ell+1)\widetilde{\delta p}\;=\;0\;.$ (68)

To find the spectrum resulting from Eq. (68), we need to impose a new boundary
condition to replace the one at the center of the star. Recall that the
center-of-the-star condition is trivially satisfied in our case because of the
regularization process in the core. The new condition is set by the scaling of
the radial dependence of the perturbations modes $V_{\ell m}$ and $W_{\ell
m}$. This will be applied for $\;r\ll R\;$, but far enough from the center so
that $r$ falls within the bulk of the frozen star, outside of the regularized
core.

Let us first consider the expansion of $V_{\ell m}$ in Eqs.(13) and (14). As
in the multipole expansion of any such quantity,

$V_{\ell m}\;\sim\;V_{\ell}\;r^{\ell}\;,$ (69)

where the relevant expansion range is $\;\ell\geq 2\;$.

Now combining Eqs. (65) and (66), we obtain

$\partial_{r}\left(\frac{\widetilde{V}_{\ell
m}}{r^{2}}\right)\;=\;\frac{2}{r^{4}}\widetilde{W}_{\ell m}\;,$ (70)

from which it follows that

$\widetilde{W}_{\ell m}\;\sim\;\widetilde{W}_{\ell}\;r^{\ell+1}$ (71)

and then

$\widetilde{W}_{\ell m}\;=\;\frac{\ell-2}{2}\widetilde{V}_{\ell m}\;.$ (72)

Additionally, it follows that $\;\widetilde{\delta p}\sim r^{\ell-2}\;$.

Using these scaling relations, we observe that

$\frac{1}{r^{3}}\partial_{r}\left(r^{5}\partial_{r}\widetilde{\delta
p}\right)\;=\;(\ell+2)(\ell-2)\widetilde{\delta p}\;,$ (73)

and then conclude from this equation and Eq. (68) that

$\widetilde{\omega}^{2}\;=\;2\ell^{2}+\ell-4\;.$ (74)

Assuming a wavelength of order $R$, we can also deduce a speed of propagation
(squared) directly from the spectrum,
$\;v^{2}\sim\gamma^{2}(2\ell^{2}+\ell-4)\;$. However, from an internal (I) or
proper-time perspective, one factor of the redshift $\gamma$ should be
dispensed with, leaving $\;v_{I}^{2}\sim\gamma(2\ell^{2}+\ell-4)\;$. The
dispersion relation (74) is in agreement with that of [21], which considered
the quasinormal-mode problem for the closely related polymer model. There, the
square of the real part of the frequencies scaled with the coupling strength
of closed strings $g_{s}^{2}$; the smallest non-classical, dimensionless
parameter in the framework and the smallest dimensionless one besides
$\;\epsilon=l_{s}/R\;$. Now recall that $\gamma$ played precisely the same
role for the deformed frozen star, as $\varepsilon$ — the redshift of the
undeformed star — is the only smaller dimensionless parameter. Another way of
seeing that $\gamma$ and $g_{s}^{2}$ play the same role in their respective
models is that they both provide a direct measure of $\frac{\Delta R}{R}$.
From the polymer point of view, this relative scaling is already clear in
[21], but see also [31, 32] for further discussion.

The solutions for the perturbations can be obtain by using Eq. (73) to solve
for $\widetilde{\delta p}$,

$\widetilde{\delta p}_{\ell}\;=\;C_{1}\leavevmode\nobreak\
r^{l-2}+C_{2}\leavevmode\nobreak\ r^{-(l+2)}\;.$ (75)

The boundary condition at $\;r\ll R\;$ requires that $\;C_{2}=0\;$, so that

$\widetilde{\delta p}_{\ell}\;=\;C_{1}\leavevmode\nobreak\ r^{l-2}\;.$ (76)

Then, from Eq. (66), it follows that

$\widetilde{V}_{\ell}\;=\;C_{1}\leavevmode\nobreak\
\frac{1}{\widetilde{\omega}^{2}}\leavevmode\nobreak\ r^{l}\;,$ (77)

whereas Eq. (65) leads to

$\widetilde{W}_{\ell}=C_{1}\leavevmode\nobreak\
\frac{l-2}{2}\leavevmode\nobreak\
\frac{1}{\widetilde{\omega}^{2}}\leavevmode\nobreak\ r^{l+1}.$ (78)

Note that $\widetilde{W}_{\ell}$ vanishes at leading order for $\;\ell=2\;$.

An interesting feature of the spectrum in [21] was that the imaginary part of
the mode frequencies (i.e., the inverse of the damping time) scaled as
$v^{2}$, and so $\;{\Im}(\omega)\ll{\Re}(\omega)\;$. This behavior is not
obvious at this point in the current analysis. But this is because we have yet
to provide a mechanism for any of the modes to couple to and then escape into
the external spacetime, which is what is required for the modes to dissipate
energy and thus experience damping. However, it is easy to see why such a
scaling is an inevitable feature of any ultra-compact object that is using
non-relativistic fluid modes to source gravitational waves. As explained in
[26], from an external observer perspective, these modes have exceptionally
long wavelengths $\;\lambda\sim R/v\;$ when compared to the size $R$ of the
object, reducing the transmission cross-section through a surface of area
$R^{2}$ by a factor of $\;\lambda^{2}/R^{2}\sim v^{2}\;$. This cross-section
determines, in turn, the power loss $\;\frac{dE}{dt}\sim\lambda^{2}/R^{2}\sim
v^{2}\;$, which then sets the damping time scale as $\;1/\tau\sim v^{2}\;$,
and so $\;{\rm Im}(\omega)\sim v^{2}\;$.

##  6 Overview

We have calculated the oscillatory mode spectrum for a deformed version of the
frozen star model: the defrosted star. The deformations were a mathematical
necessity; otherwise, the ultra-stability of the frozen star geometry would
have doomed the perturbative modes from the get go. The deformations do,
however, make sense when one considers that gravitational waves are generally
emitted as the result of some out-of-equilibrium event like the merger of a
binary system. One expects the disturbed object to eventually settle back into
its original equilibrium state. It would be interesting to understand how this
process works in the case of a frozen star.

We have found a spectrum for the star that is in agreement with that of the
polymer model. In both cases, the square of the real part of the frequency
scales with a small dimensionless parameter; the deformation parameter
$\gamma$ in the case of the defrosted version of the frozen star and the
closed-string coupling strength $g_{s}^{2}$ for the collapsed polymer. These
are seemingly very different but have two important similarities: (1) They
both represent the second-smallest dimensionless parameter in their respective
frameworks and the smallest one that is also non-classical (effectively so, in
the case of the star) and (2) they both scale as $\frac{\Delta R}{R}$. We
expect that this is really the same parameter but different observers will
assign it a different meaning depending on their vantage point. For instance,
anyone who believes that BHs are purely classical, geometric objects would
have no room for string coupling in her “story”. This reasoning falls very
much in line with what Hawking called the “Principle of Ignorance” [33].

Missing in our analysis is a direct calculation of the damping time. This
could be done by changing the boundary conditions at the surface of the star
from a vanishing wave to a standing wave, which is then matched to an outgoing
wave in the exterior. This calculation is technically quite involved and
requires the incorporation of the metric perturbations, meaning that the
Cowling approximation can no longer be applied. We do hope to perform such a
calculation in the future. As explained in the main text, it can be expected
on physical grounds that the damping time scales as $1/\gamma^{2}$, a result
that was indeed substantiated in [21]. This is well worth confirming for the
defrosted star, as the implication is a particularly long lifetime for the
emitted modes and, with it, an excellent chance of detection. Moreover, if
$\gamma$ is indeed related to the coupling strength of closed strings, the
lifetime of the fluid modes provides an unexpected window into fundamental
string theory.

## Acknowledgments

We thank Daniela Doneva and Stoytcho Yazadjiev for clarifying comments on
their paper. The research is supported by the German Research Foundation
through a German-Israeli Project Cooperation (DIP) grant “Holography and the
Swampland” and by VATAT (Israel planning and budgeting committee) grant for
supporting theoretical high energy physics. The research of AJMM received
support from an NRF Evaluation and Rating Grant 119411 and a Rhodes
Discretionary Grant SD07/2022. AJMM thanks Ben Gurion University for their
hospitality during his visit.

##  Appendix A Perpendicular projection of the perturbation

Here, it is shown in some detail how to obtain Eq. (8), the perpendicular
projection of the covariant derivative of the perturbation of the stress
tensor (3).

Using the properties of the 4-velocity (see Eq. (6)) and that the Cowling
approximation is in effect, we have

$\displaystyle\mathcal{P}_{\alpha}^{\mu}\nabla_{\nu}\delta T_{\mu}^{\nu}\;=\;$
$\displaystyle\nabla_{\nu}\left(\delta\rho+\delta
q\right)\underset{=0}{\underbrace{\left(\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\right)u_{\mu}}}u^{\nu}$
$\displaystyle+\left(\delta\rho+\delta
q\right)\left(\left(\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\right)u^{\nu}\nabla_{\nu}u_{\mu}+\underset{=0}{\underbrace{\left(\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\right)u_{\mu}}}\overset{=0}{\overbrace{\nabla_{\nu}u^{\nu}}}\right)$
$\displaystyle+\left(\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\right)\left[\nabla_{\nu}\left(\rho+q\right)\left(\delta
u_{\mu}u^{\nu}+u_{\mu}\delta
u^{\nu}\right)+\left(\rho+q\right)\nabla_{\nu}\left(\delta
u_{\mu}u^{\nu}+u_{\mu}\delta u^{\nu}\right)\right]$
$\displaystyle+\left(\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\right)\nabla_{\nu}\delta
q\delta_{\mu}^{\nu}+\mathcal{P}_{\alpha}^{\mu}\nabla_{\nu}\delta\left(\sigma
k_{\mu}k^{\nu}\right)$ $\displaystyle=$ $\displaystyle\left(\delta\rho+\delta
q\right)\left(a_{\alpha}+u_{\alpha}u^{\nu}\underset{=0}{\underbrace{u^{\mu}\nabla_{\nu}u_{\mu}}}\right)$
$\displaystyle+\nabla_{\nu}\left(\rho+q\right)\left(\left(\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\right)\delta
u_{\mu}u^{\nu}+\delta
u^{\nu}\underset{=0}{\underbrace{\left(\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\right)u_{\mu}}}\right)$
$\displaystyle+\left(\rho+q\right)\left(\left(\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\right)\nabla_{\nu}\delta
u_{\mu}u^{\nu}+\left(\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\right)\nabla_{\nu}u_{\mu}\delta
u^{\nu}\right)$
$\displaystyle+\left(\rho+q\right)\left(\delta_{\alpha}^{\mu}+u^{\mu}u_{\alpha}\right)\left(\delta
u_{\mu}\overset{=0}{\overbrace{\nabla_{\nu}u^{\nu}}}+u_{\mu}\overset{=0}{\overbrace{\nabla_{\nu}\delta
u^{\nu}}}\right)$ $\displaystyle+\nabla_{\nu}\delta
q\delta_{\alpha}^{\nu}+u^{\nu}u_{\alpha}\nabla_{\nu}\delta
q+\mathcal{P}_{\alpha}^{\mu}\nabla_{\nu}\delta\left(\sigma
k_{\mu}k^{\nu}\right)$ $\displaystyle\;=\;$
$\displaystyle\delta\left(\rho+q\right)a_{\alpha}+\nabla_{\nu}\left(\rho+q\right)\left(\delta
u_{\alpha}u^{\nu}+u_{\alpha}\underset{=0}{\underbrace{u^{\mu}\delta
u_{\mu}}}u^{\nu}\right)$
$\displaystyle+\left(\rho+q\right)\left(u^{\nu}\nabla_{\nu}\delta
u_{\alpha}+u_{\alpha}u^{\mu}u^{\nu}\nabla_{\nu}\delta
u_{\mu}+u_{\alpha}u^{\mu}\nabla_{\nu}u_{\mu}\delta
u^{\nu}+\nabla_{\nu}u_{\alpha}\delta u^{\nu}\right)$
$\displaystyle+\nabla_{\nu}\delta
q\delta_{\alpha}^{\nu}+u^{\nu}u_{\alpha}\nabla_{\nu}\delta
q+\mathcal{P}_{\alpha}^{\mu}\nabla_{\nu}\delta\left(\sigma
k_{\mu}k^{\nu}\right)\;.$ (79)

Next, we implement our choice to work in the comoving frame, so that
$\;u^{t}u_{t}=-1\;$, and recall the time independence of any background
quantity. Then,

$\displaystyle 0\;=\;\mathcal{P}_{\alpha}^{\mu}\nabla_{\nu}\delta
T_{\mu}^{\nu}$
$\displaystyle=\delta\left(\rho+q\right)a_{\alpha}+\underset{=0}{\underbrace{\partial_{t}\left(\rho+q\right)}}\delta
u_{\alpha}u^{t}$
$\displaystyle+\left(\rho+q\right)\left(u^{\nu}\nabla_{\nu}\delta
u_{\alpha}-u^{\mu}\overset{=\delta_{\alpha}^{\nu}}{\overbrace{\delta_{\alpha}^{t}\delta_{t}^{\nu}}}\nabla_{\nu}\delta
u_{\mu}\underset{=0}{\underbrace{-\overset{=\delta_{\alpha}^{\mu}}{\overbrace{\delta_{t}^{\mu}\delta_{\alpha}^{t}}}\nabla_{\nu}u_{\mu}\delta
u^{\nu}+\nabla_{\nu}u_{\alpha}\delta u^{\nu}}}\right)$
$\displaystyle+\nabla_{\nu}\delta
q\delta_{\alpha}^{\nu}+u^{\nu}u_{\alpha}\nabla_{\nu}\delta
q+\mathcal{P}_{\alpha}^{\mu}\nabla_{\nu}\delta\left(\sigma
k_{\mu}k^{\nu}\right)$
$\displaystyle\;=\;\delta\left(\rho+q\right)a_{\alpha}+\left(\rho+q\right)u^{\nu}\left(\nabla_{\nu}\delta
u_{\alpha}-\nabla_{\alpha}\delta u_{\nu}\right)$
$\displaystyle+\nabla_{\nu}\delta
q\delta_{\alpha}^{\nu}+u^{\nu}u_{\alpha}\nabla_{\nu}\delta
q+\mathcal{P}_{\alpha}^{\mu}\nabla_{\nu}\delta\left(\sigma
k_{\mu}k^{\nu}\right)\;.$ (80)

## References

  * [1] R. Ruffini and J. A. Wheeler, “Introducing the black hole,” Phys. Today 24, no.1, 30 (1971)
  * [2] R. Brustein, A. J. M. Medved and T. Simhon, “Black holes as frozen stars,” Phys. Rev. D 105, no.2, 024019 (2022) [arXiv:2109.10017 [gr-qc]].
  * [3] R. Brustein, A. J. M. Medved, T. Shindelman and T. Simhon, “Black holes as frozen stars: Regular interior geometry,” [arXiv:2301.09712 [gr-qc]].
  * [4] R. Brustein and A. J. M. Medved, “Resisting collapse: How matter inside a black hole can withstand gravity,” Phys. Rev. D 99, no.6, 064019 (2019) [arXiv:1805.11667 [hep-th]].
  * [5] R. Brustein and A. J. M. Medved, “Non-Singular Black Holes Interiors Need Physics Beyond the Standard Model,” Fortsch. Phys. 67, no.10, 1900058 (2019) [arXiv:1902.07990 [hep-th]].
  * [6] R. Penrose, “Gravitational Collapse and Space-Time Singularities,” Phys. Rev. Lett. 14, 57 (1965).
  * [7] S. W. Hawking and R. Penrose, “The singularities of gravitational collapse and cosmology,” Proc. R. Soc. Lond. A 314, 529 (1970).
  * [8] H. Buchdahl, “General Relativistic Fluid Spheres,” Phys. Rev. 116, 1027 (1959).
  * [9] S. Chandrasekhar “Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity,” Phys. Rev. Lett. 12, 114 (1964).
  * [10] S. Chandrasekhar, “The Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity,” Astrophys. J. 140, 417 (1964).
  * [11] H. Bondi, “Massive spheres in general relativity,” Proc. Roy. Soc. Lond. A 282, 303 (1964).
  * [12] P. O. Mazur and E. Mottola, “Surface tension and negative pressure interior of a non-singular ,” Class. Quant. Grav. 32, no. 21, 215024 (2015) [arXiv:1501.03806 [gr-qc]].
  * [13] V. P. Frolov and A. Zelnikov, “Quantum radiation from an evaporating nonsingular black hole,” Phys. Rev. D 95, no. 12, 124028 (2017) [arXiv:1704.03043 [hep-th]].
  * [14] R. Carballo-Rubio, F. Di Filippo, S. Liberati, C. Pacilio and M. Visser, “On the viability of regular black holes,” JHEP 1807, 023 (2018) [arXiv:1805.02675 [gr-qc]].
  * [15] P. O. Mazur and E. Mottola, “Gravitational condensate stars: An alternative to black holes,” arXiv:gr-qc/0109035.
  * [16] C. Barcelo, S. Liberati, S. Sonego and M. Visser, “Fate of gravitational collapse in semiclassical gravity,” Phys. Rev. D 77, 044032 (2008) [arXiv:0712.1130 [gr-qc]].
  * [17] R. Carballo-Rubio, “Stellar equilibrium in semiclassical gravity,” Phys. Rev. Lett. 120, no. 6, 061102 (2018) [arXiv:1706.05379 [gr-qc]].
  * [18] R. Brustein and A. J. M. Medved, “Black holes as collapsed polymers,” Fortsch. Phys. 65, no. 1, 1600114 (2017) [arXiv:1602.07706 [hep-th]].
  * [19] R. Brustein and A. J. M. Medved, “Quantum state of the black hole interior,” JHEP 1508, 082 (2015) [arXiv:1505.07131 [hep-th]].
  * [20] R. Brustein and A. J. M. Medved, “Emergent horizon, Hawking radiation and chaos in the collapsed polymer model of a black hole,” Fortsch. Phys. 65, no.2, 1600116 (2017) [arXiv:1607.03721 [hep-th]].
  * [21] R. Brustein, A. J. M. Medved and K. Yagi, “When black holes collide: Probing the interior composition by the spectrum of ringdown modes and emitted gravitational waves,” Phys. Rev. D 96, no.6, 064033 (2017) [arXiv:1704.05789 [gr-qc]].
  * [22] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], “Obser- vation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116, no. 6, 061102 (2016) [arXiv:1602.03837 [gr-qc]].
  * [23] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], “Tests of general relativity with GW150914,” Phys. Rev. Lett. 116, no. 22, 221101 (2016) [arXiv:1602.03841 [gr-qc]].
  * [24] V. Cardoso and P. Pani, “Testing the nature of dark compact objects: a status report,” Living Rev. Rel. 22, no.1, 4 (2019) [arXiv:1904.05363 [gr-qc]].
  * [25] R. Brustein, A. J. M. Medved and K. Yagi, “Discovering the interior of black holes,” Phys. Rev. D 96, no.12, 124021 (2017) [arXiv:1701.07444 [gr-qc]].
  * [26] R. Brustein and A. J. M. Medved, “Quantum hair of black holes out of equilibrium,” Phys. Rev. D 97, no.4, 044035 (2018) [arXiv:1709.03566 [hep-th]].
  * [27] T.G. Cowling, “The non-radial oscillations of polytropic stars,” MNRAS 101, 367 (1941).
  * [28] K. D. Kokkotas and B. F. Schutz, “Normal modes of a model radiating system,” Gen. Rel. and Grav. 18, 913 (1986).
  * [29] D. D. Doneva and S. S. Yazadjiev, “Gravitational wave spectrum of anisotropic neutron stars in Cowling approximation,” Phys. Rev. D 85, 124023 (2012) [arXiv:1203.3963 [gr-qc]].
  * [30] K. S. Thorne, “Multipole Expansion of Gravitational Radiation,” Rev. Mod. Phys. 52, 299 (1980).
  * [31] R. Brustein and Y. Sherf, “Quantum Love numbers,” Phys. Rev. D 105, no.2, 024043 (2022) [arXiv:2008.02738 [gr-qc]].
  * [32] R. Brustein and Y. Sherf, “Classical Love number for quantum black holes,” Phys. Rev. D 105, no.2, 024044 (2022) [arXiv:2104.06013 [gr-qc]].
  * [33] S. W. Hawking, “Particle creation by black holes,” Comm. Math. Phys. 43, 199 (1975).


Oriented unicyclic graphs with extremal skew energy

Yaoping Hou, Xiaoling Shen and Chongyan Zhang

Department of Mathematics, Hunan Normal University

Changsha, Hunan 410081, China

email<EMAIL_ADDRESS>

Abstract

Let $\overrightarrow{G}$ be an oriented graph of order $n$ and
$\lambda_{1},\lambda_{2},\cdots,\lambda_{n}$ denote all the eigenvalues of the
skew-adjacency matrix of $\overrightarrow{G}.$ The skew energy
$\displaystyle{\cal E}_{s}(\overrightarrow{G})=\sum_{i=1}^{n}|\lambda_{i}|.$
In this paper, the oriented unicyclic graphs with minimal and maximal skew
energy are determined.

_Keywords: Unicyclic graph; oriented graph; skew-adjacency matrix, skew
energy_

AMS Classification 05C50, 15A03

##  1 Introduction

An important quantum–chemical characteristic of a conjugated molecule is its
total $\pi$-electron energy. The energy of a graph has closed links to
chemistry. Since the concept of the energy of simple undirected graphs was
introduced by Gutman in [5], there have been lots of research papers on this
topic. For the extremal energy of unicyclic graph, Hou [9] showed that
$S^{3}_{n}$ is the graph with minimal energy in all unicyclic graphs; In [7],
Huo and Li showed that $P^{6}_{n}$ is the graph with maximal energy in all
unicyclic graphs. More results on the energy of unicyclic graphs see [6, 8,
10, 12, 15, 16].

There are various generalizations of energy to graph matrices [13], such as
Laplacian energy, incidence energy and distance energy. Let $G$ be a simple
undirected finite graph of order $n$ and $\overrightarrow{G}$ be an
orientation of $G,$ which assigns to each edge of $G$ a direction so that
$\overrightarrow{G}$ becomes a directed graph. All digraphs in this paper are
oriented graphs of some graphs. Let $\overrightarrow{G}$ be an oriented graph
of order $n$. The skew adjacency matrix $S(\overrightarrow{G})=(s_{i,j})$ is a
real skew symmetric matrix, where $s_{i,j}=1$ and $s_{j,i}=-1$ if
$i\rightarrow j$ is an arc of $\overrightarrow{G},$ otherwise
$s_{i,j}=s_{j,i}=0.$ The skew spectrum $Sp(\overrightarrow{G})$ of
$\overrightarrow{G}$ is defined as the spectrum of $S(\overrightarrow{G}).$
Note that $Sp(\overrightarrow{G})$ consists of only purely imaginary
eigenvalues because $S(\overrightarrow{G})$ is real skew symmetric.

Recently, the skew-energy of an oriented graph $\overrightarrow{G}$ was
defined as the energy of matrix $S(\overrightarrow{G})$ in [2], that is,

${\cal E}_{s}(\overrightarrow{G})=\sum_{\lambda\in
Sp(\overrightarrow{G})}|\lambda|.$

There are situations when chemists use digraphs rather than graphs. One such
situation is when vertices represent distinct chemical species and arcs
represent the direction in which a particular reaction takes place between the
two corresponding species. It is possible that the skew energy has similar
applications as energy in chemistry. For a graph $G$, there are any
orientations on it, it is also interesting to find what orientation has
extremal energy among all orientations of a given graph.

An unicyclic graph is the connected graph with the same number of vertices and
edges. In this paper, we are interested in studying the orientations of
unicyclic graphs with extremal skew energy. Let $G(n,\ell)$ be the set of all
connected unicyclic graphs on $n$ vertices with girth $\ell$. Denote, as
usual, the $n$-vertex path and cycle by $P_{n}$ and $C_{n}$, respectively. Let
$P^{\ell}_{n}$ be the unicyclic graph obtained by connecting a vertex of
$C_{\ell}$ with a terminal vertex of $P_{n-\ell}$, $S_{n}^{\ell}$ be the graph
obtained by connecting $n-\ell$ pendant vertices to a vertex of $C_{\ell}$
(see Fig. $1$).

The rest of the paper is organized as follows. In section $2$, a new integral
formal for ${\cal E}_{s}(\overrightarrow{G})$ is obtained and the oriented
unicyclic graph with minimal skew energy is determined. In section $3$, the
oriented unicyclic graph with maximal skew energy is determined.

##  2 Oriented unicyclic graph with minimal skew energy

Let $G$ be a graph. A _linear subgraph_ $L$ of $G$ is a disjoint union of some
edges and some cycles in $G$ ([3]). A _$k$ -matching_ $M$ in $G$ is a disjoint
union of $k$-edges. If $2k$ is the order of $G$, then a $k$-matching of $G$ is
called a _perfect matching_ of $G$. The number of $k$-matching is denoted by
$m({G},k)$.

If $C$ be any undirected even cycle of ${G},$ we say $C$ is _evenly oriented_
relative to the orientation $\overrightarrow{G}$ if it has an even number of
edges oriented in clockwise direction. Otherwise $C$ is _oddly oriented_.

We call a linear subgraph $L$ of $G$ _evenly linear_ if $L$ contains no odd
cycle and denote by $\mathcal{E}\mathcal{L}_{i}(G)$ (or
$\mathcal{E}\mathcal{L}_{i}$ for short) the set of all evenly linear subgraphs
of $G$ with $i$ vertices. For a linear subgraph
$L\in\mathcal{E}\mathcal{L}_{i}$ denote by $p_{e}(L)$(resp., $p_{o}(L)$ ) the
number of evenly (resp., oddly) oriented cycles in $L$ relative to
$\overrightarrow{G}$. Denote the characteristic polynomial of
$S(\overrightarrow{G})$ by

$P_{s}(\overrightarrow{G};x)=det(xI-S(\overrightarrow{G}))=\sum^{n}_{i=0}b_{i}x^{n-i}.$
(1)

Then (i) $b_{0}=1$, (ii) $b_{2}$ is the number of edges of $G,$ (iii) all
$b_{i}\geq 0$ and (iv) $b_{i}=0$ for all odd $i$ since the determinant of
every real skew symmetric matrix is nonnegative and is 0 if its order is odd.

We have

###### Lemma 2.1

([11]) Let $\overrightarrow{G}$ be an orientation of a graph $G$. Then

$b_{i}(\overrightarrow{G})=\sum_{L\in\mathcal{E}\mathcal{L}_{i}}(-2)^{p_{e}(L)}2^{p_{o}(L)},$
(2)

where $p_{e}(L)$ is the number of evenly oriented cycles of $L$ and $p_{o}(S)$
is the number of oddly oriented cycles of $L$ relative to
$\overrightarrow{G}$, respectively.

From Lemma 2.1, we can get the following statement.

###### Corollary 2.2

Let $G$ be a unicyclic graph with unique cycle $C$ of length $\ell$ and
$\overrightarrow{G}$ be an orientation of $G.$ Then

$b_{2k}(\overrightarrow{G})=\left\\{\begin{array}[]{cc}m({G},k),&C\ is\ an\
odd\ cycle;\\\ m({G},k)-2m({G}-C,k-\frac{\ell}{2}),&C\ is\ evenly\
oriented;\\\ m({G},k)+2m({G}-C,k-\frac{\ell}{2}),&C\ \ is\ oddly\
oriented.\end{array}\right.$

###### Lemma 2.3

([11]) Let $e=uv$ be an edge of $G$ that is on no any even cycle of $G.$ Then

$P_{s}(\overrightarrow{G};x)=P_{s}(\overrightarrow{G}-e;x)+P_{s}(\overrightarrow{G}-u-v;x).$

By Lemma 2.3, we obtain that

###### Corollary 2.4

Let $e=uv$ be an edge of $G$ that is on no even cycle of $G.$ Then

$b_{2k}(\overrightarrow{G})=b_{2k}(\overrightarrow{G}-e)+b_{2k-2}(\overrightarrow{G}-u-v).$
(3)

Furthermore, if $e=uv$ is a pendant edge with the pendant vertex $v.$ Then

$b_{2k}(\overrightarrow{G})=b_{2k}(\overrightarrow{G}-v)+b_{2k-2}(\overrightarrow{G}-u-v).$
(4)

For any orientation of a graph which does not contain any even cycle (in
particular, a tree, a unicyclic non-bipartite graph) by Theorem 2.9, we have
$b_{2k}(\overrightarrow{G})=m({G},k)$.

For the $k$-matching number of a graph $G$, the following result is well-
known.

###### Lemma 2.5

Let $e=uv$ be an edge of $G$. Then
(i) $m({G},k)=m({G}-e,k)+m({G}-u-v,k-1).$
(ii) If $G$ is a forest, then $m({G},k)\leq m(P_{n},k)$, $k\geq 1$.
(iii) If ${H}$ is a subgraph of ${G}$, then $m({H},k)\leq m({G},k)$, $k\geq
1$. Moreover, if ${H}$ is a proper subgraph of ${G}$, then the inequality is
strict.

From [2], an integral formula for skew energy was given:

${\cal
E}_{s}(\overrightarrow{G})=\frac{1}{\pi}\int^{\infty}_{-\infty}(n+x\frac{P^{\prime}_{s}(\overrightarrow{G},-x)}{P_{s}(\overrightarrow{G},-x)})dx,$
(5)

where $P^{\prime}_{s}(\overrightarrow{G},-x)$ is the derivative of
$P_{s}(\overrightarrow{G},-x).$

###### Theorem 2.6

Let $\overrightarrow{G}$ be an orientation of a graph $G$. Then

${\cal
E}_{s}(\overrightarrow{G})=\frac{1}{\pi}\int^{\infty}_{-\infty}\frac{1}{t^{2}}ln(1+\sum_{k=1}^{\lfloor\frac{n}{2}\rfloor}b_{2k}t^{2k})dt.$

Proof. From equality (5), we have

$\displaystyle{\cal E}_{s}(\overrightarrow{G})$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int^{\infty}_{-\infty}(n+x\frac{P^{\prime}_{s}(\overrightarrow{G},-x)}{P_{s}(\overrightarrow{G},-x)})dx=\frac{1}{\pi}\int^{\infty}_{-\infty}(n-\mu\frac{P^{\prime}_{s}(\mu)}{P_{s}(\mu)})d\mu$
$\displaystyle=$
$\displaystyle\frac{1}{\pi}\int^{\infty}_{0}(n-\frac{1}{t}\frac{P^{\prime}_{s}(\frac{1}{t})}{P_{s}(\frac{1}{t})})(-\frac{1}{t^{2}})dt+\frac{1}{\pi}\int^{0}_{-\infty}(n-\frac{1}{t}\frac{P^{\prime}_{s}(\frac{1}{t})}{P_{s}(\frac{1}{t})})(-\frac{1}{t^{2}})dt$
$\displaystyle=$
$\displaystyle\frac{1}{\pi}\int^{\infty}_{-\infty}(n-\frac{1}{t}\frac{P^{\prime}_{s}(\frac{1}{t})}{P_{s}(\frac{1}{t})})\frac{1}{t^{2}}dt=\frac{1}{\pi}\int^{\infty}_{-\infty}\frac{1}{t}d[ln(t^{n}P_{s}(\frac{1}{t})]$
$\displaystyle=$
$\displaystyle\frac{1}{\pi}\int^{\infty}_{-\infty}\frac{1}{t^{2}}ln(t^{n}P_{s}(\frac{1}{t}))dt.=\frac{1}{\pi}\int^{\infty}_{-\infty}\frac{1}{t^{2}}ln(1+\sum_{k=1}^{\lfloor\frac{n}{2}\rfloor}b_{2k}t^{2k})dt.\;\;\;\Box$

###### Example 2.7

Let $\overrightarrow{C_{4}}^{-}$ be an orientation of $C_{4}$ such that all
edges have the same direction. Then by Lemma 2.1, the characteristic
polynomial of $S(\overrightarrow{C_{4}}^{-})$ is $x^{4}+4x^{2}$. By Lemma 2.6,

$\displaystyle{\cal E}_{s}(\overrightarrow{C_{4}}^{-},x))$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int^{\infty}_{-\infty}\frac{1}{t^{2}}ln(1+4t^{2})dt$
$\displaystyle=$
$\displaystyle\frac{8}{\pi}\int^{\infty}_{0}\frac{1}{1+4t^{2}}dt=4.$

From Theorem 2.6, ${\cal E}_{s}(\overrightarrow{G})$ is an increasing function
of $b_{2k}(\overrightarrow{G})$, $k=0$, $1$, $\cdots$,
$\lfloor\frac{n}{2}\rfloor$. Consequently, if $\overrightarrow{G_{1}}$ and
$\overrightarrow{G_{2}}$ are oriented graphs of $G_{1}$ and $G_{2},$
respectively, for which

$b_{2k}(\overrightarrow{G_{1}})\geq b_{2k}(\overrightarrow{G_{2}})$ (6)

for all $\lfloor{\frac{n}{2}}\rfloor\geq k\geq 0$, then

${\cal E}_{s}(\overrightarrow{G_{1}})\geq{\cal
E}_{s}(\overrightarrow{G_{2}}).$ (7)

Equality in (7) is attained only if (6) is an equality for all
$\lfloor{\frac{n}{2}}\rfloor\geq k\geq 0$.

If relations (6) hold for all $k$, then we write $G_{1}\succeq G_{2}$ or
$G_{2}\preceq G_{1}$. If $G_{1}\succeq G_{2}$, but not $G_{2}\succeq G_{1}$,
then we write $G_{1}\succ G_{2}$.

Let $\overrightarrow{G}$ be an orientation of a graph $G$. Let $W$ be a subset
of $V(G)$ and $\overline{W}=V(G)\setminus W$. The orientation
$\overrightarrow{G^{\prime}}$ of $G$ obtained from $\overrightarrow{G}$ by
reversing the direction of all arcs between $\overline{W}$ and $W$ is said to
be obtained from $\overrightarrow{G}$ by switching with respect to $W$.
Moreover, two orientations $\overrightarrow{G}$ and
$\overrightarrow{G^{\prime}}$ of a graph $G$ are said to be switching-
equivalent if $\overrightarrow{G^{\prime}}$ can be obtained from
$\overrightarrow{G}$ by a sequences of switching. If two orientations
$\overrightarrow{G}$ and $\overrightarrow{G}^{\prime}$ of a graph $G$ are
switching-equivalent then their skew-adjacency matrices are similar by a
diagonal matrix whose $(i,i)$-entry is $-1$ when $i\in W$ and $1$ when
$i\not\in W$ by Lemma 3.1 in [2] and hence they have the same skew spectrum.
Thus,

###### Lemma 2.8

Let $\overrightarrow{G}$ and $\overrightarrow{G}^{\prime}$ be two orientations
of a graph $G.$ If $\overrightarrow{G}$ and $\overrightarrow{G}^{\prime}$ are
switching-equivalent, then ${\cal E}_{s}(\overrightarrow{G})={\cal
E}_{s}(\overrightarrow{G}^{\prime})$.

By Lemma 2.8 and switching-equivalence, there are only two different
orientations on a unicyclic graph $G.$ All edges on the unique cycle $C$ have
the same direction or just one edge on the cycle has the opposite direction to
the directions of other edges on the cycle regardless how the edges not on the
cycle $C$ are oriented. Denote by ${\overrightarrow{G}}^{-}$
(${\overrightarrow{G}}^{+},$ resp.) the orientation of $G$ in first (second,
resp.) case above.

###### Lemma 2.9

For any unicyclic graph $G$,
${\overrightarrow{G}}^{+}\succeq{\overrightarrow{G}}^{-}.$

Proof. By Corollary (2.2), if the girth $\ell$ of $G$ is odd, then

$b_{2k}({{\overrightarrow{G}}^{+}})=b_{2k}({{\overrightarrow{G}}^{-}})=m(G,k)\mbox{
for all $0\leq k\leq\lfloor\frac{n}{2}\rfloor$.}$

If the girth $\ell$ of $G$ is even, then

$b_{2k}({{\overrightarrow{G}}^{+}})=m(G,k)+2m(G-C_{\ell},k-\lfloor\frac{\ell}{2}\rfloor),$
$b_{2k}({{\overrightarrow{G}}^{-}})=m(G,k)-2m({G}-C_{\ell},k-\lfloor\frac{\ell}{2}\rfloor)$

for all nonnegative integer $\lfloor\frac{n}{2}\rfloor\geq
k\geq\lfloor\frac{\ell}{2}\rfloor,$ and

$b_{0}({{\overrightarrow{G}}^{+}})=b_{0}({{\overrightarrow{G}}^{-}})=1,\;\;b_{2k}({{\overrightarrow{G}}^{+}})=b_{2k}({{\overrightarrow{G}}^{-}})=m(G,k)$

for all $0<k<\lfloor\frac{\ell}{2}\rfloor$. Thus the result follows
immediately. $\Box$

###### Lemma 2.10

Let ${\overrightarrow{G}}$ be an orientation of a unicyclic graph $G\in
G(n,\ell)$, $G\neq S_{n}^{\ell}$. If unique cycle $C_{\ell}$ in
${\overrightarrow{G}}$ and $\overrightarrow{S_{n}^{\ell}}$ is the same
orientation, then ${\overrightarrow{G}}\succ\overrightarrow{S_{n}^{\ell}}$.

Proof. We prove the statement by induction on $n$. Since $G\neq S_{n}^{\ell}$,
$n\geq\ell+2$. For $n=\ell+2$, $G$ is one of the two graphs in Fig. $2$.

By Corollary 4,

$b_{2k}(\overrightarrow{P_{\ell+2}^{\ell}})=b_{2k}(\overrightarrow{P_{\ell+1}^{\ell}})+b_{2k-2}(\overrightarrow{C_{\ell}}),$

$b_{2k}({\overrightarrow{G}_{1}})=b_{2k}(\overrightarrow{P_{\ell+1}^{\ell}})+b_{2k-2}(\overrightarrow{T}),$

where $T$ is a graph obtained by connecting one of the vertices in
$P_{\ell-1}$ to a pendant vertex.

$b_{2k}(\overrightarrow{S_{\ell+2}^{\ell}})=b_{2k}(\overrightarrow{P_{\ell+1}^{\ell}})+b_{2k-2}(\overrightarrow{P_{\ell-1}}).$

By Corollary 4, if $\ell$ is odd or $\ell$ is even but $k\leq\frac{\ell}{2}$
then

$b_{2k-2}(\overrightarrow{C_{\ell}})=m({C_{\ell}},k-1);$

for $k=\frac{\ell}{2}+1.$ If $C$ is oddly oriented, then

$b_{l}(\overrightarrow{C_{\ell}})=m({C_{\ell}},\frac{\ell}{2})+2=4;$

If $C$ is evenly oriented, then

$b_{\ell}(\overrightarrow{C_{\ell}})=m({C_{\ell}},\frac{\ell}{2})-2=0,\;\;b_{\ell}(\overrightarrow{P_{\ell-1}})=0.$

Since both $T$ and $P_{\ell-1}$ are trees,

$b_{2k-2}(\overrightarrow{T})=m({T},k-1),$
$b_{2k-2}(\overrightarrow{P_{\ell-1}})=m({P_{\ell-1}},k-1).$

Since $P_{\ell-1}$ is a proper subgraph of both $C_{\ell}$ and $T$, by Lemma
2.5,

$m(\overrightarrow{P_{\ell-1}},k-1)<m({T},k-1),m({P_{\ell-1}},k-1)<m({C_{\ell}},k-1).$

The result holds immediately for $n=\ell+2$.

Suppose that ${\overrightarrow{G}}\succ\overrightarrow{S_{n^{\prime}}^{\ell}}$
for all $n^{\prime}<n$. Since ${\overrightarrow{G}}$ is a unicyclic digraph,
there is at least a pendant edge $uv$ with pendant vertex $v$ in
${\overrightarrow{G}}$, by equality (4), we get

$\displaystyle b_{2k}({\overrightarrow{G}})$ $\displaystyle=$ $\displaystyle
b_{2k}({\overrightarrow{G}}-v)+b_{2k-2}({\overrightarrow{G}}-u-v),$
$\displaystyle b_{2k}(\overrightarrow{S_{n}^{\ell}})$ $\displaystyle=$
$\displaystyle
b_{2k}(\overrightarrow{S_{n-1}^{\ell}})+b_{2k-2}(\overrightarrow{P_{\ell-1}}).$

By induction assumption, it suffices to prove that
$b_{2k-2}({\overrightarrow{G}}-u-v)\geq b_{2k-2}(\overrightarrow{P_{\ell-1}})$
for all $0\leq k\leq\lfloor\frac{n-\ell}{2}\rfloor$. For
$k>\lfloor\frac{\ell-1}{2}\rfloor$, we have
$b_{2k-2}({\overrightarrow{G}}-u-v)\geq
b_{2k-2}(\overrightarrow{P_{\ell-1}})=0$. For
$\lfloor\frac{\ell-1}{2}\rfloor\geq k\geq 0$,
$b_{2k-2}({\overrightarrow{G}}-u-v)=m({\overrightarrow{G}}-u-v,k-1)\geq
m(\overrightarrow{P_{\ell-1}},k-1)=b_{2k-2}(\overrightarrow{P_{\ell-1}})$
since ${P_{\ell-1}}$ is a subgraph of ${G}-u-v$. $\Box$

###### Lemma 2.11

Let $n\geq\ell\geq 6$ or $n>\ell=5$, then
${\overrightarrow{S_{n}^{4}}}^{-}\prec{\overrightarrow{S_{n}^{4}}}^{+}\prec{\overrightarrow{S_{n}^{\ell}}}^{-}\preceq{\overrightarrow{S_{n}^{\ell}}}^{+}$.

Proof. By Lemma 2.1, the characteristic polynomial of
$\overrightarrow{S_{n}}^{4}$ is:

$\displaystyle P_{s}({\overrightarrow{S_{n}^{4}}}^{+},x)$ $\displaystyle=$
$\displaystyle x^{n-4}(x^{4}+nx^{2}+2n-4);$ $\displaystyle
P_{s}({\overrightarrow{S_{n}^{4}}}^{-},x)$ $\displaystyle=$ $\displaystyle
x^{n-4}(x^{4}+nx^{2}+2n-8);$

Obviously,
${\overrightarrow{S_{n}^{4}}}^{-}\prec{\overrightarrow{S_{n}^{4}}}^{+}$.

So it suffices to prove that $b_{4}(\overrightarrow{S_{n}^{\ell}})>2n-4$. By
equality (2),

$\displaystyle b_{4}({\overrightarrow{S_{n}^{\ell}}}^{+})$ $\displaystyle=$
$\displaystyle b_{4}({\overrightarrow{S_{n}^{\ell}}}^{-})$ $\displaystyle=$
$\displaystyle m(P_{\ell-1},2)+(n-\ell)m(P_{\ell-1},1)+2m(P_{\ell-2},1)$
$\displaystyle=$
$\displaystyle\frac{(\ell-3)(\ell-4)}{2}+(n-\ell)(\ell-2)+2(\ell-3)$
$\displaystyle=$ $\displaystyle\frac{2n\ell-\ell^{2}+\ell-4n}{2}.$

For $n>\ell=5$, $b_{4}(\overrightarrow{S_{n}^{\ell}})=3n-10\geq 2n-4$; For
$n\geq\ell=6$, $b_{4}(\overrightarrow{S_{n}^{\ell}})=4n-15>2n-4$; for
$n\geq\ell=7$, $b_{4}(\overrightarrow{S_{n}^{\ell}})=5n-21>2n-4$; for
$\ell\geq 8$,
$b_{4}(\overrightarrow{S_{n}^{\ell}})\geq\frac{n\ell+\ell-4n}{2}\geq
2n+4>2n-4.$

By Lemma 2.9, the proof is completed. $\Box$

For $n=\ell=5$,
$b_{4}({\overrightarrow{S_{n}^{\ell}}}^{+})=\frac{2n{\ell}-{\ell}^{2}+\ell-4n}{2}=5$,
$b_{4}({\overrightarrow{S_{n}^{4}}}^{+})=10-4=6$,
$b_{4}({\overrightarrow{S_{n}^{4}}}^{-})=10-8=2$. It gives
${\overrightarrow{S_{5}^{4}}}^{-}\prec{\overrightarrow{S_{5}^{5}}}^{+}\prec{\overrightarrow{S_{5}^{4}}}^{+}$.

###### Lemma 2.12

Let $\overrightarrow{S_{n}^{3}}$ be any orientation of unicyclic graph
$S_{n}^{3}$. Then $\overrightarrow{S_{n}^{3}}\prec\overrightarrow{S_{n}^{4}}$
for $n\geq 6$;
$\overrightarrow{S_{n}^{3}}={\overrightarrow{S_{n}^{4}}}^{-}\prec{\overrightarrow{S_{n}^{4}}}^{+}$
for $n=5$;
${\overrightarrow{S_{n}^{4}}}^{-}\prec\overrightarrow{S_{n}^{3}}\prec{\overrightarrow{S_{n}^{4}}}^{+}$
for $n=4$.

Proof. By Lemma 2.1, the characteristic polynomial of
$\overrightarrow{S_{n}^{3}}$ is

$P_{s}(\overrightarrow{S_{n}^{3}},x)=x^{n-4}(x^{4}+nx^{2}+n-3).$

Since $n-3<2n-8$, the result holds by Lemmas 2.11 and 2.9.

For $n=5$, it is easy to get that for any orientation ${\cal
E}_{s}(\overrightarrow{S_{5}^{3}})=2\sqrt{5+2\sqrt{2}}$, ${\cal
E}_{s}({\overrightarrow{S_{5}^{4}}}^{-})=2\sqrt{5+2\sqrt{2}}<2(\sqrt{2}+\sqrt{3})$,
${\cal E}_{s}({\overrightarrow{S_{5}^{4}}}^{+})=2(\sqrt{2}+\sqrt{3})$. Then
$\overrightarrow{S_{n}^{3}}={\overrightarrow{S_{5}^{4}}}^{-}\prec{\overrightarrow{S_{5}^{4}}}^{+}$.

For $n=4$, ${\cal E}_{s}(\overrightarrow{S_{5}^{3}})=2\sqrt{6}$, ${\cal
E}_{s}({\overrightarrow{C_{4}}}^{+})=4\sqrt{2}$, ${\cal
E}_{s}({\overrightarrow{C_{4}}}^{-})=4$. Thus
${\overrightarrow{S_{5}^{4}}}^{-}\prec\overrightarrow{S_{n}^{3}}\prec{\overrightarrow{S_{5}^{4}}}^{+}$.
$\Box$

From Lemmas 2.10, 2.11 and 2.12, we obtain one of the main result of this
paper.

###### Theorem 2.13

Among all orientations of unicyclic graphs on $n$ vertices,
$\overrightarrow{S_{n}^{3}}$ has the minimal skew energy and
${\overrightarrow{S_{n}^{4}}}^{-}$ has the second minimal skew energy for
$n\geq 6$; both ${\overrightarrow{S_{5}^{3}}}$ and
${\overrightarrow{S_{5}^{4}}}^{-}$ have the minimal skew energy,
${\overrightarrow{S_{5}^{4}}}^{+}$ has the second minimal skew energy for
$n=5$; ${\overrightarrow{C^{4}}}^{-}$ has the minimal skew energy,
$\overrightarrow{S_{4}^{3}}$ has the second minimal skew energy for $n=4$.

##  3  Oriented unicyclic graph with maximal skew energy

By Lemma 2.9, we only need to consider ${{\overrightarrow{G}}^{+}}$ for
considering of maximum skew energy.

###### Lemma 3.1

Let ${\overrightarrow{G}}$ be an orientation of unicyclic graph $G\in
G(n,\ell)$ and $G\neq P_{n}^{\ell}$. Then
${{\overrightarrow{G}}^{+}}\prec{\overrightarrow{P_{n}^{\ell}}}^{+}$.

Proof. We prove the statement by induction on $n$. For $n=\ell+2$, there are
only two cases for $G\neq P_{n}^{l}$: one is $S_{\ell+2}^{l}$, the other is
the graph $G_{1}$ in Fig. $2$. By the proof in Lemma 2.10, we only need to
prove that $b_{2k-2}(\overrightarrow{C_{l}})\geq m({T},k-1)$. By Lemma 2.5,

$\displaystyle m({\overrightarrow{C_{\ell}}}^{+},k-1)$ $\displaystyle=$
$\displaystyle m({C_{\ell}},k-1)$ $\displaystyle=$ $\displaystyle
m({P_{\ell}},k-1)+m({P_{\ell-2}},k-2),$ $\displaystyle m({T},k-1)$
$\displaystyle=$ $\displaystyle
m({P_{\ell-1}},k-1)+m({P_{s}}\bigcup{P_{t}},k-2)$ $\displaystyle\leq$
$\displaystyle m({P_{\ell-1}},k-1)+m({P_{\ell-2}},k-2)$ $\displaystyle<$
$\displaystyle m({P_{\ell}},k-1)+m({P_{\ell-2}},k-2)$ $\displaystyle=$
$\displaystyle m({{C_{\ell}}},k-1)$ $\displaystyle\leq$ $\displaystyle
b_{2k-2}({\overrightarrow{C_{\ell}}}^{+}),\;\;\;\;\;\;\;\mbox{(}by\ Corollary\
\ \ref{bre})$

where $s+t=l-2$.

Suppose
${{\overrightarrow{G}}^{+}}\prec{\overrightarrow{P_{n^{\prime}}^{\ell}}}^{+}$
for all $n^{\prime}<n$. Then there is a pendant edge, say $uv$, with pendant
vertex $v$. By Lemma 2.3,

$\displaystyle b_{2k}({{\overrightarrow{G}}^{+}})$ $\displaystyle=$
$\displaystyle
b_{2k}({{\overrightarrow{G}}^{+}}-v)+b_{2k-2}({{\overrightarrow{G}}^{+}}-v-u);$
$\displaystyle b_{2k}({\overrightarrow{P_{n}^{\ell}}}^{+})$ $\displaystyle=$
$\displaystyle
b_{2k}({\overrightarrow{P_{n-1}^{\ell}}}^{+})+b_{2k-2}({\overrightarrow{P_{n-2}^{\ell}}}^{+}).$

By induction hypothesis, it suffices to prove that
$b_{2k-2}({{\overrightarrow{G}}^{+}}-v-u)\leq
b_{2k-2}({\overrightarrow{P_{n-2}^{\ell}}}^{+})$. If
${{\overrightarrow{G}}^{+}}-v-u$ contains a cycle, then by induction
hypothesis the inequality holds. Suppose that ${{\overrightarrow{G}}^{+}}-v-u$
is a forest, then by Lemmas 2.5 and 2.9,

$\displaystyle b_{2k-2}({{\overrightarrow{G}}^{+}}-v-u)$ $\displaystyle=$
$\displaystyle m({{G}}-v-u,k-1)$ $\displaystyle\leq$ $\displaystyle
m({P_{n-2}},k-1)\leq m({{P_{n-2}^{\ell}}},k-1)$ $\displaystyle\leq$
$\displaystyle b_{2k-2}({{P_{n-2}^{\ell}}}).\mbox{ \hskip 28.45274pt $\Box$ }$

###### Lemma 3.2

For $3\leq\ell\leq n$, $\ell\neq 4$,
${\overrightarrow{P_{n}^{\ell}}}\prec{\overrightarrow{P_{n}^{4}}}^{+}$.

Proof. We proceed the proof by induction on $n$. For $n=\ell=5$,
$b_{4}({\overrightarrow{P_{5}^{5}}}^{+})=b_{4}({\overrightarrow{C_{5}}}^{+})=5$,
$b_{4}({\overrightarrow{P_{5}^{4}}}^{+})=b_{4}({\overrightarrow{S_{5}^{4}}}^{+})=6$,
thus
$b_{4}({\overrightarrow{P_{5}^{5}}}^{+})<b_{4}({\overrightarrow{P_{5}^{4}}}^{+})$.
Suppose that
${\overrightarrow{P_{n^{\prime}}^{\ell}}}^{+}\prec{\overrightarrow{P_{n^{\prime}}^{4}}}^{+}$
for all $n^{\prime}<n$, $5\leq\ell$. By Lemma 2.1, we have

$\displaystyle b_{2k}({\overrightarrow{P_{n}^{\ell}}}^{+})$ $\displaystyle=$
$\displaystyle
b_{2k}({\overrightarrow{P_{n-1}^{\ell}}}^{+})+b_{2k-2}({\overrightarrow{P_{n-2}^{\ell}}}^{+}),$
$\displaystyle b_{2k}({\overrightarrow{P_{n}^{4}}}^{+})$ $\displaystyle=$
$\displaystyle
b_{2k}({\overrightarrow{P_{n-1}^{4}}}^{+})+b_{2k-2}({\overrightarrow{P_{n-2}^{4}}}^{+}).$

By induction hypothesis, the result follows.

Similarly, we can prove that for $\ell=3$ the inequality holds too. $\Box$

By Lemmas 3.1 and 3.2, we have

###### Theorem 3.3

Among all orientations of unicyclic graph, ${\cal
E}_{s}({\overrightarrow{P_{n}^{4}}}^{+})$ is the unique oriented graph (under
switching-equivalent) with maximal skew energy.

Acknowledgments. This project is supported by National Natural Science
Foundation of China

## References

  * [1]
  * [2] C. Adiga, R. Balakrishnan and Wasin So, The skew energy of a digraph, Linear Algebra and its Application 432 (2010) 1825-1835.
  * [3] D. M. Cvetkovic, M. Doob, H. Sachs, Spectra of Graphs, Academic Press, New York, 1979.
  * [4] M. Dehmer, F. Emmert-Streib, Analysis of complex networks: from biology to linguistics. WILEY-VCH Verlag GmbH Co. hGaA, Weinheim, 2009.
  * [5] I. Gutman, The energy of a graph, Ber. Math. Stat. Sekt. Forschungszentrum Graz 103 (1978) 1–22.
  * [6] I. Gutman, Y. Hou, Bipartite unicyclic graphs with greatest energy, MATCH Commun. Math. Comput. Chem. 43 (2001) 17–28.
  * [7] B. Huo, X. Li, Complete solution to a conjecture on the maximal energy of unicyclic graphs. European journal of combinatorics. 32 (2011) 662–673.
  * [8] H. Hua, M. Wang, Unicyclic graphs with given number of pendent vertices and minimal energy, Lin. Algebra Appl. 426 (2007) 478–489.
  * [9] Y. Hou, Unicyclic graphs with minimal energy, J. Math. Chem. 29 (2001) 163–168.
  * [10] Y. Hou, I. Gutman, C. W. Woo, Unicyclic graphs with maximal energy, Lin. Algebra Appl. 356 (2002) 27–36.
  * [11] Y. Hou, T. Lei, Characteristic polynomials of skew-adjacency matrices of oriented graphs, The Electronic Journal of Combinatorics 18(2011) # P156.
  * [12] X. Li, J. Zhang, B. Zhou, On unicyclic conjugated molecules with minimal energies, J. Math. Chem. 42 (2007) 729–740.
  * [13] V. Nikiforov, The energy of graphs and matrices, J. Math. Anal. Appl. 320,(2007) 1472–1475.
  * [14] I. Pena, J. Rada, Energy of digraphs, Linear and Multilinear Algebra 56(2008) 565-579.
  * [15] W. H. Wang, A. Chang, L. Z. Zhang, D. Q. Lu, Unicyclic H$\ddot{u}$ckel molecular graphs with minimal energy, J. Math. Chem. 39 (2006), 231–241.
  * [16] W. H. Wang, A. Chang, D. Q. Lu, Unicyclic graphs possessing Kekul structures with minimal energy, J. Math. Chem. 42 (2007) 311–320.