Post-red-giant-branch Planetary Nebulae
Jones, D.1,2,3, Hillwig, T.C.4, and Reindl, N.5
1 Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain

2 Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna,
Tenerife, Spain

3 Nordic Optical Telescope, Rambla José Ana Fernández Pérez 7, 38711, Breña
Baja, Spain

4 Department of Physics and Astronomy, Valparaiso University, Valparaiso, IN
46383, USA

5 Institute for Physics and Astronomy, University of Potsdam, Karl-Liebknecht-
Str. 24/25, D-14476 Potsdam, Germany

## Abstract

Common envelope events have been associated with the formation of a planetary
nebulae since its proposition more than forty five years ago. However, until
recently there have been doubts as to whether a common envelope while the
donor is ascending the red giant branch, rather than the subsequent asymptotic
red giant branch, would result in a planetary nebula. There is now strong
theoretical and observational evidence to suggest that some planetary nebulae
are, indeed, the products of common envelope phases which occurred while the
nebular progenitor was on the red giant branch. The characterisation of these
systems is challenging but has the potential to reveal much about the common
envelope – a critical evolutionary phase in the formation of a plethora of
interesting astrophysical phenomena.

##  1 Introduction

The common envelope (CE) phase is perhaps the most poorly understood aspect of
close-binary evolution [9], in spite of its critical nature for the formation
of a number of important phenomena – ranging from type ia supernovae
(irrespective of their formation channel) through to merging binary black
holes recently observed in gravitational waves. A CE is formed following run-
away unstable mass transfer in a close binary and ultimately leads to dramatic
in-spiral of the binary, sometimes to merger, and the ejection of the donor’s
envelope (partial ejection in the case of a merger). Planetary nebulae (PNe)
have long been seen as important probes of the CE, being identified by
Paczyński (one of the first to outline the CE scenario) as the critical proof
of its existence [16]. They are of particular interest because post-CE PNe are
the immediate after-products of the CE (where the central stars have yet to
have time to relax/evolve), and the only systems where the ejected envelope is
still observable111The envelopes surrounding (Luminous) Red Novae represent
the partial ejection of the evolution as the stars merged prior to complete
ejection [14]. (as the PN itself).

Until recently, it was thought that observable PNe would only be formed
following CE events occurring when the donor was on the asymptotic giant
branch (AGB), or close to the very tip of the red giant branch (RGB).
Otherwise, the post-RGB post-CE remnants were expected to evolve too slowly or
to too low a maximum temperature in order to ionise the ejected material
before it dissipated into the surrounding interstellar medium [8]. Updated
simulations by Hall et al. have since shown that most CE events on the RGB
should, in fact, lead to observable PNe, revising the minimum mass limit for
the post-RGB central star down to around 0.25 M⊙ (c.f. the previous estimate
of $>$0.4 M⊙). This is particularly interesting in understanding the
population of PNe as a whole given that at least 20% of all PNe are found to
be the products of a CE phase [11] and the majority of post-CE binaries are
post-RGB [17]. As such, a number of recent studies have focussed on the search
for and characterisation of post-RGB central stars of PNe, and their
implications for our understanding of the CE. In this proceedings, I present a
review of this recent work and its future prospects.

##  2 Candidates

Roughly eight strong candidate post-RGB central stars have thus far been
identified (see table 1). Although the central stars are orders of magnitude
more luminous than their low-mass (M-type) companions, these systems are
generally double-lined to the high levels of irradiation leading to the
formation of spectroscopically detectable emission lines on the day face of
the companions. Similarly, even when they are not eclipsing, the central star
systems present strong photometric variability due to the differing projection
of the irradiated face of the companion as a function of orbital phase. Here,
we will discuss further the properties of each of these candidates as
determined by detailed modelling of light and radial velocity curves, and/or
through non-local thermal equilibrium (NLTE) modelling of the spectrum of the
hot post-RGB component.

Table 1: Candidate post-RGB central stars of PNe PN | RA | Dec | Period | Reference
---|---|---|---|---
Abell 46 | 18:31:18.3 |  $+$26:56:12.9 | 0.472 d | [1]
ESO 330-9 | 16:02:18.1 |  $-$41:26:49.5 | 0.296 d | [7]
HaTr 4 | 16:45:00.2 |  $-$51:12:21.1 | 1.738 d | [6]
HaTr 7 | 17:54:09.5 |  $-$60:49:57.7 | 0.322 d | [7]
HaWe 13 | 19:31:07.2 |  $-$03:42:31.5 | – | [4, 15]
Hf 2-2 | 18:32:30.8 |  $-$28:43:20.0 | 0.399 d | [6]
Ou 5 | 21:14:20.0 |  $+$43:41:36.0 | 0.364 d | [13]
PN G283.7$-$05.1 | 09:58:32.3 |  $-$61:26:40.0 | 0.246 d | [12]

###  2.1 Abell 46

Simultaneous light and radial velocity curve modelling of the HW Vir-type
central star of Abell 46 derives a mass of $0.51\pm 0.05$ M⊙ with a luminosity
of $\log L\sim 2.2$ [1]. While the derived mass is a little on the high side
for a post-RGB central star, the uncertainties do clearly encompass post-RGB
masses (see figure 1). The low luminosity on the other hand is entirely
consistent with a post-RGB remnant. Furthermore, the nebula has been found to
have a very large abundance discrepancy factor, a property speculatively
linked with post-RGB evolution [10]. Interestingly, the surface gravity and
temperature implied by the light and radial velocity curve modelling are
consistent with a post-RGB evolutionary track for a much lower mass than that
derived by the modelling (see figure 2). An NLTE modelling of the central star
spectrum could better constrain the effective temperature and surface gravity,
placing stronger constraints on the possible post-RGB nature of Abell 46.

Figure 1:  The core mass evolution of stars for a range of initial masses. The
evolution to the tip of the RGB is shown in red with the subsequent AGB
evolution in black. This demonstrates that the maximum mass for the post-RGB
central star is approximately 0.46 M⊙.

###  2.2 ESO 330-9

NLTE modelling of the spectrum of the central star of ESO 330-9 gives a
surface gravity and effective temperature consistent with post-RGB
evolutionary tracks (see figure 2). Furthermore, an excellent fit to the light
and radial velocity curves of the binary could be found constraining the
parameters of the primary to match those derived from the NLTE modelling [7].
The light curve itself is not eclipsing so this is perhaps not a stringent
test of the system parameters but, nonetheless, ESO 330-9 is an excellent
candidate post-RGB central star.

###  2.3 HaTr 4

Multi-band modelling of the eclipsing light curve of HaTr 4 results in a range
of radii and temperatures consistent with post-RGB evolution [6].
Unfortunately, without radial velocities to properly constrain the mass or
NLTE spectral modelling, it is difficult to assess the true nature of the
central star.

###  2.4 HaTr 7

NLTE modelling of the spectrum of the central star of HaTr 7 gives parameters
which lie between post-RGB and post-AGB tracks on the Kiel diagram (figure 2).
Just as for ESO 330-9, an excellent fit to the non-eclipsing light curve and
radial velocity curves could be found with parameters which match those
predicted by the NLTE modelling [7].

Figure 2:  A Kiel diagram showing the properties of the five strongest
candidate post-RGB stars for which effective temperatures and surface
gravities have been measured (either spectroscopically or via combined light
and radial velocity curve modelling. Figure reproduced from [13].

###  2.5 HaWe 13

The central star of HaWe 13 was identified as being under-luminous based on
its gravity distance, while the nebula presents a morphology consistent with
being the product of a CE phase [4]. However, no study has as yet definitively
revealed the binary nature of the central star.

###  2.6 Hf 2-2

Just as for HaTr 4, modelling of multi-band non-eclipsing light curves
indicates a low luminosity central star. A single spectrum of the central star
was also taken allowing the mass ratio to also be constrained, indicating a
low-mass primary [6].

###  2.7 Ou 5

Simultaneous light and radial velocity curve modelling of the eclipsing binary
(HW Vir-type) at the heart of Ou 5 derives a mass of $0.50\pm 0.06$ M⊙ with a
luminosity of $\log L\sim 2.1$ [13] – very similar to the parameters derived
for Abell 46. Just as for Abell 46, the implied surface gravity and
temperature lie on a post-RGB evolutionary track for a much lower mass than
that derived from the light and radial velocity curve modelling (figure 2).

###  2.8 PN G283.7$-$05.1

Simultaneous light and radial velocity curve modelling of the eclipsing binary
(HW Vir-type) in PN G283.7$-$05.1 leads to a definitively post-RGB mass
determination of $0.34\pm 0.05$ M⊙, however the implied luminosity is
extremely high for a post-RGB star – $\log L\sim 3.5$ [12]. Similarly, the
surface gravity and temperature of the model central star reside in a
distinctly post-AGB region of the Kiel diagram (figure 2). Nonetheless,
comparison of the spectra, used to derive the radial velocities of the system,
with NLTE models indicates that the surface gravity and temperature from the
light and radial velocity curve modelling are indeed consistent with the
observed spectroscopy.

##  3 Conclusions

Theoretical studies indicate that post-RGB masses greater than approximately
0.25 M⊙ should lead to observable PNe [5]. Furthermore, comparison of the
orbital period distribution of naked post-CE binaries (i.e. those that are now
too old to reside inside a visible PNe) – a significant fraction of which are
post-RGB [17] – and the orbital period distribution of post-CE central stars
of PNe indicates that the two are consistent with being drawn from the same
population [2]. As such, there is little reason to believe that the currently
known sample of post-CE central stars (now roughly 100) contains a number of
post-RGB systems. However, as yet, no definitively post-RGB central star of a
PN has been identified, although a handful of excellent candidates have been
found. If the strongest of the candidates are indeed post-RGB then this
reveals important information about the CE – in particular, given the apparent
discrepancy between masses derived by from radial velocities and those based
on their loci in the Kiel diagram (Abell 46, Ou 5, PN G283.7$-$05.1), that
post-RGB evolutionary tracks may not be entirely representative of post-CE
evolution. Any difference between post-CE models and reality can provide an
invaluable window into the CE process. Additionally, RGB CEs have been
speculatively linked to the extreme abundance discrepancies in observed in
some PNe. As such, continued attempts to discover and characterise post-RGB
central stars may be key not only in constraining the importance of this
formation pathway among PNe but also in constraining the CE and its relation
to other long-standing astrophysical questions.

## Acknowledgments

D.J. acknowledges support from the Erasmus+ programme of the European Union
under grant number 2020-1-CZ01-KA203-078200.

## References

  * [1] Afşar, M., & Ibanoǧlu, C., 2008, MNRAS, 391, 802
  * [2] Boffin, H.M.J., & Jones, D., 2019, The Importance of Binaries in the Formation and Evolution of Planetary Nebulae (Springer Nature)
  * [3] Brown, A.J., Jones, D., Boffin, H.M.J., & Van Winckel, H., 2019, MNRAS, 482, 4951
  * [4] Frew, D.J., Parker, Q.A., & Bojičić, I.S., 2016, MNRAS, 455, 1459
  * [5] Hall, P.D., Tout, C.A., Izzard, R.G, & Keller, D., 2013, MNRAS, 435, 2048
  * [6] Hillwig, T.C., Bond, H.E., Frew, D.J., et al., 2016, AJ, 152, 34
  * [7] Hillwig, T.C., Frew, D.J., Reindl, N., et al., 2017, AJ, 153, 24
  * [8] Iben, I. Jr., & Tutukov, A.V., 1993, ApJ, 418, 343
  * [9] Ivanova, N., Justham, S., Chen, X., et al., 2013, A&ARv, 21, 59
  * [10] Jones, D., Wesson, R., García-Rojas, J., et al., 2016, MNRAS, 455, 3263
  * [11] Jones, D., & Boffin, H.M.J., 2017, Nat. Astron., 1, 0117
  * [12] Jones, D., Boffin, H.M.J., Hibbert, J., et al., 2020, A&A, 642, 108
  * [13] Jones, D., Munday, J., Corradi, R.L.M., et al., 2022, MNRAS, 510, 3102
  * [14] Kamiński, T., Steffen, W., Bujarrabal, V., et al. 2021, A&A, 646, 1
  * [15] Napiwotzki, R., 1999, A&A, 350, 101
  * [16] Paczyński, B., 1976, IAUS, 73, 75
  * [17] Rebassa-Mansergas, A., Nebot Gómez-Morán, A., Schreiber, M.R., et al., 2011, MNRAS, 413, 1121


# Light trajectory and shadow shape in the rotating naked singularity

Vishva Patel  International Center for Cosmology, Charusat University, Anand,
GUJ 388421, India  Divya Tahelyani<EMAIL_ADDRESS>International
Center for Cosmology, Charusat University, Anand, GUJ 388421, India  Ashok B.
Joshi<EMAIL_ADDRESS>International Center for Cosmology, Charusat
University, Anand, GUJ 388421, India  Dipanjan Dey<EMAIL_ADDRESS>International Center for Cosmology, Charusat University, Anand, GUJ 388421,
India  Pankaj S. Joshi<EMAIL_ADDRESS>International Centre for Space and
Cosmology, Ahmedabad University, Ahmedabad, GUJ 380009, India
International Center for Cosmology, Charusat University, Anand, GUJ 388421,
India

###### Abstract

In this paper, we investigate the light trajectories and shadow properties in
the rotating version of null naked singularity (NNS) spacetime which is
derived using the Newman- Janis algorithm without complexification method. We
discuss some of the geometrical properties and causal structure of Rotating
Naked Singularity (RNS) spacetime. The gravitational lensing in a rotating
naked singularity is analyzed, and the results are compared to those of a Kerr
black hole. In the case of a Kerr black hole, the photon sphere exists for
both prograde and retrograde photon orbits, whereas for RNS, the photon sphere
exists only for retrograde photon orbits. As a result, the naked singularity
projects an arc-shaped shadow that differs from the contour-shaped shadow cast
by a Kerr black hole.

Key words: Black hole, Naked singularity, Shadow, Newman-Janis Algorithm.

##  I Introduction

So far, Einstein’s theory of general relativity is a remarkably successful
theory of gravity. When large enough masses collapse under the pull of their
own gravity, the general theory of relativity predicts that a spacetime
singularity would develop. Under certain conditions, a singularity theorem
demonstrates the inevitability of the emergence of a spacetime singularity.
However, there was no such theorem about the existence of an event horizon
enclosing the singularity. In 1969, Roger Penrose introduced the cosmic
censorship conjecture (CCC), which asserts that horizon-less strong spacetime
singularities are not feasible Penrose:1964wq . Nonetheless, there are several
studies on the continuous gravitational collapse of the inhomogeneous matter
cloud that illustrates that spacetime singularities developed during
gravitational collapse can be seen by an outside observer joshi ; goswami ;
mosani1 ; mosani2 ; mosani3 ; mosani4 ; Deshingkar:1998ge ; Jhingan:2014gpa .
The visibility of spacetime singularity depends upon the initial conditions of
the collapsing matter Joshi:2011zm . In the regime of general relativity, such
a scenario can also arise where the continual gravitational collapse may lead
to an equilibrium static spacetime with a central visible or naked
singularity, e.g., Joshi-Malafarina-Narayan (JMN) spacetime Joshi:2011zm ; JNW
; Joshi:2020tlq .

It is evident from the above diacussion that the naked singularity can arise
as the final state of continuous gravitational collapse of an inhomogeneous
matter cloud. Therefore, if naked singularity exists in nature, it would have
distinct physical traces. Photon trajectories may extract information about
the distinguish spacetime structure near a naked singularity. Therefore,
detailed study of null geodesics near naked singularity is very much important
to predict the physical signature of naked singularity. The route is offered
by gravitational lensing by naked singularities, as described by several
authors in recent literature Virbhadra:2002ju ; Sahu:2012 ; Paul:2020fc ;
Dey:2013ya ; Dey:2020hf ; Shaikh:2019lcc . Gravitational lensing is defined as
both the deflection of light and the variation in visual brightness of the
radiating source induced by a gravitational field. The study of analyzing the
shadow properties of a naked singularity is a curious and useful technique for
identifying the spacetime structure of naked singularity. The observation of
the shadow also gives an initial way of determining the properties of the
naked singularity. In 2019, the Event Horizon Telescope (EHT) group released
the shadow image of a compact object located at the center of Messier 87 (M87)
galaxy EventHorizonTelescope:2019dse . Besides that, recently the EHT group
released the shadow image of a compact object located at the center of Milky
way galaxy Sgr $A^{*}$ EventHorizonTelescope:2022xnr ;
EventHorizonTelescope:2022tzy ; EventHorizonTelescope:2022ago ;
EventHorizonTelescope:2022vjs ; EventHorizonTelescope:2022wok ;
EventHorizonTelescope:2022urf ; EventHorizonTelescope:2022xqj ;
EventHorizonTelescope:2022xqj ; EventHorizonTelescope:2022exc ;
EventHorizonTelescope:2022gsd . They explore the shadow image cast by various
models and compared it to the shadow image of Sgr $A^{*}$, concluding that the
JMN-1 naked singularity is one of the best possible non-spinning black hole
mimickers EventHorizonTelescope:2022xqj . Furthermore, they suggest that the
metric tests cannot rule out the possibility of the naked singularities as a
Sgr $A^{*}$ compact object EventHorizonTelescope:2022xqj . These observations
have opened up new possibilities for analysing various naked singularity
spacetimes as black hole mimickers. In literature, many researchers probe into
the shadow cast by other types of compact objects such as gravastar, regular
black holes, and wormholes Gralla:2019xty ; Abdikamalov:2019ztb ;
Vagnozzi:2019apd ; Gyulchev:2019tvk ; Dey:2013yga ; Dey+15 ; Dey:2020haf ;
atamurotov_2015 ; abdujabbarov_2015b ; ohgami_2015 ; stuchlik_2019 ; Kaur ;
Sakai ; Shaikh:2021 ; Li:2021 ; Hu:2020usx . It is now known that event
horizon and photon sphere are not necessary for a shadow formation. In papers
Joshi:2020tlq , authors show that existence of upper bound of the effective
potential of null geodesic is the necessary and sufficient condition for
shadow formation.

Since every compact object in our universe possesses inherent angular
momentum, a rotating counterpart of the singular static spacetime is essential
for a more realistic scenario. In this work, we first construct a rotating
version of null naked singularity (NNS) spacetime to explore its physical
features. The rotating naked singularity spacetime could be obtained by
applying the Newmann-Janis Algorithm (NJA) to the static NNS. Furthermore,
while applying the NJA to other spacetime metrics, such as naked singularities
and other black hole solutions, the final result in Eddington-Finkelstein
coordinates (EFC) may not be transformed into the Boyer-Lindquist coordinates
(BLC) due to the complexification process mentioned in Mustapha2011 . In
Mustapha2011 , Mustapha Azreg-Anou show how the rotating spacetime obtained by
the NJA and written in the EFC could be transformed into the BLC by skipping
the complexification process. Recently, the rotating form of the polytropic
black hole and Janis-Newmann-Winicour naked singularity spacetimes have been
generated using the NJA without complexification Solanki:2021mkt . In this
paper, we apply the NJA without complexification to set the rotating
counterpart of static NNS spacetime discussed in Joshi:2020tlq ; Paul:2020fc ;
Dey:2020bgo . We analyse the properties of light geodesics and the shadow
around RNS and compare the results with those of the Kerr black hole. The
results reveal that the geometry of the rotating version of NNS spacetime has
significant different features in terms of the light trajectories and shadow
shape than the Kerr black hole, which provides an effective tool for
distinguishing the RNS from the Kerr black hole.

This paper is organised as follows. In Section (II), we obtain the rotating
naked singularity (RNS) solution by applying the NJA to static NNS spacetime.
In Section (III), we investigate whether the energy conditions are satisfied
in RNS spacetime. In Section (IV), we compare the light trajectory in the Kerr
black hole and the RNS spacetime. We then compare the shadow shape for the
Kerr black hole and the RNS spacetime. Ultimately in Section (VI), we conclude
and discuss the outcome of the work. Throughout the paper, we consider
geometrized units, where the gravitational constant(G) and the speed of
light(c) are set equal to one. The signature of the metric is (-,+,+,+).

##  II Newman-Janis Algorithm without Complexification

To apply the Newman- Janis algorithm without complexification, let us consider
a general form of spherically symmetric and static spacetime, in which, the
$g_{tt}(r)$ and $g_{rr}(r)$ are inversely proportional to each other. The
spacetime metric is,

$ds^{2}=-A(r)dt^{2}+\frac{dr^{2}}{B(r)}+C(r)d\Omega^{2}.$ (1)

where, $d\Omega^{2}=d\theta^{2}+\sin^{2}\theta d\phi^{2}$. The $g_{tt}(r)$
term is $A(r)$, $g_{rr}(r)$ term is $B(r)$ and $C(r)=r^{2}$. One can use the
null coordinates $\\{u,r,\theta,\phi\\}$ to describe the above spacetime. The
coordinate transformation can be written as,

$du=dt-\frac{dr}{\sqrt{A(r)B(r)}},$ (2)

and in that null coordinates, the spacetime metric (1) takes the following
form,

$ds^{2}=-A(r)du^{2}-2\sqrt{\frac{A(r)}{B(r)}}dudr+C(r)d\Omega^{2}.$ (3)

One can write the null tetrads of the above spacetime metric as,

$\displaystyle l^{\mu}=\delta^{\mu}_{1}\,,$ (4) $\displaystyle
n^{\mu}=\sqrt{\frac{B(r)}{A(r)}}\delta^{\mu}_{0}-\frac{1}{2}B(r)\delta^{\mu}_{1}\,,$
(5) $\displaystyle
m^{\mu}=\frac{1}{\sqrt{2C(r)}}\left(\delta^{\mu}_{2}+\frac{i}{\sin\theta}\delta^{\mu}_{3}\right)$
(6)

The complexification procedure of the NJA extrapolates the function (i.e.,
$A(r)$, $B(r)$, $C(r))$ as the real function of the radial coordinate r and
its complex conjugate $\bar{r}$. The solution is to skip this step and move on
to the next, which requires complex coordinate transformations. In the complex
coordinate transformation the spacetime metric is derived in the null
coordinates $\\{u,r,\theta,\Phi\\}$. One can transform this into the Boyer-
Lindquist coordinates (BLC) $\\{t,r,\theta,\phi\\}$ using the transformation,

$\displaystyle du=dt-\lambda(r)dr$ (7) $\displaystyle d\Phi=d\phi-\chi(r)dr$
(8)

where we get the BLC transformation function $\chi$ and $\lambda$ as,

$\displaystyle\chi(r)=\frac{a}{CB+a^{2}},$ (9)
$\displaystyle\lambda(r)=\frac{K+a^{2}}{CB+a^{2}},$ (10)

where, $K=C\sqrt{\frac{B}{A}}$. In the Boyer- Lindquist coordinates
($t,r,\theta,\phi$), the metric can be represented as,

$\displaystyle
ds^{2}=-\left(1-\frac{2f}{\rho^{2}}\right)dt^{2}+\frac{\rho^{2}}{\Delta}dr^{2}+\rho^{2}d\theta^{2}+\frac{\Sigma\sin^{2}\theta}{\rho^{2}}d\phi^{2}$
$\displaystyle-\,\frac{4af\sin^{2}\theta}{\rho^{2}}dtd\phi.\,\,\,\,\,\,\,\,\,\,$
(11)

Where,

$\displaystyle\rho^{2}=K(r)+a^{2}\cos^{2}\theta\,,$ (12) $\displaystyle
f=\frac{K(r)-C(r)B(r)}{2}\,,$ (13) $\displaystyle\Delta=C(r)B(r)+a^{2}\,,$
(14)
$\displaystyle\Sigma=\big{(}K(r)+a^{2}\big{)}^{2}-a^{2}\Delta\sin^{2}\theta.$
(15)

The rotating solution (11) reduces to the static when $a\to 0$, and
spherically symmetric spacetime metric if $A(r)$ and $B(r)$ are inversely
proportional to each other.

##  III Rotating compact object

In this section, we present a new axially symmetric, asymptotically flat RNS
spacetime. We first discuss the basic spacetime properties of this RNS
spacetime, and then we briefly discuss the Kerr black hole spacetime.

###  III.1 Rotating Naked Singularity (RNS) Spacetime

In the section (II), we develop the rotating form of static and spherically
symmetric metric using Janis- Newman algorithm without complexification. In
Joshi2020 , authors describe a null-naked singularity spacetime which
resembles with Schwarzschild spacetime far away from the center. The metric of
the null-naked singularity spacetime is,

$ds^{2}=-\frac{1}{\left(1+\frac{M}{r}\right)^{2}}dt^{2}+\left(1+\frac{M}{r}\right)^{2}dr^{2}+r^{2}d\Omega^{2}\,\,,$
(16)

where $d\Omega^{2}=d\theta^{2}+\sin^{2}\theta d\phi^{2}$ and $M$ is the
Arnowitt- Deser-Misner (ADM) mass of the above spacetime metric. One can write
down the following rotating counterpart of the above mentioned static NNS
spacetime using the NJA without complexification technique described in the
previous section,

$ds^{2}=-\left(1-\frac{2f}{\rho^{2}}\right)dt^{2}-\frac{4af\sin^{2}\theta}{\rho^{2}}dtd\phi+\frac{\Sigma\sin^{2}\theta}{\rho^{2}}d\phi^{2}+\frac{\rho^{2}}{\Delta}dr^{2}+\rho^{2}d\theta^{2},\\\
$ (17)

where,

$\displaystyle f=\frac{r^{2}M^{2}+2Mr^{3}}{2(M+r)^{2}}\,,$ (18)
$\displaystyle\rho^{2}=r^{2}+a^{2}\cos^{2}\theta\,,$ (19)
$\displaystyle\Delta=\frac{r^{2}}{\left(1+\frac{M}{r}\right)^{2}}+a^{2}\,,$
(20) $\displaystyle\Sigma=(r^{2}+a^{2})^{2}-a^{2}\Delta\sin^{2}\theta.$ (21)

Here, we note that the above rotating metric will reduce to the Eq. (16) in
the limit of $a\to 0$. Now, we must proceed to verify the validity of the
energy conditions in this rotating spacetime. As we know, the energy-momentum
tensor can be written as,

$T^{\mu\nu}=\rho_{e}u^{\mu}u^{\nu}+P_{r}e_{r}^{\mu}e_{r}^{\nu}+P_{\theta}e_{\theta}^{\mu}e_{\theta}^{\nu}+P_{\phi}e_{\phi}^{\mu}e_{\phi}^{\nu}\,\,,$
(22)

where $\rho_{e}$ is the energy density, and $P_{i}$ $(i\to r,\theta,\phi)$ are
the main pressure components. Using the Einstein field equation
$G_{\mu\nu}=T_{\mu\nu}$, we can write down the $\rho_{e}$ and $P_{i}$ as,

$\displaystyle\rho_{e}=G_{\mu\nu}\,u^{\mu}u^{\nu}\,,$ (23) $\displaystyle
P_{r}=G_{\mu\nu}\,e_{r}^{\mu}e_{r}^{\nu}=g^{rr}G_{rr}\,,$ (24) $\displaystyle
P_{\theta}=G_{\mu\nu}\,e_{\theta}^{\mu}e_{\theta}^{\nu}=g^{\theta\theta}\,G_{\theta\theta}\,,$
(25) $\displaystyle P_{\phi}=G_{\mu\nu}\,e_{\phi}^{\mu}e_{\phi}^{\nu}.$ (26)

Figure 1: Penrose diagram of rotating naked singularity

From (23), (24), (25) and (26), the components of the energy-momentum tensor
of the RNS spacetime are as follows (considering $\theta=\pi/2$),

$\displaystyle\rho_{e}=\frac{M^{2}(M+3r)}{r^{2}\left(M+r\right)^{3}}\,,$ (27)
$\displaystyle P_{r}=-\frac{M^{2}(M+3r)}{r^{2}\left(M+r\right)^{3}}\,,$ (28)
$\displaystyle P_{\theta}=P_{\phi}=\frac{3M^{2}}{(M+r)^{4}}.$ (29)

From the above expressions, it can be verified that the RNS spacetime (16)
satisfies the weak energy condition ($\rho_{e}\geq 0$, $\rho_{e}+P_{i}\geq
0$), strong energy condition ($\rho_{e}+\sum P_{i}\geq 0$, $\rho_{e}+P_{i}\geq
0$), the dominant energy condition ($\rho_{e}\geq|P_{i}|$), and the null
energy condition ($\rho_{e}+P_{i}\geq 0$) considering $\theta=\frac{\pi}{2}$.
The presence of horizons in a stationary spacetime geometry are generally
given by the condition $g_{rr}\to\infty$. Equivalently, the null values of
$\Delta(r)$ represent the locations of horizons. From $\Delta(r)=0$, we find
that, the radius of horizons becomes imaginary for RNS spacetime, which
indicates that it is a horizon-less compact object. The Kretschmann scalar
which is defined as,
$k=R_{\alpha\beta\gamma\delta}\,R^{\alpha\beta\gamma\delta}$ becomes infinite
in the RNS spacetime for $r\to 0$ and $\theta\to\frac{\pi}{2}$. This behaviour
of the Kretschmann scalar indicates the presence of a ring singularity in the
equatorial plane with a coordinate radius $a$. The causal structure of the RNS
spacetime can be examined using the Penrose diagram. To generate the Penrose
diagram, we transform the $t$ and $r$ coordinates into $T$ and $R$
coordinates, respectively, by using the following transformations,

$\displaystyle T=\frac{1}{2}(\tan^{-1}(t+r^{*})+\tan^{-1}(t-r^{*}))\,,$ (30)
$\displaystyle R=\frac{1}{2}(\tan^{-1}(t+r^{*})-\tan^{-1}(t-r^{*}))\,.$ (31)

where $r^{*}=\int\left(\frac{g_{rr}}{g_{tt}}\right)^{\frac{1}{2}}dr$. The new
coordinates $T$ and $R$ are used to describe the entire spacetime manifold in
a specific compactified diagram. We need to check the values of $T$ and $R$ at
$r\to 0$ to see if the spacetime singularity is time like or null like. For
the finite value of $t$ and $r\to 0$, if $r^{*}\to-\infty$ then one can
conclude that the singularity at $r=0$ is nulllike, since
$\lim_{r^{*}\to-\infty}R=-\pi\,\,,\lim_{r^{*}\to-\infty}T=0$. On the other
hand, for the finite value of $t$ and $r\to 0$, if the singularity exists at a
constant $r$, then the singularity is timelike. From the Penrose diagram in
Fig. (1), one can see that the singularity in the RNS spacetime is a timelike
singularity.

###  III.2 Kerr Black Hole Spacetime

In the Boyer- Lindquist coordinates $(t,r,\theta,\phi)$, the line element of a
stationary, axisymmeric, and rotating Kerr spacetime is given by Eq. (17)
with,

$\displaystyle a=\frac{J}{M}\,,$ (32) $\displaystyle f=Mr\,,$ (33)
$\displaystyle\rho^{2}=r^{2}+a^{2}\cos^{2}\theta\,,$ (34)
$\displaystyle\Delta=r^{2}-2Mr+a^{2}\,,$ (35)
$\displaystyle\Sigma=(r^{2}+a^{2})^{2}-a^{2}\Delta\sin^{2}\theta\,.$ (36)

The spin parameter $a$ has the dimensions of length in geometrized units. As
discussed above, the positions of horizons in Kerr geometry can be determined
using the condition $\Delta(r)=0$, which gives

$r_{\pm}=M\pm\sqrt{M^{2}-a^{2}}\,.$ (37)

The surface at $r=r_{-}$ represents a Cauchy horizon, while the surface at
$r=r_{+}$ represents the event horizon. In the extreme case $a=M$, both the
horizons coincide at $r=M$. For $a>M$, the horizons become imaginary, and the
singularity in the Kerr black hole becomes naked.

##  IV LIGHT TRAJECTORIES IN ROTATING NAKED SINGULARITY AND KERR BLACK HOLE

A general rotating spacetime in the Boyer-Lindquist co-ordinates can be given
by,

$ds^{2}=-g_{tt}dt^{2}+g_{rr}dr^{2}+g_{\theta\theta}d\theta^{2}+g_{\phi\phi}d\phi^{2}-2g_{t\phi}dtd\phi\,.$
(38)

In order to obtain the shadow of these compact objects, it is required to
study the null geodesics in the spacetimes. The null geodesics can be derived
from the Lagrangian given by

$2\mathcal{L}=g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}\,,$ (39)

where the overdot represents the differentiation with respect to the affine
parameter ($\lambda$) along the light geodesic. For the line element given in
Eq. (38), Eq. (39) becomes

$2\mathcal{L}=-g_{tt}\dot{t}^{2}+g_{rr}\dot{r}^{2}+g_{\theta\theta}\dot{\theta}^{2}+g_{\phi\phi}\dot{\phi}^{2}-2g_{t\phi}\dot{t}\dot{\phi}\,.$
(40)

Setting $\theta=\pi/2$ for equatorial plane in RNS (Eq. (17)) and Kerr black
hole (Eq. (36)) metric, the above equation yields,

$\displaystyle 2\mathcal{L}_{RNS}$ $\displaystyle=$
$\displaystyle-\frac{1}{(1+\frac{M}{r})^{2}}\dot{t}^{2}+\frac{r^{2}}{\Delta}\dot{r}^{2}+\frac{\left(a^{2}+r^{2}\right)^{2}-a^{2}\Delta}{r^{2}}\dot{\phi}^{2}-\frac{2aM(M+2r)}{(M+r)^{2}}\dot{t}\dot{\phi}\,,$
(41) $\displaystyle 2\mathcal{L}_{Kerr}$ $\displaystyle=$
$\displaystyle-\left(1-\frac{2M}{r}\right)\dot{t}^{2}+\frac{r^{2}}{\Delta}\dot{r}^{2}+\left(r^{2}+a^{2}+\frac{2Ma^{2}}{r}\right)\dot{\phi}^{2}-\frac{4aM}{r}\dot{t}\dot{\phi}\,.$
(42)

Since the Lagrangian is independent of $t$ and $\phi$ co-ordiantes, $p_{t}$
and $p_{\phi}$ are conserved along the light geodesic. From the above
Lagrangian, the $t$\- and $\phi$\- component of momenta can be deduced:

$\displaystyle
p_{t}=\frac{\partial\mathcal{L}}{\partial\dot{t}}=-g_{tt}\dot{t}-g_{t\phi}\dot{\phi}=-\epsilon\,,$
(43) $\displaystyle
p_{\phi}=\frac{\partial\mathcal{L}}{\partial\dot{\phi}}=-g_{t\phi}\dot{t}+g_{\phi\phi}\dot{\phi}=L\,,$
(44)

where $\epsilon$ and $L$ are, respectively, the conserved energy and angular
momentum (about the axis of symmetry) of the photon evaluated at spatial
infinity. The $r$\- component of the momentum is given by

$p_{r}=\frac{\partial\mathcal{L}}{\partial\dot{r}}=g_{rr}\dot{r}\,.$ (45)

Solving Eq. (43) and Eq. (44) for $\dot{t}$ and $\dot{\phi}$, we arrive at

$\displaystyle\dot{t}=\frac{-Lg_{t\phi}+\epsilon
g_{\phi\phi}}{g_{tt}g_{\phi\phi}+g^{2}_{t\phi}}\,,$ (46)
$\displaystyle\dot{\phi}=\frac{Lg_{tt}+\epsilon
g_{t\phi}}{g_{tt}g_{\phi\phi}+g^{2}_{t\phi}}\,.$ (47)

Four velocity $\dot{x}^{\mu}$ in Eq. (39) is normalized to zero for null
geodesics. Hence, for photons we have $\mathcal{L}=0$. Inserting the solutions
Eq. (46) and Eq. (47) in Eq. (40), we obtain

$\dot{r}^{2}=\frac{-g_{tt}L^{2}-2\epsilon
Lg_{t\phi}+\epsilon^{2}g_{\phi\phi}}{g_{rr}(g_{tt}g_{\phi\phi}+g^{2}_{\phi\phi})}\,.$
(48)

The shape of the photon orbit can be described by how $\phi$ changes with $r$.
Using Eq. (47) and Eq. (48), we find

$\frac{d\phi}{dr}=\frac{(Lg_{tt}+\epsilon
g_{t\phi})\sqrt{g_{rr}}}{\sqrt{(g_{tt}g_{\phi\phi}+g^{2}_{t\phi})(-g_{tt}L^{2}-2\epsilon
Lg_{t\phi}+\epsilon^{2}g_{\phi\phi})}}$ (49)

For the simplicity, we define a new variable in above equation such that
$u=1/r$. Now, we can write down the second order differentiation of $u$ with
respect to $\phi$ in the following way,

$\frac{d^{2}u}{d{\phi}^{2}}=\frac{2}{r^{3}}\left(\frac{dr}{d\phi}\right)^{2}-\frac{1}{2r^{2}}\frac{d}{dr}\left(\frac{dr}{d\phi}\right)^{2}\,.$
(50)

We can construct the orbit equation in the RNS and Kerr spacetimes by
employing the corresponding expression of $\frac{dr}{d\phi}$ and the
components of the metric tensor of the spacetimes in the Eq. (50). The
trajectories in the rotating equatorial geometry will differ depending on
whether the particle or photon rotates about the axis of symmetry in the same
direction (co-rotating or prograde) as the rotating compact object or in the
opposite direction (counter-rotating or retrograde). By numerically solving
the orbit equations and plotting the resultant light trajectories, we can
easily visualise how light is lensed in the RNS and Kerr spacetime. The
behaviour of light trajectories and the appearance of shadow shapes in RNS and
Kerr black hole are discussed in the next section.

##  V Equatorial Photon Orbits and Shadow Outline

In the general theory of relativity, the Hamilton-Jacobi formalism is
favourable for uncovering the nature of the shadow. The action in the
Hamilton-Jacobi equation is given by,

$2\frac{\partial S}{\partial\lambda}=-g^{\mu\nu}\frac{\partial S}{\partial
x^{\mu}}\frac{\partial S}{\partial x^{\nu}},$ (51)

where $S$ is the Jacobi action, also known as Hamilton’s principal function,
$\lambda$ is an affine parameter, and $g^{\mu\nu}$ is the metric tensor of the
geometry. If the above differential Eq. (51) has a separable solution, then in
terms of the known constants of the motion, it should take the form,

$S=\frac{1}{2}\mu^{2}\lambda-\epsilon t+L\phi+S_{r}(r)+S_{\theta}(\theta)\,.$
(52)

(a) For an rotating naked singularity

(b) For a Kerr black hole

Figure 2: Fig. (2(a)) and (2(b)) demonstrate the behaviour of $\eta$ and $\xi$
for the RNS and Kerr black hole as a function of radial coordinate $r$. Here,
we consider $a=0.8$ and $M=1$.

Here, $\mu$ is the rest mass of the test particle, and $\mu=0$ for null
geodesic. The functions $S_{r}(r)$ and $S_{\theta}(\theta)$ are variable ($r$
and $\theta$) specified functions.

Substituting Eq. (52) into Eq. (51) and using the contravariant metric tensor
components, we obtain the following equation,

$(L^{2}cosec^{2}\theta-a^{2}\epsilon^{2})cos^{2}\theta-\left(\frac{dS_{\theta}(\theta)}{d\theta}\right)^{2}=\frac{1}{\Delta}((r^{2}+a^{2})\epsilon-
aL)^{2}-(L-a\epsilon)^{2}+\Delta\left(\frac{dS_{r}(r)}{dr}\right)^{2}\,.$ (53)

Here, the equation above is the only function of $\theta$ and $r$. Therefore,
one can separate both sides of Eq. (53) using the separation constant as,

$(L^{2}cosec^{2}\theta-a^{2}\epsilon^{2})cos^{2}\theta-\left(\frac{dS_{\theta}(\theta)}{d\theta}\right)^{2}=K,$
(54)

(a) Prograde Light trajectories in the Kerr Black hole for $M=1,a=0.8$

(b) Retrograde light trajectories in the Kerr black hole for $M=1,a=0.8$

(c) Prograde Light trajectories in the rotating naked singularity for
$M=1,a=0.8$

(d) Retrograde Light trajectories in the rotating naked singularity for
$M=1,a=0.8$

Figure 3: Figs. (3(a)) and (3(b)) show the lensing of null trajectories in the
equatorial plane of a RNS and Kerr black hole, respectively. In this figure,
the spin axes of the compact objects point out of the page. Here, the red
solid circle and green dashed circle denote the critical impact parameters and
photon spheres, respectively. The scattered light rays are represented by
blue, and those plunging into the compact object are denoted by black. The
light rays denoted by pink correspond to the critical impact parameter. In
Figs. (3(a)\- 3(b)), the black shaded regions at the centre represent the
outer event horizon of the Kerr black hole. Whereas in Figs. (3(c)\- 3(d)),
the central black rings correspond to the ring singularity.

(a) shadow shape for the rotating naked singularity

(b) shadow shape for the Kerr Black hole

Figure 4: The figure shows the shadow outline cast by RNS and Kerr black hole.
For both cases, the value of mass is set to one. The shadow shape indicated by
blue to magenta color lines is formed for the spin parameters 0.9 to 0.1. The
detailed discussion of this shadow images is given in section - (VI).

$\frac{1}{\Delta}((r^{2}+a^{2})\,\epsilon-a\,L)^{2}-(L-a\,\epsilon)^{2}+\Delta\left(\frac{dS_{r}(r)}{dr}\right)^{2}=K\,.$
(55)

In above differential equation, $K$ is the separation constant. Eqs. (54) and
(55) can be considered as a part of effective potential where the tangential
part and the radial part of the effective potential can be written as,

$\displaystyle R(r)$ $\displaystyle=$
$\displaystyle((r^{2}+a^{2})-a\xi)^{2}-\Delta(\eta+(\xi-a)^{2})\,,$ (56)
$\displaystyle\Theta(\theta)$ $\displaystyle=$
$\displaystyle\eta-(\xi^{2}cosec^{2}\theta-a^{2})\,cos^{2}\theta\,.$ (57)

In Eq. (56), $\xi=\frac{L}{\epsilon}$ and $\eta=\frac{K}{\epsilon^{2}}$. In
the equatorial plane, $\xi$ is the impact parameter of light geodesics and it
can be used to distinguish the null geodesics. Note that for a motion of
photon, the $\Theta(\theta)$ and $R(r)$ should be positive. Using the
expressions of $\Theta(\theta)$ and $R(r)$, one can characterize certain
critical values of parameters $\xi~{},\eta$ and under small perturbations,
these critical values cause escape or plunge orbit for photons. On that count,
we mark that the outline of the shadow is majorly dependent on the critical
impact parameters. The standard conditions to determine the unstable circular
orbit of photons are,

$R(r_{ph})=0,\;\frac{dR(r_{ph})}{dr}=0,\;\frac{d^{2}R(r_{ph})}{dr^{2}}\leq
0\,\,.$ (58)

where, $r_{ph}$ is the radius of the photon sphere. Using the conditions (58),
one can obtain the critical impact parameters corresponding to the maxima of
the $R(r)$ for RNS and Kerr Black hole spacetimes as,

$\displaystyle\eta_{RNS}=\frac{-r^{6}+4a^{2}M(M+r)^{3}}{a^{2}(r+2M)^{2}},$
$\displaystyle\xi_{RNS}=\frac{-2a^{2}M^{3}-r(6a^{2}M^{2}+4a^{2}Mr+a^{2}r^{2}+r^{4})}{ar^{2}(r+2M)}\,.$
(60)

and

$\displaystyle\eta_{Kerr}=\frac{r^{3}\left(4a^{2}M-r(r-3M)^{2}\right)}{a^{2}(M-r)^{2}},$
(61) $\displaystyle\xi_{Kerr}=\frac{a^{2}M+a^{2}r-3Mr^{2}+r^{3}}{a(M-r)}\,.$
(62)

Since $\eta=0$ at the equatorial plane, the roots of Eqs. (V) and (61)
determine the critical orbits of the photons ($r_{ph}$) in RNS and Kerr
spacetimes, respectively at that plane. In Figs. (2(a)) and (2(b)), the
behaviour of $\eta$ and $\xi$ as a function of $r$ is shown, where we consider
$a=0.8$ and $M=1$. Orbits with a positive $\xi$ value represent the prograde
equatorial motion of photons, and those with a negative $\xi$ value depict the
retrograde equatorial orbit. For RNS, $\eta_{RNS}$ vanishes at $r=2.0389$, as
shown in Fig. (2(a)). The circular orbit corresponding to this radius is
retrograde, since $\xi_{RNS}$ is negative at $r_{ph}=2.0389$, and we observe
that no photon sphere exists for prograde equatorial motion of photons. Hence,
there is only one photon sphere at $r_{ph}=2.0389$ for the retrograde motion
of null geodesics in RNS spacetime, and the corresponding critical impact
parameter is $\xi_{RNS}(r_{ph})=4.49759$. On the other hand, in Fig. (2(b)),
we see that there are three roots of $\eta_{Kerr}$. As we are only interested
in orbits outside the event horizon, the first root $r_{1}$ must be discarded,
since $r_{1}<r_{+}$. Here, $r_{ph+}$ is associated with the photon orbit
having $\xi_{Kerr}>0$. Therefore, it represents the prograde equatorial orbit
of a photon. In the same way, $r_{ph-}$ corresponds to the retrograde
equatorial orbit, as it is associated with $\xi_{Kerr}<0$. At the equatorial
plane, there is a critical impact parameter ($\xi_{ph}$) corresponding to the
unstable circular photon orbits $(r_{ph})$ below which the incoming null
geodesics will plunge into the central compact object. The null geodesics with
impact parameters larger than $\xi_{ph}$ will be scattered away, while those
with impact parameters equal to $\xi_{ph}$ will remain trapped in the photon
sphere at $r=r_{ph}$. Since the photon sphere corresponds to the unstable
photon orbits, any small perturbation in the $r$-direction would cause the
photons at $r=r_{ph}$ to fall into the singularity or scatter to infinity.
Only the scattered photons from the compact object can reach the observer at
infinity and create bright patches in the observer’s sky. On the other hand,
photons falling into a compact object do not reach the observer at infinity,
resulting in a dark region. The shadow of the compact object is created by the
combination of light and dark spots in the observer’s sky.

Figs. (3(a), 3(b)) and Figs. (3(c), 3(d)) depict the deflection of light
trajectories in the equatorial plane of Kerr black hole and a RNS respectively
with $a=0.8$ and $M=1$. In these figures, both the compact objects have spins
in counter-clockwise direction. As shown in Fig. (3(c)), there is no photon
sphere for prograde light trajectories; therefore, the light rays are
dispersed away from the singularity. However, due to the finite size of the
ring singularity in the equatorial plane, if the impact parameter of the light
geodesics is smaller than the coordinate radius $a$ of the singularity, the
light geodesics will become discontinuous as they approach the singularity and
never reach the observer. For retrograde orbits of light geodesics, there
exists a photon sphere in the equatorial plane, which is shown in Fig. (3(d))
by a green dashed circle. Due to the presence of the photon sphere, light rays
with impact parameters smaller than the critical impact parameter are trapped
inside the photon sphere, creating dark patches in the observer’s sky. Here,
$\xi_{RNS}(r_{ph})=4.49759$ represents the radius of the shadow in the
equatorial plane.

As explained earlier, the parameters $\eta$ and $\xi$ are essential to finding
the shape of the shadow. If we consider the celestial coordinates $\alpha$ and
$\beta$, which are the coordinates of the asymptotic observer sky, then using
these celestial coordinates, one can obtain the apparent shape of the shadow
seen by the asymptotic observer. The general expression to derive celestial
coordinates $\alpha$ and $\beta$ are,

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\lim_{r_{0}\to\infty}\left(-r_{0}^{2}\sin{\theta_{0}}\frac{d\phi}{dr}\Big{|}_{(r_{o},\theta_{o})}\right),$
(63) $\displaystyle\beta$ $\displaystyle=$
$\displaystyle\lim_{r_{0}\to\infty}\left(r_{0}^{2}\frac{d\theta}{dr}\Big{|}_{(r_{o},\theta_{o})}\right),$
(64)

where $r_{0}$ and $\theta_{0}$ are the coordinates of the asymptotic observer.
Here, we consider asymptotically flat rotating metric. Therefore, in the RNS
and Kerr black hole spacetimes, the expressions (63) and (64) will take the
form as,

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\xi\,cosec\,\theta_{0}\,\,,$ (65) $\displaystyle\beta$
$\displaystyle=$
$\displaystyle\pm\sqrt{\eta+a^{2}\cos^{2}{\theta_{0}}-\xi^{2}\cot^{2}{\theta_{0}}}\,\,,$
(66)

By substituting Eqs. (61) and (62) for Kerr black hole, and Eqs. (V) and (60)
for RNS spacetime in above expressions, one can deduce the corresponding
celestial coordinates $\alpha$ and $\beta$ in Kerr and RNS spacetimes. We
construct the outline of the shadows using the radius of unstable circular
orbit in the expressions of celestial coordinates for RNS and Kerr black hole
spacetimes, where we have considered the inclination angle of the observer
$\theta=\pi/2$. Fig. (4) represents the apparent size of the shadow cast by
RNS and Kerr black hole for different values of spin parameters. The analysis
of shadow shape is discussed in the next section.

##  VI Discussion and Conclusion

In this paper, we carry out a comparative study of the properties of the
lightlike geodesics and shadow in the rotating naked singularity (RNS) and
Kerr black hole spacetimes. First, we obtain the stationary and axisymmetric
rotating null singularity spacetime by applying the Newman-Janis algorithm
without the complexification process from the static and spherically symmetric
null naked singularity (NNS) spacetime. The orbit equation for null geodesics
is then derived in order to plot the light trajectories around RNS and Kerr
black holes. From the Hamilton-Jacobi formalism, we obtain the separable
solution of radial $(R(r))$ and azimuthal $(\Theta(\theta))$ parts of the
effective potential. Using the conditions of unstable circular orbits, we have
studied the equatorial prograde and retrograde photon orbits. The following
are the outcomes of this paper:

  * •

In the Kerr black hole, there are two horizons and two stationary limit
surfaces. However, no such horizons or stationary limit surfaces are present
in the RNS spacetime. Therefore, the ergoregion does not exist in an RNS.

  * •

We have studied the properties like nature of singularity, light trajectory,
and shadow shape due to the photon sphere of RNS. As per the given spacetime
in Joshi2020 , the static null naked singularity spacetime contains the null
like singularity. Whereas, due to the inclusion of non-zero spin, the RNS
spacetime has the time like singularity. Furthermore, there is no photon
sphere for static null singularity Joshi2020 . In the case of RNS spacetime,
however, a photon sphere exists for the retrograde motion of the null
geodesics.

  * •

In Fig. (2(b)), we can see that there are two radii of the photon sphere in
the equatorial plane of the Kerr black hole. One radius corresponds to the
prograde motion of light geodesics, which has a positive value of $\xi$ when
$\eta$ = 0. Another radius is for a retrograde motion of light, where the
value of $\xi$ is negative when $\eta$ = 0. Besides that, for an RNS in Fig.
(2(a)), there exists only one photon sphere in the equatorial plane, and that
is for retrograde motion of light geodesics. As a result, in the case of an
RNS, all equatorial light rays travelling in prograde motion will be scattered
by the central singularity. While for the retrograde light orbits, light
geodesics with impact parameters smaller than the critical impact parameter
will fall into the singularity. This can be understood from the behavior of
light trajectories shown in Fig. (3).

  * •

We know that the rotating gravitational sources drag the surrounding
spacetime. This dragging effect can be observed from the deflection of light
geodesics in Fig. (3). From Fig. (3(b)), it is apparent that the dragging
effect is largest near the outer event horizon of the Kerr black hole, causing
the retrograde light orbits to change their direction of motion in the
direction of the spin of the Kerr black hole. Whereas for RNS spacetime, the
dragging effect is most significant around the ring singularity, as seen in
Fig. (3).

  * •

The RNS casts an arc shape shadow which can be seen in Fig. (4(a)). That is
because of the particular occurrence of the photon sphere around RNS. The
shadow in Fig. (4(b)) is cast by Kerr black hole. In both cases, shadow shape
is shown for different spin parameters indicated by blue to magenta color
lines are formed for the spin parameters 0.9 to 0.1, respectively. The shadow
cast by RNS is almost half of that of the Kerr black hole. With the increasing
spin, the shadow size in RNS increases, whereas in Kerr black hole, the
deformation in the shadow shape increases.

In this work, we restricted ourselves to examining the shadow shape due to the
presence of the photon sphere in RNS spacetime. In general, one may show that
in the absence of photon sphere, the shadow can be cast by a rotating compact
object due to the existence of the ring singularity.

###### Acknowledgements.

V. P. would like to acknowledge the support of the SHODH fellowship (ScHeme Of
Developing High quality research - MYSY).

## References

  * (1) R. Penrose, Phys. Rev. Lett. 14, 57-59 (1965).
  * (2) P. S. Joshi and I. H. Dwivedi, Phys. Rev. D 47, 5357 (1993).
  * (3) R. Goswami and P. S. Joshi, Classical Quantum Gravity 21, 3645 (2004).
  * (4) K. Mosani, D. Dey, and P. S. Joshi, Phys. Rev. D 101, 044052 (2020).
  * (5) K. Mosani, D. Dey and P. S. Joshi, Phys. Rev. D 102, no.4, 044037 (2020).
  * (6) K. Mosani, D. Dey and P. S. Joshi, Monthly Notices of the Royal Astronomical Society,  504, 4, 4743–4750 (2021).
  * (7) D. Dey, P. S. Joshi, K. Mosani and V. Vertogradov, arXiv:2103.07190 [gr-qc].
  * (8) S. S. Deshingkar, S. Jhingan and P. S. Joshi, Gen. Rel. Grav. 30, 1477-1499 (1998).
  * (9) S. Jhingan and S. Kaushik, Phys. Rev. D 90, no.2, 024009 (2014).
  * (10) P. S. Joshi, D. Malafarina and R. Narayan, Class. Quant. Grav. 28, 235018 (2011).
  * (11) Allen I. Janis, Ezra T. Newman, and Jeffrey Winicour Phys. Rev. Lett. 20, 878, (1968).
  * (12) A. B. Joshi, D. Dey, P. S. Joshi and P. Bambhaniya, Phys. Rev. D 102, no.2, 024022 (2020).
  * (13) S. Paul, Phys. Rev. D 102, no.6, 064045 (2020).
  * (14) K. S. Virbhadra and G. F. R. Ellis, Phys. Rev. D 65, 103004 (2002).
  * (15) S. Sahu, M. Patil, D. Narasimha and P. S. Joshi, Phys. Rev. D 86, 063010 (2012).
  * (16) D. Dey, K. Bhattacharya and T. Sarkar, Phys. Rev. D 88, 083532 (2013).
  * (17) D. Dey, R. Shaikh and P. S. Joshi, Phys. Rev. D 102, no.4, 044042 (2020).
  * (18) R. Shaikh, P. Kocherlakota, R. Narayan and P. S. Joshi, Mon. Not. Roy. Astron. Soc. 482, no.1, 52-64 (2019).
  * (19) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 875, L1 (2019).
  * (20) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 no.2, L12 (2022).
  * (21) A. E. Broderick et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 no.2, L21 (2022).
  * (22) M. Wielgus et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 no.2, L19 (2022).
  * (23) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 no.2, L13 (2022).
  * (24) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 no.2, L14 (2022).
  * (25) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 no.2, L16 (2022).
  * (26) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 no.2, L17 (2022).
  * (27) K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 no.2, L15 (2022).
  * (28) J. Farah et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 no.2, L18 (2022).
  * (29) S. E. Gralla, D. E. Holz and R. M. Wald, Phys. Rev. D 100, no. 2, 024018 (2019).
  * (30) A. B. Abdikamalov, A. A. Abdujabbarov, D. Ayzenberg, D. Malafarina, C. Bambi and B. Ahmedov, Phys. Rev. D 100, no. 2, 024014 (2019).
  * (31) S. Vagnozzi and L. Visinelli, Phys. Rev. D 100, no. 2, 024020 (2019).
  * (32) Y. Chen, R. Roy, S. Vagnozzi and L. Visinelli, arXiv:2205.06238 [astro-ph.HE] (2022).
  * (33) G. Gyulchev, P. Nedkova, T. Vetsov and S. Yazadjiev, Phys. Rev. D 100, no. 2, 024055 (2019).
  * (34) D. Dey, K. Bhattacharya and T. Sarkar, Phys. Rev. D 88, 083532 (2013).
  * (35) D. Dey, K. Bhattacharya and T. Sarkar, Phys. Rev. D 87, 10, 103505 (2013).
  * (36) D. Dey, R. Shaikh and P. S. Joshi, Phys. Rev. D 102, no.4, 044042 (2020).
  * (37) F. Atamurotov, B. Ahmedov, and A. Abdujabbarov, Phys. Rev. D 92, 084005 (2015).
  * (38) A. A. Abdujabbarov, L. Rezzolla, and B. J. Ahmedov, Mon. Not. R. Astron. Soc. 454, 2423 (2015).
  * (39) T. Ohgami and N. Sakai, Phys. Rev. D 91, 124020 (2015).
  * (40) Z. Stuchlik and J. Schee, Eur. Phys. J. C 79, 44 (2019).
  * (41) K. P. Kaur, P. S. Joshi, D. Dey, A. B. Joshi and R. P. Desai, arXiv:2106.13175 [gr-qc].
  * (42) N. Sakai, H. Saida and T. Tamaki, Phys. Rev. D 90, no.10, 104013 (2014).
  * (43) R. Shaikh, S. Paul, P. Banerjee and T. Sarkar, arXiv:2105.12057 [gr-qc].
  * (44) G. P. Li and K. J. He, JCAP 06, 037 (2021).
  * (45) Z. Hu, Z. Zhong, P. C. Li, M. Guo and B. Chen, Phys. Rev. D 103, no.4, 044057 (2021).
  * (46) R. Shaikh, P. Kocherlakota, R. Narayan and P. S. Joshi,  Mon. Not. Roy. Astron. Soc. 482, 52 (2019).
  * (47) A. B. Joshi, D. Dey, P. S. Joshi and P. Bambhaniya, Phys. Rev. D 102, no.2, 024022 (2020).
  * (48) E. T. Newman, and A. I. Janis, J. Math. Phys. 6, 915 (1965).
  * (49) E. T. Newman, E. Couch, K. Chinnapared, A. Exton, A. Prakash, and R. Torrence, J. Math. Phys. 6, 918 (1965).
  * (50) S.P. Drake and R. Turolla, Class. Quantum Grav. 14, 1883 (1997).
  * (51) Drake, S.P., Szekeres, P., General Relativity and Gravitation 32, 445–457 (2000).
  * (52) Mustapha Azreg-Aïnou, Class.Quantum Grav.28, 148001, (2011).
  * (53) D. N. Solanki, P. Bambhaniya, D. Dey, P. S. Joshi and K. N. Pathak, Eur. Phys. J. C 82, no.1, 77 (2022).
  * (54) Azreg-Aïnou, M., Eur. Phys. J. C 74, 2865 (2014).
  * (55) Mustapha Azreg-Aïnou, Phys. Lett. B 730, 95-98 (2014).
  * (56) Azreg-Aïnou, Mustapha, Phys. Rev. D 90, 064041 (2014).
  * (57) C. W. Misner, K. S. Thorne and J. A. Wheeler, Pinceton University Press, Princeton and Oxford,(2017).
  * (58) M. P. Hobson, G. P. Efstathiou and A. N. Lasenby, Cambridge University Press, Cambridge, (2006).
  * (59) D. Dey, R. Shaikh and P. S. Joshi, Phys. Rev. D 103, no.2, 024015 (2021).
  * (60) P. Bambhaniya, A. B. Joshi, D. Dey and P. S. Joshi, Phys. Rev. D 100, no.12, 124020 (2019).
  * (61) P. Bambhaniya, D. N. Solanki, D. Dey, A. B. Joshi, P. S. Joshi and V. Patel, Eur. Phys. J. C 81, no.3, 205 (2021).
  * (62) D. Tahelyani, A. B. Joshi, D. Dey and P. S. Joshi, arXiv:2205.04055 [gr-qc], (2022).
  * (63) X. M. Deng, Phys. Dark Univ. 30, 100629 (2020).
  * (64) B. Gao and X. M. Deng, Annals Phys. 418, 168194 (2020).


Accepted by PESGM 2021 [gray]0.1 0.65cm 90 20.5cm

# A Cooperative Game Theory-based Approach to Under-frequency Load Shedding
Control

Mukesh Gautam, _Student Member, IEEE_ , Narayan Bhusal, _Student Member, IEEE_
,
and Mohammed Benidris _Member, IEEE_ ,
Department of Electrical & Biomedical Engineering,
University of Nevada, Reno, Reno, NV 89557
(emails<EMAIL_ADDRESS><EMAIL_ADDRESS>and
<EMAIL_ADDRESS>
###### Abstract

This paper proposes a cooperative game theory-based under-frequency load
shedding (UFLS) approach for frequency stability and control in power systems.
UFLS is a crucial factor for frequency stability and control especially in
power grids with high penetration of renewable energy sources and restructured
power systems. Conventional UFLS methods, most of which are off-line, usually
shed fixed amounts of predetermined loads based on a predetermined schedule
which can lead to over or under curtailment of load. This paper presents a co-
operative game theory-based two-stage strategy to effectively and precisely
determine locations and amounts of loads to be shed for UFLS control. In the
first stage, the total amount of loads to be shed, also referred to as deficit
in generation or the disturbance power, is computed using the initial rate of
change of frequency (ROCOF) referred to the equivalent inertial center. In the
second stage, the Shapley value, one of the solution concepts of cooperative
game theory, is used to determine load shedding amounts and locations. The
proposed method is implemented on the reduced 9-bus 3-machine Western
Electricity Coordinating Council (WECC) system and simulated on Real-time
Digital Simulators (RTDS). The results show that the proposed UFLS approach
can effectively return the system to normal state after disturbances.

###### Index Terms:

Adaptive load shedding, co-operative game theory, frequency control, Shapley
Value.

##  I Introduction

Electric utility restructuring and increased penetration of distributed energy
resources has contributed in frequency stability and control challenges in
modern power systems. When existing resources and reserve margin cannot meet
system load, frequency drop occurs in power systems. If the drop in system
frequency reaches a predefined threshold, under-frequency load shedding (UFLS)
schemes are usually applied. Since load shedding is the ultimate way to
prevent system collapse, it should be prompt, able to conserve essential and
critical loads, and able to prevent the system from reaching an underfrequency
condition. Also, control and protection systems of UFLS should be reliable to
avoid unnecessary load shedding to keep the system within its stability
limits. When designing and employing UFLS schemes, the following factors must
be considered: minimum allowable frequency for secure system operation,
frequency threshold for load shedding, number of load shedding locations, and
load shedding amount [1]. Conventionally, under-frequency relays have been
used to gradually shed a fixed amount of predetermined load based on
predetermined schedule to retain balance between load and the available
generation [2].

Adaptive and semi-adaptive UFLS schemes have been known to be more efficient
and effective than shedding fixed amounts of predetermined loads based on
predetermined schedules [1]. Out of these schemes, adaptive UFLS has the
advantage of determining the minimum amount of loading for given system
conditions. These schemes calculate the deficit in generation and then shed
loads in adaptive steps. The computation and implementation of load shedding
should be very quick to avoid frequency instability. Therefore, a reliable and
quick methodology that can simultaneously determine minimum amounts of load
curtailments and best locations for load shedding is pivotal for under-
frequency load shedding schemes.

Various methods have been presented in the literature for frequency control
based on UFLS. An adaptive methodology to determine the minimum amount of load
shedding by underfrequency relays based on the initial rate of change of
frequency (ROCOF) has been proposed in [3]. A two-stage approach for load
shedding has been presented in [4] to address severe under frequency
conditions resulting from islanding. In [5], a coordinated load shedding
strategy based on four variable/dimensional analysis has been proposed to
adaptively specify amount, location, and speed of load shedding. A two-stage
UFLS scheme has been proposed in [6] in which the amount of load shedding is
estimated based on the magnitude of disturbance. A swing equation-based load-
shedding scheme that maintains frequency stability of the grid with minimal
load shedding has been presented in [7]. A centralized three-stage adaptive
UFLS scheme based on frequency-voltage stability has been proposed in [8]. A
sensitivity-based adaptive UFLS approach using Lagrange multipliers has been
proposed in [9].

Conventional UFLS schemes based on analytical modeling pose modeling and
calculation complexities and their computation time increases with the
increase in system size [10]. Also, they may over or under shed loads in
response to frequency drop, which may lead to system instability. Cooperative
game theory-based approaches have the potential to determine minimum load
shedding while being computationally efficient and precise.

Game theory-based approaches (both cooperative and non-cooperative) have been
successfully applied in various fields of power systems. These applications
include power system reliability enhancement, loss allocation, and
transmission expansion planning, to name a few. Game theory-based approaches
have been used for comparing various methods of loss allocation in power
systems [11] and for using vehicle-to-grid in frequency regulation [12]. A
cooperative game theory-based approach is implemented in [13] for loss
reduction allocation of distributed generations using Shapley values (a
solution concept in cooperative game theory).

This paper proposes a cooperative game theory-based UFLS approach that
determines optimal load shedding amounts and locations. The proposed approach
distributes load shedding among various load points based on their Shapley
values. The proposed approach is implemented in two-stages. In the first
stage, deficit in generation is calculated based on the ROCOF in case of
disturbances. In the second stage, the deficit in generation is distributed
based on the equivalent Shapley values of load points. The proposed UFLS
approach is applied on the reduced 9-bus 3-machine Western Electricity
Coordinating Council (WECC) system and simulated on Real-time Digital
Simulators (RTDS). The results show that the proposed UFLS approach can
effectively return the system to normal state after contingencies.

The remainder of the paper is organized as follows. Section II describes the
cooperative game theory including the Shapley value. Section III presents the
formulation of the load cooperative model, which is essential for optimal load
shedding. Section IV explains the proposed UFLS scheme for frequency control.
Section V describes case studies on WECC-9 bus system. Section VI provides
concluding remarks.

##  II Cooperative Game Theory and Shapley Value

In game theory, a cooperative game (or coalitional game) refers to a special
class of games in which players are allowed to compete with each other by
forming alliances (coalitions) among themselves. Mathematically, a cooperative
game is defined through specifying a certain value for each coalition. A
cooperative game has the following components:

  * •

a finite set $\mathscr{N}$, and

  * •

a real-valued set function $V$, called characteristic function, defined on all
sub-sets of $\mathscr{N}$ and satisfies $V(\phi)=0$.

In the game theory terms, $\mathscr{N}$ is defined as the player set, and
$V(S):2^{\mathscr{N}}\rightarrow\mathbb{R}$ is defined as the “worth” or
“value” of coalition $S$, i.e., the total utility that members of $S$ can
acquire if a coalition is formed among themselves and the game is played
without assistance from other players.

###  II-A The Core of a Cooperative Game

In game theory, the core refers to the set of feasible allocations that cannot
be further improved through any other coalitions. Generally, outcomes of a
cooperative game are expressed as n-tuples of utility:
$\alpha=\\{\alpha^{i}:i\in\mathscr{N}\\}$, called payoff vectors that are
measured in some common monetary unit [14]. The core is the set of imputations
under which all sets of coalitions have values less than or equal to the sum
of its members’ payoffs. Thus, $\alpha$ is core if and only if,

$\displaystyle\alpha.e^{S}\geq V(S),\forall S\subset\mathscr{N}$ (1)
$\displaystyle\alpha.e^{\mathscr{N}}=V(\mathscr{N})$ (2)

where $e^{S}$ denotes the n-vector having $e_{i}^{S}=1$ if $i\in S$ and
$e_{i}^{S}=0$ if $i\in\mathscr{N}-S$. Equation (1), also referred to as
stability or coalitional rationality, guarantees that the sum of payoffs of
all players of a coalition is greater than or equal to the worth of the
coalition. Equation (2), also referred to as efficiency, ensures that the
whole amount of the grand coalition is allocated among players.

###  II-B The Shapley Value

The Shapley value is one solution concept of cooperative game theory. The
Shapley value allocates the payoffs in such a way that is fair for cooperative
solutions. The Shapley value of a cooperative game is given as follows [15].

$\psi_{j}(V)=\sum_{S\in 2^{\mathscr{N}},j\in S}\frac{(\lvert
S\rvert-1)!(n-\lvert S\rvert)!}{n!}[V(S)-V(S\backslash\\{j\\})]$ (3)

where $n=\lvert\mathscr{N}\rvert$ is the total number of players.

The Shapley value satisfies the following axioms:

  1. 1.

Efficiency: $\sum_{j\in\mathscr{N}}\psi_{j}(V)=V(\mathscr{N})$.

  2. 2.

Individual Rationality: $\psi_{j}(V)\geq V(\\{j\\}),\forall j\in\mathscr{N}$.

  3. 3.

Symmetry: If $j$ and $k$ are such that $V(S\cup\\{j\\})=V(S\cup\\{k\\})$ for
every coalition $S$ not containing $j$ and $k$, then
$\psi_{j}(V)=\psi_{k}(V)$.

  4. 4.

Dummy Axiom: If $j$ is such that $V(S)=V(S\cup\\{j\\})$ for every coalition
$S$ not containing $j$, then $\psi_{j}(V)=0$.

  5. 5.

Additivity: If $V$ and $W$ are characteristic functions, then
$\psi(V+W)=\psi(V)+\psi(W)$.

##  III Load Cooperative Model Formulation

A load cooperative model is formulated to determine the Shapley value of each
load participating in the load shedding scheme. In this paper, two types of
characteristic functions are implemented for the cooperative model to capture
both transient and steady-state parts of the frequency response. The first
characteristic function utilizes the steady state rise in the system frequency
(referred to the equivalent inertial center) resulting due to participation of
loads in a particular coalition, which captures the steady state frequency
response of the system. The second characteristic function utilizes the
initial ROCOF for computing the worth of loads and their possible coalitions.
The second characteristic function is employed in addition to the first
characteristic function in order to capture the transient part of the system
frequency response. Using each of these characteristic functions, two Shapley
values are computed for each load using (3) and the equivalent Shapley value
of each load is determined by taking their average.

The load cooperative model formulation for the proposed approach can be
enumerated as follows.

  1. 1.

Collect system data including generation data, transmission line data, load
data, etc., which serve as input to the load cooperative model.

  2. 2.

Identify loads that participate in load shedding based on their priority or
any other optimization techniques.

  3. 3.

Generate the list of all possible coalitions of loads. For example, if three
loads ($L_{1}$, $L_{2}$, and $L_{3}$) are participating in load shedding, the
set of all possible coalitions, denoted by $2^{\mathscr{N}}$, is as follows.

$2^{\mathscr{N}}=\\{\phi,\\{L_{1}\\},\\{L_{2}\\},\\{L_{3}\\},\\{L_{1},L_{2}\\},\\{L_{1},L_{3}\\},\\\
\\{L_{2},L_{3}\\},\\{L_{1},L_{2},L_{3}\\}\\}$,
where $\phi$ denotes an empty set.

  4. 4.

For each load and its possible coalitions, compute the steady state rise (due
to load shedding) in system frequency with respect to the equivalent center of
inertia (COI) and the initial ROCOF. The steady state rise in system
frequency, $\Delta f$, and the initial ROCOF of each coalition are computed by
making all the loads of a given coalition participate in ULFS. These values
serve as worth of each load and their coalitions.

  5. 5.

Compute two Shapley values, $\psi_{i}^{1}$ and $\psi_{i}^{2}$, of each load,
$L_{i}$, using the characteristic functions determined in step 4 using (3).

  6. 6.

Determine the equivalent Shapley value, $\psi_{i}^{eqv}$, of each load,
$L_{i}$, taking the average of two Shapley values computed in step 5. The
equivalent Shapley values are computed by taking the average of two Shapley
values to give equal weights to both transient and steady state parts of the
frequency response.

$\psi_{i}^{eqv}=\frac{\psi_{i}^{1}+\psi_{i}^{2}}{2}$ (4)

The load cooperative model formulated by utilizing the procedure explained in
the above steps is essential for determining locations and amounts of loads to
be shed for UFLS control.

##  IV The Proposed UFLS Approach

The load shedding is the fastest and ultimate way of preventing power systems
from blackouts and damages initiated due to frequency drops. Load shedding is
done when spinning reserves are exhausted. In this work, the proposed two-
stage under-frequency load-shedding scheme is implemented in the following
steps.

  1. 1.

Calculation of disturbance power, $P_{d}$, based on observed ROCOF; and

  2. 2.

Determination of locations and amounts of load shedding using the Shapley
values.

In order to calculate disturbance power $P_{d}$, the initial ROCOF referred to
equivalent COI and the swing equation of the $i$th machine with inertia
constant $H_{i}$ are utilized as follows [3, 16].

$(2H_{i}/f_{n})\times(df_{i}/dt)=\Delta p_{i}$ (5) $P_{d}=\sum_{i=1}^{M}\Delta
p_{i}=2\sum_{i=1}^{M}H_{i}\times(df_{c}/dt)/f_{n}$ (6)

where $f_{n}$ is the rated system frequency, $f_{i}$ is the frequency of the
$i$th machine, $\Delta p_{i}$ is the power deviation of the $i$th machine, $M$
is the total number of machines, and $f_{c}$ is the frequency of the
equivalent COI.

The frequency of the equivalent COI describes the average system frequency at
the time of electromechanical transients when individual machine frequencies
are not the same, which can be calculated using (7).

$f_{c}=\sum_{i=1}^{M}H_{i}f_{i}\Big{/}\sum_{i=1}^{M}H_{i}$ (7)

The disturbance power, $P_{d}$, obtained using (7) is distributed among
candidate loads based on their equivalent Shapley values. The locations with
priority or critical loads can be removed from the list of candidate load
shedding locations. Load-shedding steps are distributed among the $n$
candidate loads using distribution factor, $D_{k}$, defined for the $k$th
location as in (8).

$D_{k}=\psi^{eqv}_{k}\Big{/}\sum_{k=1}^{n}\psi^{eqv}_{k}$ (8)

where $\psi^{eqv}_{k}$ is the equivalent Shapley value of the $k$th load.

Now, since the sum of distribution factors of all participating loads is equal
to unity, the amount of load to be shed at the $k$th location can be
calculated as follows.

$p_{k}=D_{k}\times P_{d}$ (9)

The load-shedding amount, $p_{k}$, at the $k$th location obtained from (9) may
not be practical to implement. In such cases, load shedding values are rounded
off in such a manner to make their sum still equal to the total disturbance
power, $P_{d}$.

The proposed methodology or the solution algorithm to determine the location
and amount of load shedding can be summarized follows.

  1. 1.

Provide system data related to lines, loads, transformers, and generators.

  2. 2.

Calculate the disturbance power using initial ROCOF of equivalent COI using
(6).

  3. 3.

Determine the number of load-shedding locations and enumerate all possible
coalitions of loads that participate in load shedding and compute their
characteristic functions.

  4. 4.

For each load-shedding location, compute the equivalent Shapley value of each
load based on (4) and the respective amount of load to be shed using (9). For
example, if load-shedding is to be distributed among three loads, the
equivalent Shapley values and the amounts of loads to be shed are determined
for each of those three loads.

The flowchart of the proposed solution algorithm is shown in Fig. 1. The steps
for calculating characteristic functions of each coalition (shown inside
dashed rectangle in the flowchart) can be done offline. Since the calculation
of characteristic functions can be done off-line, the proposed UFLS approach
can be applied for primary frequency control as it only requires the
calculation of initial ROCOF based on the magnitude of the disturbance.

Figure 1: Flow chart of proposed UFLS scheme

##  V Case Study and Discussions

The proposed method is implemented on the reduced WECC 9-bus system. This
system consists of three generators and three load points with total loading
of 315 MW and 115 MVar. System data and configuration are provided in [17] and
inertia constants are given in Table I. This system has been extensively used
in several power system stability studies.

TABLE I: Machine inertia data for WECC-9 bus system Machine number |  Inertia constant (H)
---|---
1 | 23.64
2 | 6.4
3 | 3.01

In order to compute the worth (or characteristic function) for each load and
their possible coalitions in the tested system, the steady state rise in the
system frequency referred to the equivalent COI, $\Delta f$, and the initial
ROCOF are computed through dynamic simulations using RTDS. When the load of
bus $5$ is allowed to participate in UFLS, the steady-state rise in system
frequency and initial ROCOF are, respectively, $1.2757$ Hz and $1.1189$ Hz/s.
Similarly, if both loads at buses $5$ and $6$ are allowed to participate in
UFLS, the steady-state rise in system frequency and the initial ROCOF are,
respectively, $2.1869$ Hz and $1.9169$ Hz/s. Table II shows the characteristic
functions of possible coalitions in terms of $\Delta f$ and initial ROCOF for
all possible coalitions for this system.

TABLE II: Characteristic Function of Possible Coalitions in Terms of $\Delta f$ and Initial $ROCOF$ Load bus coalitions |  $\Delta f$ (Hz) | Initial $ROCOF$ (Hz/s)
---|---|---
5 | 1.2757 | 1.1189
6 | 0.9196 | 0.7990
8 | 0.9887 | 0.8890
5 & 6 | 2.1869 | 1.9169
5 & 8 | 2.2727 | 2.0092
6 & 8 | 1.9144 | 1.6866
5, 6 & 8 | 3.1883 | 2.8071

For the reduced WECC 9-bus system, the outage of machine-3 is simulated. As
there are three load points in the system, system operators have to determine
a certain number (or steps) of load shedding locations and amounts. Therefore,
in this paper, cases of two and three load-shedding locations are considered.
For the disturbance power of 85 MW, if two locations are chosen for load
shedding, the results show that $49$ MW is to be shed from bus $5$ and $36$ MW
from bus $8$. The load-shedding amount along with the values of equivalent
Shapley values for different loads are given in Table III.

TABLE III: Locations and Sizes of Two Step Load Shedding Load-shedding | Equivalent Shapley | Load-Shedding
---|---|---
Locations | values | Steps
Bus $5$ | $1.2222$ |  $49$ MW
Bus $8$ | $0.9187$ |  $36$ MW

If three load shedding locations are chosen for the same disturbance power of
$85$ MW, the following amounts of loads are shed from buses $5$, $6$, and $8$
respectively: $34$ MW, $24$ MW, and $27$ MW. The load-shedding amounts and the
values of equivalent Shapley values for a three-step load-shedding are as
shown in Table IV.

TABLE IV: Locations and Sizes of Three Step Load-Shedding Load-shedding | Equivalent Shapley | Load-Shedding
---|---|---
Locations | values | Steps
Bus $5$ | $1.1973$ |  $34$ MW
Bus $6$ | $0.8581$ |  $24$ MW
Bus $8$ | $0.9424$ |  $27$ MW

In order to study system dynamics with and without load-shedding, the WECC
9-bus system is simulated on RTDS considering the outage of machine-2 at
second 5 and load-shedding after two seconds of the outage. The three-step
load-shedding amounts as shown in Table IV are used for dynamic simulations on
the RTDS. The plot of system frequency (without and with load-shedding) of the
reduced WECC 9-bus system, which are simulated on the RTDS, is shown in Fig.
2. From the frequency response we can see that without load shedding, the
system frequency drops to nearly $57.7$ Hz at second 10.5, where steady-state
frequency being nearly $58.6$ Hz. However, when load shedding is performed at
second 7 (i.e., two seconds after the outage), the system frequency recovers
to the normal frequency level. The frequency plots in Fig. 2 show that the
load-shedding is effective in returning the system to normal frequency.

Figure 2: Plot of system frequency

The proposed UFLS approach is compared with the sensitivity-based approach
presented in [9]. When the study is performed on a PC with 64-bit Intel i5
core, 3.15 GHz processor, 8 GB RAM, and Windows OS, the sensitivity-based
approach presented in [9] takes $0.23$ seconds for the reduced WECC 9-bus
system and $180.36$ seconds for the New England 39-bus system to distribute
the load shedding amounts. On the other hand, the proposed UFLS approach can
distribute the load shedding amounts within a fraction of milliseconds
irrespective of the system size.

##  VI Conclusion

In this paper, a cooperative game theory-based two-stage UFLS scheme has been
proposed. In the first stage, the computation of disturbance power was
performed using initial ROCOF. In the second stage, equivalent Shapley values
were determined for the loads that are participating in load shedding. In
order to calculate the equivalent Shapley values of the loads, the
characteristic functions of loads and their possible coalitions were
calculated using steady-state rise in system frequency and inital ROCOF
evaluated through dynamic simulation of the system. Using the equivalent
Shapley values of loads, the disturbance power was distributed among different
loads. Through the dynamic simulations of WECC-9 bus system using RTDS, it was
shown that the proposed methodology can be effectively implemented for
frequency control of power systems.

## Acknowledgement

This work was supported in Part by the U.S. Department of Energy (DOE) under
Grant DE-EE0009022. This work was prepared as an account of work sponsored by
an agency of the United States Government. Neither the United States
Government nor any agency thereof, nor any of their employees, makes any
warranty, express or implied, or assumes any legal liability or responsibility
for the accuracy, completeness, or usefulness of any information, apparatus,
product, or process disclosed, or represents that its use would not infringe
privately owned rights.

## References

  * [1] B. Delfino, S. Massucco, A. Morini, P. Scalera, and F. Silvestro, “Implementation and comparison of different under frequency load-shedding schemes,” in _2001 Power Engineering Society Summer Meeting. Conference Proceedings (Cat. No.01CH37262)_ , vol. 1, Vancouver, BC, Canada, Canada, July 2001, pp. 307–312.
  * [2] W. C. New and W. C. New, _Load shedding, load restoration and generator protection using solid-state and electromechanical underfrequency relays_. General Electric Company, Power Systems Management, 1974.
  * [3] P. Anderson and M. Mirheydar, “An adaptive method for setting underfrequency load shedding relays,” _IEEE Transactions on Power Systems_ , vol. 7, no. 2, pp. 647–655, May 1992.
  * [4] Q. Zhou, Z. Li, Q. Wu, and M. Shahidehpour, “Two-stage load shedding for secondary control in hierarchical operation of islanded microgrids,” _IEEE Trans. on Smart Grid_ , vol. 10, no. 3, pp. 3103–3111, May 2019.
  * [5] S. Nourollah and G. B. Gharehpetian, “Coordinated load shedding strategy to restore voltage and frequency of microgrid to secure region,” _IEEE Trans. on Smart Grid_ , vol. 10, no. 4, pp. 4360–4368, July 2019.
  * [6] V. V. Terzija, “Adaptive underfrequency load shedding based on the magnitude of the disturbance estimation,” _IEEE Transactions on Power Systems_ , vol. 21, no. 3, pp. 1260–1266, Aug 2006.
  * [7] Y. Tofis, S. Timotheou, and E. Kyriakides, “Minimal load shedding using the swing equation,” _IEEE Transactions on Power Systems_ , vol. 32, no. 3, pp. 2466–2467, May 2017.
  * [8] T. Shekari, A. Gholami, F. Aminifar, and M. Sanaye-Pasand, “An adaptive wide-area load shedding scheme incorporating power system real-time limitations,” _IEEE Systems Journal_ , vol. 12, no. 1, pp. 759–767, March 2018.
  * [9] M. Gautam, N. Bhusal, and M. Benidris, “A sensitivity-based approach to adaptive under-frequency load shedding,” in _2020 IEEE Texas Power and Energy Conference (TPEC)_. IEEE, 2020, pp. 1–5.
  * [10] H. Haes Alhelou, M. E. Hamedani-Golshan, T. C. Njenda, and P. Siano, “Wide-area measurement system-based optimal multi-stage under-frequency load-shedding in interconnected smart power systems using evolutionary computing techniques,” _Applied Sciences_ , vol. 9, no. 3, p. 508, 2019.
  * [11] D. A. Lima, J. Contreras, and A. Padilha-Feltrin, “A cooperative game theory analysis for transmission loss allocation,” _Electric Power Systems Research_ , vol. 78, no. 2, pp. 264–275, 2008.
  * [12] J. Tan and L. Wang, “A game-theoretic framework for vehicle-to-grid frequency regulation considering smart charging mechanism,” _IEEE Transactions on Smart Grid_ , vol. 8, no. 5, pp. 2358–2369, 2017.
  * [13] K. Shaloudegi, N. Madinehi, S. Hosseinian, and H. A. Abyaneh, “A novel policy for locational marginal price calculation in distribution systems based on loss reduction allocation using game theory,” _IEEE transactions on power systems_ , vol. 27, no. 2, pp. 811–820, 2012.
  * [14] L. S. Shapley and M. Shubik, “Competitive outcomes in the cores of market games,” _International Journal of Game Theory_ , vol. 4, no. 4, pp. 229–237, 1975.
  * [15] I. Curiel, _Cooperative game theory and applications: cooperative games arising from combinatorial optimization problems_. Springer Science & Business Media, 2013, vol. 16.
  * [16] S. K., S. N. Singh, and S. C. Srivastava, “A synchrophasor assisted frequency and voltage stability based load shedding scheme for self-healing of power system,” _IEEE Transactions on Smart Grid_ , vol. 2, no. 2, pp. 221–230, June 2011.
  * [17] P. Anderson and A. A. Faud, _Power System Control and Stability_. IA Ames: Iowa State University, 1977.