Circular spectropolarimetric sensing of vegetation in the field; possibilities
for the remote detection of extraterrestrial life

C.H. Lucas Patty1,2,3,*, Inge Loes ten Kate4, Wybren Jan Buma5, Rob J.M. van
Spanning1, Gábor Steinbach6, Freek Ariese7, Frans Snik8

1 Amsterdam Institute for Molecules, Medicine and Systems (AIMMS), VU
Amsterdam, De Boelelaan 1108, 1081 HZ Amsterdam, The Netherlands
2 Institute of Plant Biology, Biological Research Centre, Hungarian Academy of
Sciences, P.O. Box 521, H-6701 Szeged, Hungary
3 Biofotonika R&D Ltd, Szeged, Hungary
4 Faculty of Earth Sciences, Utrecht University, Budapestlaan 4, 3584 CD
Utrecht, The Netherlands
5 HIMS, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The
Netherlands
6 Institute of Biophysics, Biological Research Centre, Hungarian Academy of
Sciences, P.O. Box 521, H-6701 Szeged, Hungary
7 LaserLaB, VU Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The
Netherlands
8 Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The
Netherlands

<EMAIL_ADDRESS>
## Abstract

Homochirality is a generic and unique property of all biochemical life and the
fractional circular polarization of light it induces therefore constitutes a
potentially unambiguous biosignature. However, while high-quality circular
polarimetric spectra can be easily and quickly obtained in the laboratory,
accurate measurements in the field are much more challenging due to large
changes in illumination and target movement. In this study we have measured
various targets in the field, up to distances of a few kilometers, using the
dedicated circular spectropolarimeter TreePol. We show how photosynthetic life
can readily be distinguished from abiotic matter. We underline the potential
of circular polarization signals as a remotely accessible means to
characterize and monitor terrestrial vegetation, e.g. for agriculture and
forestry. Additionally, we discuss the potential of circular polarization for
the remote detection of extraterrestrial life.

##  1 Introduction

Exoplanetary science has been advancing rapidly over the last few decades.
Estimates are that on average with every star of the 100-400 billion stars in
our galaxy there is at least one planet [1]. Additionally, estimates on the
occurrence of rocky exoplanets in the habitable zone range from 2 % - 20 % per
stellar system ([2] and references therein). One of the primary motivations of
exoplanetary science is to find signs of life beyond Earth and therefore
substantial research is devoted towards the identification, verification and
validation of remotely detectable spectral characteristics of exoplanets as a
diagnostic of the presence of life [3, 4, 5, 6, 7, 8, 9].

Research on remotely detectable biosignatures has focused on detecting
particular atmospheric constituents such as liquid H2O, which indicates
possible planetary habitability. Additionally, the simultaneous presence of
gases in thermodynamic disequilibrium, such as O2 and CH4, has been
investigated [10, 11]. However, detection of these gases is not free of false-
positive scenarios [12, 13, 14, 15]. Other suggested remotely detectable
(surface) biosignatures include the well-known red edge effect, resulting from
terrestrial vegetation [16, 17], and pigment signatures resulting from other
organisms [18], with the risk of possible false-positives by mineral
reflectance.

Chiral molecules in their simplest form exist in a left-handed (L-) and a
right-handed (D-) version that are not superimposable. Unlike abiotic
chemistry, different classes of biological molecules tend to select
exclusively only one of these configurations (homochirality). For instance,
all living organisms mainly synthesize amino acids in the L-configuration
while sugars are predominantly synthesized in the D-configuration.
Homochirality also manifests itself within biological macromolecules and
biomolecular architectures. The $\alpha$-helix, for example, a common
secondary structure of proteins, is exclusively right-hand-coiled.
Homochirality is required for processes ranging from proper enzymatic
functioning to self replication. The latter is of importance as it implies the
prerequisite of homochirality for life [19, 20, 21]. It is likely that
homochirality is a universal feature of life and therefore may serve as a
unique and unambiguous biosignature.

When interacting with light, chiral molecules cause optical rotation (which is
the rotation of the orientation of the linear polarization plane upon
interaction with a sample) and exhibit circular dichroism (which is the
differential absorption of left- or right-handed circularly polarized incident
light). These phenomena are most evident upon interaction with polarized
light, but can also lead to the circular polarizance of unpolarized light,
such as emitted from a star [22]. Starlight reflecting off an extraterrestrial
body with homochiral molecules and molecular systems on its surface will thus
carry this information, which, in principle, can be sensed remotely [23, 24,
25, 26, 9, 27].

Polarimetry is in general advantageous for both the detection and the
characterization of exoplanets. Polarimetry enables one to enhance the
contrast between the starlight and the very dim light reflected off its
planets (which often is very strongly linearly polarized) and thus offers the
potential to characterize the atmosphere and surface of an exoplanet [28, 29,
30]. Induced linear polarization can potentially be a biosignature. Both
biotic and abiotic matter, however, can create large amounts of linear
polarization [31, 32]. As such, linear polarization might offer complementary
information in addition to scalar reflectance signatures, but caution should
be exercised when using it for distinguishing biotic from abiotic components.

As homochirality is exclusive to life, the circular polarization it induces
constitutes a strong and potentially unambiguous biosignature. While also
abiotic matter can create circular polarization (e.g. through multiple
scattering), these signals generally are orders of magnitude smaller than
those created by biological homochirality and have a much smoother and broader
spectral shape [26, 24, 30].

Circular polarizance is detected by measuring the induced fractional circular
polarization of unpolarized incident light. For the simplest molecules, the
polarizance will be strongest where there is higher spectral absorbance and
the polarizance will be opposite in handedness to that of the absorbing
biological molecule. Excitonic interactions further enhance the magnitude of
these signals, but a particularly interesting phenomena arises when observing
large and dense supramolecular systems and complexes. These systems can induce
anomalously large circular polarization signals with optical activity even
outside of the absorbance bands, and has been labeled Polymer and Salt Induced
(PSI) circular polarization [33, 34, 35]. While the psi-type circular
polarizance can thus extend beyond the absorption bands, the polarization
signal will still primarily occur where the total absorption is very high,
such that it must be detected within a relatively weak signal, one of the
primary challenges for remote observation.

The amplitude of the signal is also strongly wavelength dependent. Amino
acids, for instance, are strongly polarizing in the vacuum-ultraviolet region
and thin films (500 nm) of (homochiral) alanine can readily reach a
polarizance of 0.6$\%$ at 180 nm [36]. Outside the water absorption band
($>$190 nm), the maximum polarizance of the same film is, however, only
0.07$\%$. While abiotic materials can also create circular polarization, often
through multiple scattering [37], the risk of a false positive scenario is
smaller. Various minerals consistently show a much weaker and spectrally
different circular polarization signal [26, 24]. Also simulations of clouds on
Earth or Venus show circular polarizance, but this is small and has a much
broader spectral shape [30]. Additionally, circular spectropolarimetric
sensing of the surface of Mars using ground-based telescopes did not reveal
any significant signals, confirming a general lack of false positives [38].
Needless to say, chiral molecules will need to have a large enough circular
polarizance of the same sign and a large enough abundance at a planetary
surface in order to be detectable.

On Earth, algae and phototrophic bacteria have the potential to be detected
locally, whereas vegetation has a cover widespread enough to display remotely
detectable circularly polarizing features from afar. Photosynthesis is one of
the most important hallmarks of life on Earth and is the major life process
underlying global primary productivity. Photosynthesis evolved soon after the
emergence of life itself [39, 40, 41]. This early evolution on Earth and the
advantages of being able to utilize the radiation from the host star as an
energy source support expectations that photosynthesis will likely evolve on
other planets as well. In terms of productivity, surface features, and
evolutionary drive, photosynthesis would thus constitute a likely target.

The circular polarization spectra of terrestrial vegetation relates to the
absorbance of its pigments [35]. The circular polarization features with by
far the largest magnitude are found around the chlorophyll a Q absorption band
($\sim$680 nm). Typically, a split signal is observed, with a negative (left-
handed) band at $\sim$670 nm and a positive (right-handed) band at $\sim$690
nm that are relatively independent. This split circular polarization signal is
the result of the superposition of bands of opposite sign originating from
different chiral macrodomains [42, 43, 44, 45]. It has furthermore been
suggested that the local alignments of the chloroplasts might affect the
spatial variation in circular polarization and could thus affect the overall
signal on the leaf and canopy scale [45]. Beside its use as a biosignature for
the remote detection of extraterrestrial life, the circular polarizance of
vegetation may also be informative for vegetation physiology due to its
dependency on the molecular architecture [27].

While high-quality circular polarimetric spectra can be obtained in the
laboratory [27, 46, 45], the next key step is to take dedicated polarimeters
into the field. Under these circumstances the instruments will have to cope
with very dynamic scenes involving large changes and variability in
illumination (i.e. in direction, diffuseness and intensity). Additionally, due
to reflection or scattering a target may be linearly polarized at levels
$>10\%$ and, depending on the modulation approach, mitigation of linear-to-
circular polarization crosstalk is crucial.

In the present study we investigate the circular polarizance of vegetation in
the field using TreePol, a dedicated circular spectropolarimeter based on a
fast ferro-liquid-crystal (FLC) modulation with dual-beam implementation [27].
We report the detection of circularly polarized biosignatures of various
plants.  While previous studies have only measured transmission and reflection
of single leaves, we show to the best of our knowledge for the first time the
results of measurements on whole plants and even canopies measured at large
distances.

##  2 Materials and methods

###  2.1 Polarization

Polarization in general is described in terms of the four parameters of the
Stokes vector S. With the electric field vectors $E_{x}$ in the x direction
($0^{\circ}$) and $E_{y}$ in the y direction ($90^{\circ}$), $z=0$,
$i=\sqrt{-1}$ and with ∗ representing the complex conjugate, the Stokes vector
is given by:

$\textbf{S}=\begin{pmatrix}I\\\ Q\\\ U\\\ V\\\
\end{pmatrix}=\begin{pmatrix}\left\langle
E_{x}E^{*}_{x}+E_{y}E^{*}_{y}\right\rangle\\\ \left\langle
E_{x}E^{*}_{x}-E_{y}E^{*}_{y}\right\rangle\\\ \left\langle
E_{x}E^{*}_{y}-E_{y}E^{*}_{x}\right\rangle\\\ i\left\langle
E_{x}E^{*}_{y}-E_{y}E^{*}_{x}\right\rangle\\\
\end{pmatrix}=\begin{pmatrix}I_{0^{\circ}}+I_{90^{\circ}}\\\
I_{0^{\circ}}-I_{90^{\circ}}\\\ I_{45^{\circ}}-I_{-45^{\circ}}\\\
I_{RHC}-I_{LHC}\\\ \end{pmatrix}$ (1)

The Stokes parameters $I$, $Q$, $U$ and $V$ refer to intensities thus relating
to remotely measurable quantities. The absolute intensity is given by Stokes
$I$. Stokes $Q$ and $U$ denote the differences in intensity after filtering
linear polarization at perpendicular directions, where Q denotes the
difference between horizontal and vertical polarization and U the difference
in linear polarization but with a 45∘ offset. Stokes $V$ denotes the
difference between right-handed and left-handed circularly polarized light.
$I_{0^{\circ}},I_{90^{\circ}},I_{45^{\circ}}$ and $I_{-45^{\circ}}$ are the
intensities in the respective planes perpendicular to the propagation axis
while $I_{LHC}$ and $I_{RHC}$ are the intensities of left- and right-handed
circularly polarized light, respectively. If the absolute intensity $I$ is
known, the polarization state can thus be completely described by the
normalized quantities $Q/I$, $U/I$ and $V/I$.

###  2.2 Spectropolarimetry

Circular polarization measurements, in the lab and in the field, were carried
out using TreePol (see also Figure 1 for a schematic) [27]. TreePol is a
dedicated spectropolarimetric instrument developed by the Astronomical
Instrumentation Group at the Leiden Observatory (Leiden University), the
Netherlands. The instrument was specifically designed to measure the
fractionally induced circular polarization ($V/I$) as a function of wavelength
(400 nm to 900 nm) and is capable of fast measurements with a sensitivity of
$\sim 1*10^{-4}$. Treepol measures the fractional induced circular
polarization of light after interaction of the sample with unpolarized ambient
light. The polarimetric sensitivity is obtained through ferro-liquid-crystal
(FLC) modulation, which is synchronized with a dual spectrometer. An
achromatic Fresnel rhomb (i.e. a $\lambda/4$ retarder) converts the circular
polarization induced by a target into linear polarization which is modulated
by the FLC, a $\lambda/2$ retarder with a $45^{\circ}$ switching fast-axis.
Additionally, TreePol applies spectral multiplexing with the implementation of
a dual-beam approach in which a polarizing beam-splitter feeds the two
spectrographs with orthogonally polarized light.

The combination of temporal polarization modulation (i.e. the fast FLC
combined with a high-speed spectrometer) with spatial modulation (i.e.
simultaneous recording of orthogonal polarization states using two
synchronized spectrographs) ensures that systematic differential effects are
canceled out to the first order. As such, no spurious polarization signals
down to the $10^{-5}$ level are induced [47]. In order to further mitigate
possible linear polarization crosstalk, the original design of TreePol (see
[27]) was upgraded with a fast continuously spinning (5 Hz) half-wave plate in
front of the Fresnel Rhomb. The angle of view of TreePol is approximately
1.14∘, which allows for accurate target selection. TreePol was targeted using
a calibrated targeting scope mounted on top of the instrument.

The measurements were carried out with varying integration times, depending on
the illumination of the target. Measurements outside were performed using
ambient light (i.e. direct sunlight or diffuse light on days with overcast).
All measurements were carried out on and around the campus of the Vrije
Universiteit Amsterdam, the Netherlands, the Hortus Botanicus Vrije
Universiteit Amsterdam, the Netherlands and around the Biological Research
Center of the Hungarian Academy of Sciences, Hungary.

Figure 1: Schematic of TreePol.

##  3 Results

###  3.1 Circular Spectropolarimetric measurements

In previous studies we demonstrated the use of sensitive circular
spectropolarimetry in transmission in the laboratory. In the current study we
build on the knowledge previously acquired. We will first demonstrate the use
on leaves in reflection in the laboratory. We then compare measurements of the
same plant measured in the laboratory and outside while trying to mimic the
same conditions. Subsequently, we demonstrate the effect (or lack thereof) of
ambient light conditions on the circular polarimetric spectra of nearby
canopy. Finally, we will show the potential to discriminate between abiotic
and biotic matter, even for targets several kilometers away.

###  3.2 Laboratory measurements

In a previous study we used transmission spectropolarimetry to follow the
chiroptical signature in leaves that had been cut from their stems and were
decaying over time in the dark or under daylight conditions [27]. These
results showed a strong decrease of $V/I$ in time while this decrease was much
lower for the chlorophyll a concentration. In Figure 2 we show a small set
(same leaves as [27]) of spectra obtained in _reflection_ for leaves stored in
the light. These reflection measurements are spectrally very similar to those
in transmission, albeit 1 order of magnitude smaller (with the bands for
healthy plants showing a maximum amplitude of $+9*10^{-4}$ and $-7*10^{-4}$).

Figure 2: The circular polarimetric spectrum of Hedera helix leaves measured
in the laboratory using a halogen light source. The same set of leaves (n=3)
was measured after 9 days, showing a remarkable reduction of the chiral
macrodomains (see also [27]). Shaded areas denote the standard error.

###  3.3 Laboratory versus in the field measurements

To further investigate possible differences between laboratory measurements
and in-the-field measurements, the circular polarimetric spectrum of an
ornamental house plant was measured outside under cloudless conditions and in
the laboratory using a halogen light source. The target was measured at a
distance of approximately 2 meters (which gives a measurement of a single
leaf) while the illumination angle in the laboratory was kept at an
approximately equal angle as that of the solar irradiation during the
measurements of the same plant outside. Additionally, during the outside
measurements any movement due to wind was minimized by a windscreen. Figure 3
shows that inside and outside similar spectral signatures were found, but with
some distinct differences: the outside measurements had a wavelength-
independent shift to negative values and exhibited a band at 690 nm that was
less pronounced.

Figure 3: The circular polarimetric spectrum of a Ficus elastica measured in
the laboratory using a halogen light source (left) and measured under a
cloudless overhead sky (right), n=3 for the same area. The red circle (picture
not taken at measurement distance) gives the approximate measuring area for
the two sets of measurements. The shaded areas denote the standard error.

###  3.4 Ambient light conditions

The circularly polarizance of vegetation is directly dependent on the
photosynthetic machinery located within the cell’s chloroplasts, but
illumination conditions are expected to have an effect on the magnitude of the
signals. We have therefore investigated the effect of light conditions on a
single tree. Shown in Figure 4 are the circular spectropolarimetric
measurements of the same tree taken under a cloudless sky (left) and under
overcast conditions (right). These measurements were taken at noon and eight
days apart from each other. Interestingly, the results are almost identical,
indicating that at least for this particular measurement the angle of
incidence has no effect. On both days, wind conditions were comparable and
below 2 bft.

Figure 4: Circular polarimetric spectrum of an unidentified tree measured
under a cloudless overhead sky (left) and under overcast conditions (right),
n=3 for different parts of the same tree. The measured radius was
approximately 0.5 meter. Shaded areas denote the standard error.

###  3.5 Biotic versus abiotic matter

We have measured several distant targets from the roof of one of the
laboratory buildings on the campus of the Vrije Universiteit Amsterdam. Two of
these measurements are shown in Figure 5. The lower left panel shows that the
circular spectropolarimetric measurements of a local sports field, consisting
of artificial turf/grass, yield no net polarization. The slight deviation
($<5*10^{-4}$) around the zero point shows leftover fringes resulting from the
FLC (see also the discussion). Measurements on the tree canopies of the nearby
park, ’het Amsterdamse Bos’ shown in the right panels, however, show a clear
non-zero polarization signal ($\sim 3*10^{-3}$).

Figure 5: The circular polarimetric spectrum of artificial turf (left) versus
that of distant trees (right). The red circle (with a difference between the
two photos due to cropping) provides the approximate measuring area. The same
area was measured three times, shaded areas denote the standard error.

###  3.6 Field measurements in general

An ensemble of various (n=30) measurements of very diverse shrubs and trees in
the field and their mean (in red) is shown in Figure 6. Included in these
results are trees native to the Netherlands and Hungary measured from the
rooftops, but also various more exotic shrubs measured at the Hortus Botanicus
Vrije Universiteit. Especially the positive band of these reflection
measurements is generally smaller than measured in transmission in the
laboratory. In transmission, vegetation shows a negative circular polarizance
band (at $\sim$ 660 nm) of $\sim-3*10^{-3}$ and a positive band at $\sim$ 680
nm of $\sim 8*10^{-3}$.

Figure 6: An ensemble of randomly selected outside circular polarimetric
spectra of highly diverse vegetation targets, n=30. Measurements were taken of
different targets under different ambient conditions. The red graph denotes
the mean, the shaded area denotes the standard error.

##  4 Discussion and conclusions

We show here, as far as we know for the first time in the refereed literature,
the systematic investigation of circular spectropolarimetric sensing of
vegetation in the field. Using our current setup, the circular polarizance of
biotic matter can be readily distinguished from that of abiotic matter over a
distance of up to several kilometers. This demonstrates the potential of
circular polarization spectroscopy for detection of signatures of life. While
edges in scalar reflectance spectra and signatures in linear polarization can
readily emerge in the spectrometry of abiotic matter, circular
spectropolarimetry does not suffer from such drawbacks and is thus able to
provide an unambiguous means to detect biotic matter. Consequently, although
the circular polarization signals are small, they are a much more exclusive
biosignature as they are directly linked to the molecular complexity of life.

The TreePol instrument has been designed especially for highly sensitive and
accurate circular spectropolarimetric measurements in the field. One of the
limitations to its polarimetric performance is the temporal variation of the
scene (i.e., variability in illumination, moving targets, etc.). To overcome
these effects, a fast ferroelectric liquid crystal (FLC) based modulation was
implemented. Both the switch angle and the birefringence of FLC’s are,
however, temperature dependent. As such, FLC’s can introduce systematic
effects like polarized spectral fringing. Usually, these effects can be
accounted for by recalibration, but they can not completely be removed (also
due to temperature dynamics). They are therefore, for instance, still visible
in the measurement of the artificial turf (Figure 5) and limit the
polarimetric sensitivity.

Importantly, although linear polarization cross-talk mitigation strategies are
employed, it may not be possible to eliminate these effects completely,
especially in a highly dynamic environment. Moreover, while most environmental
linear polarization features are relatively spectrally broad, those produced
by vegetation (see for example [48, 49, 50, 51, 52]) are essentially inversely
related to the red-edge and can thus show steep features around the
chlorophyll absorbance band. Obtaining information on the linear polarization
properties of the target will be beneficial, even if only for calibration
purposes.

The circular polarization signals from light scattering are smaller than those
obtained in transmission spectropolarimetry. While the negative band is
similar (on average $\sim 3*10^{-3}$), we found that the positive band is on
average one order of magnitude smaller ($\sim 7*10^{-4}$). Importantly, while
in transmission most leaves have a very similar magnitude of the circular
polarizance normalized by the light intensity, we found much more variation in
signals between plants measured under the same conditions in reflection (see
also Figure 6). A possible explanation may involve the variation in leaf
optical properties and the orientation of the leaves.

The spectral characteristics in some of the measurements appear to have a
different spectral shape compared to the measurements taken in transmission
(cf. [27]). One of the most noticeable features is the much lower intensity of
the largest positive $V/I$ band usually observed. This positive band was more
pronounced in the Ficus target, both measured in the laboratory and measured
outside, than in the deciduous forest canopy. In these measurements, however,
the former target had been positioned relatively close to the polarimeter
(approximately 1.5 meter), which allowed the measurements to be conducted on
single leaves (large leaves; approximately 20$\times$12 cm). As such, rather
than observing a canopy with more random and variable orientations, these
measurements monitor only one static orientation. At this point we have not
yet systematically investigated the influence of solar zenith and azimuth
relative to the leaf surface.

We expected large spectropolarimetric differences depending on the light
conditions, because it has been demonstrated that the angle of incidence
influences the scalar reflectance properties [53]. Even if the circular
component would not be influenced directly, the total intensity might be
variable which would lead to a different value of V/I. We expected that this
phenomenon would be most prominent when comparing measurements on sunny days
and on days with overcast. On days with overcast conditions, most of the
incident light is diffuse and scattered by the clouds as opposed to the
conditions on a sunny day. Our results, however, show that the measurements
taken under a cloudless sky are virtually identical (see Figure 4) to those
taken with overcast sky. Possibly, the results might vary with leaf glossiness
and the angle of measurements. Additionally, as linear polarization is heavily
influenced by light conditions and the angle of incidence, the similar results
between the measurements on sunny days and days under overcast conditions
attest to the effectiveness of our crosstalk mitigation strategies.

Additionally, we have shown circular spectropolarimetric signals for decaying
leaves in reflection spectroscopy (Figure 2). Similar measurements were
already reported in a transmission setup [27]. The results presented in the
present study demonstrate how also in reflection healthy leaves can be
distinguished from unhealthy/dying leaves by the change in $V/I$ spectra. This
signal, which is unique to vegetation, rapidly decreases over time and shows a
more than 3-fold decrease in magnitude for the positive band after 9 days. We
argue that this decrease arises because the macrostructures are not actively
maintained due to, for example, drought stress, since the water supply to the
leaves has been cut off [27]. This highlights a possible significance of
circular polarization spectroscopy as a remote-sensing tool for vegetation and
crop monitoring on Earth.

In the context of astrobiology, the circular polarization of biomolecules is a
powerful biosignature [24, 26, 25, 27, 9, 23]. Compared to other surface
biosignatures there are no significant signals produced by abiotic matter, and
therefore no false positives. The results in this study show a vegetation
signal level in terms of circular polarizance on the order of
$10^{-4}-10^{-3}$. While promising in terms of robustness, the magnitude of
the signal is quite low and therefore detection of these signals from
exoplanets will be challenging. This is especially the case when the signal is
further diluted such as can be expected in a planetary disk average, i.e., by
surfaces with a higher reflection than those creating the circular
polarization. To be able to observe a potentially habitable planet and measure
the signals in circular polarization with a high enough signal-to-noise ratio,
an extremely large space-based telescope is required. At the moment, only the
proposed Large UV/Optical/Infrared Surveyor A (LUVOIR-A) [54] meets this
requirement. Additionally, such a telescope will need adaptive optics and
advanced coronagraphy to suppress the light of the host star to provide a high
enough contrast. In view of the scientific return, such an instrument would,
however, certainly be worth its investment.

One of the other factors that should be taken into account is that while the
circular polarization signals are more or less similar for various types of
terrestrial vegetation, it is unknown what these levels are for the dominant
photosynthesic organisms on other planets. Certain terrestrial brown algae,
for instance, display signals varying up to three orders of magnitude in
strength and displayed signals up to $2*10^{-2}$ [46]. The signals can also
result from smaller molecular organisations, leading to excitonic
polarization, such as displayed by various phototropic bacteria [26]. While we
cannot predict the molecular and structural organisation of such biomolecules
beyond Earth, the here observed lack of false positives in our approach is
encouraging and unique. Any significant signal that is observed is therefore
very likely to originate from a highly organized molecular assembly that is
prevalent enough to be detectable in the first place, thus suggesting the
abundance of something organized in terms of homochiral polymers, and hence
probably life.

We have successfully demonstrated the use of circular spectropolarimetry in
the field and our results underline the potential significance of circular
polarization both as a remotely accessible means of detecting the presence of
extraterrestrial life, and as a valuable remotely applicable tool for
vegetation monitoring on Earth. An important next step will be to use these
results in (exo)planetary models with realistic components such as different
surfaces and clouds (e.g. [30, 55, 29]), while future laboratory and field
studies (such as the utilization of canopy models to account for leaf angle
relative to solar zenith angle) should continue to explore the versatility and
potential of this technique.

## Acknowledgments

We dedicate this work to the memory of our dear friend, colleague and the
initiator of this project, Dr. Wilfred Röling, who unexpectedly passed away on
Friday the 25th of September, 2015.

This work was supported by the Planetary and Exoplanetary Science Programme
(PEPSci), grant 648.001.004, of the Netherlands Organisation for Scientific
Research (NWO) and partly by the Economic Development and Innovation Operative
Programme (GINOP), grants GINOP-2.1.7-15-2016-00713 and
GINOP-2.3.2-15-2016-00001 from the Hungarian Ministry for National Economy.

## References

  * [1] A. Cassan, D. Kubas, J.-P. Beaulieu, M. Dominik, K. Horne, J. Greenhill, J. Wambsganss, J. Menzies, A. Williams, U. G. Jørgensen, A. Udalski, D. P. Bennett, M. D. Albrow, V. Batista, S. Brillant, J. A. R. Caldwell, A. Cole, C. Coutures, K. H. Cook, S. Dieters, D. D. Prester, J. Donatowicz, P. Fouqué, K. Hill, N. Kains, S. Kane, J.-B. Marquette, R. Martin, K. R. Pollard, K. C. Sahu, C. Vinter, D. Warren, B. Watson, M. Zub, T. Sumi, M. K. Szymański, M. Kubiak, R. Poleski, I. Soszynski, K. Ulaczyk, G. Pietrzyński, and Ł. Wyrzykowski, “One or more bound planets per milky way star from microlensing observations,” Nature, vol. 481, no. 7380, pp. 167–169, 2012.
  * [2] S. Seager, “The search for habitable planets with biosignature gases framed by a ‘biosignature drake equation’,” International Journal of Astrobiology, pp. 1–9, 2017.
  * [3] J. L. Grenfell, “A review of exoplanetary biosignatures,” Physics Reports, 2017.
  * [4] E. W. Schwieterman, N. Y. Kiang, M. N. Parenteau, C. E. Harman, S. DasSarma, T. M. Fisher, G. N. Arney, H. E. Hartnett, C. T. Reinhard, S. L. Olson, V. S. Meadows, C. S. Cockell, S. I. Walker, J. L. Grenfell, S. Hegde, S. Rugheimer, R. Hu, and T. W. Lyons, “Exoplanet biosignatures: A review of remotely detectable signs of life,” Astrobiology, 2018.
  * [5] Y. Fujii, D. Angerhausen, R. Deitrick, S. Domagal-Goldman, J. L. Grenfell, Y. Hori, E. Palle, N. Siegler, K. Stapelfeldt, and H. Rauer, “Exoplanet biosignatures: Observational prospects,” arXiv preprint arXiv:1705.07098, 2017.
  * [6] E. W. Schwieterman, “Surface and temporal biosignatures,” arXiv preprint arXiv:1803.05065, 2018.
  * [7] S. I. Walker, W. Bains, L. Cronin, S. DasSarma, S. Danielache, S. Domagal-Goldman, B. Kacar, N. Y. Kiang, A. Lenardic, C. T. Reinhard, et al., “Exoplanet biosignatures: Future directions,” arXiv preprint arXiv:1705.08071, 2017.
  * [8] S. Seager, W. Bains, and J. Petkowski, “Toward a list of molecules as potential biosignature gases for the search for life on exoplanets and applications to terrestrial biochemistry,” Astrobiology, vol. 16, no. 6, pp. 465–485, 2016.
  * [9] C. H. L. Patty, I. L. ten Kate, W. B. Sparks, and F. Snik, “Remote sensing of homochirality: a proxy for the detection of extraterrestrial life,” in Chiral Analysis, Second Edition: Advances in Spectroscopy, Chromatography and Emerging Methods (P. Polavarapu, ed.), ch. 2, pp. 29–69, Elsevier Science, 2018\.
  * [10] L. Kaltenegger, W. A. Traub, and K. W. Jucks, “Spectral evolution of an Earth-like planet,” The Astrophysical Journal, vol. 658, no. 1, p. 598, 2007.
  * [11] D. J. Des Marais, M. O. Harwit, K. W. Jucks, J. F. Kasting, D. N. Lin, J. I. Lunine, J. Schneider, S. Seager, W. A. Traub, and N. J. Woolf, “Remote sensing of planetary properties and biosignatures on extrasolar terrestrial planets,” Astrobiology, vol. 2, no. 2, pp. 153–181, 2002.
  * [12] S. D. Domagal-Goldman, A. Segura, M. W. Claire, T. D. Robinson, and V. S. Meadows, “Abiotic ozone and oxygen in atmospheres similar to prebiotic Earth,” The Astrophysical Journal, vol. 792, no. 2, p. 90, 2014.
  * [13] E. W. Schwieterman, V. S. Meadows, S. D. Domagal-Goldman, D. Deming, G. N. Arney, R. Luger, C. E. Harman, A. Misra, and R. Barnes, “Identifying planetary biosignature impostors: Spectral features of ${CO}$ and ${O_{4}}$ resulting from abiotic ${O_{2}/O_{3}}$ production,” The Astrophysical Journal Letters, vol. 819, no. 1, p. L13, 2016.
  * [14] C. E. Harman, E. W. Schwieterman, J. C. Schottelkotte, and J. F. Kasting, “Abiotic O2 levels on planets around F, G, K, and M stars: Possible false positives for life?,” The Astrophysical Journal, vol. 812, no. 2, p. 137, 2015.
  * [15] R. Wordsworth and R. Pierrehumbert, “Abiotic oxygen-dominated atmospheres on terrestrial habitable zone planets,” The Astrophysical Journal Letters, vol. 785, no. 2, p. L20, 2014.
  * [16] S. Seager, E. L. Turner, J. Schafer, and E. B. Ford, “Vegetation’s red edge: a possible spectroscopic biosignature of extraterrestrial plants,” Astrobiology, vol. 5, no. 3, pp. 372–390, 2005.
  * [17] N. Y. Kiang, J. Siefert, and R. E. Blankenship, “Spectral signatures of photosynthesis. I. review of Earth organisms,” Astrobiology, vol. 7, no. 1, pp. 222–251, 2007.
  * [18] E. W. Schwieterman, C. S. Cockell, and V. S. Meadows, “Nonphotosynthetic pigments as potential biosignatures,” Astrobiology, vol. 15, no. 5, pp. 341–361, 2015.
  * [19] R. Popa, Between necessity and probability: searching for the definition and origin of life.  Springer Science & Business Media, 2004.
  * [20] W. A. Bonner, “Chirality and life,” Origins of Life and Evolution of the Biosphere, vol. 25, no. 1-3, pp. 175–190, 1995.
  * [21] F. Jafarpour, T. Biancalani, and N. Goldenfeld, “Noise-induced mechanism for biological homochirality of early life self-replicators,” Physical Review Letters, vol. 115, no. 15, p. 158101, 2015.
  * [22] J. C. Kemp, G. D. Henson, C. T. Steiner, and E. R. Powell, “The optical polarization of the sun measured at a sensitivity of parts in ten million,” Nature, vol. 326, no. 6110, pp. 270–273, 1987.
  * [23] R. D. Wolstencroft, G. E. Tranter, and D. D. Le Pevelen, “Diffuse reflectance circular dichroism for the detection of molecular chirality: An application in remote sensing of flora,” in Bioastronomy 2002: Life Among the Stars, vol. 213, p. 149, 2004.
  * [24] M. Pospergelis, “Spectroscopic measurements of the four stokes parameters for light scattered by natural objects.,” Soviet Astronomy, vol. 12, p. 973, 1969.
  * [25] W. Sparks, J. Hough, L. Kolokolova, T. Germer, F. Chen, S. DasSarma, P. DasSarma, F. Robb, N. Manset, I. Reid, F. Macchetto, and W. Martin, “Circular polarization in scattered light as a possible biomarker,” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 110, no. 14-16, pp. 1771–1779, 2009.
  * [26] W. B. Sparks, J. Hough, T. A. Germer, F. Chen, S. DasSarma, P. DasSarma, F. T. Robb, N. Manset, L. Kolokolova, N. Reid, et al., “Detection of circular polarization in light scattered from photosynthetic microbes,” Proceedings of the National Academy of Sciences, vol. 106, no. 19, pp. 7816–7821, 2009.
  * [27] C. H. L. Patty, L. J. J. Visser, F. Ariese, W. J. Buma, W. B. Sparks, R. J. M. van Spanning, W. F. M. Röling, and F. Snik, “Circular spectropolarimetric sensing of chiral photosystems in decaying leaves,” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 189, pp. 303–311, 2017\.
  * [28] D. M. Stam, W. A. D. Rooij, G. Cornet, and J. W. Hovenier, “Integrating polarized light over a planetary disk applied to starlight reflected by extrasolar planets,” Astronomy & Astrophysics, vol. 452, no. 2, pp. 669–683, 2006.
  * [29] D. M. Stam, “Spectropolarimetric signatures of earth-like extrasolar planets,” Astronomy & Astrophysics, vol. 482, no. 3, pp. 989–1007, 2008\.
  * [30] L. Rossi and D. M. Stam, “Circular polarization signals of cloudy (exo) planets,” arXiv preprint arXiv:1805.08686, 2018.
  * [31] Y. Shkuratov, S. Bondarenko, A. Ovcharenko, C. Pieters, T. Hiroi, H. Volten, O. Muñoz, and G. Videen, “Comparative studies of the reflectance and degree of linear polarization of particulate surfaces and independently scattering particles,” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 100, no. 1, pp. 340–358, 2006.
  * [32] R. A. West, L. R. Doose, A. M. Eibl, M. G. Tomasko, and M. I. Mishchenko, “Laboratory measurements of mineral dust scattering phase function and linear polarization,” Journal of Geophysical Research: Atmospheres, vol. 102, no. D14, pp. 16871–16881, 1997.
  * [33] C. Bustamante, M. F. Maestre, and I. Tinoco Jr, “Circular intensity differential scattering of light by helical structures. i. theory,” The Journal of Chemical Physics, vol. 73, no. 9, pp. 4273–4281, 1980.
  * [34] D. Keller and C. Bustamante, “Theory of the interaction of light with large inhomogeneous molecular aggregates. II. psi-type circular dichroism,” The Journal of Chemical Physics, vol. 84, no. 6, pp. 2972–2980, 1986.
  * [35] G. Garab and H. van Amerongen, “Linear dichroism and circular dichroism in photosynthesis research,” Photosynthesis Research, vol. 101, no. 2-3, pp. 135–146, 2009.
  * [36] U. J. Meierhenrich, J.-J. Filippi, C. Meinert, J. H. Bredehöft, J. ichi Takahashi, L. Nahon, N. C. Jones, and S. V. Hoffmann, “Circular dichroism of amino acids in the vacuum-ultraviolet region,” Angewandte Chemie International Edition, vol. 49, no. 42, pp. 7799–7802, 2010.
  * [37] W. Martin, E. Hesse, J. Hough, and T. Gledhill, “High-sensitivity stokes spectropolarimetry on cyanobacteria,” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 170, pp. 131–141, 2016.
  * [38] W. B. Sparks, J. H. Hough, and L. E. Bergeron, “A search for chiral signatures on Mars,” Astrobiology, vol. 5, no. 6, pp. 737–748, 2005.
  * [39] D. J. Des Marais, “When did photosynthesis emerge on earth?,” Science, vol. 289, no. 5485, pp. 1703–1705, 2000.
  * [40] J. Xiong and C. E. Bauer, “Complex evolution of photosynthesis,” Annual Review of Plant Biology, vol. 53, no. 1, pp. 503–521, 2002.
  * [41] M. F. Hohmann-Marriott and R. E. Blankenship, “Evolution of photosynthesis,” Annual Review of Plant Biology, vol. 62, no. 1, pp. 515–548, 2011.
  * [42] L. Finzi, C. Bustamante, G. Garab, and C.-B. Juang, “Direct observation of large chiral domains in chloroplast thylakoid membranes by differential polarization microscopy,” Proceedings of the National Academy of Sciences, vol. 86, no. 22, pp. 8748–8752, 1989.
  * [43] G. Garab, L. Finzi, and C. Bustamante, “Differential polarization imaging of chloroplasts: Microscopic and macroscopic linear and circular dichroism,” Light in Biology and Medicine, vol. 2, 1991.
  * [44] G. Garab, S. Wells, L. Finzi, and C. Bustamante, “Helically organized macroaggregates of pigment-protein complexes in chloroplasts: evidence from circular intensity differential scattering,” Biochemistry, vol. 27, no. 16, pp. 5839–5843, 1988.
  * [45] C. H. L. Patty, D. A. Luo, F. Snik, F. Ariese, W. J. Buma, I. L. ten Kate, R. J. M. van Spanning, W. B. Sparks, T. A. Germer, G. Garab, and M. W. Kudenov, “Imaging linear and circular polarization features in leaves with complete Mueller matrix polarimetry,” Biochimica et Biophysica Acta (BBA) - General Subjects, vol. 1862, no. 6, 2018.
  * [46] C. H. L. Patty, F. Ariese, W. J. Buma, I. L. ten Kate, R. J. M. van Spanning, and F. Snik, “Circular spectropolarimetric sensing of higher plant and algal chloroplast structural variations,” Photosynthesis Research, pp. 1–11, 2018.
  * [47] F. Snik and C. U. Keller, “Astronomical polarimetry: Polarized views of stars and planets,” in Planets, Stars and Stellar Systems, pp. 175–221, Springer, 2013.
  * [48] J. I. Peltoniemi, M. Gritsevich, and E. Puttonen, Reflectance and polarization characteristics of various vegetation types.  Springer Nature, 2015.
  * [49] V. C. Vanderbilt, L. Grant, L. Biehl, and B. Robinson, “Specular, diffuse, and polarized light scattered by two wheat canopies,” Applied optics, vol. 24, no. 15, pp. 2408–2418, 1985.
  * [50] V. C. Vanderbilt, L. Grant, and C. Daughtry, “Polarization of light scattered by vegetation,” Proceedings of the IEEE, vol. 73, no. 6, pp. 1012–1024, 1985.
  * [51] V. C. Vanderbilt, C. S. T. Daughtry, M. K. Kupinski, C. L. Bradley, A. N. French, K. Bronson, R. A. Chipman, and R. P. Dahlgren, “Estimating the relative water content of leaves in a cotton canopy,” in Polarization Science and Remote Sensing VIII (F. Snik and J. A. Shaw, eds.), vol. 10407, p. 104070Z, International Society for Optics and Photonics, SPIE, 2017.
  * [52] L. Grant, C. Daughtry, and V. Vanderbilt, “Polarized and specular reflectance variation with leaf surface features,” Physiologia Plantarum, vol. 88, no. 1, pp. 1–9, 1993.
  * [53] S. Kaasalainen, O. Nevalainen, T. Hakala, and K. Anttila, “Incidence angle dependency of leaf vegetation indices from hyperspectral lidar measurements,” Photogrammetrie-Fernerkundung-Geoinformation, vol. 2016, no. 2, pp. 75–84, 2016.
  * [54] LUVOIR Team, “The luvoir mission concept study interim report,” arXiv preprint arXiv:1809.09668, 2018.
  * [55] T. Karalidi and D. Stam, “Modeled flux and polarization signals of horizontally inhomogeneous exoplanets applied to earth-like planets,” Astronomy & Astrophysics, vol. 546, p. A56, 2012.


${\mathbf{T}}_{k}$ in (110) accordingly, we obtain on the left hand side
$\boldsymbol{\beta}_{k}{\mathbf{T}}_{k}=\tilde{\boldsymbol{\beta}}_{k}{\mathbf{T}}^{\star}$,
where ${\mathbf{T}}^{\star}={\mathbf{P}}{\mathbf{T}}_{k}$. Split
${\mathbf{T}}^{\star}$ in the following way:

$\displaystyle{\mathbf{T}}^{\star}=\left(\begin{array}[]{cc}{\mathbf{T}}^{\star}_{1}&{\mathbf{T}}^{\star}_{3}\\\
({\mathbf{T}}^{\star}_{3})^{\top}&{\mathbf{T}}^{\star}_{2}\\\
\end{array}\right).$

From (A.27) we obtain:

$\displaystyle\boldsymbol{\beta}^{\star}_{k}{\mathbf{T}}_{k}=\tilde{\boldsymbol{\beta}}^{\star}_{k}{\mathbf{T}}^{\star}=\left(\begin{array}[]{cc}\boldsymbol{\beta}^{\star}_{r}{\mathbf{T}}^{\star}_{1}&\boldsymbol{\beta}^{\star}_{r}{\mathbf{T}}^{\star}_{3}\end{array}\right)=\left(\begin{array}[]{cc}\boldsymbol{\Lambda}^{\star}&{\mathbf{O}}\end{array}\right).$

Since $\boldsymbol{\beta}^{\star}_{r}$ has full column rank, we obtain from
$\boldsymbol{\beta}^{\star}_{r}{\mathbf{T}}^{\star}_{3}={\mathbf{O}}$ by left
multiplication with
$((\boldsymbol{\beta}^{\star}_{r})^{\top}\boldsymbol{\beta}^{\star}_{r})^{-1}(\boldsymbol{\beta}^{\star}_{r})^{\top}$
that ${\mathbf{T}}^{\star}_{3}={\mathbf{O}}$. Furthermore,
$\boldsymbol{\beta}^{\star}_{r}{\mathbf{T}}^{\star}_{1}=\boldsymbol{\Lambda}^{\star}$.
Application of Theorem 1 to the GLT matrix $\boldsymbol{\beta}^{\star}_{r}$,
which satisfies
$\boldsymbol{\beta}^{\star}_{r}(\boldsymbol{\beta}^{\star}_{r})^{\top}=\boldsymbol{\beta}_{k}^{\star}(\boldsymbol{\beta}_{k}^{\star})^{\top}=\boldsymbol{\Lambda}^{\star}(\boldsymbol{\Lambda}^{\star})^{\top}$,
yields $\boldsymbol{\beta}_{r}^{\star}=\boldsymbol{\Lambda}^{\star}$.
Therefore, the pivots $\tilde{l}_{1},\ldots,\tilde{l}_{r}$ are equal to the
pivots ${l}_{1},\ldots,{l}_{r}$ of $\boldsymbol{\Lambda}$. Furthermore, we
obtain ${\mathbf{T}}^{\star}_{1}={{\mathbf{I}}}_{r}$.

Consider now the $s$ rows $\boldsymbol{\beta}_{k}^{d}$ that were deleted from
$\boldsymbol{\beta}_{k}$ and denote the corresponding rows in the reordered
matrix $\tilde{\boldsymbol{\beta}}_{k}$ as
$\tilde{\boldsymbol{\beta}}_{k}^{d}$. The first $r$ columns of
$\tilde{\boldsymbol{\beta}}_{k}^{d}$ are the rows $\boldsymbol{\beta}_{r}^{d}$
deleted from $\boldsymbol{\beta}_{r}$ to define
$\boldsymbol{\beta}^{\star}_{r}$ and the last $s$ columns are a lower
triangular matrix ${\mathbf{L}}$ with pivot elements $n_{1},\ldots,n_{s}$ on
the main diagonal. From
$\boldsymbol{\beta}_{k}^{d}{\mathbf{T}}_{k}=\tilde{\boldsymbol{\beta}}_{k}^{d}{\mathbf{T}}^{\star}$,
we obtain

$\displaystyle\left(\begin{array}[]{cc}\boldsymbol{\beta}_{r}^{d}&{\mathbf{L}}{\mathbf{T}}^{\star}_{2}\end{array}\right)=\left(\begin{array}[]{cc}\boldsymbol{\Lambda}^{d}&{\mathbf{M}}_{0}\end{array}\right).$

Hence, the deleted rows $\boldsymbol{\beta}_{r}^{d}$ are identical to the rows
$\boldsymbol{\Lambda}^{d}$ deleted from $\boldsymbol{\Lambda}$ and the first
$r$ columns of $\tilde{\boldsymbol{\beta}}_{k}$ are equal to
$\boldsymbol{\Lambda}$. Since
${\mathbf{L}}{\mathbf{T}}^{\star}_{2}({\mathbf{T}}^{\star}_{2})^{\top}{\mathbf{L}}^{\top}={\mathbf{L}}{\mathbf{L}}^{\top}={\mathbf{M}}_{0}{\mathbf{M}}_{0}^{\top}={\mathbf{D}}_{0}^{2}$,
where ${\mathbf{D}}_{0}$ is a diagonal matrix, it follows that ${\mathbf{L}}$
is a diagonal matrix and equal to ${\mathbf{D}}_{0}$. Hence, in the GLT
framework, the last $s$ columns of $\tilde{\boldsymbol{\beta}}_{k}$,
corresponding the spurious factor matrix, can be represented as a spurious
ordered GLT matrix ${\mathbf{M}}^{\Lambda}_{s}$, where the pivot rows are
equal to ${\mathbf{L}}$. This completes the proof:

$\displaystyle\left(\begin{array}[]{cc}\boldsymbol{\Lambda}&{\mathbf{M}}^{\Lambda}_{s}\end{array}\right)=\tilde{\boldsymbol{\beta}}_{k}=\boldsymbol{\beta}_{k}{{\mathbf{P}}}=\boldsymbol{\beta}_{k}{\mathbf{P}}_{\pm}{\mathbf{P}}_{\rho}.$


# New Taub-NUT Black Holes with Massive Spin-2 Hair

Yu-Qi Chen<EMAIL_ADDRESS>Hai-Shan Liu<EMAIL_ADDRESS>Center for
Joint Quantum Studies and Department of Physics,
School of Science, Tianjin University, Tianjin 300350, China

###### Abstract

We consider Einstein gravity extended with quadratic curvature invariants,
where the well-known Ricci-flat Taub-NUT black hole remains a solution. An
analysis of the unstable Lichnerowicz modes in the Taub-NUT background enables
us to identify the mass and NUT parameters $(m,n)$ where new Taub-NUT black
holes can emerge. We then adopt numerical technique to construct these new
Taub-NUT black holes that bifurcate away from the Ricci-flat ones. Unlike the
Ricci-flat Taub-NUT, there can exist two new black holes for an appropriate
given temperature, making it a total of three if we include the Ricci-flat
one.

Introduction.\- In the beginning of last century, Einstein discovered General
Relativity (GR), which has two important predictions. One is the black hole
and the other is gravitational wave. In this century, people have directly
detected gravitational waves LIGOScientific:2017vwq ; LIGOScientific:2017zic
and taken images of black holesEventHorizonTelescope:2019dse ;
EventHorizonTelescope:2019uob ; EventHorizonTelescope:2019jan ;
EventHorizonTelescope:2019ths ; EventHorizonTelescope:2019pgp ;
EventHorizonTelescope:2019ggy . These great successes validate GR as a
classical theory. However, GR is not a perfect theory; it is not
renormalizable and can have have an infinite possibility of higher order
corrections tHooft:1974toh . It was pointed out that one can obtain a
renormalizable theory by including also the quadratic curvature terms to the
Einstein-Hilbert action. The price to pay is that these higher-derivative
terms that make the theory renormalizable will inevitably introduce ghost
modes Stelle:1976gc . Many works have been done on addressing the ghost
phenomena of higher derivative gravity Horava:2009uw ; Lu:2009em ; Li:2008dq ;
Lu:2010ct ; Lu:2010cg ; Lu:2011zk  and it was pointed out that the ghost
problem might not be severe to the theory Smilga:2013vba .

There are three quadratic curvature terms in general, Ricci scalar square,
Ricci tensor square and Riemann tensor
square($R^{2},R_{\mu\nu}R^{\mu\nu},R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$).
We can regroup them as the Gauss-Bonnet combination
$R^{2}+4R_{\mu\nu}R^{\mu\nu}+R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$, Weyl-
squared term $C^{\mu\nu\rho\sigma}C_{\mu\nu\rho\sigma}$ and $R^{2}$. In four
dimensions, the Gauss-Bonnet term is a total derivative and hence doesn’t
contribute to the equations of motion. The remaining Weyl-squared term and
$R^{2}$ will introduce massive spin-2 and massive scalar modes respectively.
The absence of the massive scalar mode can lead to critical gravityLu:2011zk .
One important property of the quadratic curvature extension in four dimensions
is that Ricci-flat metrics continue to be the solutions of the extended
theory. This raises a question: does the theory admit new black holes beyond
the Schwarzschild black hole? It turns out that the Ricci scalar must be zero
owing to the application of the no scalar theorem in higher-derivative gravity
Lu:2015cqa . On the other hand, black holes beyond the Schwarzschild black
hole can carry massive spin-2 hair and they were constructed in Lu:2015cqa .

Taub-NUT spacetime is another one of the simplest Ricci-flat solutions in
Einstein gravity which was constructed in 1960s Taub:1950ez ; Newman:1963yy .
Compared to the famous Schwarszchild black hole, Taub-NUT spacetime has an
additional integration constant $n$, which is called NUT parameter. Due to
this NUT parameter, it has a string-like singularity along the polar axis,
which is now known as Misner string singularityMisner:1963fr . Surprisingly,
the Taub-NUT metric has no curvature singularity though it possesses an event
horizon and two other Killing horizons. Those peculiar properties of Taub-NUT
spacetime have been attracting many attentions for years Clement:2015cxa ;
Clement:2015aka ; Hennigar:2019ive ; Bordo:2019slw ; Durka:2019ajz ;
Awad:2022jgn ; Wu:2019pzr ; Rodriguez:2021hks ; Liu:2022wku ; Liu:2023uqf ;
Jiang:2019yzs ; Chen:2023eio ; Siahaan:2022jrl ; Yang:2023hll .

It is clear the Ricci-flat Taub-NUT black hole will continue to be the
solution in the extended gravity theory. In this paper, we would like to
investigate whether new black holes carrying the NUT parameter beyond the
Ricci-flat solution can arise. As in the case of Lu:2015cqa , the no-hair
theory will rule out the possibility of having the hair associated with the
massive scalar mode. We focus on the construction involving the massive spin-2
mode. The general solution should contain two parameters $(m,n)$. In the
Ricci-flat case, any $(m,n\neq 0)$ can lead to an event horizon, but this is
no longer true when higher-derivative terms are involved. The parameters
$(m,n)$ has to be fine tuned in order to form an event horizon. This fine-
tuning leads to difficulties in constructing the numerical solutions, since we
have to start at some fiducial point of the parameter space.

Since the massive spin-2 mode of the quadratic extension is relevant for the
new solution, we perform a spin-2 tensor mode analysis upon the Ricci-flat
Taub-NUT background, which is called Lichnerowicz mode analysis. We find that
there exists a Licnerowicz mode for some specific values of $(m,n)$ with a
negative eigenvalue, which implies that the Taub-NUT background of this
$(m,n)$ is not stable and a more stable and new black hole solution should
emerge at this point of the parameter space. With this guidance we construct
new Taub-NUT like black hole solutions carrying spin-2 hair in the theory of
Einstein gravity added with quadratic curvature terms.

Theroy.\- We consider Einstein gravity extended with the most general
quadratic terms in four dimensional spacetime

$I=\int d^{4}x\sqrt{-g}(\kappa R-\alpha
C_{\mu\nu\rho\sigma}C^{\mu\nu\rho\sigma}+\beta R^{2})\,,$ (1)

where $C_{\mu\nu\rho\sigma}$ is Weyl tensor and $\kappa,~{}\alpha$ and $\beta$
are constants. For simplicity we shall set $\kappa=1$. The equations of motion
can be obtained through the variation of the metric and written in a simple
form

$\begin{split}&R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R-4\alpha
B_{\mu\nu}+2\beta(R_{\mu\nu}\\\
&-\frac{1}{2}g_{\mu\nu}R+g_{\mu\nu}\nabla_{\rho}\nabla^{\rho}R-\nabla_{\mu}\nabla_{\nu}R)=0\end{split}$
(2)

where
$B_{\mu\nu}=(\nabla^{\rho}\nabla^{\sigma}+\frac{1}{2}R^{\rho\sigma})C_{\mu\nu\rho\sigma}$
is the trace free Bach tensor.

It is worth pointing out that we don’t use the Gauss-Bonnet terms to get rid
off the Riemann square but to group the quadratic curvature terms into Weyl
square. The advantage is that the Einstein gravity theory with quadratic
curvature terms contains three modes in the linear spectrum, while the three
terms in the theory (1) correspond to the three modes respectively, a usual
massless spin-2 graviton, a massive spin-2 mode with mass
$m_{2}=1/\sqrt{2\alpha}$ and a massive scalar mode with mass
$m_{0}=1/\sqrt{6\beta}$. It was proved that the scalar hair is forbidden for
static black hole solutionLu:2015cqa , since the trace part of the equation of
motion leads to the vanishing of the Ricci scalar, $R=0$. However, the massive
spin-2 mode can survive, and the theory can admit black hole solutions
carrying spin-2 hair . In this paper, we want to explore the situation in
Taub-NUT spacetimes.

The Taub-NUT solution was constructed in pure Einstein gravity,

$ds^{2}_{4}=-h(dt+2n\cos\theta
d\phi)^{2}+\frac{dr^{2}}{f}+(r^{2}+n^{2})d\Omega^{2}\,,$ (3)

with

$h=f=\frac{r^{2}-2mr-n^{2}}{r^{2}+n^{2}}$ (4)

The solution is a generalization of the famous Schwarzschild black hole
solution. In add to the the mass parameter $m$, the Taub-NUT has one more
integration constant $n$ which is usually called NUT parameter. Setting $n=0$,
the Taub-NUT solution reduces to Schwarzschild solution.

As is mensioned, the Taub-NUT solution is naturally a solution of theory (1),
since the solution is Ricci flat, $R_{\mu\nu}=0$. It is easily to verify that
the Ricci scalar also vanishes under the Taub-NUT like spacetime (3) with
arbitrary $h(r)$ and $f(r)$ satisfying condition of $h$ and $f$ vanish on the
event horizon. The trace of the equations of motion (2) is $\beta(\Box-
m_{0}^{2})R=0$. Multiplying $R$ and integrating over the horizon to infinity
gives $\int dx^{4}\sqrt{g}R(\Box-m_{0}^{2})R=\int
dx^{4}\Big{[}\partial_{r}(\sqrt{g}g^{rr}R\partial_{r}R)-\sqrt{g}(g^{rr}\partial_{r}R\partial_{r}R+m_{0}^{2}R^{2})\Big{]}=0$.
Since $g^{rr}$ vanishes on the horizon and $\partial_{r}R$ goes to 0 at
spatial infinity, the total derivative term gives no contribution, and the
remaining two terms are non positive, thus the Ricci scalar has to be 0. Any
Taub-NUT like solution of metric ansatz (3), with undetermined $h$ and $f$, in
the theory (1) will have a vanishing Ricci scalar. This property of vanishing
Ricci scalar makes a great simplification to the equations of motion.

Lichnerowicz Mode.\- Since Taub-NUT black hole is a solution of the theory
(1), we want to analyze the linear exited modes under the Taub-NUT background.
The Taub-NUT spacetime is locally flat, $R_{\mu\nu}=0$, making a perturbation
of Ricci tensor, we find

$\begin{split}\delta R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\delta
R+(2\beta-\frac{2}{3}\alpha)(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu})\delta
R\\\ -2\alpha\Box(\delta R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\delta R)-4\alpha
R_{\mu\rho\nu\sigma\delta}\delta R^{\rho\sigma}=0.\end{split}$ (5)

The trace of above equation is$(6\beta\Box-1)\delta R=0$. With the same
analysis in the last section, we can deduce that the $\delta R$ must vanish.
It means that any Taub-NUT like black hole which is perturbatively around
Taub-NUT spacetime must have a vanishing Ricci scalar at the linearised level.
Substituting $\delta R=0$ back to equation (5) gives

$(\Delta_{L}+\frac{1}{2\alpha})\delta R_{\mu\nu}=0,$ (6)

where

$\Delta_{L}\delta R_{\mu\nu}=-\Box\delta
R_{\mu\nu}-2R_{\mu\nu\rho\sigma}\delta R^{\rho\sigma}$ (7)

is called Lichnerowicz operator. For Ricci flat spacetime($R_{\mu\nu}=0$) such
as Taub-NUT black hole, firstly, we can deduce that $g^{\mu\nu}\delta
R_{\mu\nu}=0$ from $\delta R=0$. Secondly, taking variation of the Bianchi
identity $\nabla^{\mu}R_{\mu\nu}-(1/2)\nabla_{\nu}R=0$, we can have
$\nabla^{\mu}\delta R_{\mu\nu}=0$. Thus, $\delta R_{\mu\nu}$ is a transverse
traceless tensor on the Ricci flat background.

This implies that once the Lichnerowicz operator has a transverse traceless
eigenfunction $\psi_{\mu\nu}$ with eigenvalue $\lambda=-1/(2\alpha)$,
$\Delta_{L}\psi_{\mu\nu}=\lambda\psi_{\mu\nu}$, there exists a linearized
perturbation away from the Taub-NUT solution.

In the large $r$ limit, the Lichnerowicz equations can be solved and two
Yukawa like modes, $e^{\pm\sqrt{-\lambda r}}$, emerge, which is under
expectation. In general, the Lichnerowicz equation can not be solved
analytically, and to solve the equation numerically, one need to use shooting
method to seek proper values of the parameters to elude the diverging mode
$e^{+\sqrt{\lambda r}}$, left only the normalisable mode $e^{-\sqrt{\lambda
r}}$. Without loss of generality, we set $\alpha=1/2$, and find that there
exists one mode corresponding to each value of $n$, and relation between the
starting point $r_{0}$ of the linear perburbation and $n$ are shown in Fig 1.
The technical details can be found in the Appendix. We can see that when
$n=0$, $r_{0}=0.876$, which covers the result for Schwarszchild case
Lu:2015cqa .

Figure 1:  The bifurcation point ($r_{0}$, n) between linear perturbation and
the analytic Taub-NUT solution with the fixed parameter $\alpha=1/2$.

New Taub-NUT Black Hole Solutions.\- The Lichnerowicz mode analysis shows that
a spin-2 perturbation exists under the analytic Taub-NUT background, which
implies that there must exist black hole solutions carrying spin-2 hair. At
this stage, we carry out a numerical search of Taub-NUT like black hole
solutions which carry spin-$2$ hair in the theory (1). We start with Taub-NUT
like metric ansatz (3) and suppose there exists a black hole with event
horizon located at radius $r_{0}$, where metric functions $h$ and $f$ both
vanish. Then, we take a Taylor expansion around the horizon for $h$ and $f$,

$\begin{split}f(r)=f_{1}(r-r_{0})+f_{2}(r-r_{0})^{2}+f_{3}(r-r_{0})^{3}+...\\\
h(r)=h_{1}(r-r_{0})+h_{2}(r-r_{0})^{2}+h_{3}(r-r_{0})^{3}+...\end{split}$ (8)

Substituting these back into the equations of motion, the coefficients $h_{i}$
and $f_{i}$ for $i>1$ can be solved and expressed in terms of $f_{1}\,,h_{1}$
and $r_{0}$, the first two coefficients are

$\begin{split}f_{2}=\frac{-3n^{2}[1+4\alpha
f_{1}(f_{1}-h_{1})-f_{1}r_{0}]}{8\alpha f_{1}r_{0}(r_{0}^{2}+n^{2})}\\\
+\frac{8\alpha f_{1}(1-2f_{1}r_{0})+3r_{0}(f_{1}r_{0}-1)}{8\alpha
f_{1}(r_{0}^{2}+n^{2})}\\\ h_{2}=\frac{h_{1}n^{2}[1+4\alpha
f_{1}(f_{1}-h_{1})-f_{1}r_{0}]}{8\alpha f_{1}r_{0}(r_{0}^{2}+n^{2})}\\\
+\frac{8\alpha f_{1}(-1+2f_{1}r_{0})+r_{0}(f_{1}r_{0}-1)}{8\alpha
f_{1}(r_{0}^{2}+n^{2})}\end{split}$ (9)

We use these near-horizon expansion to set initial data and then integrate the
solution from horizon to large $r$. Since Taub-NUT solution corresponds to
$f_{1}=1/r_{0}=h_{1}$, we set
$f_{1}=\frac{1}{r_{0}}+\delta,~{}h_{1}=\frac{1}{r_{0}}$. And $\delta$
characterizes the deviation of our newly constructed numerical solution from
the analytical Taub-NUT solution.

Figure 2: The metric functions $f(r)$ and $h(r)$ for the asymptotically-flat
spacetimes with $n=0.5$. The left and right plots show the balk holes with
$r_{0}=1$ and $r_{0}=1.5$ respectively.

We set the theory parameter $\alpha=1/2$ as we did in the last section, and
$n=1/2$, then start the investigation of new solution around $r_{0}=0.852$ as
suggested by the result of the Lichnerowicz mode analysis. We find a newly
Taub-NUT like black holes carrying spin-2 hair under expectation. We show the
metric functions $h$ and $f$ in the Fig.2, it is obvious that the two metric
functions $h$ and $f$ are different, unlike that of the analytical Taub-NUT
solution in which $h$ and $f$ are the same. The profile of metric functions
are qualitatively different for different horizon radius. The Fig.2 shows the
plot of metric function with $r_{0}=1$ and $r_{0}=1.5$, corresponding to black
hole with a positive mass parameter and black hole with a negative mass
parameter, respectively. We read of the mass parameter $m$ from the fall off
of metric function $f\sim 1-2m/r+\dots$ in the large $r$. As is for analytic
Taub-NUT black hole solution, there is no objection for the mass parameter
taking negative value, and the mass of the black hole can have contributions
from both mass parameter $m$ and the NUT parameter $n$ Liu:2022wku ;
Liu:2023uqf . We also show the pattern of the profile of metric function with
negative or positive mass parameter for various NUT parameter $n$, whilst keep
$r_{0}$ fixed in Fig.3.

Figure 3: The metric function $f(r)$ with $r_{0}=1$ and variable NUT
parameters $n$. The function $f(r)$ with $n=0.1,0.5,0.6,1,2$ are shown as
orange, red, blue, pink, and green lines.

In fact, this pattern also emerges in the analytical Taub-NUT solution, where
one can easily derive that the condition for critical value of mass parameter
$m=0$ is $r_{0}=n$. Here, for our numerical solution, we can not give an
analytical relation of ($r_{0}\,,n$) for $m=0$, but we can plot the curve of
the relation ($r_{0}\,,n$) in the $r_{0}-n$ space, which is displayed in
Fig.4. We also plot the cure $r_{0}=n$ of analytic Taub-NUT solution as a
comparison, the critical curve of the newly found Taub-NUT like solution
crosses with the critical curve of analytical Taub-NUT at $n=0.83$. It is
worth pointing out that when $n=0$, the critical value of $m=0$ emerges at
$r_{0}=1.143$, which is consistent with the result in Lu:2015cqa .

Figure 4: The $r_{0}-n$ curve of $m=0$ in Taub-NUT (violet line) and the new
Taub (orange line). The three regions with different $m$ are shown as
$Region~{}A,B,C$. The intersection of Taub-NUT and new Taub is tagged at
$n\approx 0.83$.

We also plot the mass parameter as a function of horizon radius in Fig.5. The
red line is the new found Taub-NUT like black holes, while the blue line
represents the analytic Taub-NUT black hole. One can see that the mass
parameter decreases as the event horizon $r_{0}$ increases. For $n=0.5$, the
mass parameter can be negative when the black hole is larger enough. However,
for larger n (like n=1.5), the mass parameter is always negative for the newly
found Taub-NUT like black holes. It is consistent with Fig.2, mass parameter
decreases as NUT parameter grows. Compared with the Schwarzschild case in
Lu:2015cqa , our new Taub-NUT solutions have two different properties. The
first is that the new solutions have two branches and the smooth connected.
The other is that, for fixed NUT parameter n, there is a minimal value of
event horizon radius,$r_{min}=0.721$ for $n=0.5$ and $r_{min}=0.67$ for
$n=1.5$.

Figure 5: The $m-r_{0}$ curve of analytic Taub-NUT (blue line) and the new
Taub-NUT black holes (red line) for different NUT parameter (n=0.5 for left
plot and n=1.5 for right plot).

We now turn to the thermodynamical properties of the newly found Taub-NUT
black hole solution. The Hawking temperature can be calculated through the
standard method

$T=\frac{\sqrt{f^{\prime}(r_{0})h^{\prime}(r_{0})}}{4\pi}=\frac{\sqrt{1+r_{0}\delta}}{4\pi
r_{0}}\,.$ (10)

Figure 6: The Hawking temperature and Wald entropy in Taub-NUT and the new
Taub-NUT solution with fixed $n=0.5$ are shown in the upleft and upright
plots. And the homologous plots with $n=1.5$ are shown in the bottom. The
Taub-NUT black holes is shown as blue full line, the first new Taub black
holes are shown as red dash lines .

The temperature of the solution as a function of event horizon radius is shown
in Fig.6 for $n=0.5$ and $n=1.5$. As a comparison, we also plot the
temperature cure of Taub-NUT solution, the two curves cross at $r_{0}=0.852$
for $n=0.5$ and at $r_{0}=0.716$ for $n=1.5$. It is consistent with the
Lichnerowicz mode analysis in the previous section. There is a maximum value
of temperature, which is more obvious for $n=1.5$. There can exist two black
holes in the same temperature. Furthermore, with $r_{0}$ decreasing, a notable
inflection point will occur at $r_{0}=0.721$ in the $n=0.5$ case. As is seen
in Fig.6, the two branches of the new solution are smoothly connected in the
minimal radius. And the temperature of one branch is always higher than the
other one.

The entropy of the newly found black hole is not equal to one quarter of the
event horizon area due to the higher derivative terms. According to Wald’s
entropy formula Wald:1993nt ; Iyer:1994ys , the expression of entropy with
higher curvature correction is

$S=-\frac{1}{8}\int\sqrt{h}d\Omega\epsilon_{ab}\epsilon_{cd}\frac{\partial
L}{\partial R_{abcd}}=\pi(r_{0}^{2}+n^{2})-4\pi\alpha r_{0}\delta.$ (11)

As is mentioned in the beginning of this paper, the Gauss-Bonnet terms is a
total derivative term in four dimensional spacetime and do not contribute to
the equations of motion. However, the Gauss-Bonnet terms can have a constant
contribution to the entropy. We use this freedom to make the entropy is equal
to one quarter of the event horizon area for analytic Taub-NUT black hole
solution. We plot the $S-r_{0}$ curve for $n=0.5$ and $n=1.5$ in Fig.6, the
solid line is the cure for analytic Taub-NUT solution. Here, the two branches
of the new solution meet at the minimal radius, and the entropy of one branch
is always larger than the other branch. However, it is worth pointing out that
the branch with lower temperature has larger entropy, and vice versa.

Conculsion.\- We constructed new Taub-NUT-like black hole solutions in
Einstein gravity plus quadratic curvature terms. A Lichnerowicz mode analysis
was carried out on the Ricci-flat Taub-NUT background, which gives us values
of $(m,n)$ where there could exist new black holes carrying the massive spin-2
hair. Then new Taub-NUT like black hole solutions which carry spin-2 hair were
constructed numerically, enlarging the black hole family in higher derivative
gravity theories. The newly found solutions bifurcate from the Ricci-flat
Taub-NUT solutions. Compared with the analytic Taub-NUT solutions, the
temperature of one branch new solution has a maximum for a large $n$, and
there exists two black holes with the same temperature. Since the Ricci-flat
Taub-NUT is a solution, there can be a total of three black holes carrying NUT
parameter for an appropriate given temperature. This implies that there could
be phase transitions between these black holes and the phase diagram becomes
much richer.

As pointed out in Liu:2022wku  the mass of the analytic Taub-NUT black hole
has extra contribution from the Minser strings which makes the total mass of
the spacetime non-negative, though the mass parameter of the Taub-NUT
solutions can be negative. The situation here is similar, the mass parameter
can be negative, but negativity of the total mass could be avoided due to the
extra contribution from the polar strings. However, the spin-2 modes are ghost
like, the solution may still suffer from negative mass problem. The definition
of mass is a difficult problem, but is very crucial at the same time, and so
is the NUT charge. Both the definitions of mass and NUT charge are worth
further exploring.

Acknowledgement.\- We thanks Ze Li and Hong Lu for useful discussion. We are
also grateful for Ling-Yu Zhang’s help in figures. This work is supported in
part by NSFC (National Natural Science Foundation of China) Grants No.
12075166.

## References

  * (1) B. P. Abbott et al. [LIGO Scientific and Virgo], GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, no.16, 161101 (2017),arXiv:1710.05832 [gr-qc]
  * (2) B. P. Abbott et al. [LIGO Scientific, Virgo, Fermi-GBM and INTEGRAL], Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A, Astrophys. J. Lett. 848, no.2, L13 (2017), arXiv:1710.05834 [astro-ph.HE].
  * (3) K. Akiyama et al. [Event Horizon Telescope], First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole, Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA].
  * (4) K. Akiyama et al. [Event Horizon Telescope], First M87 Event Horizon Telescope Results. II. Array and Instrumentation, Astrophys. J. Lett. 875, no.1, L2 (2019), [arXiv:1906.11239 [astro-ph.IM].
  * (5) K. Akiyama et al. [Event Horizon Telescope], First M87 Event Horizon Telescope Results. III. Data Processing and Calibration, Astrophys. J. Lett. 875, no.1, L3 (2019), arXiv:1906.11240 [astro-ph.GA].
  * (6) K. Akiyama et al. [Event Horizon Telescope], First M87 Event Horizon Telescope Results. IV. Imaging the Central Supermassive Black Hole, Astrophys. J. Lett. 875, no.1, L4 (2019), arXiv:1906.11241 [astro-ph.GA].
  * (7) K. Akiyama et al. [Event Horizon Telescope], First M87 Event Horizon Telescope Results. V. Physical Origin of the Asymmetric Ring, Astrophys. J. Lett. 875, no.1, L5 (2019), arXiv:1906.11242 [astro-ph.GA].
  * (8) K. Akiyama et al. [Event Horizon Telescope], First M87 Event Horizon Telescope Results. VI. The Shadow and Mass of the Central Black Hole, Astrophys. J. Lett. 875, no.1, L6 (2019) ,arXiv:1906.11243 [astro-ph.GA].
  * (9) G. ’t Hooft and M. J. G. Veltman, One loop divergencies in the theory of gravitation, Ann. Inst. H. Poincare Phys. Theor. A 20, 69-94 (1974)
  * (10) K. S. Stelle, Renormalization of Higher Derivative Quantum Gravity, Phys. Rev. D 16, 953-969 (1977)
  * (11) P. Horava, Quantum Gravity at a Lifshitz Point, Phys. Rev. D 79, 084008 (2009), arXiv:0901.3775 [hep-th].
  * (12) H. Lu, J. Mei and C. N. Pope, Solutions to Horava Gravity, Phys. Rev. Lett. 103, 091301 (2009) , arXiv:0904.1595 [hep-th].
  * (13) W. Li, W. Song and A. Strominger, Chiral Gravity in Three Dimensions JHEP 04, 082 (2008) , arXiv:0801.4566 [hep-th].
  * (14) H. Lu, C. N. Pope and E. Sezgin, Massive Three-Dimensional Supergravity From $R+R^{2}$ Action in Six Dimensions, JHEP 10, 016 (2010) , arXiv:1007.0173 [hep-th].
  * (15) H. Lu and Y. Pang, On Hybrid (Topologically) Massive Supergravity in Three Dimensions, JHEP 03, 050 (2011) , arXiv:1011.6212 [hep-th].
  * (16) H. Lu and C. N. Pope, Critical Gravity in Four Dimensions, Phys. Rev. Lett. 106, 181302 (2011) , arXiv:1101.1971 [hep-th].
  * (17) A. V. Smilga, Supersymmetric field theory with benign ghosts, J. Phys. A 47, no.5, 052001 (2014) , arXiv:1306.6066 [hep-th].
  * (18) H. Lu, A. Perkins, C. N. Pope and K. S. Stelle, Black Holes in Higher-Derivative Gravity, Phys. Rev. Lett. 114, no.17, 171601 (2015), arXiv:1502.01028 [hep-th].
  * (19) A. H. Taub, Empty space-times admitting a three parameter group of motions, Annals Math. 53, 472-490 (1951)
  * (20) E. Newman, L. Tamburino and T. Unti, Empty space generalization of the Schwarzschild metric, J. Math. Phys. 4, 915 (1963)
  * (21) C. W. Misner, The Flatter regions of Newman, Unti and Tamburino’s generalized Schwarzschild space, J. Math. Phys. 4, 924-938 (1963)
  * (22) G. Clément, D. Gal’tsov and M. Guenouche, Rehabilitating space-times with NUTs, Phys. Lett. B 750, 591-594 (2015) , arXiv:1508.07622 [hep-th].
  * (23) G. Clément, D. Gal’tsov and M. Guenouche, NUT wormholes, Phys. Rev. D 93, no.2, 024048 (2016) , arXiv:1509.07854 [hep-th].
  * (24) R. A. Hennigar, D. Kubizňák and R. B. Mann, Thermodynamics of Lorentzian Taub-NUT spacetimes, Phys. Rev. D 100, no.6, 064055 (2019) , arXiv:1903.08668 [hep-th].
  * (25) A. B. Bordo, F. Gray and D. Kubizňák, Thermodynamics and Phase Transitions of NUTty Dyons, JHEP 07, 119 (2019) , arXiv:1904.00030 [hep-th].
  * (26) R. Durka, The first law of black hole thermodynamics for Taub–NUT spacetime, Int. J. Mod. Phys. D 31, no.04, 2250021 (2022) ,arXiv:1908.04238 [gr-qc].
  * (27) A. Awad and S. Eissa, Lorentzian Taub-NUT spacetimes: Misner string charges and the first law, Phys. Rev. D 105, no.12, 124034 (2022) , arXiv:2206.09124 [hep-th].
  * (28) S. Q. Wu and D. Wu, Thermodynamical hairs of the four-dimensional Taub-Newman-Unti-Tamburino spacetimes, Phys. Rev. D 100, no.10, 101501 (2019) , arXiv:1909.07776 [hep-th].
  * (29) N. H. Rodríguez and M. J. Rodriguez, First law for Kerr Taub-NUT AdS black holes, JHEP 10, 044 (2022) , arXiv:2112.00780 [hep-th].
  * (30) H. S. Liu, H. Lu and L. Ma, Thermodynamics of Taub-NUT and Plebanski solutions, JHEP 10, 174 (2022) , arXiv:2208.05494 [gr-qc].
  * (31) J. F. Liu and H. S. Liu, Thermodynamics of Taub-NUT-AdS Spactimes, arXiv:2309.01609 [hep-th].
  * (32) J. Jiang, B. Deng and X. W. Li, Holographic complexity of charged Taub-NUT-AdS black holes, Phys. Rev. D 100, no.6, 066007 (2019), arXiv:1908.06565 [hep-th].
  * (33) S. Chen, Y. Pei, L. Li and T. Yang, Charged Taub-NUT-AdS Black Holes in $f(R)$ Gravity and Holographic Complexity, Int. J. Theor. Phys. 62, no.2, 16 (2023)
  * (34) H. M. Siahaan, Hidden conformal symmetry and pair production near the cosmological horizon in Kerr-Newman-Taub-NUT-de Sitter spacetime, Chin. Phys. C 47, no.3, 035104 (2023), arXiv:2212.03058 [gr-qc].
  * (35) S. J. Yang, W. D. Guo, S. W. Wei and Y. X. Liu, First law of black hole thermodynamics and the weak cosmic censorship conjecture for Kerr-Newman Taub-NUT black holes, arXiv:2306.05266 [gr-qc].
  * (36) H. Lü, A. Perkins, C. N. Pope and K. S. Stelle, Lichnerowicz Modes and Black Hole Families in Ricci Quadratic Gravity, Phys. Rev. D 96, no.4, 046006 (2017), arXiv:1704.05493 [hep-th].
  * (37) R. M. Wald, Black hole entropy is the Noether charge, Phys. Rev. D 48, no.8, R3427-R3431 (1993), arXiv:gr-qc/9307038 [gr-qc].
  * (38) V. Iyer and R. M. Wald, Some properties of Noether charge and a proposal for dynamical black hole entropy, Phys. Rev. D 50, 846-864 (1994), arXiv:gr-qc/9403028 [gr-qc].

##  Appendix A Full Analysis of Linechrowicz modes

In section 3, it was shown that once the Lichnerowicz equation
$\Delta_{L}\psi_{\mu\nu}=\lambda\psi_{\mu\nu}$, with eigenvalue
$\lambda=-1/(2\alpha)$, has a a transverse traceless eigenfunction
$\psi_{\mu\nu}$, there exists a linearized perturbation away from the Taub-NUT
solution. Here, we want to give a comprehensive analysis of the Lichnerowicz
equation. Since the perturbation is under the Taub-NUT black hole background
(3), we consider a symmetric tensor whose non-zero components are

$\begin{split}\psi_{00}=h\psi_{0}(r),~{}\psi_{03}=\psi_{30}=2nh\cos\theta\psi_{0}(r)\,,\\\
\psi_{11}=h^{-1}\psi_{1}(r),~{}\psi_{ij}=(r^{2}+n^{2})\gamma_{ij}\psi_{2}(r),\end{split}$
(12)

where $\psi_{0}$, $\psi_{1}$ and $\psi_{2}$ are unknown functions of radial
coordinate $r$ only and $\gamma_{ij}$ are the induced metric of $d\Omega^{2}$.
Imposing the tracefree condition $g^{\mu\nu}\psi_{\mu\nu}=0$ and the
transverse condition $\nabla^{\mu}\psi_{\mu\nu}=0$, we can obtain

$-\psi_{0}+\psi_{1}+2\psi_{2}=0,$ (13)

and

$\psi_{1}^{\prime}+\frac{h^{\prime}}{2h}(\psi_{0}+\psi_{1})+\frac{2r}{r^{2}+n^{2}}(\psi_{1}-\psi_{2})=0.$
(14)

With these equations, we can express $\psi_{0}$ and $\psi_{2}$ in terms of
$\psi_{1}$ and $\psi_{1}^{\prime}$. Then the $rr$ direction of the
Lichnerowicz equations $\Delta_{L}\psi_{\mu\nu}=\lambda\psi_{\mu\nu}$ gives a
equation of second derivative of $\psi_{1}$

$\begin{split}&\frac{2\left(2m^{2}\left(7n^{2}r^{2}+n^{4}-6r^{4}\right)+m\left(-30n^{2}r^{3}+7n^{4}r+11r^{5}\right)-6n^{4}r^{2}+13n^{2}r^{4}+3n^{6}-2r^{6}\right)\psi_{1}^{\prime}}{\left(n^{2}+r^{2}\right)^{2}\left(-mn^{2}+3mr^{2}+3n^{2}r-r^{3}\right)}\\\
&+\frac{\left(r^{2}-2mr-n^{2}\right)\psi_{1}^{\prime\prime}}{n^{2}+r^{2}}+\left(\lambda-\frac{8\left(4m^{2}n^{2}r+m\left(-6n^{2}r^{2}+n^{4}+r^{4}\right)+4n^{2}r^{3}\right)}{\left(n^{2}+r^{2}\right)^{2}\left(m\left(n^{2}-3r^{2}\right)-3n^{2}r+r^{3}\right)}\right)\psi_{1}=0.\end{split}$
(15)

In the large $r$, the equation turns out to be

$\psi_{1}^{\prime\prime}+\frac{4}{r}\psi_{1}^{\prime}+\lambda\psi_{1}=0\,.$
(16)

This equation can be easily solved and two Yukawa like
modes,$\psi_{1}~{}e^{\pm\sqrt{-\lambda r}}$, emerge. Generally, the full
equation can not be solved analytically, and for numerical solution, one need
to explore the appropriate values of the parameters to eliminate the diverging
mode $e^{+\sqrt{-\lambda r}}$ through the shooting method, and only the
normalisable mode $e^{-\sqrt{-\lambda r}}$ survives.

To numerically solve the quation (15), we first do a near-horizon Taylor
expansion for $\psi_{1}$ and then use this as initial data at the horizon to
perform numerical integration out to large $r$. Without loss of generality, we
set $\alpha=1/2$, and find that there exists one mode corresponding to each
value of $n$, and the relation between the starting point $r_{0}$ of the
linear perturbation and $n$ are shown in Fig 1.