# A mechanistic model of the neural entropy increase elicited by psychedelic
drugs

Rubén Herzog  Centro Interdisciplinario de Neurociencia de Valparaíso,
Universidad de Valparaíso, Pje Harrington 287, 2360103 Valparaíso, Chile.
<EMAIL_ADDRESS>Pedro A.M. Mediano  Department of Psychology,
University of Cambridge, Cambridge CB2 3EB, UK  Fernando E. Rosas  Centre for
Psychedelic Research, Department of Brain Science, Imperial College London,
London SW7 2DD, UK.  Data Science Institute, Imperial College London, London
SW7 2AZ, UK  Centre for Complexity Science, Imperial College London, London
SW7 2AZ, UK  Robin Carhart-Harris  Centre for Psychedelic Research, Department
of Brain Science, Imperial College London, London SW7 2DD, UK.  Yonatan Sanz
Perl  Buenos Aires Physics Institute and Physics Department, University of
Buenos Aires, Buenos Aires, Argentina  National Scientific and Technical
Research Council, Buenos Aires, Argentina  Universidad de San Andres, Buenos
Aires, Argentina  Enzo Tagliazucchi  Buenos Aires Physics Institute and
Physics Department, University of Buenos Aires, Buenos Aires, Argentina
National Scientific and Technical Research Council, Buenos Aires, Argentina
Rodrigo Cofre  CIMFAV-Ingemat, Facultad de Ingeniería, Universidad de
Valparaíso, Valparaíso, Chile.

###### Abstract

Psychedelic drugs, including lysergic acid diethylamide (LSD) and other
agonists of the serotonin 2A receptor (5HT2A-R), induce drastic changes in
subjective experience, and provide a unique opportunity to study the
neurobiological basis of consciousness. One of the most notable
neurophysiological signatures of psychedelics, increased entropy in
spontaneous neural activity, is thought to be of relevance to the psychedelic
experience, mediating both acute alterations in consciousness and long-term
effects. However, no clear mechanistic explanation for this entropy increase
has been put forward so far. We sought to do this here by building upon a
recent whole-brain model of serotonergic neuromodulation, to study the
entropic effects of 5HT2A-R activation. Our results reproduce the overall
entropy increase observed in previous experiments in vivo, providing the first
model-based explanation for this phenomenon. We also found that entropy
changes were not uniform across the brain: entropy increased in some regions
and decreased in others, suggesting a topographical reconfiguration mediated
by 5HT2A-R activation. Interestingly, at the whole-brain level, this
reconfiguration was not well explained by 5HT2A-R density, but related closely
to the topological properties of the brain’s anatomical connectivity. These
results help us understand the mechanisms underlying the psychedelic state
and, more generally, the pharmacological modulation of whole-brain activity.

## Introduction

Psychedelic drugs provide a privileged opportunity to study the mind-brain
relationship, and promise to revolutionise some of our current mental health
treatments.[Carhart-Harris2014, Scott2019, Carhart-Harris2016a] However, while
some aspects of the neurochemical action of psychedelics at the neuronal and
sub-neuronal level are well known,[Grasso2016, Wacker2017] our current
understanding of their action at the whole-brain level is still very limited.
A deeper understanding of the mechanisms that trigger the changes in conscious
experience produced by psychedelics would greatly advance our knowledge of
human consciousness, and medical development of psychedelics.

At a subjective level, serotonergic psychedelics (including LSD,
dimethyltryptamine [DMT] and psilocybin, among others) potentially modulate
mood, cognition, perception and self-awareness. At a neurophysiological level,
recent research has identified (among many) two particularly prominent
signatures of the psychedelic state: an overall dysregulation of neural
population activity, most clearly seen as suppression of spectral power in the
alpha (8-$12\text{\,}\mathrm{Hz}$) band;[Carhart-Harris2012, Carhart-
Harris2016, Timmermann2019] and an increase in the signal diversity of the
neural activity, measured through the information-theoretic notion of
entropy.[Schartner2017] In particular, this acute entropy increase has been
linked to both short- and long-term effects of the psychedelic experience,
including certain aspects of the reported subjective experience[Schartner2017]
and subsequent personality changes.[Lebedev2016]

Interestingly, the opposite effects have been reported for states of loss of
consciousness, where a strong decrease in brain entropy has been repeatedly
observed. This seems to be a core feature of loss of consciousness,
generalising across states such as deep sleep,[Schartner2017a] general
anaesthesia,[Zhang2001] and disorders of consciousness.[Casali2013] Together,
these facts are yielding converging evidence that entropy and related measures
offer simple and powerful indices of consciousness.

Based on these findings, Carhart-Harris and colleagues have put forward the
Entropic Brain Hypothesis (EBH) an entropy-based theory of conscious
states.[Carhart-Harris2014, Carhart-Harris2018] The EBH proposes the simple,
yet powerful idea that the variety of states of consciousness can be indexed
by the entropy of key parameters of brain activity; or, in other words, that
the richness of subjective experience is directly related to the richness of
on-going neural activity, where richness can be defined most simply as
diversity.  Note that even when different states can be indexed by their
corresponding diversity, we are not claiming that one state is more conscious
than other. Therefore, investigating the neurobiological origins of such
changes in brain activity is a key step in the study of altered states of
consciousness.

Multiple experiments in humans and animal models have established that the
mind-altering effects of psychedelics depend on agonism specifically at the
serotonin 2A receptor (5HT2A-R).[Preller2018, Kraehenmann2017, Quednow2011] A
recent simulation study involving whole-brain computational modelling
confirmed that the topographic distribution of 5HT2A-R in the human brain is
critical to reproduce the functional connectivity dynamics (FCD) of human fMRI
data recorded under the acute effects of LSD.[Deco2018] Here we build upon
this model to characterise the interplay between the entropy of brain signals,
the distribution of 5HT2A receptors, and structural connectivity properties of
the brain, with the overarching goal of explaining the sources of entropic
effects, and thus altered consciousness, elicited by psychedelic drugs.

## Results

We simulated whole-brain activity using the Dynamic Mean-Field (DMF) model by
Deco et al.,[Deco2014] using parameter values fit to reproduce the dynamics of
fMRI recordings in humans during wakeful rest, as well as under the acute
effects of LSD.[Deco2018] The model consists of interacting pools of
excitatory and inhibitory neural populations, coupled via long-range
excitatory connections informed by the anatomical connectivity of the brain.
The DMF model combines a theoretical model of neural and synaptic dynamics
with two empirical sources of information about the human brain: the human
connectome, i.e. DTI-estimated connectivity between the 90 regions of the
AAL[Tzourio-Mazoyer2002] atlas; and average 5HT2A-R expression across the
human brain obtained with Positron Emission Tomography (PET)
scans.[Beliveau2017]

Each simulation of the DMF model generates 90 time series of excitatory firing
rates, one for each region of the AAL atlas (Fig. 1A). These excitatory firing
rates have a non-linear dependency on the local excitatory inputs, determined
by a frequency-current (F-I) curve.  Given that 5HT2A-R seems to modulate the
sensitivity of neurons rather than driving them towards excitation or
inhibition [Aghajanian1997, Licata1993, Grasso2016], we followed Deco et al.
[Deco2018] and modelled 5HT2A-R activation as a response-gain modulation
[Abbott2005] of this F-I curve dependent on the receptor density at a given
region (Fig. 1B).  We simulated the model in two conditions, with and without
5HT2A-R activation, which, in an analogy with neuroscientific terminology, we
refer to as placebo (PLA) and 5HT2A-R conditions, respectively. In this way,
we obtain 90 time series for each condition, which we keep for further
analysis.

Finally, using these time series, we estimate Shannon’s differential entropy
for each region under both conditions, yielding a topographical distribution
of entropy values (Fig. 1C) that constitutes the main subject of analysis in
this study. Further details about model specification and entropy estimation,
as well as other methodological caveats, are presented in the Methods section
at the end of this article.

Figure 1: Modelling the effect of 5HT2A-R activation on the whole-brain
topographical distribution of entropy. (A) Resting state activity is simulated
using the Dynamic Mean-Field (DMF) model, in which each region’s activity is
represented by a time series of excitatory firing rates  (constrained to 0-15
Hz for visualisation).  The probability density function (PDF) and
differential entropy (h(X)) of each region is then estimated, obtaining a
topographical distribution of entropy values. (B) 5HT2A-R agonism is modelled
as a receptor-density-dependent response gain modulation. Black line is the
frequency-current (F-I) curve of a population without 5HT2A-R agonism, and
coloured curves show the resulting F-I curves of regions with increasing
5HT2A-R agonism. (C) 5HT2A-R activation changes the topographical distribution
of entropy with respect to resting state activity, which constitutes the main
subject of analysis in this study.

### 5HT2A-R activation causes a heterogeneous, non-linear increase in the
entropy of simulated brain activity

In our first analysis, we used the DMF model to test the main prediction of
the EBH: that 5HT2A-R activation causes an increase in the overall entropy of
neural signals (Fig. 2). In line with previous experimental findings with
psychedelic drugs,[Schartner2017] the model shows a significant increase in
the brain’s entropy $h$ due to 5HT2A-R activation, with an average entropy of
$h^{\text{PLA}}=2.20\,\mathrm{nat}$ in the placebo condition, and
$h^{\text{5HT2A}}=2.30\,\mathrm{nat}$ in the 5HT2A-R condition (dashed
vertical lines in Fig. 2B, Wilcoxon signed-rank test $p<10^{-6}$, Cohen’s
$d=0.262\pm 0.011$).  The average firing rates of regions, in turn, remained
unchanged. A closer look at the distribution of entropy changes, however,
reveals a more nuanced picture, with some regions increasing and some
decreasing their entropy as a result of 5HT2A-R activation (Fig. 2D). This
suggests that, according to the model, 5HT2A-R agonism might trigger a complex
reconfiguration of the topographically specific distribution of entropy, and
not a mere homogeneous overall increase.

Figure 2: Non-linear heterogeneous increase of entropy following 5HT2A-R
activation. (A) Effect of 5HT2A-R agonism on the local entropy each of region
in the AAL atlas.  See Supplementary Table 1 for abbreviations. Bars indicate
the (bilateral) average relative change in local entropy, $\Delta h_{n}$, and
error bars indicate 1 standard deviation across 1000 simulations. (B)
Histograms of local entropy values for the condition with (red) and without
(blue) 5HT2A-R activation. 5HT2A-R activation increased both the average and
the spread of the local entropy distribution. (C) Topographical map of entropy
changes. Brain regions are coloured according to their $\Delta h_{n}$ values.
(D) 5HT2A-R agonism changed the topographical distribution of entropy in a
heavily non-linear manner. Each circle indicates the averages of each region
across 1000 simulations.

To investigate this reconfiguration, we analysed the local changes in entropy
by plotting the entropy of the $n$th region in both PLA and 5HT2A conditions,
denoted by $h^{\text{5HT2A}}_{n}$ and $h^{\text{PLA}}_{n}$ respectively (Fig.
2D). Our results show that 5HT2A-R activation affects local entropy in a
highly non-linear manner, especially in regions with low baseline resting
state entropy. In particular, in such regions the effect of 5HT2A-R activation
could either increase or decrease entropy, such that the local entropy in the
5HT2A-R condition could not always be determined by the region’s baseline
entropy. This finding hints towards a more general theme: that local dynamical
properties alone are not able to explain the local changes in activity induced
by 5HT2A-R agonists like psychedelic drugs. We will explore this phenomenon in
depth in the following sections.

Next, we studied the effect of 5HT2A-R agonism on local entropy by considering
the relative change scores,

$\displaystyle\Delta
h_{n}=\frac{h^{\text{5HT2A}}_{n}-h^{\text{PLA}}_{n}}{h^{\text{PLA}}_{n}}\leavevmode\nobreak\
.$

Based on recent in vivo experiments with serotonergic psychedelics we expect
to find localized entropy increases on occipital, cingulate and parietal
regions [Lebedev2016, Schartner2017] as well as several changes on regions
belonging to Resting State Networks (RSNs).[Carhart-Harris2016,
Tagliazucchi2014, Lord2019] Accordingly, we split brain regions by standard
anatomical and functional groupings, and found that occipital and cingulate
areas tend to show a strong relative increase in their entropy, while parietal
and subcortical areas tend to show a decrease (Fig. 2A, Supp. Fig. 1).
Additionally, regions participating in the primary visual and Default Mode
(DMN) RSNs – including occipital and cingulate areas, respectively – showed a
marked tendency to increase their entropy, while regions participating in the
Frontoparietal Executive Network (FPN) showed the opposite behaviour. Further
discussion about the relation between these results and recent in vivo
experiments with psychedelic drugs is presented in the Discussion section.

Figure 3: Changes in local entropy are explained best by connectivity
strength, then receptor density. (A) Changes in entropy were overall
independent from receptor density, although (B) they were well predicted by
the connectivity strength of each region. We split into groups of low (blue),
intermediate (grey) and higher (red) connectivity strength, only those regions
with low and high connectivity strength were well predicted by their receptor
density. Regions with intermediate strength show no significant relationship
with receptor density, but are even more accurately predicted by their
strength. (C) Topographical localisation of the three groups, following the
same colour code. Low-strength regions are mainly located in the parietal
area, while high-strength ones are in occipital and cingulate areas.

### The topography of entropy changes is explained by local connectivity
strength and 5HT2A-R density

Given the results above, and given the potential clinical and neuroscientific
relevance of entropy reconfiguration in the psychedelic state, our next task
is to elucidate what neurobiological factors underlie such phenomenon. In this
section we investigate which structural and dynamical features of the model
are able to predict local entropy changes due to 5HT2A-R activation, and what
we can learn about the large-scale action mechanism of psychedelic drugs and
other 5HT2A-R agonists.

To explain the effects of 5HT2A-R activation, the first natural step is to
factor in the density of 5HT2A-R at each region. Somewhat surprisingly, at a
whole-brain level, receptor density was a very poor predictor of the entropy
change due to 5HT2A-R activation (Fig. 3A, $R^{2}=0.078\pm 0.004$). In
contrast, we found that the connectivity strength (based on the DTI human
connectome), defined as the sum of all the weighted links connecting a given
region, exhibits a strong correlation with local entropy changes (Fig. 3B,
$R^{2}=0.801\pm 0.006$).

By visual inspection, however, it is clear that the relationship between DTI-
informed connectivity strength and $\Delta h_{n}$ is linear for intermediate
values of strength, but it departs from linearity in both extremes of the
distribution. To investigate this phenomenon, we implemented a simple
optimisation algorithm to find the set of regions with the strongest linear
relationship between strength and $\Delta h_{n}$. This yielded a set of
regions with a higher strength-$\Delta h_{n}$ correlation ($R^{2}=0.886\pm
0.007$, as opposed to $R^{2}=0.801\pm 0.006$ for the whole-brain), and induced
a partition of all brain regions into three groups: those with intermediate
strength (IS, $76\pm 4$ brain regions, grey dots in Fig. 3A), high strength
(HS, $12\pm 3$ regions, red dots), and low strength (LS, $2\pm 2$ regions,
blue dots).

Building up on this partition, we studied the effect of receptor density on
those regions where entropy change could not be predicted from connectivity
strength – the HS and LS groups. We found a strong relationship between
receptor density and entropy changes for both the HS and LS groups, which
resulted in positive ($R^{2}=0.95\pm 0.02$) and negative ($R^{2}=0.999\pm
0.001$) correlation, respectively. On the contrary, receptor density does not
predict entropy changes for areas in the IS category. This shows a
complementary role of density and strength in LS, IS, and HS regions: entropy
changes in HS and LS regions strongly depend on the receptor density, but not
connectivity strength; while changes in IS regions depend on connectivity
strength, but not receptor density.

Overall, to assess the predictive power of connectivity strength and the
receptor density on local entropy changes, we built a linear mixed model for
$\Delta h_{n}$ using connectivity strength, 5HT2A-R density, and the
aforementioned three-way separation of brain regions as the predictor
variables. Together these variables explain $95.91\pm 0.01\%$ of the variance
of $\Delta h_{n}$, confirming that they provide an accurate model for
predicting the entropic effects of 5HT2A-R activation. This suggests that
psychedelic drugs and other 5HT2A-R agonists do not have a simple, localised
effect on brain activity, but instead amplify the fundamentally collective,
emergent properties of the brain as a complex system of interacting elements.

### The specific connectivity strength distribution explains relative changes
in entropy

As a final analysis, we set out to investigate exactly which topological
properties of the brain’s structural connectivity explain the observed changes
in entropy. With this exercise, we aim to answer two questions: whether any
property simpler than connectivity strength can explain the results; and
whether any property more complicated than connectivity strength is needed to
do so.

To this end, we ran further simulations of the DMF model using suitable null
network models of the human connectome (Fig. 4A), with the 5HT2A-R density map
held fixed. In particular, we used three null models designed to test
increasingly strict null hypotheses that preserved different network
attributes of the original connectome: i) the overall density and strength
(RAND, Fig. 4B), ii) the degree distribution (degree-preserving randomisation
[DPR], Fig. 4C); and iii) the strength distribution (strength-preserving
randomisation, [DSPR], Fig. 4D). We simulated the DMF model in these surrogate
networks, with and without neuromodulation, and computed the resulting entropy
changes $\Delta h_{n}$.

Our first result is that random and degree-preserving surrogate networks are
unable to reproduce the entropy changes observed in the human connectome: when
compared against the $\Delta h_{n}$ values obtained from the unperturbed
network, neither of them produce values close to the original (Figs. 4E, 4F).
This result asserts the findings in the previous section, showing that indeed
node strength plays an important role in shaping the global pattern of entropy
change associated with the action of psychedelics and other 5HT2A-R agonists.
Simpler topological features, like the degree distribution, are not enough to
explain such changes.

Perhaps more interestingly, the connectivity strength-preserving surrogate
networks not only reproduced the entropy changes, but did so almost exactly
(Fig. 4G). Furthermore, we repeated the analyses in the previous section
trying to predict $\Delta h_{n}$ using other local connectivity measures (e.g.
betweenness, eigenvector centrality), and none of them produced a model as
accurate as strength (Supp. Fig. 2). Together, these results show that, once
the receptor distribution is fixed, the network strength distribution is both
necessary and sufficient to explain the entropic effects of 5HT2A-R
activation. Of course, this is not to say that strength explains every aspect
of 5HT2A-R action – other topological network features are known to mediate
transitions of consciousness in other contexts,[Chennu2016] and investigating
which network properties explain high-order dynamical signatures[Rosas2019] of
psychedelics remains an exciting avenue of future work.

Figure 4: Relative changes in entropy are reproduced by a strength-preserving
null model of the connectome. (A)–(D) Connectivity matrices used to control
the role of local properties of the connectome on $\Delta h_{n}$. See main
text for the description of the matrices and randomisation algorithm. (E)–(G).
Scatter plots of $\Delta h_{n}$ for the human connectome against the three
null models. DSPR yielded almost the same results than the human connectome,
showing that only local network properties of human connectome are sufficient
to capture the effect of 5HT2A-R activation.

## Discussion

In this study we investigated the brain entropy changes induced by
serotonergic psychedelics by simulating whole-brain resting state activity
with and without 5HT2A-R activation. In contrast to empirical studies, which
usually only have access to coarse-grained fMRI or M/EEG data, our approach
allows us to study firing rates of brain regions and hence to directly assess
the effect of 5HT2A-R agonism at the level of neural population activity.

In agreement with the Entropic Brain Hypothesis,[Carhart-Harris2018] the model
shows a significant increase of global brain entropy, albeit in a highly
inhomogeneous and non-linear reconfiguration. The diversity of effects
stresses the importance of extending the scope of the EBH from a simple level-
based approach towards a multi-dimensional perspective, which might better
characterise the richness of 5HT2A-R activation effect both on the brain and
consciousness.

We propose a possible mechanism in the sense that the average increase in
neural entropy induced by serotonergic psychedelic drugs can be explained (in
statistical sense) in terms of quantities that provide a neurobiological
interpretation about how this phenomena could occur in the brain. These
results have important consequences for our understanding of consciousness,
neurodynamics, and the psychedelic state. In the following, we highlight some
of those implications, and propose possibilities for future work.

### 5HT2A-R-induced entropy changes are regionally heterogeneous

Our main result in this study is a mechanistic explanation of the main
prediction of the EBH that serotonergic psychedelics increase the entropy of
brain activity. However, one of the main takeaways is that this overall
increase is tied to a spatially heterogeneous reconfiguration, rather than a
globally consistent increase, in entropy.

Our simulation results show that 5HT2A-R activation triggers an entropy
increase in sensory areas, as observed in fusiform and olfactory cortices, as
well as the primary visual RSN and occipital regions (Supp. Fig. 1). This
agrees with the increased perceptual ‘bandwidth’ that characterises the
psychedelic state,[Bayne2018] as higher entropy might be related to richer
perceptual experience in a given moment of consciousness – potentially related
to reduced gating. As a matter of fact, all the primary visual RSN regions
(with the exception of the lingual area) are part of the high-strength group,
suggesting that the effect of psychedelics on perceptual experience might be
directly related to 5HT2A-R density in those regions. The localisation of the
entropy increases may relate to domain-specific changes in consciousness,
which could be interpreted as consistent with a recent dimensions of
consciousness characterisation of the psychedelic experience.[Bayne2016]

### Comparison to in vivo experiments with psychedelic drugs

Throughout the paper, we have focused on one particular signature of 5HT2A-R
agonists on the brain – a global increase in average entropy. But much more is
known about the effect of psychedelics on the brain, and studying these more
nuanced effects is key to understanding the rich phenomenology of the
psychedelic state. In this section we deepen this connection by providing a
more detailed comparison between the behaviour of the model and experimental
results with psychedelic drugs.

Our findings are in agreement with earlier studies where the effect of
psychedelic drugs on the topographical distribution of entropy-related
measures was correlated with subjective effects. For example, Schartner et
al.[Schartner2017] studied the effect of LSD, psilocybin, and ketamine on the
entropy rate of binarized MEG signals, and found localised increases in
entropy rate on occipital and parietal regions. Similarly, Lebedev et
al.[Lebedev2016] analysed the sample entropy of fMRI recordings of subjects
under the effect of LSD, finding localized increases on frontoparietal, medial
occipital, posterior and dorsal cingulate regions. Many of those regions
showed a consistent increase of their entropy with 5HT2A-R agonism in our
study (Fig. 2A).

Moreover, many of those regions actually belong to the high connectivity
strength group (c.f. Fig. 3), which suggests that their entropy increase in
experimental data might be directly related to the high 5HT2A-R density in
those regions.  We speculate that regions of the high-strength group may
systematically exhibit significant differences in local entropy on future
experiments. Together, these findings support the conclusion that the DMF
model, once optimised, can reproduce not only functional
connectivity,[Deco2018] but also some of the most salient localised entropy
increases found on in vivo human studies up to date.

Also, on a more fundamental level, our findings suggest that psychedelics
disrupt the functional organisation of the brain with an especially focal and
pronounced action on highly anatomically connected brain regions. In fact,
experimental and modelling evidence points towards a key role of connector
hubs on psychedelic and psychopatological states. [Muthukumaraswamy2013,
Muller2018] Perturbing the activity in these specific regions (via 5HT2A-R
agonism) could have particularly profound implications for the regularity of
brain activity and the quality of conscious experience.  However, sub cortical
regions may also play an important role on the modulation of brain activity,
conscious experience, and the psychedelic state [Barrett2020].

On a separate line of inquiry, there is strong evidence associating the DMN to
high-level cognitive functions, including the sense of space, time and
recognition of (self) boundaries. Disruptions to the DMN have been linked to
fundamental changes in experience, such as ego dissolution.[Tagliazucchi2016,
Carhart-Harris2018] Our simulation results show that 5HT2A-R activation
increases the entropy of all DMN regions (with the exception of the angular
cortex), which is consistent with the reported decrease in the DMN network
integrity[Muthukumaraswamy2013] induced by psychedelic drugs. In contrast,
low-level motor functions such as motor regulation remain largely unaffected
during the psychedelic state,[Carhart-Harris2019] which is consistent with the
modest entropy changes observed in the lingual and superior motor areas.

Finally, it is worth noting that both angular regions showed particularly
important decreases in local entropy, specially in the left hemisphere (Fig.
2A). Damage to this region is associated with impairments on language
processing, and electrical stimulation can induce out-of-body
experiences[Seghier2013] (both experiences can feature within psychedelic
states[Griffiths2006]). Since both structural damage and electrical
stimulation can be related to entropy reduction (e.g. by stimulus-driven
synchronisation), these experimental findings are consistent with the strong
entropy decrease in angular regions predicted by the model. We speculate that
the reduced entropy observed on the angular region induced by 5HT2A-R agonism
could be related to a smaller local bandwidth, which in turn might impair the
multi-modal and integrative functions of this region.[Seghier2013]

### Current limitations and future research

The approach employed in this work presents certain limitations related to
several aspects of the simulation and analysis. Acknowledging and
understanding these limitations can help us extend and improve our approach,
while introducing new questions in the field of psychedelic computational
neuroscience.

To our knowledge, the DMF model is the only whole-brain model that implements
neuromodulation and is capable of reproducing neuroimaging data in the placebo
and psychedelic states. Nonetheless, it makes some important simplifications
that are worth discussing. At the network level, the DTI-based connectome used
here is known to be incomplete, thus improvements could be made to the model
parameters of brain connectivity.[Markov2014] At the dynamical level, the DMF
model models neuronal populations as perfectly homogeneous within a given
region, and it is known that finer-grained local structure of certain brain
regions is likely to be key to explaining certain subjective effects of the
psychedelic state (e.g. lattice structure in the visual cortex and geometric
visual hallucinations[Ermentrout1979]). Additionally, the version of the DMF
model used here only considers 5HT2A-R agonism, while classic serotonergic
psychedelic drugs also have high binding affinities for other receptors (e.g.
in the case of LSD, the D1 and D2 dopamine receptors[Carhart-Harris2017]).

These simplifications do not prevent the model from reproducing statistical
features of brain signals under the placebo and 5HT2A-R conditions, but could
result in an inability to reproduce finer aspects of the dynamics of the
whole-brain activity in these conditions. Extending the model to reproduce
other dynamical signatures of psychedelics (like alpha suppression
[Timmermann2019] or reduced directed functional connectivity [Barnett2020])
constitutes a natural extension of this work,  with the recent example of
Kringelbach et al., [Kringelbach2020] where an extension of the DMF can
reproduce dynamical features of the brain in placebo and in the psychedelic
state

Another potentially fruitful line of future work involves making more detailed
comparisons with in vivo psychedelic neuroimaging data, and, potentially,
subjective experience reports. For example, one natural option would be to use
forward models of fMRI[Stephan2007] or M/EEG[Mosher1999] to bridge between the
firing rates produced by the DMF model and other data modalities, to produce
simulated data that is more directly comparable with available empirical data.

Another exciting possibility is to explore model parameters to examine
potential non-linearities and their implications for different relevant
aspects of brain function. For example, there are reasons to believe that the
dose-response relationship is non-linear for psychedelics and that over a
certain threshold dosage (level of 5HT2A-R stimulation) new subjective and
global brain function properties can appear.[Madsen2019]

Most interestingly, a potentially very useful extension of this work would be
to include individual subject-level connectome and receptor data to build
personalised models of response to psychedelic drugs. This would enable a much
more comprehensive modelling framework, capable of correlating structural
brain features with subjective experience reports. Such a framework could
potentially make individualised predictions of the action of serotonergic
psychedelics on specific individuals, aiding patient stratification and
treatment customisation.[Scott2019]

Finally, it is worth noting that all our analyses here are based on the
univariate statistics of individual brain regions, not including any
correlation or information flow between them. However, it is known that some
high-level subjective effects of psychedelics (such as complex
imagery[Carhart-Harris2016] and ego dissolution[Tagliazucchi2016]) are related
to network, as opposed to single region, dynamics. Therefore, building a
richer statistical description of the brain’s dynamics using recent
information-theoretic tools (such as multivariate extensions of mutual
information[Rosas2019] or Integrated Information Theory, IIT
[balduzzi2008integrated]) remains an exciting open problem.  In fact, recent
attempts of unifying IIT and the EBH in a single framework to understand the
effects of psychedelic drugs [Swanson2018] are helping to bridge the gap
between the univariate EBH and the multivariate IIT.

### Final remarks

In this paper we have provided the first mechanistic explanation of the neural
entropy increase elicited by psychedelic drugs, using a whole-brain dynamical
model with 5HT2A-R neuromodulation. Furthermore, we built a simple linear
model able to predict a region’s relative change in entropy from its local
5HT2A-R density and topological properties, showing that, somewhat
paradoxically, at a whole-brain level receptor density is a poor predictor of
5HT2A-R activation effect.

Key to developing this predictive model was a three-way partition of brain
regions according to their connectivity strength, suggesting a differentiated
action mechanism of 5HT2A-R agonists that depends on the local topology of
brain regions. In summary, our results suggest that the local changes in
entropy, as well as the global entropy increase, induced by 5HT2A-R activation
can be explained from a region-specific interplay between structural
connectivity and receptor density. Finally, controlled experiments with null
network models confirm that receptor density and connectivity strength are not
only necessary, but also sufficient, to explain the entropic effects of
5HT2A-R activation.

The spatially heterogeneous, complex nature of the observed effects of 5HT2A-R
activation opens a challenging problem for understanding the clinical and
scientific relevance of psychedelic drugs and their entropic effect.
Furthermore, it stresses the necessity of moving beyond the current
unidimensional approach to consciousness to a multi-dimensional one, that
better captures the phenomenological and neurodynamical richness of
psychedelic state.

## Methods

### Dynamic mean-field model with 5HT2A-R neuromodulation

The main computational tool used in this study is the Dynamic Mean-Field (DMF)
model by Deco et al.,[Deco2014, Deco2018] which consists of a set of coupled
differential equations modelling the average activity of one or more
interacting brain regions. In this model, each brain region $n$ is modelled as
two reciprocally coupled neuronal masses, one excitatory and one inhibitory,
and the excitatory and inhibitory synaptic currents $I^{(E)}$ and $I^{(I)}$
are mediated by NMDA and GABAA receptors respectively. Different brain regions
are coupled via their excitatory populations only, and the structural
connectivity is given by the matrix $C$.

The full model, including the neuromodulatory effect described below, is given
by

$\displaystyle I_{n}^{(E)}$
$\displaystyle=W_{E}I_{0}+w_{+}J_{NMDA}S_{n}^{(E)}+GJ_{NMDA}\sum_{p}^{N}C_{np}S_{p}^{(E)}-J^{FIC}_{n}S_{n}^{(I)}$
$\displaystyle I_{n}^{(I)}$
$\displaystyle=W_{I}I_{0}+J_{NMDA}S_{n}^{(E)}-S_{n}^{(I)}$ $\displaystyle
r_{n}^{(E)}$
$\displaystyle=F\left(I_{n}^{(E)}\right)=\frac{g_{NM}\;g_{E}\left(I_{n}^{(E)}-I_{thr}^{(E)}\right)}{1-\exp\left(-d_{E}\;g_{NM}\;g_{E}\left(I_{n}^{(E)}-I_{thr}^{(E)}\right)\right)}$
$\displaystyle r_{n}^{(I)}$
$\displaystyle=F\left(I_{n}^{(I)}\right)=\frac{g_{I}\left(I_{n}^{(I)}-I_{thr}^{(I)}\right)}{1-\exp\left(-d_{I}\;g_{I}\left(I_{n}^{(I)}-I_{thr}^{(I)}\right)\right)}$
$\displaystyle\frac{dS_{n}^{(E)}(t)}{dt}$
$\displaystyle=-\frac{S_{n}^{(E)}}{\tau_{NMDA}}+\left(1-S_{n}^{(E)}\right)\gamma
r_{n}^{(E)}+\sigma v_{n}(t)$ $\displaystyle\frac{dS_{n}^{(I)}(t)}{dt}$
$\displaystyle=-\frac{S_{n}^{(I)}}{\tau_{GABA_{A}}}+r_{n}^{(I)}+\sigma
v_{n}(t)$ $\displaystyle g^{NM}_{n}$
$\displaystyle=1+s_{E}d_{n}^{rec}\leavevmode\nobreak\ .$

Above, for each excitatory $(E)$ and inhibitory $(I)$ neural mass, the
quantities $I_{n}^{(E,I)}$, $r_{n}^{(E,I)}$, and $S_{n}^{(E,I)}$ represent its
total input current (nA), firing rate (Hz) and synaptic gating variable,
respectively. The function $F(\cdot)$ is the transfer function (or F-I curve),
representing the non-linear relationship between the input current and the
output firing rate of a neural population. Finally, $J^{FIC}_{n}$ is the local
feedback inhibitory control of region $n$, which is optimized[Deco2014] to
keep its average firing rate at approximately $3\,\mathrm{Hz}$, and $v_{n}$ is
uncorrelated Gaussian noise injected to region $n$. The interested reader is
referred to the original publication for further details.[Deco2014]

Key to this study is the recent extension of this model including
neuromodulatory effects.[Deco2018] In the equations above, $g^{NM}_{n}$ is a
neuromodulatory scaling factor modulating the F-I curve of all brain regions
in the model, which affects region $n$ depending on its density of the
receptor of interest, $d_{n}^{rec}$. Including neuromodulation, the free
parameters of the model are $G$, the global coupling parameter, and $s_{E}$,
the excitatory neuromodulatory gain, that are set to 2 and 0.2 respectively
following Deco et al.[Deco2018] The model was simulated using a standard
Euler-Maruyama integration method,[Kloeden2013] using the parameter values
shown in Table 1.

Table 1: Dynamic Mean Field (DMF) model parameters Parameter | Symbol | Value
---|---|---
External current | $I_{0}$ | 0.382 nA
Excitatory scaling factor for $I_{0}$ | $W_{E}$ | 1
Inhibitory scaling factor for $I_{0}$ | $W_{I}$ | 0.7
Local excitatory recurrence | $w_{+}$ | 1.4
Excitatory synaptic coupling | $J_{NMDA}$ | 0.15 nA
Threshold for $F(I_{n}^{(E)})$ | $I_{thr}^{(E)}$ | 0.403 nA
Threshold for $F(I_{n}^{(I)})$ | $I_{thr}^{(I)}$ | 0.288 nA
Gain factor of $F(I_{n}^{(E)})$ | $g_{E}$ | 310 nC-1
Gain factor of $F(I_{n}^{(I)})$ | $g_{I}$ | 615 nC-1
Shape of $F(I_{n}^{(E)})$ around $I_{thr}^{(E)}$ | $d_{E}$ | 0.16 s
Shape of $F(I_{n}^{(I)})$ around $I_{thr}^{(I)}$ | $d_{I}$ | 0.087 s
Excitatory kinetic parameter | $\gamma$ | 0.641
Amplitude of uncorrelated Gaussian noise $v_{n}$ | $\sigma$ | 0.01 nA
Time constant of NMDA  receptor | $\tau_{NMDA}$ | 100 ms
Time constant of GABA -A receptor | $\tau_{GABA_{A}}$ | 10 ms

### Parcellation, connectome, and 5HT2A receptor maps

In addition to the DMF equations above, to simulate whole-brain activity we
need three more ingredients: a suitable parcellation of the cortex into well-
defined regions, a structural connectivity matrix between those regions, and
the density map of a receptor of interest – in our case, the serotonin 2A
receptor.  We clarify that all the empirical inputs of the model were produced
elsewhere. For a detailed description of the procedure, please follow the
corresponding references cited below.

As a basis for our simulation, we used the Automated Anatomical Labelling
(AAL) atlas,[Tzourio-Mazoyer2002] a sulci-based parcellation of brain volume
registered in MNI space. The AAL parcellation specifies a partition of the
brain into 90 Regions Of Interest (ROIs), that provides sufficient level of
detail to obtain a picture of the topographical distribution of entropy, while
being suitable for computationally intensive simulations.

Structural connectivity data between the 90 AAL ROIs was obtained from 16
healthy subjects (5 female) using diffusion Magnetic Resonance Imaging (dMRI),
registered to the MNI space and parcellated according to the AAL atlas on the
subject’s native space. Then, for each subject the histogram of fiber
directions at each voxel was obtained, yielding an estimate of the number of
fibers passing trough each voxel. The weight of the connection between regions
$i$ and $j$ was defined as the proportion of fibers passing through region $i$
that reach region $j$, which (since the directionality of connections cannot
be determined using dMRI) yields a symmetric $90\\!\times\\!90$ Structural
Connectivity (SC) matrix. Finally, the SC matrix of each subject was
thresholded, pruning any connection lower than 0.1%, and SC matrices of all
subjects were averaged to obtain a single SC matrix used in all simulations.
This average matrix was kindly provided by Deco et al.[Deco2018] For a
detailed description of these methods see or refer to the Supplementary
Subsection 7.

For the density map of 5HT2A receptors, we used the Positron Emission
Tomography (PET) data made public by Beliveau et al.,[Beliveau2017] which can
be combined with the AAL parcellation to obtain an estimate of receptor
density in every AAL region. Further details can be found in the original
publications.[Deco2018, Beliveau2017]

Together with the DMF equations, these three elements fully specify our whole-
brain model and allow us to run simulations and obtain time series of
excitatory firing rates, $r_{n}^{(E)}$, which we save for further analysis.

### Estimating the differential entropy of firing rates

The result of each simulation is a set of 90 time series representing the
firing rate of each excitatory population ($r_{n}^{(E)}$), which we analysed
to measure the entropy of the activity of each ROI. For a continuous random
variable $X$ with associated probability distribution $p(x)$ and support
$\mathcal{X}$, its differential entropy is defined as[Cover2012]

$\displaystyle h(X)=-\int_{\mathcal{X}}p(x)\ln
p(x)\;\text{d}x\leavevmode\nobreak\ .$

Estimating differential entropy from data is in general a hard problem,
stemming (for the most part) from the difficulties in estimating a probability
density from samples [Lizier2014]. We found that the probability distributions
of firing rates $r_{n}^{(E)}$ were well approximated by a gamma distribution:
goodness-of-fit, measured with a standard Kolmogorov-Smirnov statistic, was
satisfactory and comparable across brain regions and conditions (see Supp.
Fig. 3).

This observation greatly simplifies the estimation of the differential entropy
of each region, which is now reduced to estimating the shape and scale
parameters of the gamma distribution, $k$ and $\theta$. Once these parameters
are estimated (in our case by standard maximum likelihood estimation), the
differential entropy of a gamma distribution can be computed analytically in
closed form as

$\displaystyle h(X)=k+\ln\theta+\ln\Gamma(k)+(1-k)\psi(k),$

where $\Gamma(\cdot)$ and $\psi(\cdot)$ are the standard gamma and digamma
functions, respectively.

### Linear models for $\Delta h_{n}$ and three-way separation of AAL regions

To explain the changes in differential entropy across regions and conditions,
we trained several linear mixed-effects models[McLean1991] using different
sets of covariates designed to test specific hypotheses about the relation
between entropy, 5HT2A-R density, and network topology.

First, we implemented a simple algorithm to find the group of regions where
the change in entropy depends linearly on region strength. To do so, we
designed an optimisation procedure to find the optimal linear fit between
$\Delta h_{n}$ and connectivity strength after the removal of some regions
from the fit by setting a threshold on $h^{\text{PLA}}_{n}$. The procedure
consisted in iteratively i) setting a threshold on $h^{\text{PLA}}_{n}$; ii)
fitting a linear model with strength as a regressor and $\Delta h_{n}$ as
target for those regions below the threshold; and iii) computing the goodness
of fit measured by $R^{2}$. This process is repeated for the whole range of
entropy values, yielding one $R^{2}$ value per threshold. The threshold that
maximises goodness of fit gives the optimal separation between regions,
separating those regions whose change in entropy linearly depends on
connectivity strength from those that it does not.

With this three-way grouping fixed, we investigated the power of several
dynamical and topological features to explain the changes in local entropy
($\Delta h_{n}$) by including them as covariates in the linear model. The
three-way separation of brain regions induces a model with 7 terms: a constant
term and 2 fixed-effect terms for each group of regions: one for the
connectivity strength and one for the receptor density. The model was fit by
maximum likelihood using off-the-shelf software, and explanatory power
measured using $R^{2}$.

### Null network models of the human connectome

Null network models of the human structural connectivity were used to evaluate
the role of the local connectivity properties[Rubinov2011] on the local
changes in entropy induced by 5HT2A-R activation. To this end, we applied
three different randomisation schemes to the structural connectivity in order
to produce suitable surrogate networks. These surrogates were designed to
preserve different network attributes of the original connectome: i) the
overall density and strength (RAND), ii) the degree distribution (degree-
preserving randomisation [DPR]); and iii) the strength distribution (strength-
preserving randomisation [DSPR]) After randomisation, the DMF model was run
with and without 5HT2A-R activation, and entropies estimated and analysed
following the same procedure as in the rest of the article. Every surrogate
model was run 120 times, and the results averaged across runs.

In addition to the network surrogates, we computed several topological
measures to include in the linear model reported in Fig. 3. These measures
included betweenness, eigenvector centrality, closeness, communicability,
PageRank index, and sub-graph centrality. All of these computations, including
the surrogate networks and the topology measures, were performed using the
Brain Connectivity Toolbox[Rubinov2010].

## Acknowledgements

R.H. is funded by CONICYT scholarship CONICYT-PFCHA/Doctorado
Nacional/2018-21180428. P.M. is funded by the Wellcome Trust (grant no.
210920/Z/18/Z). F.R. is supported by the Ad Astra Chandaria Foundation. R.C.
is supported by CONICYT-PAI Inserción 79160120, Proyecto REDES ETAPA INICIAL,
Convocatoria 2017 REDI170457 and Fondecyt Iniciación 2018 Proyecto 11181072.

## Author contributions statement

R.H. performed the numerical analysis and run the simulations. All the authors
of this article conceptualized this research, design the methodology. All the
authors carefully analysed the results, drafted and reviewed the manuscript.


#  Rotating Lee-Wick Black Hole and Thermodynamics

Dharm Veer Singha<EMAIL_ADDRESS>Sudhaker Upadhyayb,c,d,e
<EMAIL_ADDRESS><EMAIL_ADDRESS>Md Sabir Alif
<EMAIL_ADDRESS>aDepartment of Physics, GLA University, Mathura 281406
India  bDepartment of Physics, K. L. S. College, Nawada, Bihar 805110, India
cDepartment of Physics, Magadh University, Bodh Gaya, Bihar 824234, India
dInter-University Centre for Astronomy and Astrophysics (IUCAA), Pune,
Maharashtra 411007, India  eSchool of Physics, Damghan University, P.O. Box
3671641167, Damghan, Iran  fCenter for Theoretical Physics, Jamia Millia
Islamia, New Delhi, India

###### Abstract

We derive a singular solution for the rotating counterpart of Lee-Wick gravity
having a point source in a higher-derivative theory. We critically analyze the
thermodynamics of such a thermal system by evaluating mass parameters, angular
velocity, and Hawking temperature. The system follows the first law of
thermodynamics and leads to the expression of entropy. We further discuss the
stability and phase transition of the theory by evaluating heat capacity and
free energy. The phase transition occurs at the point of divergence and the
temperature is maximum. Remarkably, the black hole is unstable for a small
horizon radius and stable for a large horizon radius.

##  I Introduction

Quantizing the theory of gravity is one of the most challenging tasks of
theoretical physics. It is known that the gravity theory with forth derivative
terms is renormalizable which in fact violate the unitarity of the theory [1].
However, the fourth derivative terms lead to the negative-norm states in the
physical spectrum of the theory. This justifies instabilities in the classical
solutions. A consistent theory of quantum gravity theory which agrees with
both unitarity and renormalizability requires a modification in Einstein
theory of gravity. One of the possible solutions to this problem is to
consider a local gravitational theory with higher than fourth derivative
terms. Such theories with higher derivatives terms have very interesting
properties as they are superrenormalizable [2]. Also, these theories are
unitary in the Lee-Wick formalism for the massive complex poles [3, 4]. Hence,
these theories remove the conflict between unitarity and renormalizability in
quantum gravity [5, 6, 7, 8, 9, 10, 11]. Also, the proof of unitarity of these
models do not require any fine-tuning and therefore prove themselves as a
better alternative to the non-local theories [12, 13, 14, 15, 16].

The singular solutions (black holes) are the most fundamental objects in
gravity theory as they provide powerful probes to the study of various aspects
of the theory. Therefore, the search of black-hole solutions in the higher-
derivative gravity theories are of considerable interests. A static vacuum
solution of an approximate equation of motion of the Lee-Wick gravity is
studied recently where the solution is generated by a point-like massive
source [17]. The authors of Ref. [17] discussed different scenarios by
comparing the mass of the source with a critical value. They also shed light
on the thermodynamics of these black holes and their evaporation process. The
thermodynamics of the black hole plays a vital role in understanding the
quantum theory of gravity [18, 19, 20, 21, 22]. Hawking and Page had found
that the black hole solutions in asymptotically $AdS$ space have thermodynamic
properties [23]. The rotational counterpart of these gravity theory remains
unexplored. This provides us an opportunity to generalize the work and this is
the motivation of present study.

In this paper, we obtain a rotating Lee-Wick black hole solution by using a
modified Newman-Janis algorithm [24, 25, 26, 27]. To do so, we first convert
static spherically symmetric Lee-Wick black hole solution into its rotational
counterpart, we transform the coordinates of the Lee-Wick metric from the
Boyer-Lindquist to the Eddington-Finkelstein. The numerical analysis leads to
the possibilities to find non-vanishing values of rotation parameter which
corresponds minimum metric function [28, 29]. We also investigate the basic
thermodynamics properties of rotating Lee-Wick black hole. We first derive
mass parameter and angular velocity. Moreover, we compute the Hawking
temperature of rotating lee-Wick black hole in terms of temperature of Kerr
black hole. A comparative analysis suggests that rotating Lee-Wick black holes
with small rotation parameter are comparatively hotter. Also, the rotation
parameter plays more significant role for small black holes. In order to study
the stability of such black holes, we need to calculate heat capacity of the
black hole and check the sign of heat capacity and free energy. Here, we find
three phase transitions. For the stationary case when rotation parameter is
switched-off, the phase transition occurs only once. In contrast, when
rotation parameter acquires non-zero values, the phase transition occurs.
Also, in stationary case, the larger black holes are stable.

The paper is organized as follows. In section II, we obtain a rotational
counterpart of Lee-Wick black hole solution. In section III, we discussed the
thermodynamics of rotational Lee-Wick black hole and studies their local and
global stability. Finally, we conclude the work in the section V.

##  II Generating rotating Lee-Wick black hole solutions

. We start with the action of the minimal theory having unitarity with the
Lee-Wick prescription [17]

$S=\frac{1}{8\pi G_{N}}\int
d^{4}x\sqrt{-g}\left[R+{\Lambda^{-4}}G_{\mu\nu}\Box R^{\mu\nu}\right],$ (1)

where $R$ is the Ricci scalar, $\Lambda$ represents the UV scale of the theory
and $G_{\mu\nu}$ is the Einstein tensor. Varying the action (1), we obtain the
equations of motion (EoM)

$\left(1+\frac{\Box^{2}}{\Lambda^{4}}\right)G_{\mu\nu}+O(R^{2}_{\mu\nu})=8\,\pi\,G_{N}T_{\mu\nu}$
(2)

where energy momentum tensor (EMT)
$T_{\mu\nu}=\left(1+{\Box^{2}}/{\Lambda^{4}}\right){\tilde{T}}_{\mu\nu}$.
Ignoring higher-order ($O(R^{2}_{\mu\nu})$) terms, this EoM takes following
form [17]:

$G_{\mu\nu}=8\,\pi\,G\,\left(1+\frac{\Box^{2}}{\Lambda^{4}}\right)^{-1}{T}_{\mu\nu}.$
(3)

The only non-vanishing component of EMT is $T_{t}^{t}=-M\delta(x)$ and, for
this , the effective EMT is given by

${\tilde{T}}^{\mu}_{\nu}=\text{diag}\left(-{\tilde{\rho}},{\tilde{P}_{r}},{\tilde{P}_{\theta}},{\tilde{P}_{\theta}}\right),$
(4)

where ${\tilde{\rho}}$, ${\tilde{P}_{r}}$ and ${\tilde{P}_{\theta}}$ are the
effective energy density, effective radial pressure and effective tengential
pressure, repectively. The effective energy density is calculated as [17]

$\displaystyle\tilde{\rho}=\left(1+\frac{\Box^{2}}{\Lambda^{4}}\right)^{-1}$
$\displaystyle=$
$\displaystyle-M\delta(x)=M\int\frac{dk^{3}}{(2\pi)^{3}}\frac{e^{ikx}}{1+k^{4}/\Lambda^{4}}$
(5) $\displaystyle=$ $\displaystyle\frac{M\Lambda^{2}}{4\pi
r}\sin\frac{r\Lambda}{\sqrt{2}}e^{-\frac{r\Lambda}{\sqrt{2}}}.$

Now, we write the static and spherically symmetric Lee-Wick solution which
matches with the form of Schwarzschild metric as follows [17],

$\displaystyle
ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}d\Omega_{2}^{2}\quad\text{with}\quad
f(r)=1-\frac{2m(r)}{r}.$ (6)

Here $m(r)$ refers to effective mass and depends on the radial coordinate
($r$) only due to the spherical symmetry. The expression for mass function
$m(r)$ is given by [17]

$\displaystyle m(r)=M\left[1-e^{-\frac{\Lambda
r}{\sqrt{2}}}\left(1+\frac{\Lambda r}{\sqrt{2}}\right)\cos\left(\frac{\Lambda
r}{\sqrt{2}}\right)+\frac{\Lambda r}{\sqrt{2}}\sin\left(\frac{\Lambda
r}{\sqrt{2}}\right)\right],$ (7)

where $\Lambda$ is the UV scale and $M$ corresponds to the mass of a static
point-like source. It should be noted that this is not an exact solution of
the equations of motion of the model (1) but rather a solution to the
‘approximate’ equations of motion (2) in [17] when the quadratic and higher-
order terms in the Ricci tensor are neglected.

Upon effective mass expansion around the centre ($r=0$), one finds de Sitter
core [17]

$\displaystyle
ds^{2}\approx-\left(1-\frac{\Lambda_{eff}r^{2}}{3}\right)dt^{2}+\frac{1}{\left(1-\frac{\Lambda_{eff}r^{2}}{3}\right)}dr^{2}+r^{2}d\Omega^{2}_{2}$
(8)

where $\Lambda_{eff}=\sqrt{2}\Lambda^{3}M$ is the effective cosmological
constant. The metric (8) describes a regular spacetime and thus curvature is
singularity free.

Now, we proceed to obtain the rotating counterpart of the metric (6). In order
to do so, we follow the method proposed originally by Newman and Janis [24]
and further modified by Azreg-Ainou [26]. In order to convert static
spherically symmetric Lee-Wick black hole solution into its rotational
counterpart, we first transform the coordinates of the spacetime metric (6)
from the Boyer-Lindquist ($t,r,\theta,\phi$) to the Eddington-Finkelstein
($u,r,\theta,\phi$) as following:

$\displaystyle du=dt-\frac{dr}{f(r)}.$ (9)

This leads to the final spacetime metric in the following form:

$\displaystyle ds^{2}=-fdu^{2}-2dudr+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta
d\phi^{2}.$ (10)

It is known that contravariant components of the metric tensor in the advanced
null Eddington-Finkelstein coordinates can be expressed by the null tetrad of
the form

$\displaystyle
g^{\mu\nu}=-l^{\mu}n^{\nu}-l^{\nu}n^{\mu}+m^{\mu}\bar{m}^{\nu}+m^{\nu}\bar{m}^{\mu},$
(11)

Following the method as discussed in Ref. [27], we first write the metric
coordinates in null tetrad and do complex coordinate transformations in the
$u-r$ plane as

$\displaystyle u\rightarrow u-ia\cos\theta,\quad r\rightarrow r-ia\cos\theta.$
(12)

This converts the metric functions into a new form: $f(r)\rightarrow
F(r,a,\theta)$, $r^{2}\rightarrow\Sigma(r,a,\theta)$. Furthermore, null
tetrads also take the following new form:

$\displaystyle l^{\mu}=\delta_{r}^{\mu}\ ,\quad
m^{\mu}=\frac{1}{\sqrt{2\Sigma}}\left[\delta_{\theta}^{\mu}+ia\sin\theta(\delta_{u}^{\mu}-\delta_{r}^{\mu})+\frac{i}{\sin\theta}\delta_{\phi}^{\mu}\right]\
,$ $\displaystyle
n^{\mu}=\sqrt{\frac{G}{F}}\delta_{u}^{\mu}-\frac{1}{2}F\delta_{r}^{\mu}\
,\quad\bar{m}^{\mu}=\frac{1}{\sqrt{2\Sigma}}\left[\delta_{\theta}^{\mu}-ia\sin\theta(\delta_{u}^{\mu}-\delta_{r}^{\mu})-\frac{i}{\sin\theta}\delta_{\phi}^{\mu}\right].$
(13)

Then we can rewrite the contravariant components of the metric tensor
$g^{\mu\nu}$ by using (11) as

$\displaystyle g^{\mu\nu}=\left(\begin{array}[]{c c c
c}a^{2}\sin^{2}\theta/\Sigma&-\sqrt{G/F}-a^{2}\sin^{2}\theta/\Sigma&0&a/\Sigma\\\
-\sqrt{G/F}-a^{2}\sin^{2}\theta/\Sigma&G+a^{2}\sin^{2}\theta/\Sigma&0&-a/\Sigma\\\
0&0&1/\Sigma&0\\\ a/\Sigma&-a/\Sigma&0&1/(\Sigma\sin^{2}\theta)\\\
\end{array}\right).$ (18)

The covariant components read

$\displaystyle g_{\mu\nu}=\left(\begin{array}[]{c c c
c}-F&-\sqrt{\frac{F}{G}}&0&a\left(F-\sqrt{\frac{F}{G}}\right)\sin^{2}\theta\\\
-\sqrt{\frac{F}{G}}&0&0&a\sqrt{\frac{F}{G}}\sin^{2}\theta\\\ 0&0&\Sigma&0\\\
a\left(F-\sqrt{\frac{F}{G}}\right)\sin^{2}\theta&a\sqrt{\frac{F}{G}}\sin^{2}\theta&0&\sin^{2}\theta\left[\Sigma+a^{2}\left(2\sqrt{\frac{F}{G}}-F\right)\sin^{2}\theta\right]\\\
\end{array}\right).$ (23)

The last step of the algorithm is to turn back from the Eddington-Finkelstein
coordinates to the Boyer-Lindquist coordinates by using the following
coordinate transformations:

$\displaystyle du=dt+\lambda(r)dr\ \qquad\text{and}\qquad\quad
d\phi=d\phi+\chi(r)dr\ .$ (24)

The transformation functions $\lambda(r)$ and $\chi(r)$ are found due to the
requirement that, except the coefficient $g_{t\phi}$ ($g_{\phi t}$), all the
non-diagonal components of the metric tensor are equal to zero. Thus

$\displaystyle\lambda(r)=-\frac{r^{2}+a^{2}}{f(r)^{2}+a^{2}}\
,\qquad\qquad\qquad\qquad\chi(r)=-\frac{a}{f(r)r^{2}+a^{2}},$ $\displaystyle
F(r,\theta)=\frac{(f^{2}+a^{2}\cos^{2}\theta)\Sigma}{(r^{2}+a^{2}\cos^{2}\theta)^{2}},\qquad\text{and}\qquad
G(r,\theta)=\frac{f^{2}+a^{2}\cos^{2}\theta}{\Sigma}.$ (25)

Consequently, we obtain the following line-element for the rotational
counterpart of the Lee-Wick black hole solution

$\displaystyle ds^{2}$ $\displaystyle=$
$\displaystyle-\left(1-\frac{2Mr-2Me^{-\frac{\Lambda
r}{\sqrt{2}}}\left[\left(1+\frac{\Lambda r}{\sqrt{2}}\right)\cos\frac{\Lambda
r}{\sqrt{2}}+\sin\frac{\Lambda
r}{\sqrt{2}}\right]}{r^{2}+a^{2}\cos^{2}\theta}\right)dt^{2}+\frac{r^{2}+a^{2}\cos^{2}\theta}{\Delta}dr^{2}$
(26) $\displaystyle-$ $\displaystyle
2a\sin^{2}\theta\left(\frac{2Mr-2Me^{-\frac{\Lambda
r}{\sqrt{2}}}\left[\left(1+\frac{\Lambda r}{\sqrt{2}}\right)\cos\frac{\Lambda
r}{\sqrt{2}}+\sin\frac{\Lambda
r}{\sqrt{2}}\right]}{r^{2}+a^{2}\cos^{2}\theta}\right)d\phi
dt+(r^{2}+a^{2}\cos^{2}\theta)d\theta^{2}$ $\displaystyle+$
$\displaystyle\sin^{2}\theta\left[r^{2}+a^{2}+a^{2}\sin^{2}\theta\left(\frac{2Mr-2Me^{-\frac{\Lambda
r}{\sqrt{2}}}\left[\left(1+\frac{\Lambda r}{\sqrt{2}}\right)\cos\frac{\Lambda
r}{\sqrt{2}}+\sin\frac{\Lambda
r}{\sqrt{2}}\right]}{r^{2}+a^{2}\cos^{2}\theta}\right)\right]d\phi^{2},$

with

$\displaystyle\Delta=r^{2}-2Mr+a^{2}+2Me^{-\frac{\Lambda
r}{\sqrt{2}}}\left[\left(1+\frac{\Lambda r}{\sqrt{2}}\right)\cos\frac{\Lambda
r}{\sqrt{2}}+\sin\frac{\Lambda r}{\sqrt{2}}\right],$ (27)

where $a$ is the rotation parameter of black hole. Arising from the UV scale
theory, $\Lambda$ measures the potential deviation from the Kerr metric. The
obtained black hole solution (26) is independent of $t,\phi$, which implies
that it admits two Killing vectors given by $\eta^{\mu}=\delta^{\mu}_{t}$ and
$\xi^{\mu}=\delta^{\mu}_{\phi}$. The horizons of rotating Lee-Wick black holes
are solution of equation
$(\eta^{\mu}\xi_{\mu})^{2}-(\eta_{\mu}\eta^{\mu})(\xi_{\nu}\xi^{\nu})=0$,
which gives

$\displaystyle g_{t\phi}^{2}-g_{tt}g_{\phi\phi}=0.$ (28)

Here, this yields

$\displaystyle r^{2}-2Mr+a^{2}+2Me^{-\frac{\Lambda
r}{\sqrt{2}}}\left[\left(1+\frac{\Lambda r}{\sqrt{2}}\right)\cos\frac{\Lambda
r}{\sqrt{2}}+\sin\frac{\Lambda r}{\sqrt{2}}\right]=0.$ (29)

This is a transcendental equation which cannot be solved analytically. The
numerical analysis of the $\Delta=0$ on varying the angular momentum $a$ with
fixed mass $M=1$ is depicted in the Fig. 1. The numerical analysis of
$\Delta=0$ reveals that it is possible to find non-vanishing value of angular
momentum ($a$) and UV scale parameter $x=\Lambda r/\sqrt{2}$ for which metric
function $\Delta$ is minimum, i.e, $\Delta(x_{+})=0$ and this will give two
real roots $x_{+}$ and $x_{-}$which correspond to the Cauchy horizon and event
horizon.

---
Figure 1: The plot of metric function $\Delta(r)$ as the function of dimensionless parameter $x=\Lambda r/\sqrt{2}$ for different values of angular momentum $a$ with fixed value of $\Lambda=4\sqrt{2}\pi$.  $a$ | $x_{+}$ | $x_{-}$ | $\delta=x_{+}-x_{-}$
---|---|---|---
0.10 | 1.85 | 2.31. | 0.46
0.20 | 1.92 | 2.24 | 0.32
0.272 | 2.096 | 2.096 | 0
Table 1: The table for horizon radius for different value of rotation
parameter $a$ with fixed value of $\Lambda=4\sqrt{2}\pi$.

From the Fig. (1) and Table 1, we can see that the horizon of the black hole
decreases with increasing rotation parameter and at $a=0.268$ two horizons
reduces to one, i.e., $r_{\pm}=r_{E}$ such that
$\Delta(r_{E})=\Delta^{\prime}{r_{E}}=0$.

##  III Thermodynamics

In this section, we analyse the thermodynamics of the rotating Lee-Wick black
hole described by metric (26). Let us begin by deriving the mass by $g^{rr}=0$
as following:

$\displaystyle
M_{+}=\frac{\left(a^{2}+r_{+}^{2}\right)}{2r_{+}\left(1-e^{-\frac{\Lambda
r_{+}}{\sqrt{2}}}\left(\left(1+\frac{\Lambda
r_{+}}{\sqrt{2}}\right)\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]+\frac{\Lambda r_{+}}{\sqrt{2}}\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)\right)}.$ (30)

In the absence of rotation parameter ($a$), the mass (30) reduces to the mass
of Lee-Wick black hole [17]; however this reduces to mass of Kerr black hole
in the limit of $\Lambda\to 0$.

For the stationary axially symmetric metric like (26), we have two associated
Killing vectors $\eta^{\mu}=\partial_{t}$ corresponding to the time
translational symmetry along the $t$-axis and $\xi^{\mu}=\partial_{\phi}$
corresponding to the rotational symmetry about $\phi$-axis. In therms of these
killing vectors a Killing field $\chi{{}^{\mu}}$ can be expressed as
$\chi{{}^{\mu}}=\eta^{\mu}+\Omega\xi^{\mu}$, where $\Omega$ refers to the
angular velocity of the metric. Being a null vector at the event horizon (i.e.
$\chi{{}^{\mu}}\chi_{\mu}=0$), the Killing field decides the angular velocity
and leads to [30, 31]

$\displaystyle g_{tt}+2\Omega g_{t\phi}+\Omega^{2}g_{\phi\phi}=0,$ (31)

where angular velocity is given by

$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle-\frac{g_{t\phi}}{g_{\phi\phi}}\pm\sqrt{\left(\frac{g_{t\phi}}{g_{\phi\phi}}\right)^{2}-\frac{g_{tt}}{g_{\phi\phi}}},$
(32) $\displaystyle=$
$\displaystyle\frac{-a\left[\Sigma-(\Delta-a^{2}\sin^{2}\theta)\right]}{\Sigma^{2}+a^{2}\sin^{2}\theta\left[2\Sigma-(\Delta-a^{2}\sin^{2}\theta)^{2}\right]}$
$\displaystyle\pm$
$\displaystyle\frac{\Sigma\Delta^{\frac{1}{2}}}{\sin\theta\left[(a^{2}\sin^{2}\theta+\Sigma)^{2}-\Delta
a^{2}\sin^{2}\theta\right]}.$

In terms of horizon radius, the angular velocity reads

$\displaystyle\Omega_{+}=-\left.\frac{g_{t\phi}}{g_{\phi\phi}}\right|_{r=r_{+}}=\frac{a}{r_{+}^{2}+a^{2}}.$
(33)

Before computing temperature of the black holes, it is worth to first work out
the surface gravity at the event horizon for the obtained solution (26). This
gives

$\displaystyle\kappa=\sqrt{-\frac{\nabla_{\mu}\chi{{}_{\nu}}\nabla{{}^{\mu}}\chi{{}^{\nu}}}{2}}=\frac{\Delta^{\prime}(r_{+})}{2(r_{+}^{2}+a^{2})}.$
(34)

Here ′ denotes differentiation with respect to horizon radius. Utilizing the
standard relation of Surface gravity and Hawking temperature
$T_{+}=\kappa/2\pi$, we obtain the Hawking temperature of the rotating Lee-
Wick black hole metric (26) as

$\displaystyle
T_{+}=T_{+}^{Kerr}\left[\frac{e^{-\frac{r_{+}\Lambda}{\sqrt{2}}}\left(\left(1+\frac{\Lambda
r_{+}}{\sqrt{2}}\right)\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]+r_{+}\left(\frac{\Lambda
r_{+}}{\sqrt{2}}+\frac{\Lambda}{2}\left(\frac{a^{2}+r_{+}^{2}}{a^{2}-r_{+}^{2}}\right)\right)\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)-1}{e^{-\frac{r_{+}\Lambda}{\sqrt{2}}}\left(\left(1+\frac{\Lambda
r_{+}}{\sqrt{2}}\right)\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]+\frac{\Lambda r_{+}}{\sqrt{2}}\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]-r_{+}\right)}\right],$ (35)

where $T_{+}^{Kerr}={\left(r_{+}^{2}-a^{2}\right)}/{\left(4\pi
r_{+}(r_{+}^{2}+a^{2})\right)}$ is the temperature of the Kerr black hole. The
temperature (35) reduces to the temperature of the Kerr black hole in the
limit of $\Lambda\to 0$.

---
Figure 2: The plot of Hawking temperature $T_{+}$ versus horizon radius $x_{+}$ for different values of rotation parameter $a$ with fixed Value of $\Lambda=4\sqrt{2}\pi$.  $a$ | $0.30$ | $0.50$ | $0.70$ | $0.90$
---|---|---|---|---
$T_{max}$ | 0.0279 | 0.0275. | 0.0268 | 0.0260
$x_{c}$ | 3.750 | 3.797 | 3.816 | 3.872
Table 2: The table for temperature for different value of rotation paramete
$a$ with fixed value of $\Lambda=4\sqrt{2}\pi$.

From the Fig. 2 and table (2), it is clear that the black hole becomes hotter
for small rotation parameter and the critical radius increases with the
rotation parameter $(a)$. The maximum Hawking temperature decreases for higher
values of the critical radius $(r_{c})$. Thus, This signifies that the effect
of rotation parameter is more significant for small black holes.

This black hole can be considered as a thermodynamic system only if the
quantities associated with it must satisfies the first-law of thermodynamics

$\displaystyle dM_{+}=T_{+}dS_{+}+\Omega_{+}da.$ (36)

The validity of above relation leads to the expression for entropy

$S=-4\pi\int
dr_{+}\frac{(a^{2}+r_{+}^{2})e^{-r_{+}\Lambda/\sqrt{2}}}{e^{-\frac{\Lambda
r_{+}}{\sqrt{2}}}\left(\left(2+{\sqrt{2}}\,r_{+}\Lambda\right)\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]+{\sqrt{2}}\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)}.$ (37)

From this this expression it is clear that the entropy does not depends upon
roation parameter ($a$). The entropy is the function of horizon radius which
is depicted in Fig. (3)

$\begin{array}[]{c c c c}\includegraphics[width=325.215pt]{s.eps}\end{array}$

Figure 3: The plot of entropy $S_{+}$ vs $x_{+}$ with fixed value of
$\Lambda=4\sqrt{2}\pi$.

The definition of first-law of thermodynamics is worth mentioning to check the
validity of these thermodynamic quantities. It is easy to verify the complete
form of the first-law of thermodynamics as

$dM_{+}=T_{+}dS_{+}+\Omega_{+}dJ$ (38)

Substituting the value of mass, temperature and entropy from the Eqs. (30),
(35) and (37), respectively, given in Eq. (38), we can see that the rotating
Lee-Wick black hole follows the first-law of thermodynamics.

##  IV Local and Global Stability

Finally, we study the thermodynamic stability of rotating Lee-Wick black hole
by estimating the heat capacity. The stability of the black hole can be
estimated from sign of the heat capacity ($C_{+}$) as positive heat capacity
reflects the stable black hole and negative reflects the unstable one. The
heat capacity of the black hole can be calculated by [30, 31, 32, 33]

$\displaystyle
C_{+}=\frac{dM_{+}}{dT_{+}}=\frac{dM_{+}}{dr_{+}}\frac{dr_{+}}{dT_{+}}$ (39)

substituting the mass and temperature from Eq. (30) and Eq. (35) into Eq. (39)
we get

$\displaystyle
C_{+}=2\pi\frac{\left(r^{2}-a^{2}\right)\left(r^{2}+a^{2}\right)^{2}Y}{X},$
(40)

with

$\displaystyle Y$ $\displaystyle=$ $\displaystyle e^{-\frac{\Lambda
r_{+}}{\sqrt{2}}}\left(\left(1+\frac{\Lambda
r_{+}}{\sqrt{2}}\right)\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]+r_{+}\left(\frac{\Lambda
r_{+}}{\sqrt{2}}+\frac{\Lambda}{2}\left(\frac{a^{2}+r_{+}^{2}}{a^{2}-r_{+}^{2}}\right)\right)\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)-1,$ $\displaystyle X$ $\displaystyle=$
$\displaystyle
2\left(a^{4}+4a^{2}r_{+}^{2}-r_{+}^{4}\right)\left[\left(e^{\frac{{\Lambda
r_{+}}}{\sqrt{2}}}-\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)^{2}+\frac{{2\Lambda
r_{+}}}{\sqrt{2}}\left(e^{\frac{{\Lambda
r_{+}}}{\sqrt{2}}}-\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)\right.$ (41)
$\displaystyle\left.\left(\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]+\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)-2\left(\frac{\Lambda
r_{+}}{\sqrt{2}}\right)^{2}\right]+4\left(r_{+}^{2}+a^{2}\right)^{2}\left[-\left(\frac{\Lambda
r_{+}}{\sqrt{2}}\right)^{2}\right.$ $\displaystyle+$
$\displaystyle\left.\left(\frac{\Lambda
r_{+}}{\sqrt{2}}\right)^{7}\left(e^{-\frac{{\Lambda
r_{+}}}{\sqrt{2}}}-\left(\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]-\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)\right)\right].$

$\begin{array}[]{c c c c}\includegraphics[width=325.215pt]{4.eps}\end{array}$

Figure 4: The plot of specific heat $C_{+}$ as the function of dimensionless
parameter $x=\Lambda r/\sqrt{2}$ for different values of angular momentum $a$
with fixed value of $\Lambda=4\sqrt{2}\pi$.

In order to study the stability of rotating Lee-Wick black hole, we plot the
heat capacity with respect to dimensionless parameter $x$ in Fig 5. From Eq.
(40), it can be observed that the specific heat diverges at $r_{+}=3.75$ this
discontinuity of the heat capacity represents a point of phase transition.

$\begin{array}[]{c c c
c}\includegraphics[width=325.215pt]{free.eps}\end{array}$

Figure 5: The plot of free energy $C_{+}$ as the function of dimensionless
parameter $x=\Lambda r/\sqrt{2}$ for different values of angular momentum $a$
with fixed value of $\Lambda=4\sqrt{2}\pi$.

Now, we study the behaviour of Gibbs free energy for the global stability of
black hole thermodynamics. The definition of Gibbs free energy is given by

$\displaystyle G_{+}=M_{+}-T_{+}S_{+}.$ (42)

Substituting the mass, temperature and entropy from Eqs. (30), (35) and (37)
into (42) we get

$\displaystyle G_{+}$ $\displaystyle=$
$\displaystyle\frac{\left(a^{2}+r_{+}^{2}\right)}{2r_{+}\left(1-e^{-\frac{\Lambda
r_{+}}{\sqrt{2}}}\left(\left(1+\frac{\Lambda
r_{+}}{\sqrt{2}}\right)\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]+\frac{\Lambda r_{+}}{\sqrt{2}}\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)\right)}+\left(\frac{r_{+}^{2}-a^{2}}{4\pi(r_{+}^{2}+a^{2})}\right)$
(43)
$\displaystyle\left(\frac{e^{-\frac{r_{+}\Lambda}{\sqrt{2}}}\left(\left(1+\frac{\Lambda
r_{+}}{\sqrt{2}}\right)\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]+r_{+}\left(\frac{\Lambda
r_{+}}{\sqrt{2}}+\frac{\Lambda}{2}\left(\frac{a^{2}+r_{+}^{2}}{a^{2}-r_{+}^{2}}\right)\right)\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)-1}{e^{-\frac{r_{+}\Lambda}{\sqrt{2}}}\left(\left(1+\frac{\Lambda
r_{+}}{\sqrt{2}}\right)\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]+\frac{\Lambda r_{+}}{\sqrt{2}}\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]-r_{+}\right)}\right)$ $\displaystyle\int
dr_{+}\frac{(a^{2}+r_{+}^{2})e^{-r_{+}\Lambda/\sqrt{2}}}{e^{-\frac{\Lambda
r_{+}}{\sqrt{2}}}\left(\left(2+{\sqrt{2}}\,r_{+}\Lambda\right)\cos\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]+{\sqrt{2}}\sin\left[\frac{\Lambda
r_{+}}{\sqrt{2}}\right]\right)}.$

We analyse the stability of the black hole by studying the nature of free
energy which is plotted for the different values of rotation parameter $(a)$
in Fig. 5. Here we see that the free energy has a maximum value where the
specific heat diverges (see Fig. 4) which can be identified as the extremal
points of the Hawking temperature (see Fig. 2). We also notice that the
rotating black hole is unstable for small values of $r_{+}$ and stable for
large horizon radius. It means that the large black hole is stable.

##  V Conclusions

We have extended the investigation of static and spherically symmetric Lee-
Wick solution. In fact, we have obtained a rotating Lee-Wick black hole
solution by exploiting a corrected Newman-Janis algorithm. Specifically, in
order to derive a rotational counterpart, we first considered static
spherically symmetric Lee-Wick black hole solution and made a coordinate
transformation from the Boyer-Lindquist coordinates to the Eddington-
Finkelstei coordinates. Moreover, we have written the contravariant components
of the metric tensor in the advanced null Eddington-Finkelstei coordinates in
the null tetrad form. After that we have performed a complex coordinate
transformations followed by the coordinate transformation in order to retain
the metric in Boyer-Lindquist coordinates. Since we are left with a
transcendental equation which can not be solved analytically. Hence, we have
tried to do a numerical analysis by comparing the metric function $\Delta$
equal to zero which reveals the possibilities for finding the non-zero values
of angular momentum and cosmological constant.

Also, we have discussed the basic thermodynamics properties of this rotating
Lee-Wick black hole. In order to do so, first we have evaluated a mass
parameter and angular velocity. Following standard black hole chemistry from
areal-law, we have derived the Hawking temperature of rotating Lee-Wick black
hole and expressed in terms of temperature of Kerr black hole. In order to
study the behaviour of black hole with horizon radius, we have plotted a graph
and found that the temperature of rotating Lee-Wick black holes is increasing
for smaller rotation parameter. We also found that the effects of rotation
parameter become more significant for relatively smaller black holes. The
stability of this black holes is studied by evaluating heat capacity of the
black hole and checking the sign of heat capacity. Remarkably, the system
exhibits various phase transitions. A comparative analysis is done and found
that for the stationary case (when rotation parameter is turned-off) the phase
transition occurs only once for small sized black holes. In contrast, when
rotation parameter acquires non-zero values, the phase transition occurred at
various points. The graph signified that in contrast to the stationary case,
we have found a stable region for the larger rotating black holes also. The
results of the paper are quite interesting and yields new insight for future
investigations.

## References

  * [1] R. Utiyama and B.S. DeWitt, J. Math. Phys. 3 (1962) 608.
  * [2] M. Asorey, J. L. Lopez, and I.L. Shapiro, Int. J. Mod. Phys. A 12, 5711 (1997).
  * [3] L. Modesto and I.L. Shapiro, Phys. Lett. B 755 (2016) 279.
  * [4] L. Modesto, Nucl. Phys. B 909 (2016) 584.
  * [5] T. D. Lee and G. C. Wick, Nucl. Phys. B 9, 209 (1969).
  * [6] T. D. Lee and G. C. Wick, Phys. Rev. D 2, 1033 (1970).
  * [7] R. E. Cutkosky, P. V. Landshoff, D. I. Olive and J. C. Polkinghorne, Nucl. Phys. B 12, 281 (1969).
  * [8] L. Modesto, Nucl. Phys. B 909, 584 (2016) [arXiv:1602.02421 [hep-th]].
  * [9] B. L. Giacchini, arXiv:1609.05432 [hep-th].
  * [10] G. P. de Brito, P. I. C. Caneda, Y. M. P. Gomes, J. T. G. Junior and V. Nikoofard, arXiv:1610.01480 [hep-th].
  * [11] A. Accioly, B. L. Giacchini and I. L. Shapiro, arXiv:1610.05260 [gr-qc].
  * [12] I. L. Shapiro, Phys. Lett. B 744, 67 (2015).
  * [13] L. Modesto and I. L. Shapiro, Phys. Lett. B 755, 279 (2016).
  * [14] T. D. Lee and G. C. Wick, Nucl. Phys. B 9, 209 (1969).
  * [15] T. D. Lee and G. C. Wick, Phys. Rev. D 2, 1033 (1970).
  * [16] R. E. Cutkosky, P. V. Landshoff, D. I. Olive and J. C. Polkinghorne, Nucl. Phys. B 12, 281 (1969).
  * [17] C. Bambi, L. Modesto and Y. Wang, Phys. Lett. B 764, 306 (2017).
  * [18] S. Carlip, Int. J. Mod. Phys. D. 23, 1430023 (2014).
  * [19] A. Bekenstein, Nuovo Cimento Lett. 4, 99 (1972).
  * [20] S. W Hawking, Comm. in Math. Phys. 43, 199 (1975)
  * [21] J. D. Bekenstein, Physical Review D. 9, 3292 (1974).
  * [22] R. M Wald, Living Reviews in Relativity. 4,1 (2001)
  * [23] S. W . Hawking and D. N. Page, Commun, Math. Phys. 87, 577 (1983).
  * [24] E. Newman and A. Janis, J. Math. Phys. 6, 915 (1965).
  * [25] D. Hansen and N. Yunes, Phys. Rev. D 88, 104020 (2013).
  * [26] M. Azreg-Aınou, Phys. Rev. D 90 (2014) 064041.
  * [27] B. Toshmatov, Z. Stuchlík, B. Ahmedov, Eur. Phys. J. Plus 132 (2017) 2.
  * [28] F. Ahmed, D. V. Singh and S. G. Ghosh, arXiv:2008.10241 [gr-qc].
  * [29] F. Ahmed, D. V. Singh and S. G. Ghosh, arXiv:2002.12031 [gr-qc].
  * [30] M. S. Ali and S. G. Ghosh, arXiv:1906.11284 [gr-qc].
  * [31] M. S. Ali and S. G. Ghosh, Phys. Rev. D 99 (2019) 024015.
  * [32] D. V. Singh, M. S. Ali and S. G. Ghosh, Int. J. Mod. Phys. D 27 (2018) no.12, 1850108.
  * [33] P. Chaturvedi, N. K. Singh and D. V. Singh, Int. J. Mod. Phys. D 26 (2017) no.08, 1750082.
  * [34] R. M. Wald, Phys. Rev. D 43, R3427 (1993).

##  Appendix A Verification of black hole solution by equation of motion

In order to compute the energy momentum tensor, we use the following tetrad

$\displaystyle e^{\mu}_{t}=\frac{1}{\Sigma\sqrt{\Delta}}(r^{2}+a^{2},0,0,a),$
(44) $\displaystyle e^{\mu}_{r}=\frac{\sqrt{\Delta}}{\Sigma}(0,1,0,0),$ (45)
$\displaystyle e^{\mu}_{\theta}=\frac{1}{\Sigma}(0,0,1,0),$ (46)
$\displaystyle e^{\mu}_{\phi}=-\frac{1}{\Sigma\sin\theta}(r^{2}+a^{2},0,0,a).$
(47)

The non-zero component of energy momentum tensor are given by

$\displaystyle{\tilde{\rho}(r)}=-{\tilde{P}}_{r}=\frac{M\Lambda^{2}}{4\pi
r}\sin\frac{r\Lambda}{\sqrt{2}}e^{-\frac{r\Lambda}{\sqrt{2}}},$ (48)
$\displaystyle{\tilde{P}}_{\theta}={\tilde{P}}_{\phi}=\frac{M\Lambda^{2}e^{-\frac{r\Lambda}{\sqrt{2}}}}{16\pi
r}\left(\sqrt{2}\cos\frac{r\Lambda}{\sqrt{2}}+(2-\sqrt{2}r)\sin\frac{r\Lambda}{\sqrt{2}}-8\sin\frac{r\Lambda}{\sqrt{2}}\frac{M-rM^{\prime}}{2M-r}\right).$
(49)

The equation of motion of the action (1) can be written as

$\displaystyle G_{t}^{t}$ $\displaystyle=$
$\displaystyle-\frac{1}{4\Sigma^{2}(f(r,\theta)E(r,\theta)+B(r,\theta))^{2}}\Big{[}-E^{2}\,f\,\frac{\Sigma}{\Delta}\Delta^{\prime}-E^{2}f\left(\frac{\Sigma}{\Delta}\right)^{\prime
2}\Delta$ (50) $\displaystyle-$ $\displaystyle EB^{2}\Delta^{\prime
2}\frac{\Sigma}{\Delta}-E\Delta
B^{2}\left(\frac{\Sigma}{\Delta}\right)^{\prime 2}+2\Sigma\Delta
E^{\prime\prime}fE+2E\left(\frac{\Sigma}{\Delta}\right)^{\prime\prime}\Sigma
B^{2}+2E\Delta^{\prime\prime}\Sigma B^{2}$ $\displaystyle+$ $\displaystyle
2E^{2}f\left(\frac{\Sigma}{\Delta}\right)^{\prime\prime}\Sigma+2\left(\frac{\Sigma}{\Delta}\right)^{2}\Delta
E^{\prime\prime}B^{2}+2E^{2}f\Delta^{\prime\prime}\Sigma-E^{2}f\Delta^{\prime}\left(\frac{\Sigma}{\Delta}\right)^{\prime}\frac{\Sigma}{\Delta}$
$\displaystyle-$ $\displaystyle
E^{2}f\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta^{\prime}\Delta-E\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta^{\prime}\Delta
B^{2}-E\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta^{\prime}\frac{\Sigma}{\Delta}B^{2}+\left(\frac{\Sigma}{\Delta}\right)\Delta
B^{\prime}E$ $\displaystyle+$
$\displaystyle\left(\frac{\Sigma}{\Delta}\right)^{2}\Delta
B^{\prime}E-\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta^{\prime}E^{\prime}B^{2}+2\Sigma\Delta
E^{\prime\prime}fE+\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Sigma
E^{\prime}fE$ $\displaystyle+$ $\displaystyle\Delta^{\prime}\Sigma
E^{\prime}fE-2\Sigma\Delta
BB^{\prime}E^{\prime}-2\left(\frac{\Sigma}{\Delta}\right)^{2}\Delta
BB^{\prime}E^{\prime}-\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta^{2}E^{\prime}fE$
$\displaystyle-$
$\displaystyle\left(\frac{\Sigma}{\Delta}\right)^{\prime}\left(\frac{\Sigma}{\Delta}\right)\Delta
E^{\prime}B^{2}-\Delta^{\prime}\left(\frac{\Sigma}{\Delta}\right)^{2}E^{\prime}fE-\Delta^{\prime}\Sigma
E^{\prime}B^{2}-E^{\prime 2}\Sigma\Delta f$ $\displaystyle-$
$\displaystyle\Delta^{\prime}\left(\frac{\Sigma}{\Delta}\right)^{2}E^{\prime}B^{2}-E^{\prime
2}\left(\frac{\Sigma}{\Delta}\right)^{2}\Delta f\Big{]},$ $\displaystyle
G^{t}_{\phi}$ $\displaystyle=$
$\displaystyle-\frac{1}{4\Sigma^{2}(f(r,\theta)E(r,\theta)+B(r,\theta))^{2}}\Big{[}\Delta^{\prime
2}\frac{\Sigma}{\Delta}B^{3}+\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)^{2}\Delta
B^{3}$ (51) $\displaystyle+$ $\displaystyle
EfB\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)^{2}\Delta+EfB\Delta^{\prime
2}\frac{\Sigma}{\Delta}-2EfB\left(\frac{d^{2}}{d\theta^{2}}\frac{\Sigma}{\Delta}\right)\Sigma-2B^{\prime\prime}\Sigma\Delta
fE$ $\displaystyle-$ $\displaystyle
2B^{\prime\prime}\frac{\Sigma^{2}}{\Delta}fE-\frac{dB}{d\theta}\frac{dC}{d\theta}\Sigma
fE-2EfB\Delta^{\prime\prime}\Sigma-B^{\prime}\Delta^{\prime}\Sigma
fE-2EfB\Delta^{\prime\prime}\Sigma$ $\displaystyle-$ $\displaystyle
f^{\prime}\Sigma\Delta
BE^{\prime}+\frac{dA}{d\theta}\frac{dB}{d\theta}\frac{\Sigma^{2}}{\Delta}E-\frac{dA}{d\theta}\frac{dE}{d\theta}\frac{\Sigma^{2}}{\Delta}B+B^{\prime}C^{\prime}\Delta^{2}fE$
$\displaystyle-$ $\displaystyle\frac{dB}{d\theta}\frac{dC}{d\theta}\Sigma
B^{2}+E^{\prime}\Sigma\Delta
fB^{\prime}+\frac{dB}{d\theta}\frac{d\Delta}{d\theta}\frac{\Sigma^{2}}{\Delta^{2}}fE-B^{\prime}\Delta^{\prime}\Sigma
B^{2}+\frac{dE}{d\theta}\frac{dB}{f}\frac{\Sigma^{2}}{\Delta}{d\theta}$
$\displaystyle-$ $\displaystyle 2B^{\prime\prime}\Sigma\Delta
B^{2}-2\frac{d^{2}B}{d\theta^{2}}B^{2}\frac{\Sigma^{2}}{\Delta}+EfBC^{\prime}\Delta^{\prime}\Delta+EfB\frac{dC}{d\theta}\frac{d\Delta}{d\theta}\frac{\Sigma}{\Delta}-2\Delta^{\prime\prime}\Sigma
B^{3}$ $\displaystyle+$ $\displaystyle C^{\prime}\Delta^{\prime}\Delta
B^{3}+B\Sigma\Delta
B^{\prime\prime}+B\frac{\Sigma^{2}}{\Delta}\left(\frac{dB}{d\theta}\right)^{2}+B^{\prime}C^{\prime}\Delta^{2}B^{2}+\frac{dB}{d\theta}\frac{d\Delta}{d\theta}\frac{\Sigma^{2}}{\Delta^{2}}$
$\displaystyle+$ $\displaystyle A^{\prime}\Sigma\Delta
EB^{\prime}-2\frac{d^{C}}{d\theta^{2}}\Delta B^{3}\Big{]},$ $\displaystyle
G^{r}_{r}$ $\displaystyle=$
$\displaystyle\frac{1}{4\Sigma^{2}(f(r,\theta)E(r,\theta)+B(r,\theta))^{2}}\Big{[}\Delta^{2}B^{\prime
2}B^{2}+2\Delta^{\prime}\Delta B^{3}B^{\prime}-\Sigma
B^{2}\left(\frac{dB}{d\theta}\right)^{2}$ (52) $\displaystyle-$
$\displaystyle\Sigma
E^{2}\left(\frac{dA}{d\theta}\right)^{2}-2\left(\frac{d\Delta}{d\theta}\right)\frac{\Sigma}{\Delta}B^{3}\left(\frac{dB}{d\theta}\right)-\Sigma
f^{2}\left(\frac{dE}{d\theta}\right)^{2}+\Delta^{2}B^{\prime 2}fE$
$\displaystyle+$
$\displaystyle\Delta^{2}E^{\prime}f^{\prime}B^{2}+2\Delta^{\prime}\Delta
BB^{\prime}fE+\Sigma\left(\frac{dE}{d\theta}\right)\left(\frac{df}{d\theta}\right)fE-4\Sigma
EB\left(\frac{df}{d\theta}\right)\left(\frac{dB}{d\theta}\right)$
$\displaystyle-$ $\displaystyle 4\Sigma
Bf\left(\frac{dB}{d\theta}\right)\left(\frac{dE}{d\theta}\right)-2\left(\frac{d\Delta}{d\theta}\right)\left(\frac{dB}{d\theta}\right)\frac{\Sigma}{\Delta}BfE+3\Sigma\left(\frac{dE}{d\theta}\right)\left(\frac{df}{d\theta}\right)B^{2}$
$\displaystyle+$ $\displaystyle
3\Sigma\left(\frac{dE}{d\theta}\right)\left(\frac{dA}{d\theta}\right)B^{2}+3\Sigma\left(\frac{dB}{d\theta}\right)^{2}fE-\left(\frac{d\Delta}{d\theta}\right)\frac{\Sigma}{\Delta}E^{2}f\left(\frac{df}{d\theta}\right)$
$\displaystyle-$
$\displaystyle\left(\frac{d\Delta}{d\theta}\right)\frac{\Sigma}{\Delta}EB^{2}\left(\frac{dA}{d\theta}\right)-\left(\frac{d\Delta}{d\theta}\right)\frac{\Sigma}{\Delta}A^{2}E\left(\frac{dE}{d\theta}\right)-\left(\frac{d\Delta}{d\theta}\right)\frac{\Sigma}{\Delta}fB^{2}\left(\frac{dE}{d\theta}\right)$
$\displaystyle+$
$\displaystyle\Delta^{2}E^{\prime}f^{\prime}fE+\Delta^{\prime}\Delta
E^{2}f^{\prime}f+\Delta^{\prime}\Delta Ef^{\prime}B^{2}+\Delta^{\prime}\Delta
fE^{\prime}E+\Delta^{\prime}\Delta fE^{\prime}B^{2}$ $\displaystyle+$
$\displaystyle 4\Sigma B\left(\frac{dB}{d\theta^{2}}\right)fE+4\Sigma
E^{2}\left(\frac{d^{2}f}{d\theta^{2}}\right)f+2\Sigma
fB^{2}\left(\frac{d^{E}}{d\theta^{2}}\right)+2\Sigma
f^{2}\left(\frac{d^{2}E}{d\theta^{2}}\right)E$ $\displaystyle+$ $\displaystyle
2\Sigma E\left(\frac{d^{2}A}{d\theta^{2}}\right)B^{2}+4\Sigma
B^{2}\left(\frac{d^{2}B}{d\theta^{2}}\right)\Big{]},$ $\displaystyle
G_{\theta}^{\theta}$ $\displaystyle=$
$\displaystyle\frac{1}{4\Sigma^{2}(f(r,\theta)E(r,\theta)+B(r,\theta))^{2}}\Big{[}\Sigma
B^{\prime 2}B^{2}-\Sigma E^{2}A^{\prime 2}-2C^{\prime}\Delta
B^{3}B^{\prime}-\Sigma f^{2}E^{\prime 2}$ (53) $\displaystyle+$ $\displaystyle
2\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)\frac{\Sigma}{\Delta}B^{3}\left(\frac{dB}{d\theta}\right)+\left(\frac{dE}{d\theta}\right)\left(\frac{df}{d\theta}\right)B^{2}\left(\frac{\Sigma}{\Delta}\right)^{2}\left(\frac{dB}{d\theta}\right)^{2}fE$
$\displaystyle+$
$\displaystyle\left(\frac{\Sigma}{\Delta}\right)BfE\left(\frac{dB}{d\theta}\right)\left(\frac{dC}{d\theta}\right)+\Sigma
E^{\prime}f^{\prime}fE-4\Sigma Ef^{\prime}BB^{\prime}-4\Sigma
BB^{\prime}fE^{\prime}$ $\displaystyle-$ $\displaystyle
2\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta
BB^{\prime}fE+\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)\frac{\Sigma}{\Delta}E^{2}\frac{df}{d\theta}f+\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)\frac{\Sigma}{\Delta}E\frac{df}{d\theta}B^{2}$
$\displaystyle+$
$\displaystyle\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)\frac{\Sigma}{\Delta}A^{2}\frac{dE}{d\theta}E+\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)\frac{\Sigma}{\Delta}\frac{dE}{d\theta}fB^{2}+\left(\frac{\Sigma}{\Delta}\right)^{2}\frac{dE}{d\theta}\frac{df}{d\theta}+3\Sigma
E^{\prime}f^{\prime}B^{2}$ $\displaystyle+$ $\displaystyle 3\Sigma B^{\prime
2}fE-\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta
E^{2}ff^{\prime}-\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta
Ef^{\prime}B^{2}-\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta f^{\prime
2}EE^{\prime}$ $\displaystyle-$
$\displaystyle\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta
fE^{\prime}B^{2}+2\Sigma fE^{\prime\prime}B^{2}+2\Sigma
f^{2}E^{\prime\prime}E+2\Sigma E^{2}f^{\prime\prime}f$ $\displaystyle+$
$\displaystyle 4\Sigma B^{3}B^{\prime\prime}+4\Sigma
BB^{\prime\prime}fE+\frac{\Sigma^{2}}{\Delta^{2}}B^{2}\left(\frac{dB}{d\theta}\right)^{2}\Big{]},$
$\displaystyle G_{\phi}^{\phi}$ $\displaystyle=$
$\displaystyle\frac{1}{4\Sigma^{2}(f(r,\theta)E(r,\theta)+B(r,\theta))^{2}}\Big{[}2\Delta^{\prime\prime}\Sigma
fB^{2}+2E\left(\frac{d^{2}}{d\theta^{2}}\frac{\Sigma}{\Delta}\right)f^{2}+2E\Delta^{\prime\prime}\Sigma
f^{2}$ (54) $\displaystyle-$ $\displaystyle\Delta^{\prime
2}\frac{\Sigma}{\Delta}fB^{2}-\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)^{2}\Delta
fB^{2}-E\Delta^{\prime
2}\frac{\Sigma}{\Delta}f^{2}-E\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)^{2}\Delta
f^{2}$ $\displaystyle-$ $\displaystyle
f^{\prime}\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta^{2}fE-2f^{\prime}\Sigma\Delta
BB^{\prime}+\frac{df}{d\theta}\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)\Sigma
B^{2}-\frac{df}{d\theta}\frac{d\Delta}{d\theta}f\left(\frac{\Sigma}{\Delta}\right)^{2}\frac{dB}{d\theta}$
$\displaystyle+$
$\displaystyle\frac{df}{d\theta}\left(\frac{\Sigma}{\Delta}\right)^{2}\Delta
B^{2}+f^{\prime}\Delta^{\prime}\Sigma
B^{2}+2\left(\frac{\Sigma}{\Delta}\right)^{2}\frac{d^{2}f}{d\theta^{2}}\Delta
B^{2}+2\frac{d^{2}f}{d\theta^{2}}\Sigma\Delta
B^{2}+2f^{\prime\prime}\Sigma\Delta fB^{2}$ $\displaystyle+$ $\displaystyle
2\left(\frac{d^{2}}{d\theta^{2}}\frac{\Sigma}{\Delta}\right)\Sigma
fB^{2}+\frac{d^{2}f}{d\theta^{2}}\frac{\Sigma^{2}}{\delta}fE+f^{\prime\prime}\Sigma\Delta
fE+\frac{df}{d\theta}\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)\Sigma
fE+f^{\prime}\Delta^{\prime}\Sigma fE$ $\displaystyle-$ $\displaystyle
f^{\prime}\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta^{2}B^{2}-f^{\prime
2}\Sigma\Delta E+B^{\prime 2}\Sigma\Delta
f-\frac{df}{d\theta}\frac{d\Delta}{d\theta}\left(\frac{\Sigma}{\Delta}\right)^{2}B^{2}-\left(\frac{d}{d\theta}\frac{\Sigma}{\Delta}\right)\frac{d\Delta}{d\theta}\frac{\Sigma}{\Delta}fB^{2}$
$\displaystyle-$
$\displaystyle\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta^{\prime}\Delta
fB^{2}-E\frac{d\Delta}{d\theta}\frac{\Sigma}{\Delta}f^{2}-E\left(\frac{\Sigma}{\Delta}\right)^{\prime}\Delta^{\prime}\Delta
f^{2}-\left(\frac{df}{d\theta}\right)^{2}\frac{\Sigma^{2}}{\Delta}E+$
$\displaystyle\left(\frac{dB}{d\theta}\right)^{2}\frac{\Sigma^{2}}{\Delta}f\Big{]},$

where

$\displaystyle f=1-\frac{2rm(r)}{\Sigma},$ (55) $\displaystyle
B=\frac{r\,a\,m(r)\sin^{2}\theta}{\Sigma},$ (56)
$\displaystyle\Sigma=r^{2}+a^{2}\,\cos^{2}\theta,$ (57)
$\displaystyle\Delta=r^{2}-2Mr+a^{2}+2Me^{-\frac{\Lambda
r}{\sqrt{2}}}\left[\left(1+\frac{\Lambda r}{\sqrt{2}}\right)\cos\frac{\Lambda
r}{\sqrt{2}}+\sin\frac{\Lambda r}{\sqrt{2}}\right],$ (58) $\displaystyle
E=r^{2}+a^{2}+\frac{2m(r)a^{2}\sin^{2}\theta}{\Sigma}.$ (59)

By substituting the above value of $f,B,\Sigma,\Delta$ and $E$ in Eqs. (50),
(51), (52), (53) and (54) and then substituting the values of
$G_{tt},G_{rr},G_{\theta\theta}\,\text{and}\,G_{\phi\phi}$ with
$T_{tt},T_{rr},T_{\theta\theta}$ and $T_{\phi\phi}$ in Eq. (3), Einstein field
equation gets satisfied.

##  Appendix B Black Hole Entropy

Next, let us focus to an important thermodynamic quantity, namely, entropy
$S_{+}$ of the black hole in term of horizon radius. The entropy of the
rotating black hole black hole is calculated by using the Wald entropy
relation [34] as

$\displaystyle S_{+}=\int\frac{\delta I}{\delta
R_{ab\alpha\beta}}\epsilon^{a\alpha}\epsilon^{b\beta}\sqrt{h}d^{2}\Omega_{2},\qquad\text{with}\qquad
I=\int d^{4}x\sqrt{-g}\mathcal{L},$ (60)

where $S_{+}$, $I$, $\epsilon^{ab}$, $h$ and $R_{ab\alpha\beta}$ are the
entropy, action, bi-normal to the horizon, induced metric on the horizon and
the Riemann tensor, respectively. To calculate the entropy an antisymmetric
second rank tensor $\epsilon_{ab}$ is constructed so that
$\epsilon_{rt}=\epsilon_{tr}=1$. The specific form of Lagrangian and its
derivative with respect to Riemann tensor are given by

$\displaystyle\mathcal{L}=\frac{1}{16\pi}R_{ab\alpha\beta}g^{b\alpha}g^{a\beta}\qquad\mbox{and}\qquad\frac{\partial\mathcal{L}}{\partial
R_{ab\alpha\beta}}=\frac{1}{16\pi}\frac{1}{2}(g^{a\alpha}g^{b\beta}-g^{b\alpha}g^{a\beta}).$
(61)

The definition (60) corresponding to the above expressions leads to the Wald
entropy as

$\displaystyle
S_{+}=\frac{1}{8}\int\frac{1}{2}(g^{a\alpha}g^{b\beta}-g^{b\alpha}g^{a\beta})\epsilon_{ab}\epsilon^{\alpha\beta}\sqrt{h}d\theta
d\phi,$ $\displaystyle S_{+}=\frac{1}{8}g^{tt}g^{rr}\,\sqrt{h}\,d\theta
d\phi=\frac{1}{4}\int r^{2}\,d\theta d\phi,$ $\displaystyle S_{+}=\pi
r_{+}^{2}=\frac{A}{4}.$ (62)

This expression of entropy resembles with the standard Bekenstein-Hawking
area-law.