# The role of electron inertia in Laser-Plasma Interaction

James R. Young  Pierre-Alexandre Gourdain  Department of Physics and
Astronomy, University of Rochester, Rochester, New York 14627, USA
Laboratory for Laser Energetics, University of Rochester, Rochester, New York
14627, USA
.

###### Abstract

Incorporating lasers in macroscopic plasma models is made challenging by the
fast temporal and small spatial scales when electrons interact with laser
fields. When ion motion can be neglected, like in laser wakefield
acceleration, particle-in-cell models capture most of features present in the
plasma dynamics. But when the ion dynamics takes central stage, only multiple-
fluid models can be used, as particle-in-cell simulations become too expensive
to run. In this letter, we show that extended magnetohydrodynamics can be used
instead, as it captures exactly the macroscopic coupling between laser fields,
electron flows and the ion dynamics, on the ion time scale. While it
replicates several features usually associated with explicit electron physics,
such as cut-off density, or electromagnetic wave reflection and absorption, it
does so with little computational overhead compared to standard
magnetohydrodynamics. Since this approach computes the electron dynamics
implicitly using generalized Ohm’s law, rather than electron momentum
conservation, the spatial and temporal scales of extended magnetohydrodynamics
fit well between multi-fluid and standard magnetohydrodynamics scales,
allowing to study a new class of problems that would be too expensive to run
otherwise.

††preprint: APS/123-QED

## introduction

Laser-plasma interactions (LPI) are key features when a high-power laser
drives materials into ionized states. For example they can cause electron and
ion density modulations that grow and steepen[1, 2, 3] and whose sharp nature
has been used to form plasma gratings[4, 5, 6]. Additionally, stimulated
Brillouin [7] combined with ponderomotively driven plasma response has long
been a major concern for inertial confinement fusion[8, 9], especially when
cross-beam energy transfer (CBET)[10, 11] is involved. In addition to
providing a source of magnetization in plasmas[12, 13], the very nature of
LPIs can change in the presence of strong magnetic fields[14], like those
generated around magnetars[15] and white-dwarves [16], and prior research even
links fast radio bursts in the relativistic winds of magnetars[17, 18] to
ponderomotive self-focusing of electromagnetic radiation[19, 20].

In one regime, when the ions are mostly frozen in place, a common situation in
laser wakefield acceleration simulations, particle-in-cell models [21] capture
LPIs extremely well. In the other regime, when the ion and electron motion is
ambipolar, the essence of LPIs is deeply rooted in inverse Bremsstrahlung
(IB)[7, 22], where ray-tracing (RT) [23, 24] can be used to track the energy
budget left inside the electric field as the laser beam travels inside the
plasma. While the first approach is collisionless and uses the Vlasov
approximation, the second is collisional and uses conservation laws to track
the ion dynamics[25]. State-of-the-art MHD codes (e.g. FLASH [26], LILAC [27],
DRACO [28], HYDRA [29]) often drop some or all of these electron terms. They
split the LPI into a radiation component (eikonal approximation of the laser
field) and an absorption component (i.e. IB) while only coupling the electron
and rays via the plasma dispersion relation [30]. Unfortunately, problems that
fall in between regimes are more difficult to model. Hybrid codes [31] have
been used recently to bridge this scale gap, running a multi-fluid
approximation when possible and switching to kinetic when the fluid approach
breaks down.

However when LPIs occur in a fully formed plasma [32], such as those produced
in modern pulsed-power machines [33, 34, 35, 36], the ray approximation
usually breaks down. The large density gradients and complex geometrical
structures present before the laser is fired will generate situations where
the paraxial assumption is violated and ray-tracing cannot be used. Up to now,
such situations relied on blending multiple computational approaches to study
these regimes[37, 38]. While hybrid codes or multi-fluid codes could be used
instead, this letter shows that the physics necessary to capture LPIs is fully
accounted for by the generalized Ohm’s law (GOL), where electron terms have
been retained. We show this model correctly captures IB, the primary delivery
mechanism for low energy lasers, and transitions away from IB when the laser
power is increased. While we restrict our discussion to non-relativistic
plasma interaction using high power (i.e. long pulse) lasers, this framework
can be extended to relativistic plasmas by modifying the conservation laws and
the GOL [39].

The first section uses dimensional analysis to show which terms in the GOL are
important to LPIs in the extended magnetohydrodynamics (XMHD) framework. Then,
we show quantitatively how XMHD handles LPIs in slab geometries using PERSEUS
[40, 41]. Finally, we highlight how energy transfer compares to basic analytic
theory.

##  I The Generalized Ohm’s law scaled for laser-plasma interactions

###  I.1  Redimensionalization of the generalized Ohm’s law using LPI
characteristic scales

Dimensional analysis comes into play to make important parameters appear in
XMHD equations [40]. Removing all units in conservation laws and Maxwell’s
equation can be done with taking three independent characteristics scales such
as time $t_{0}$, speed $v_{0}$ and mass density $\rho_{0}$ [42]. One key step
is necessary to arrive to a self consistent result. Conservation of momentum
dictates that the characteristic speed $v_{0}$ should also be the
characteristic Alfvén speed. In the framework of laser plasma interactions, we
can take the electromagnetic wave speed $c$ as our characteristic speed,
$v_{0}=c$. Following Faraday’s law, we get relate the electric and magnetic
fields characteristic scales as $E_{0}=cB_{0}$. As we remove dimensions from
each XMHD equation, we find the following relationships between the different
characteristic scales: the thermal pressure scale is $p_{0}=\rho_{0}c^{2}$,
the resistivity scale $\eta_{0}=t_{0}/\varepsilon_{0}$, the magnetic field
scale is $B_{0}=\sqrt{\rho_{0}/\varepsilon_{0}}$, and the current density
scale is $j_{0}=\sqrt{\rho_{0}/\mu_{0}}/t_{0}$.

We assume here quasi-neutrality and separate electron and ion energies.
However the electron momentum is computed implicitly, using the time evolution
of the electrical current density via the GOL, given by

$\displaystyle\mathbf{E}=-\mathbf{v}\times\mathbf{B}+\eta\textbf{j}+\frac{1}{\omega_{pe}t_{0}}\sqrt{\frac{m}{n_{e}}}(\textbf{j}\times\textbf{B}-\vec{\nabla}p_{e})$
$\displaystyle+\left(\frac{1}{\omega_{pe}t_{0}}\right)^{2}\left(\frac{\partial\mathbf{j}}{\partial
t}+\vec{\nabla}\cdot(\mathbf{vj+jv})\right)$ (1)
$\displaystyle-\left(\frac{1}{\omega_{pe}t_{0}}\right)^{3}\sqrt{mn_{e}}\,\vec{\nabla}\cdot\left(\frac{1}{n_{e}}\mathbf{jj}\right).$

Here E is the electric field, v is the flow speed, B is the magnetic field, j
is the current density, $n_{e}$ is the electron number density, $m$ is the
ratio of the ion mass to the electron mass, $\rho$ is the mass density, $p$ is
the thermal pressure, $\eta$ is the electric resistivity and $\epsilon$ is the
energy density, all dimensionless.

The dimensional electron plasma frequency
$\omega_{pe}=\sqrt{N_{e}e^{2}/(m_{e}\varepsilon_{0})}$ was highlighted by the
nondimensionalization process. And it is now in front of all the terms
connected to electron physics inside Eq. (I.1). Further, $N_{e}$ is electron
number density, e is the elementary charge and $m_{e}$ is the electron mass,
all dimensional.

###  I.2  Characterization of LPIs in XMHD

To understand what type of LPI takes place in XMHD, we can look at the
interaction of a plane electromagnetic wave traveling to the left with a
plasma slab located in the region $x<0$. At early times, the ion fluid is
static and has no spatial dependence in the region $x<0$. We can equate the
laser electric field at the plasma-vacuum interface, the LHS of Eq. (I.1), to
the electric field inside the plasma, the RHS of Eq. (I.1). In this case, Eq.
(I.1) simplifies to

$\mathbf{E}\approx\eta\textbf{j}+\frac{1}{\omega_{pe}^{2}t_{0}^{2}}\frac{\partial\mathbf{j}}{\partial
t},$ (2)

dropping the third order terms in $1/\omega_{pe}t_{0}$ since
$t_{0}>>1/\omega_{pe}$ because the characteristic time scale of the
interaction $t_{0}$ is much larger than the electron time scale (i.e.
$1/\omega_{pe}$). In other terms, we look at effects that are much slower than
the electron plasma frequency. Note that we also dropped the first order term
in $1/\omega_{pe}t_{0}$ only because

$\mathbf{j}\times\mathbf{B}=\vec{\nabla}p_{e},$

from conservation of momentum. At this early point in time the ions have not
moved yet and the Lorentz force equilibrates exactly the electron pressure
gradient.

Modifying Maxwell-Ampère’s law using Eq. (2) while taking $\eta$ and
$\omega_{pe}$ to have small time variations compared to $t_{0}$, we get

$\frac{1}{\omega_{pe}^{2}t_{0}^{2}}\frac{\partial^{2}\mathbf{j}}{\partial
t^{2}}+\eta\frac{\partial\textbf{j}}{\partial
t}+\textbf{j}\approx\nabla\times\textbf{B}$ (3)

Eq. (3) is the equation of a damped harmonic oscillator with a resonant
frequency, $\omega_{pe}$, driven by the laser magnetic field B.

This dimensional analysis shows that LPIs are mediated by electron inertia and
Maxwell-Ampère’s law in XMHD. So, the laser energy is transferred to the
electron fluid via harmonic oscillation. This is possible because we can now
store mechanical energy inside the electron fluid, a transformation that
requires the electrons having mass (hence inertia). Resonance allows this
energy to become large enough so it can interact non-linearly with the
incoming electromagnetic wave. Consequently, XMHD recovers some common
features found in the two-fluid models[25, 43], while resolving the electron
motion on the ion time scale only.

A fully reflecting laser beam becomes especially important for the formation
of standing waves, which leads to a ponderomotive steepening [3] first
introduced by Estabrook [1] as a formation of ion and electron density peaks
due to a laser pulse reflection from overdense material. We will demonstrate
in a future work that the ion-scale ponderomotive action does occur in XMHD
through the Hall term of Eq. (I.1). And so by including electron inertia, an
important ingredient for CBET is also naturally present.

Although electron inertia has been important in studies of collisionless
magnetic reconnection [44, 45, 46, 47], often the components are broken to
include only the the convective acceleration term
($\mathbf{v_{e}}\cdot\nabla\mathbf{j}$) rather than considering temporal
variations of current directly [48], which are key in the LPIs energy transfer
as we just saw.

##  II XMHD simulations

Figure 1:  The electric field $\mathrm{E_{z}}$ output from the simulation
compared to the components computed from Generalized Ohm’s law. The electron
inertial component ($\mathrm{\frac{\partial j}{\partial t}}$) dominates in the
higher temperature upstream region, but the ohmic term becomes important when
looking at the downstream region. Also plotted is the $\mathrm{N_{ion}}$ to
show the validity across a plasma-vacuum interface.

We now extend these results to much longer time scales so the ions have the
time to move. A non-relativistic ($\mathrm{a_{o}<0.01}$) s-polarized laser was
propagated into a simulation domain by setting an oscillatory electromagnetic
field as a boundary condition. The magnetic field component, $B_{y}$, was
defined on the vertical boundary so the Poynting flux was directed to the
left. The beam has a uniform spatial profile,
$\mathrm{B_{y}=B_{max}\cos{\left(kx+\omega t\right)}}$. Here,
$\mathrm{E_{max}=\sqrt{2I_{laser}/(c\epsilon_{0})}}$,
$\mathrm{k=2\pi/\lambda}$, and $\mathrm{\omega=ck}$ with $\mathrm{I_{laser}}$
varying from $\mathrm{10^{16}\rightarrow 10^{18}\text{Wm}^{-2}}$
($a_{0}\approx 0.001\rightarrow 0.01$) and $\mathrm{\lambda=527\text{nm}}$.

This work will mostly focus on an under-dense regime as this allows us to
investigate both resonance and collisional absorption with a single domain.
Additionally, since the plasma response can lead to additional electromagnetic
waves, results are shown from both a continuous laser and one with 6-periods
only. The latter allows us to observe re-emitted waves without resorting to a
separate vacuum solution.

Although we restrict this work to fluid-like resistivity and distribution
functions, the Spitzer model generally needs to be modified for high- and low-
frequency lasers [22], as well as for varying ionization states [49] and
densities [50]. The deviation from a Maxwellian towards a super-Gaussian can
be quantified with the Langdon parameter
($\alpha=Zv_{osc}^{2}/v_{thermal}^{2}$) in the so-called Langdon effect [51,
52]. Although corrections can be made to account for super-Gaussian
distribution functions, this work restricts itself to $\alpha\approx 1$ as
previous work has shown the collisional absorption error to be relatively
small [53]. This restriction allows us to focus on more fundamental
electromagnetic coupling.

In XMHD, the coupling occurs directly between the laser fields and plasma
itself. Fig. 1 shows the electric field ($\mathrm{E_{z}}$) computed by the
code in an initially isotropic under-dense step-density plasma. This field is
primarily composed of electron inertia
($\mathrm{\frac{\partial\mathbf{j}}{\partial t}}$) and the ohmic terms of the
GOL. Initially, displacement currents dominate inside the low-density/vacuum
portion of the plasma domain. Then physical currents take over as the density
is large enough to supply electrons easily. Although this would not be a
surprising result for a purely explicit 2-fluid code, we see that our semi-
implicit 1-fluid (2-energy) XMHD formulation also captures this interaction.

Figure 2:  The main cumulative contributions of the energy equation for
electron fluids computed using $\mathrm{\mathbf{E}\cdot\mathbf{j}}$. Also
shown is the instantaneous partition of energy for electron fluid is as a sum
of thermal ($\mathrm{k_{b}T}$) and directed kinetic ($\mathrm{mv^{2}}$). The
inset plot shows a zoomed in section at early time to demonstrate resonant
absorption increases prior to collisional heating. Energy has been normalized
for this plot. The geometry here is a wedge-density profile pre-heated to 1 eV
and with constant $\mathrm{ln\Lambda}=1$. The laser power is $S=3.5\times
10^{17}$ $\mathrm{W/m^{2}}$ ($a_{0}\approx 0.003$)

.

We now look at how the energy is transferred from the laser field to the
electrons. This analysis again used a step-density profile for simplicity.
However, we will use results of both a continuous and 6-period laser pulse.

While there is collisional momentum transfer between ions and electrons, a
detailed analysis found it relatively insignificant from the electron
perspective. The major contributions of energy exchange in LPIs is best
described by considering Poynting’s theorem,

$\mathrm{\frac{d}{dt}(u_{E}+u_{B})+\nabla\cdot\mathbf{S}=-\mathbf{j}\cdot\mathbf{E}}.$
(4)

where $u_{E}$ and $u_{B}$ are the electric and magnetic energy densities,
respectively. The heating term on the RHS only gets turned into
$\mathbf{E}\cdot\mathbf{j}\approx\eta j^{2}$ when
$\mathbf{E}\approx\mathbf{\eta j}$. However, the ohmic term in Eq. (I.1) does
not dominate the electric field in an under-dense plasma (see Fig. 1).
Instead, Fig. 2 shows the laser field heats the plasma via resonant
absorption, forcing oscillatory kinetic energy into the electron fluid. This
energy is then quickly thermalized through resistive collisions.

The inset plot shows the plasma first gains kinetic energy before resistivity
initiates energy transfer. Under these assumptions the plasma behaves as an
RLC circuit. The oscillatory motion stores reactive energy, while the
collisional energy corresponds to dissipation. When the 6-period laser turns
off, the remaining reactive oscillatory energy smoothly turns into thermal
energy.
$\mathrm{[\mathbf{E}\cdot\mathbf{j}]_{reactive}\approx\frac{m_{e}}{e^{2}N_{e}}j\frac{\partial
j}{\partial t}=\frac{\partial{\frac{1}{2}\rho_{e}v_{e}^{2}}}{\partial t}}$ is
just the kinetic energy given to the electron fluid by the laser. If
$\mathrm{ln\Lambda}$ is allowed to vary, then the plasma can store much more
oscillatory energy as the increase in temperature lowers $\mathrm{ln\Lambda}$
and thereby reduces $\eta$ and increases the current density.

Figure 3:  The scaling factor between a basic inverse Bremsstrahlung model and
the predictions from an XMHD model for two different laser powers with a
plasma pre-heated to 1 eV (to ensure ionization), (a)
$\mathbf{S}=\mathrm{3.5\times 10^{16}}$ $\mathrm{W/m^{2}}$ and (b)
$\mathbf{S}=\mathrm{3.5\times 10^{17}}$ $\mathrm{W/m^{2}}$. $\mathbf{S}$
refers to the Poynting flux, with the green curve computed from XMHD outputs,
the dotted lines computed from analytic IB theory at two different
temperatures, whose location in the plasma is shown in the filled-in circles.
The vertical axes are all unitless and refer to the fraction of maximum value
for each quantity, except the density, which is the fraction of critical
density. Additionally, the red curve shows the fraction of directed kinetic
energy ($\mathrm{\frac{E_{k}}{E_{T}}}$) in the system relative to the total
energy (including thermal).

##  III Connection to Inverse Bremsstrahlung

Figure 4:  Comparison between the final incident edge electron-temperature as
predicted by a simplified inverse Bremsstrahlung calculation from a pulsed
laser and the measurement from simulation using a wedge-density profile pre-
heated to 1 eV. The dotted red line shows a fit of the data to a polynomial
form predicted by IB theory from Eq. (7), $\mathrm{(ax+b)^{c}}$. The black
dots represent the final temperature for individual simulation predictions for
a laser power density specified on the x-axis. The inset plot shows an example
edge temperature evolution for the simulation in black vs. theory in red.

To better define the role IB plays in this type of model we compare it against
some basic analytic theory [7, 54, 55]. In this context, the dispersion
relation,

$\mathrm{\omega^{2}=k^{2}c^{2}+\omega_{pe}^{2}-i\frac{\omega\omega_{pe}^{2}}{\nu_{ei}}},$
(5)

shows lasers have their electromagnetic waves damped at a rate
$\mathrm{k_{ib}=\nu_{ei}\omega_{pe}^{2}/\mu c\omega^{2}}$, which leads to the
laser intensity (S) decaying as $\mathrm{dS/dx=-k_{ib}S}$ with a corresponding
solution, $\mathrm{\mathbf{S}=\mathbf{S}_{0}e^{-ax}sin(bx+c)}$. $\mu$ is the
real refractive index and $\mathrm{k_{ib}}$ is the inverse Bremsstrahlung
coefficient.

Using a continuous laser, Fig. 3 shows the Poynting flux,
$\mathrm{\mathbf{S}_{xmhd}}$, computed directly from outputs of PERSEUS along
with the damped exponential function, $\mathrm{\mathbf{S}_{T_{min/max}}}$,
with constants matching the IB theory from Eq. (5). However, a choice of what
rapidly varying $T_{e}$ to use must be made when computing $\nu_{ei}$. The
color and spatial location of the solid circles show values and location used
for corresponding curves, $\mathrm{\mathbf{S}_{T_{min}}}$ and
$\mathrm{\mathbf{S}_{T_{max}}}$.

While $\mathrm{\mathbf{S}_{T_{max}}}$ in Fig. 3a matches
$\mathrm{\mathbf{S}_{xmhd}}$ well, Fig. 3b clearly does not show the same. In
fact, the decay does not match well the functional form. From Fig. 1 we see
the electric field in the plasma is primarily composed of electron inertia for
the first few wavelengths. This along with Eq. (4), shows the absorbed
electromagnetic energy is partitioned into reactive and thermal kinetic
energy, a clear departure from IB theory. $T_{e}$ rapidly increases in the
more intense laser fields. And so $\eta\propto T^{-3/2}$ becomes low enough to
nearly stop collisional absorption altogether. Instead, energy is pumped into
the directed kinetic energy of the plasma, where it shows up as an electron
current, $\mathrm{\mathbf{j}_{L}}$. Fig. 3b shows this effect as the
$\mathrm{E_{k}/E_{T}}$ (ratio of directed-kinetic to total energy) curve
increases relative to Fig. 3a. Since Eq. (5) requires immediate thermalization
of laser-deposited energy, then its solution is no longer compatible with the
system. Since the energy is retained as directed kinetic energy, dynamics
associated with $\mathrm{\mathbf{j}_{L}\times\mathbf{B}}$ can now potentially
takeover.

To isolate incoming laser from the reflected/excited waves we now turn to
results using only a 6-period input laser. Now there is a definite ”laser-off”
portion, which shows how energy is partitioned among the electron fluid.
Similar to an approach in an earlier work [56], we computed the edge-
temperature using a 1D method. After integrating the electron energy equation,
Eq. (6), we get $\mathrm{T_{e}}$ as a function of the laser fluence,
$\mathrm{\Phi}$, shown in Eq. (7).

$\mathrm{\frac{3}{2}n_{e}k_{b}\frac{\partial T_{e}}{\partial
t}=\frac{\kappa}{T_{e}^{3/2}}S},$ (6)

$\mathrm{T_{e}(x=0,t)=(\frac{5}{3}\frac{\kappa}{n_{e}k_{b}}\Phi_{0}+T_{e0}^{5/2})^{2/5}},$
(7)

$\mathrm{\Phi_{0}=\int_{0}^{t}S(z=0,t^{\prime})dt^{\prime}}.$ (8)

This method assumes constant Coulomb logarithm ($\mathrm{ln\Lambda}$), Spitzer
collisions, no thermal conduction and immobile ions. This first assumption may
not be particularly valid, since when temperature increases,
$\mathrm{ln\Lambda}$ can also vary significantly. In our simulations
$\mathrm{ln\Lambda}$ falls between 0 and 3, so we chose a constant value of 1.

Fig. 4 shows there is a good agreement of the edge-temperature at the plasma-
vacuum region when plotting the value from Eq. (7) (using the laser poynting
flux, S) along with the actual rise in temperature produced by the simulation,
indicating XMHD captures IB relatively well.

The spatial shape of $\mathrm{T_{e}}$ is different from the analytic theory
because electron inertia implies a time lag. This along with Spitzer
resistivity, which is ultimately velocity-based, implies there is an inherent
time required to pump energy into electron oscillations before they can
thermalize it.

##  IV  Conclusions

This letter explores how XMHD can be used to study the intricacies of LPIs
across a broad range of parameters, that remains elusive to particle-in-cell
and MHD models. While this work explores fundamental effects that are well
understood, it does it by using a completely different approach. Leveraging
the GOL, XMHD can capture the oscillatory nature of the electrons and shows
how and where resonance comes into play in the proposed framework. We have
shown that LPI is a mix of reactance, via resonance, and dissipation, via
resistivity. The reactive component stems from harmonic oscillations which
stores mechanical energy inside the electron motion. This is only possible
because electron have a non-zero mass in our model. The model can recover many
fundamental aspects of LPI usually captured by particle-in-cell or multi-fluid
MHD codes, such as reflection, absorption and even a critical density. Energy
transfer and its subsequent decay match well analytic theory for simple
configurations, and at low energies. Because of its accuracy, XMHD is
appropriate for regimes that can be too expensive to simulate with multi-fluid
codes and that MHD cannot capture.

Since XMHD accounts for the entire $d\mathbf{j}/dt=\partial\mathbf{j}/\partial
t+\mathbf{v}\cdot\nabla\mathbf{j}$, there is a potential richness of spatio-
temporal conditions. This framework naturally bridges the gap between the
kinetic and hydrodynamics scales, it also ensures computational continuity
between the short scales, dominated by electron physics, and the large scale,
dominated by MHD. As a result, we can compute the effects of microscopic
physics and study how they impact the microscopic scales by using the same
code, especially when using adaptive mesh refinement[57]. Under such
assumptions, using one code to compute LPI and how they impact plasma dynamics
on larger scales is clearly possible and will be explore in future work.

## Funding

This research was supported by the NSF CAREER Award PHY-1943939 and by the
Laboratory for Laser Energetics Horton Fellowships.

## References

  * Estabrook  _et al._ [1975] K. G. Estabrook, E. J. Valeo, and W. L. Kruer, Two-dimensional relativistic simulations of resonance absorption, Physics of Fluids 18, 1151 (1975).
  * Dragila and Krepelka [1978] R. Dragila and J. Krepelka, Laser plasma density profile modification by ponderomotive force, Journal de Physique 39, 617 (1978).
  * Smith  _et al._ [2019] J. R. Smith, C. Orban, G. K. Ngirmang, J. T. Morrison, K. M. George, E. A. Chowdhury, and W. M. Roquemore, Particle-in-cell simulations of density peak formation and ion heating from short pulse laser-driven ponderomotive steepening, Physics of Plasmas 26, 123103 (2019).
  * Lehmann and Spatschek [2019] G. Lehmann and K. H. Spatschek, Plasma photonic crystal growth in the trapping regime, Physics of Plasmas 26, 013106 (2019).
  * Sheng  _et al._ [2003] Z.-M. Sheng, J. Zhang, and D. Umstadter, Plasma density gratings induced by intersecting laser pulses in underdense plasmas, Applied Physics B 77, 673 (2003).
  * Plaja and Roso [1997] L. Plaja and L. Roso, Analytical description of a plasma diffraction grating induced by two crossed laser beams, Physical Review E 56, 7142 (1997).
  * Kruer [2003] W. L. Kruer,  _The physics of laser plasma interactions_ , Frontiers in physics (Boulder, Colo. [u.a.] Westview, 2003).
  * Bezzerides  _et al._ [1996] B. Bezzerides, H. X. Vu, and J. M. Wallace, Convective gain of stimulated Brillouin scattering in long‐scale length, two‐ion‐component plasmas, Physics of Plasmas 3, 1073 (1996).
  * Eliseev  _et al._ [1995] V. V. Eliseev, W. Rozmus, V. T. Tikhonchuk, and C. E. Capjack, Stimulated Brillouin scattering and ponderomotive self‐focusing from a single laser hot spot, Physics of Plasmas 2, 1712 (1995).
  * Hüller  _et al._ [2020] S. Hüller, G. Raj, W. Rozmus, and D. Pesme, Crossed beam energy transfer in the presence of laser speckle ponderomotive self-focusing and nonlinear sound waves, Physics of Plasmas 27, 022703 (2020).
  * Ruocco [2020] A. Ruocco,  _Modelling of ponderomotive laser self-focusing in a plasma with a hydrodynamic code in the context of direct-drive inertial confinement fusion_ , phdthesis, Université de Bordeaux (2020).
  * Gradov and Stenflo [1983] O. Gradov and L. Stenflo, Magnetic-field generation by a finite-radius electromagnetic beam, Physics Letters A 95, 233 (1983).
  * Shukla  _et al._ [2009] N. Shukla, P. K. Shukla, and L. Stenflo, Magnetization of a warm plasma by the nonstationary ponderomotive force of an electromagnetic wave, Physical Review E 80, 027401 (2009).
  * Ivanov  _et al._ [2021] V. V. Ivanov, A. V. Maximov, R. Betti, L. S. Leal, J. D. Moody, K. J. Swanson, and N. A. Huerta, Generation of strong magnetic fields for magnetized plasma experiments at the 1-MA pulsed power machine, Matter and Radiation at Extremes 6, 046901 (2021).
  * Palmer  _et al._ [2005] D. M. Palmer, S. Barthelmy, N. Gehrels, R. M. Kippen, T. Cayton, C. Kouveliotou, D. Eichler, R. a. M. J. Wijers, P. M. Woods, J. Granot, Y. E. Lyubarsky, E. Ramirez-Ruiz, L. Barbier, M. Chester, J. Cummings, E. E. Fenimore, M. H. Finger, B. M. Gaensler, D. Hullinger, H. Krimm, C. B. Markwardt, J. A. Nousek, A. Parsons, S. Patel, T. Sakamoto, G. Sato, M. Suzuki, and J. Tueller, A giant $\gamma$-ray flare from the magnetar sgr 1806-20, Nature 434, 1107 (2005).
  * Jordan  _et al._ [2007] S. Jordan, R. Aznar Cuadrado, R. Napiwotzki, H. M. Schmid, and S. K. Solanki, The fraction of DA white dwarfs with kilo-Gauss magnetic fields, Astronomy & Astrophysics 462, 1097 (2007).
  * Sobacchi  _et al._ [2023] E. Sobacchi, Y. Lyubarsky, A. M. Beloborodov, L. Sironi, and M. Iwamoto, Saturation of the Filamentation Instability and Dispersion Measure of Fast Radio Bursts, The Astrophysical Journal Letters 943, L21 (2023).
  * Ghosh  _et al._ [2022] A. Ghosh, D. Kagan, U. Keshet, and Y. Lyubarsky, Nonlinear Electromagnetic-wave Interactions in Pair Plasma. I. Nonrelativistic Regime, The Astrophysical Journal 930, 106 (2022).
  * Kaw [1973] P. Kaw, Filamentation and trapping of electromagnetic radiation in plasmas, Physics of Fluids 16, 1522 (1973).
  * Del Pizzo  _et al._ [1979] V. Del Pizzo, B. Luther-Davies, and M. R. Siegrist, Self-focussing of a laser beam in a multiply ionized, absorbing plasma, Applied physics 18, 199 (1979).
  * Hockney and Eastwood [1981] R. Hockney and J. Eastwood,  _Computer simulations using particles_ (McGraw-Hill, 1981).
  * Dawson and Oberman [1962] J. Dawson and C. Oberman, High‐Frequency Conductivity and the Emission and Absorption Coefficients of a Fully Ionized Plasma, The Physics of Fluids 5, 517 (1962).
  * Kaiser [2000] T. B. Kaiser, Laser ray tracing and power deposition on an unstructured three-dimensional grid, Physical Review E 61, 895 (2000).
  * Colaïtis  _et al._ [2021] A. Colaïtis, I. Igumenshchev, J. Mathiaud, and V. Goncharov, Inverse ray tracing on icosahedral tetrahedron grids for non-linear laser plasma interaction coupled to 3D radiation hydrodynamics, Journal of Computational Physics 443, 110537 (2021).
  * Braginskii [1965] S. Braginskii, Transport properties in a plasma, Review of plasma physics 1, 205 (1965).
  * Fryxell  _et al._ [2000] B. Fryxell, K. Olson, P. Ricker, F. X. Timmes, M. Zingale, D. Q. Lamb, P. MacNeice, R. Rosner, J. W. Truran, and H. Tufo, FLASH: An Adaptive Mesh Hydrodynamics Code for Modeling Astrophysical Thermonuclear Flashes, The Astrophysical Journal Supplement Series 131, 273 (2000).
  * Delettrez  _et al._ [1987] J. Delettrez, R. Epstein, M. C. Richardson, P. A. Jaanimagi, and B. L. Henke, Effect of laser illumination nonuniformity on the analysis of time-resolved x-ray measurements in uv spherical transport experiments, Physical Review A 36, 3926 (1987).
  * Keller  _et al._ [1999] D. Keller, T. J. B. Collins, J. A. Delettrez, P. W. McKenty, P. B. Radha, R. P. J. Town, B. Whitney, and G. A. Moses, Draco - a new multidimensional hydrocode, Bull. Am. Phys. Soc. 41, 39 (1999).
  * Marinak  _et al._ [1996] M. M. Marinak, R. E. Tipton, O. L. Landen, T. J. Murphy, P. Amendt, S. W. Haan, S. P. Hatchett, C. J. Keane, R. McEachern, and R. Wallace, Three‐dimensional simulations of Nova high growth factor capsule implosion experiments, Physics of Plasmas 3, 2070 (1996).
  * Stix [1992] T. H. Stix,  _Waves in plasmas_ (Springer Science & Business Media, 1992).
  * Holderied  _et al._ [2021] F. Holderied, S. Possanner, and X. Wang, Mhd-kinetic hybrid code based on structure-preserving finite elements with particles-in-cell, Journal of Computational Physics 433, 110143 (2021).
  * Young  _et al._ [2021] J. R. Young, M. B. Adams, H. Hasson, I. West-Abdallah, M. Evans, and P.-A. Gourdain, Using extended MHD to explore lasers as a trigger for x-pinches, Physics of Plasmas 28, 102703 (2021).
  * Shapovalov  _et al._ [2017] R. V. Shapovalov, R. B. Spielman, and G. R. Imel, An oil-free compact X-pinch plasma radiation source: Design and radiation performance, Review of Scientific Instruments 88, 063504 (2017).
  * Shapovalov  _et al._ [2019] R. Shapovalov, M. Adams, M. Evans, H. Hasson, J. Young, I. West-Abdallah, and P.-A. Gourdain, Low-Inductance Load Test of a New 250-Ka, 150-Ns Pulser for Fast X-Pinch Sources, in  _2019 IEEE Pulsed Power & Plasma Science (PPPS)_ (IEEE, Orlando, FL, USA, 2019) pp. 1–4.
  * Zucchini  _et al._ [2015] F. Zucchini, S. N. Bland, C. Chauvin, P. Combes, D. Sol, A. Loyen, B. Roques, and J. Grunenwald, Characteristics of a molybdenum x-pinch x-ray source as a probe source for x-ray diffraction studies, Review of Scientific Instruments 86, 033507 (2015).
  * Strucka  _et al._ [2022] J. Strucka, J. Halliday, T. Gheorghiu, H. Horton, B. Krawczyk, P. Moloney, S. Parker, G. Rowland, N. Schwartz, S. Stanislaus, _et al._ , A portable x-pinch design for x-ray diagnostics of warm dense matter, Matter and Radiation at Extremes 7, 016901 (2022).
  * Davies  _et al._ [2017] J. R. Davies, D. H. Barnak, R. Betti, E. M. Campbell, P.-Y. Chang, A. B. Sefkow, K. J. Peterson, D. B. Sinars, and M. R. Weis, Laser-driven magnetized liner inertial fusion, Physics of Plasmas 24, 062701 (2017).
  * Hansen  _et al._ [2018] E. C. Hansen, D. H. Barnak, R. Betti, E. M. Campbell, P.-Y. Chang, J. R. Davies, V. Y. Glebov, J. P. Knauer, J. Peebles, S. P. Regan, and A. B. Sefkow, Measuring implosion velocities in experiments and simulations of laser-driven cylindrical implosions on the OMEGA laser, Plasma Physics and Controlled Fusion 60, 054014 (2018).
  * Hamlin and Seyler [2016] N. D. Hamlin and C. E. Seyler, Relativistic modeling capabilities in perseus extended-MHD simulation code for HED plasmas, in  _2016 IEEE International Conference on Plasma Science (ICOPS)_ (IEEE, Banff, AB, Canada, 2016) pp. 1–1.
  * Seyler and Martin [2011] C. E. Seyler and M. R. Martin, Relaxation model for extended magnetohydrodynamics: Comparison to magnetohydrodynamics for dense Z-pinches, Physics of Plasmas 18, 012703 (2011).
  * Martin [2010] M. R. Martin, Generalized ohms law at the plasma-vacuum interface, Ph. D. thesis, Cornell University (2010).
  * Gourdain [2017] P.-A. Gourdain, The impact of the Hall term on tokamak plasmas (2017), arXiv:1703.00987 [physics].
  * Chapman and Cowling [1990] S. Chapman and T. G. Cowling,  _The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases_ (Cambridge university press, 1990).
  * Birn and Priest [2007] J. Birn and E. R. Priest, eds.,  _Reconnection of Magnetic Fields: Magnetohydrodynamics and Collisionless Theory and Observations_ (Cambridge University Press, Cambridge, 2007).
  * Drake  _et al._ [1997] J. F. Drake, D. Biskamp, and A. Zeiler, Breakup of the electron current layer during 3-D collisionless magnetic reconnection, Geophysical Research Letters 24, 2921 (1997).
  * Jones and Parker [2003] S. T. Jones and S. E. Parker, Including electron inertia without advancing electron flow, Journal of Computational Physics 191, 322 (2003).
  * Muñoz  _et al._ [2018] P. Muñoz, N. Jain, P. Kilian, and J. Büchner, A new hybrid code (CHIEF) implementing the inertial electron fluid equation without approximation, Computer Physics Communications 224, 245 (2018).
  * Kuznetsova  _et al._ [1998] M. M. Kuznetsova, M. Hesse, and D. Winske, Kinetic quasi-viscous and bulk flow inertia effects in collisionless magnetotail reconnection, Journal of Geophysical Research: Space Physics 103, 199 (1998).
  * Skupsky [1987] S. Skupsky, “Coulomb logarithm” for inverse-bremsstrahlung laser absorption, Physical Review A 36, 5701 (1987).
  * Lee and More [1984] Y. T. Lee and R. M. More, An electron conductivity model for dense plasmas, Physics of Fluids 27, 1273 (1984).
  * Langdon [1980] A. B. Langdon, Nonlinear Inverse Bremsstrahlung and Heated-Electron Distributions, Physical Review Letters 44, 575 (1980).
  * Turnbull  _et al._ [2020] D. Turnbull, A. Colaïtis, A. M. Hansen, A. L. Milder, J. P. Palastro, J. Katz, C. Dorrer, B. E. Kruschwitz, D. J. Strozzi, and D. H. Froula, Impact of the Langdon effect on crossed-beam energy transfer, Nature Physics 16, 181 (2020).
  * Decker  _et al._ [1994] C. D. Decker, W. B. Mori, J. M. Dawson, and T. Katsouleas, Nonlinear collisional absorption in laser‐driven plasmas, Physics of Plasmas 1, 4043 (1994).
  * Schiavi [2004] A. M. M. Schiavi,  _Study of Laser Produced Plasmas by X-ray and Proton Radiography_ , Ph.D. thesis, Imperial College London (2004).
  * Johnston [1973] T. W. Johnston, Correct values for high-frequency power absorption by inverse bremsstrahlung in plasmas, Physics of Fluids 16, 722 (1973).
  * Denavit and Phillion [1994] J. Denavit and D. W. Phillion, Laser ionization and heating of gas targets for long‐scale‐length instability experiments, Physics of Plasmas 1, 1971 (1994).
  * Berger and Colella [1989] M. J. Berger and P. Colella, Local adaptive mesh refinement for shock hydrodynamics, Journal of computational Physics 82, 64 (1989).


Scalarization of Taub-NUT Black Holes in Extended scalar-tensor-Gauss-Bonnet
Theory

Hai-Shan Liu and Lei Zhang

Center for Joint Quantum Studies and Department of Physics,
School of Science, Tianjin University, Tianjin 300350, China

ABSTRACT

Recently, scalarization of Schwarzschild black hole are extensively studied.
In this work, we explore the scalarization of Taub-NUT black hole. The theory
we consider is the extended scalar-tensor-Gauss-Bonnet theory, which admits
Ricci-flat Taub-NUT black hole as a solution. An analysis of probe scalar
field is carried out to identify the mass parameter and NUT parameter
$(m\,,n)$ where the hairy black holes start to emerge. Then, we use shooting
method to construct the scalarized Taub-NUT black hole numerically. Being
different from the Schwarzschild case, there exists two branches of new hairy
black holes which are smoothly connected to each other. We calculate the
entropy of scalarized black holes and compare it with the entropy of scalar-
free Taub-NUT black holes, it turns out that the entropy of the new hairy
black holes are larger than that of scalar-free black holes. A novel phenomena
emerges in this system that the entropy of the black holes at the bifurcation
point is constant for positive mass parameter. We then conjecture a maximal
entropy bound for all the scalarized black hole whose mass parameter at the
bifurcation point is greater than zero.

<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Whether there exists black holes carrying scalar hair has been attracting much
attentions for years. And a no-hair theorem was proposed, which states that
Einstein’s general relativity can not support static black holes with a non-
trivial scalar[1, 2]. Recently, there is a big breakthrough in this direction.
Many black holes with non-trivial scalar field have been successfully
constructed in the extended scalar-tensor theories[3, 4, 5, 6, 7]. These black
holes are resulting from the scalarization procedure and are often called
scalarized black holes, more details can be found in[9, 8]. It is found that a
wide class of extended scalar-tensor theories can allow spontaneous
scalarization, which means that the theories can accommodate scalarized black
holes as well as the standard vacuum solutions of general relativity.

Schwarzschild black hole is the simplest black hole solutions of general
relativity, and the scalarization processes which are connected to
Schwarzschild black hole have been extensively studied. It is well-known that
Taub-NUT spacetime is also a solution of general relativity[10, 11]. The Taub-
NUT solution has two integration constants, the mass parameter and NUT
parameter, when the NUT parameter approaches zero, the solution returns back
to Schwarzschild black hole solution. Though the Taub-NUT black hole solution
was discovered in the middle of last century, the thermodynamics of this
solution is still controversial, especially, the definition and physical
interpretation of its mass and NUT charge. There are some new progress on this
subject in the last few years [12, 13, 14], triggering a new surge of interest
in this area. In this work, we shall study the scalarization of black holes
which are connected to Taub-NUT black holes, providing a wider platform to
explore the thermodynamics of Taub-NUT like black holes.

We try to obtain Taub-NUT black holes with non-trivial scalar fields through
coupling the scalar field directly to spacetime curvature invariants. In
particular, we shall consider the extended scalar-tensor-Gauss-Bonnet theories
in four dimensions. It is shown that this theory can support spontaneous
scalarization of Schwarzschild black holes for certain classes of coupling
functions, the Schwarzschild black hole will become unstable below a critical
value of mass, then new scalarized black holes will emerge. We shall choose
one coupling function such that the theory can admit Taub-NUT black hole as a
solution. A probe scalar field analysis will be carried out to identify the
bifurcation points where scalarized black holes start to appear. Then we shall
apply shooting method to construct the scalarized Taub-NUT black holes
numerically.

The structure of this paper is as follows. In section 2, we present the
extended scalar-tensor-Gauss-Bonnet gravity theory with a particular coupling
function which admit Taub-NUT black hole as a solution and take a probe scalar
analysis to identify parameter range where the scalarized black holes start
emerge. In section 3, we analysis the behavior of Taub-NUT like black holes
near the horizon and in the infinity, respectively. We numerically construct
scalarized Taub-NUT black holes and study their properties in section 4.
Finally, a conclusion of our results are present in Section 5.

##  2 Extended scalar-tensor-Gauss-Bonnet theory

In this section, we consider extended scalar-tensor-Gauss-Bonnet theory, which
is given by

$S=\frac{1}{16\pi}\int
d^{4}x\sqrt{-g}\left(R-2\nabla_{\mu}\varphi\nabla^{\mu}\varphi-V\left(\varphi\right)+\lambda^{2}f\left(\varphi\right){R_{GB}}^{2}\right)\,,$
(1)

where $R$ is the Ricci scalar with respect of the metric $g_{\mu\nu}$,
$\varphi$ is the scalar field with the potential $V(\varphi)$, $f(\varphi)$ is
the coupling function with coupling constant $\lambda$ and $R_{GB}^{2}$ is the
Gauss-Bonnet invariant.

$R_{GB}^{2}=R^{2}+R_{\mu\nu\rho\lambda}R^{\mu\nu\rho\lambda}-4R_{\mu\nu}R^{\mu\nu}$
(2)

The equations of motion of the theory can be obtained from the variation of
the action(1) with respect to the metric $g_{\mu\nu}$ and the scalar field
$\varphi$,

$R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}+\varGamma_{\mu\nu}=2\nabla_{\mu}\varphi\nabla_{\nu}\varphi-
g_{\mu\nu}\nabla_{\alpha}\varphi\nabla^{\alpha}\varphi-\frac{1}{2}g_{\mu\nu}V\left(\varphi\right)\,,$
(3)
$\nabla_{\alpha}\nabla^{\alpha}\varphi=\frac{1}{4}\frac{dV(\varphi)}{d\varphi}-\frac{\lambda^{2}}{4}\frac{dF(\varphi)}{d\varphi}R_{GB}^{2}\,.$
(4)

Here $\nabla_{\mu}$ is the covariant derivative with respect to the metric
$g_{\mu\nu}$ and the $\varGamma_{\mu\nu}$ are defined by

$\begin{split}\varGamma_{\mu\nu}=&-2R\nabla_{\mu}\varPsi_{\nu}-4\nabla^{\alpha}\varPsi_{\alpha}(R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu})+8R_{\mu\alpha}\nabla^{\alpha}\varPsi_{\nu}\\\
&-4g_{\mu\nu}R^{\alpha\beta}\nabla_{\alpha}\varPsi_{\beta}+4R^{\beta}_{\mu\alpha\nu}\nabla^{\alpha}\varPsi_{\beta}\,.\end{split}$
(5)

with $\varPsi_{\mu}$ given by

$\varPsi_{\mu}=\lambda^{2}\frac{dF(\varphi)}{d\varphi}\nabla_{\mu}\varphi\,.$
(6)

For simplicity, we set the scalar potential $V(\varphi)$ to zero. It was shown
that the theory can admit Schwarzschild black hole solution, when the coupling
function satisfies ${\frac{\delta f}{\delta\varphi}}|_{\varphi=0}=0$[3].
Further, if ${\frac{\delta^{2}f}{\delta\varphi^{2}}}|_{\varphi=0}>0$,
scalarized black holes can emerge. Various coupling functions are studied for
scalarization of Schwarzschild black holes[3, 4] In this article, we choose
the coupling function as that in [3]

$F\left(\varphi\right)=\frac{1}{12}\left[1-\exp(-6\varphi^{2})\right]\,.$ (7)

It is easy to check that the above two conditions related to coupling function
$F$ are satisfied. If the scalar field is trivial $\varphi=0$, the theory will
return to general relativity, and the equations of motion admit Ricci-flat
Schwarzschild black hole as a solution.

As we mentioned in the introduction, the Ricci-flat Taub-NUT black hole
solution is also a solution of general relativity. Thus, Taub-NUT black hole
with trivial scalar $\varphi=0$ is also a solution for the extended scalar-
tensor-Gauss-Bonnet theory. The Taub-NUT metric has the form of[15]

$ds^{2}=-f\left(dt+2n\cos\theta
d\phi\right)^{2}+\frac{dr^{2}}{h}+\left(r^{2}+n^{2}\right)\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)\,$ (8)

with

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

The Taub-NUT solution has two integration constants, one is the mass parameter
$m$ and the other is $n$ which is usually called NUT parameter. As one can
see, when $n=0$ the Taub-NUT solution goes back to Schwarzschild solution. In
this article, we want to explore scalarization of black hole under Taub-NUT
background.

Following the strategy in [16], we first consider a test scalar field
$\delta\varphi$, and the scalar equation on Taub-NUT background is given by

$\nabla_{\alpha}\nabla^{\alpha}\delta\varphi+\frac{\lambda^{2}}{4}R_{GB}^{2}\delta\varphi=0\,,$
(10)

the Gauss-Bonnet combination $R_{GB}^{2}$ on the Taub-NUT background is

$\displaystyle R_{GB}^{2}=\frac{48}{\left(n^{2}+r^{2}\right)^{6}}$
$\displaystyle\left[m^{2}\left(-n^{6}+15n^{4}r^{2}-15n^{2}r^{4}+r^{6}\right)+4m\left(3n^{6}r-10n^{4}r^{3}+3n^{2}r^{5}\right)\right.$
$\displaystyle\left.+n^{8}-15n^{6}r^{2}+15n^{4}r^{4}-n^{2}r^{6}\right]\,.$
(11)

Then, it becomes an eigenvalue problem, under appropriate boundary conditions,
fixing $(n\,,\lambda)$ the equation can admit solutions for a discrete mass
parameter $m$. These points are the bifurcation points where the scalarized
black hole start to emerge. Specifically, performing a spherical harmonics
decomposition of the scalar field[17]
$\delta\varphi\left(r,\theta,\phi\right)=\sum_{lm}Y_{lm}\left(\theta,\phi\right)U_{l}\left(r\right)$,
the radial equation turns out to be

$\frac{1}{\left(r^{2}+n^{2}\right)}\frac{d}{dr}\left[\left(r^{2}-2mr-n^{2}\right)\frac{dU_{l}}{dr}\right]-\left[\frac{l\left(l+1\right)}{n^{2}+r^{2}}-\frac{\lambda^{2}}{4}Rgb^{2}\right]U_{l}=0\,.$
(12)

Here, we only consider the zero mode $l=0$. We require the probe scalar field
is regular on the horizon and have the form

$U_{l}\left(r\right)=u_{0}+\frac{3u_{0}\left(n^{2}-r_{h}^{2}\right)\lambda^{2}}{r_{h}\left(n^{2}+r_{h}^{2}\right)^{2}}\left(r-r_{h}\right)+{\cal
O}(r-r_{h})^{2}+\dots\,.$ (13)

The probe scalar equation can not be solved analytically, but can be solved
numerically. Given the value of $(n\,,\lambda)$, one can numerically obtain
the value of $m$. We plot these bifurcation points in Fig.1, which shows the
existence line of the scalarized black hole.

Figure 1: The line of bifurcation point where the new hairy black hole will
start to emerge, the point of ($n=0$, $m/\lambda=0.5864$) covers the result of
Schwarzschild case.

When $n=0$, $m/\lambda=0.5864$, it is consistent with the Schwarzschild
case[3]. Starting from these bifurcation points, we shall construct
numerically the hairy black hole in the next section.

##  3 Scalarization of Taub-NUT black hole

With the guidance in the previous section, we are now in a stage to obtain new
hairy Taub-NUT black holes in extended scalar-tensor-Gauss-Bonnet theory. We
consider the Taub-NUT like metric

$ds^{2}=-f\left(r\right)\left(dt+2n\cos\theta
d\phi\right)^{2}+\frac{dr^{2}}{h\left(r\right)}+\left(r^{2}+n^{2}\right)\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)\,,$ (14)

with $h(r)$ and $f(r)$ are undetermined functions of $r$. And the scalar field
is static and spherical $\varphi=\varphi(r)$. Substituting these into the
equations of motion, we can get four independent field equations

$\displaystyle\
4fh^{2}r\left[2n^{2}\varPsi_{r}+r(n^{2}+r^{2})\varPsi_{r}^{{}^{\prime}}\right]-3n^{2}f^{2}(n^{2}+r^{2})\left(-1+2\varPsi_{r}h^{{}^{\prime}}+4h\varPsi_{r}^{{}^{\prime}}\right)-$
$\displaystyle\
fh(n^{2}+r^{2})^{3}\varphi^{{}^{\prime}2}-fh(n^{2}+r^{2})\left[2n^{2}+r^{2}-6r^{2}\varPsi_{r}h^{{}^{\prime}}+4(n^{2}+r^{2})\varPsi_{r}^{{}^{\prime}}\right]+$
$\displaystyle\
f(n^{2}+r^{2})^{2}\left[1-\left(r+2\varPsi_{r}\right)h^{{}^{\prime}}\right]=0,$
(15) $\displaystyle\
4f^{2}h\left[n^{2}+(n^{2}+r^{2})^{2}\varphi^{{}^{\prime}2}\right]-h(n^{2}+r^{2})(n^{2}+r^{2}-4hr\varPsi_{r})f^{{}^{\prime}2}+2f^{2}r(n^{2}+r^{2})h^{{}^{\prime}}-$
$\displaystyle\
4n^{2}f^{3}\left(2\varPsi_{r}h^{{}^{\prime}}+4h\varPsi_{r}^{{}^{\prime}}-1\right)+f(n^{2}+r^{2})^{2}f^{{}^{\prime}}h^{{}^{\prime}}+2fhr(n^{2}+r^{2})\left(1-6\varPsi_{r}^{{}^{\prime}}h^{{}^{\prime}}\right)f^{{}^{\prime}}+$
$\displaystyle\
2fh(n^{2}+r^{2})^{2}f^{{}^{\prime\prime}}-8fh^{2}r(n^{2}+r^{2})\varPsi_{r}f^{{}^{\prime\prime}}-8fh^{2}\left[n^{2}\varPsi_{r}+r(n^{2}+r^{2})\varPsi_{r}^{{}^{\prime}}\right]f^{{}^{\prime}}=0,$
(16) $\displaystyle\
fh(n^{2}+r^{2})^{4}\varphi^{{}^{\prime}}f^{{}^{\prime}}+f^{2}\left(n^{2}+r^{2}\right)^{4}\varphi^{{}^{\prime}}h^{{}^{\prime}}+2f^{2}h\left(n^{2}+r^{2}\right)^{2}\left[2r(n^{2}+r^{2})\varphi^{{}^{\prime}}+(n^{2}+r^{2})^{2}\varphi^{{}^{\prime\prime}}\right]+$
$\displaystyle\
\left\\{fh(n^{2}+r^{2})^{2}\left[3r^{2}\lambda^{2}h^{{}^{\prime}}f^{{}^{\prime}}-2\lambda^{2}(n^{2}+r^{2})f^{{}^{\prime\prime}}\right]-3n^{2}\lambda^{2}ff^{{}^{\prime}}\left(n^{2}+r^{2}\right)^{2}\left(f^{{}^{\prime}}h+fh^{{}^{\prime}}\right)+\right.$
$\displaystyle\
\phantom{=\;\;}\left.4n^{2}\lambda^{2}f^{3}\left[2h\left(n^{2}-3r^{2}\right)+r\left(n^{2}+r^{2}\right)h^{{}^{\prime}}\right]+\lambda^{2}h(n^{2}+r^{2})(n^{2}+r^{2}-r^{2}h)f^{{}^{\prime}}+\right.$
$\displaystyle\
\phantom{=\;\;}\left.2\lambda^{2}fh^{2}r(n^{2}+r^{2})\left[2n^{2}f^{{}^{\prime}}+r(n^{2}+r^{2})f^{{}^{\prime\prime}}\right]-\lambda^{2}f\left(n^{2}+r^{2}\right)^{3}f^{{}^{\prime}}h^{{}^{\prime}}+\right.$
$\displaystyle\
\phantom{=\>\>}\left.24n^{2}rf^{2}h(n^{2}+r^{2})^{2}f^{{}^{\prime}}-6n^{2}\lambda^{2}f^{2}h\left(n^{2}+r^{2}\right)^{2}f^{{}^{\prime\prime}}\right\\}\frac{dF(\varphi)}{d\varphi}=0,$
(17) $\displaystyle\
(n^{2}+r^{2})h\left[r(n^{2}+r^{2})+2(n^{2}+r^{2}-3r^{2}h)\varPsi_{r}\right]f^{{}^{\prime}}-n^{2}f^{2}(n^{2}+r^{2}+8rh\varPsi_{r})-(n^{2}+r^{2})^{2}f-$
$\displaystyle\
(n^{2}+r^{2})fh\left[-r^{2}-6n^{2}f^{{}^{\prime}}\varPsi_{r}+(n^{2}+r^{2})^{2}\varphi^{{}^{\prime}2}\right]=0,$
(18)

where ′ is derivative with respect to radial coordinate $r$ and $\varPsi_{r}$
is defined by

$\varPsi_{r}=\lambda^{2}\frac{dF(\varphi)}{d\varphi}\frac{d\varphi}{dr}\,.$
(19)

It’s difficult to obtain analytic solutions for field equations (15)-(18) with
a nontrivial scalar field, thus, we will turn to numerical method and use
shooting method to construct new hairy black hole solutions. We shall explore
the behavior of the fields on both boundaries.

On the horizon $r_{h}$, we require the metric functions $h$ and $f$ vanish
simultaneously and impose regularity conditions for $\varphi$

$f\left(r\right)|_{r\rightarrow\\\ r_{h}}\rightarrow 0,\quad
h\left(r\right)|_{r\rightarrow\\\ r_{h}}\rightarrow
0,\quad\varphi\left(r\right)|_{r\rightarrow\\\ rh}\rightarrow\varphi_{h}\,.$
(20)

Taking Taylor expansion near the event horizon $r_{h}$, the equations of
motion tell us that the fields have the form

$\displaystyle\
h\left(r\right)=\frac{1+\delta}{r_{h}}\left(r-r_{h}\right)+{\cal
O}(r-r_{h})^{2}\,,$ (21) $\displaystyle
f\left(r\right)=\frac{r_{h}\left[e^{\left(-12\varphi_{h}^{2}\right)}\left(n^{2}+r_{h}^{2}\right)\delta+6\varphi_{h}^{2}\left(1+\delta\right)^{2}\lambda^{4}\right]}{6n^{2}\varphi_{h}^{2}\left(1+\delta\right)\lambda^{4}}\left(r-r_{h}\right)+{\cal
O}(r-r_{h})^{2}\,,$ (22)
$\displaystyle\varphi\left(r\right)=\varphi_{h}-\frac{e^{6\varphi_{h}^{2}}r_{h}\delta}{2\varphi_{h}\left(1+\delta\right)\lambda^{2}}\left(r-r_{h}\right)+{\cal
O}(r-r_{h})^{2}\,.$ (23)

There are three parameters on the horizon, $r_{h}$, $\delta$ and
$\varphi_{h}$. It is worth pointing out that the effect of the Gauss-Bonnet
coupling term appears in the leading term of the fields.

In the infinity, we take a large $r$ expansion and let the fields be

$\displaystyle\
h\left(r\right)=h_{0}+\frac{h_{1}}{r}+\frac{h_{2}}{r^{2}}+\frac{h_{3}}{r^{3}}+\frac{h_{4}}{r^{4}}+\frac{h_{5}}{r^{5}}+...$
(24) $\displaystyle
f\left(r\right)=f_{0}+\frac{f_{1}}{r}+\frac{f_{2}}{r^{2}}+\frac{f_{3}}{r^{3}}+\frac{f_{4}}{r^{4}}+\frac{f_{5}}{r^{5}}+...$
(25)
$\displaystyle\varphi\left(r\right)=\varphi_{0}+\frac{\varphi_{1}}{r}+\frac{\varphi_{2}}{r^{2}}+\frac{\varphi_{3}}{r^{3}}+\frac{\varphi_{4}}{r^{4}}+\frac{\varphi_{5}}{r^{5}}+...$
(26)

Substituting into the field equations (15-18), and solving the equations of
motion order by order in $1/r$, actually, we calculate up to order $1/r^{6}$.
It is turns out that $h_{1}=f_{1}$ and the coefficients
($h_{i}\,,f_{i}\,,\varphi_{i}$) can be expressed in terms of
($f_{1}\,,\varphi_{1}$) for $i\geq 2$. Usually, $f_{1}$ is related to black
hole mass while $\varphi_{1}$ is related to scalar charge, thus we rename them
to be $f_{1}=-2m$ and $\varphi_{1}=D$. We show the coefficients to fifth order
below

$\displaystyle\
h\left(r\right)=1-\frac{2M}{r}+\frac{D^{2}-2n^{2}}{r^{2}}+\frac{M\left(D^{2}+2n^{2}\right)}{r^{3}}+\frac{\left(4D^{2}M^{2}+6n^{4}+2n^{2}D^{2}\right)}{3r^{4}}$
$\displaystyle\quad\quad\quad+\frac{M}{12r^{5}}\left\\{63D^{4}-4D^{2}\left(142M^{2}+65n^{2}-4\lambda^{2}\right)+8\left(80M^{4}+12n^{2}-51n^{2}\lambda^{2}\right)\right.\phantom{=\;\;}$
$\displaystyle\left.\quad\quad\quad+8M^{2}\left(80n^{2}-21\lambda^{2}\right)\right\\}+...,$
(27) $\displaystyle
f\left(r\right)=1-\frac{2M}{r}-\frac{2n^{2}}{r^{2}}+\frac{M\left(D^{2}+6n^{2}\right)}{3r^{3}}+\frac{2D^{2}M^{2}+2n^{2}\left(D^{2}+3n^{2}\right)}{3r^{4}}$
$\displaystyle\quad\quad\quad+\frac{M}{60r^{5}}\left\\{55D^{4}-4D^{2}\left(130M^{2}+50n^{2}-28\lambda^{2}\right)+8\left(80M^{4}-51n^{2}\lambda^{2}\right)\right.\phantom{=\;\;}$
$\displaystyle\left.\quad\quad\quad+8M^{2}\left(80n^{2}-2\lambda^{2}\right)\right\\}+...,$
(28)
$\displaystyle\varphi\left(r\right)=\frac{D}{r}+\frac{DM}{r^{2}}-\frac{D\left[D^{2}-2\left(4M^{2}+n^{2}\right)\right]}{6r^{3}}+\frac{DM\left(6M^{2}-2D^{2}+3n^{2}\right)}{3r^{4}}$
$\displaystyle\quad\quad\quad+\frac{D}{120r^{5}}\left\\{9D^{4}+24\left[16M^{4}+n^{4}+3n^{2}\lambda^{2}+3M^{2}\left(4n^{2}-\lambda^{2}\right)\right]\right.\phantom{=\;\;}$
$\displaystyle\left.\quad\quad\quad-4D^{2}\left(58M^{2}+13n^{2}\right)\right\\}+...$
(29)

Here, we can see that the coupling constant $\lambda$ appears in order
$1/r^{5}$, which implies the Gauss-Bonnet term has minor effect in the
infinity.

We are now in a stage to construct new black holes numerically, we shall apply
numerical method to connect the horizon data to the ones in the infinity. One
choice is to integrate the solution from horizon to infinity. We can perform a
Taylor expansions on the horizon and use it as initial data, then integrate
the solution from it to the infinity. However, one can see from the complexity
of $h_{1}$ and $f_{1}$, we can not push the Taylor expansion to much higher
orders to guarantee the accuracy. Thus, we take a alternative direction of
integrating the solution from infinity to the middle, since the series
expansion in the infinity can be carried to a much higher order (We actually
calculate up to order of $1/r^{6}$, and higher than this order is reachable.).
The emergence of a black hole solution is signified by the simultaneously
vanishing of $f(r)$ and $h(r)$.

##  4 Numerical Results

With the guidance of the probe scalar field analysis under the Taub-NUT
solution, we start to search new Taub-NUT black hole around the bifurcation
point and do find a new branch Taub-NUT like black hole solutions with scalar
hair. In Fig.2, we show the metric functions $f\left(r\right)$,
$h\left(r\right)$ and scalar field $\varphi\left(r\right)$ for positive and
negative mass parameter respectively. It is obvious that the metric functions
have a peak higher than 1 for negative mass parameter, whilst there is no peak
for positive mass parameter. In Fig.3, we show the metric function $f$ for
fixed NUT parameter $n=1$ and various value of mass parameter. As can one see,
a peak emerges on the profile of metric function from bottom to top as the
mass decreases from a positive value to a negative value. The event horizon
radius becomes smaller as the mass parameter decreases, which is under
expectation.

(a) $m>0$ (b) $m<0$

Figure 2: The newly constructed Taub-NUT black hole solutions for two
different type of mass parameters. The left represents $m>0$, the right one is
for $m<0$, the profile (h,f) of which have a peak.  Figure 3: The metric
function $f\left(r\right)$ for fixed NUT parameter n and various mass
parameter $m$, the one with negative mass parameter has a peak.

We plot the mass parameter $m$ as a function of black hole radius $r_{h}$ in
Fig.4. Since the theory admits analytic Taub-NUT black hole as a solution, we
also include the $m-r_{h}$ curve of Taub-NUT black hole as a comparison. The
solid line represents Taub-NUT black hole, while the dashed line is the new
scalarized black hole. The new black holes start from the bifurcation point
which is consistent with the probe scalar analysis. There exists two
particular properties for our new constructed black holes compared with the
Schwarzschild case. One is that there is a minimal value of $r_{h}$ for the
new black holes. Another is there are two branches of new solutions and these
two new black holes smoothly connected to each other at the minimal radius.
The situation changes when we increase the value of NUT parameter $n$. There
is just one branch of new solution for large $n$, however, there still exists
a minimal value of black hole radius. We present the $m-r_{h}$ curve for
various NUT parameters in the bottom of Fig.4, including the $n=0$ case.

(a)  (b)

(c)

Figure 4: The mass parameter $m$ as a function of the horizon radius $r_{h}$
for different values of $n/\lambda$. (a) shows ($m-r_{h}$) curves for
$n=0,andn/\lambda=5/7$, the $n=0$ is for comparison. The dashed line
represents the scalarized black hole, whilst the solid line represents the
scalar-free black hole. The scalarized black hole for $n/\lambda=5/7$ has two
branches of solution, while the Schwarzschild case of $n=0$ has only one
branch. (b) shows the curve for $n/\lambda=5/3$, which has only one branch new
solution, the notation is the same, the dashed line represents the scalarized
black hole, the other is for scalar-free black hole. (c) shows ($m-r_{h}$)
curves for various NUT parameter.

Next, we want to explore the properties of the scalar hair of the new black
hole. It is obvious that the theory is invariant under
$\varphi\rightarrow-\varphi$, thus the scalar charge parameter $D$ is
symmetric. We first plot the scalar charge parameter $D$ as a function of mass
parameter $m$ for various NUT parameter in the top of Fig.5. We also plot the
case for $n=0$, which is closed, the scalar charge $D$ increases from zero to
maximal value and then decreases back to zero as mass parameter decreases.
Whilst, the curve is not closed for our new found black holes. The scalar
charge $D$ of our new black hole increases from zero at the bifurcation point
then terminates at certain value of $m$. From these curves, we can see that
every scalarized black hole has a maximal value of mass parameter, which is at
the bifurcation point. From this $D-m$ curve, we can not distinguish the two
branches of the new black hole solution, since they are smoothly connected.
Then, we plot the scalar charge parameter $D$ as a function of black hole
radius $r_{h}$ in the bottom of Fig.5, where we can see clearly that there
exists two branches of new black hole solutions with a minimal $r_{h}$ where
these two are connected. One can see that the second branch is short and not
obvious when $n$ is small (n=0 corresponding to Schwarzschild case, where
there is only one branch of hairy black hole solution.), and the second branch
is more obvious when $n$ is large. However, there exists only one branch black
hole solution when n is large enough.

(a)  (b)

(c)  (d)

Figure 5: The scalar charge $D$ as functions of the mass parameter $m$ is
shown in the top two figures, all the scalarized black hole have a maximal
value of mass parameter. The scalar charge $D$ as functions of the horizon
radius $r_{h}$ for different values of $n/\lambda$ in bottom two figures, one
can see the scalarized black holes have two branches which of smoothly
connected to each other for small value of n, while have only one branch when
the NUT parameter is large.

After exploring the basic properties of the new black holes, we turn to the
thermodynamic properties of the solutions, the entropy and temperature. Since
the presence of higher curvature Gauss-Bonnet term in the action, the entropy
is no longer equal to one quarter of the event horizon area $A_{h}$ anymore.
The correct entropy of black hole can be calculated through Wald entropy
formula[18, 19]

$S_{H}=2\pi\int_{r_{h}}\left(r^{2}+n^{2}\right)^{2}\epsilon_{ab}\epsilon_{cd}\frac{\partial
L}{\partial R_{abcd}}\,,$ (30)

where L is the Lagrangian and $\epsilon_{ab}$ is the volume form on the
2-dimensional cross section $r_{h}$ of the horizon. For the extended scalar-
tensor-Gauss-Bonnet theory we consider, the explicit expression of black hole
entropy is given by

$S_{H}=\pi(r_{h}^{2}+n^{2})+4\pi\lambda^{2}F\left(\varphi_{h}\right)\,.$ (31)

It is obvious that the entropy is not equal to one quarter of event horizon
area in general, and reduces back to one quarter of event horizon when the
Guass-Bonnet term is turned off by setting $\lambda=0$.

Before presenting result for entropy, we want to show the relations between
the event horizon area $A_{h}$ and the mass parameter $m$ of the new found
black holes in Fig.6. As one can see, the event horizon areas of the
scalarized black hole (dashed lines) are smaller than that of scalar-free
black hole(solid lines) for small value of NUT parameter $n$. As the NUT
parameter $n$ is increasing, the situation changes, and finally the areas of
scalarized black holes are large than that of scalar-free black holes. We
notice another novel phenomena that once the mass parameter $m$ is positive
$(m>0)$ at the bifurcation points where the scalarized black holes begin, the
areas of event horizon are almost the same for different NUT parameters at the
bifurcation points, see the horizontal line in Fig.6.

Figure 6: The area of the black hole horizon $A_{H}$ as functions of the mass
parameter $m$ for different values of $n/\lambda$, the value of area is almost
the same at the starting point for positive mass parameter.

At this stage, we present the results of entropy in Fig.7. We plot entropy $S$
as a function of mass parameter $m$ in the left panel of Fig.7. There is one
big difference between the area-mass relation $A_{h}-m$ in Fig.6. Here, the
entropy of scalarized black hole is always larger than that of scalar-free
black holes, the reason is that the Gauss-Bonnet term has an additional
contribution to entropy for scalarized black holes. The hairy black holes with
larger entropy implies that they are more thermodynamically stable and more
favorable than the scalar-free black holes. The right panel of Fig.7 shows the
entropy as a function of event horizon radius. From the right panel, we can
see two branches of new hairy black hole solutions smoothly connected, again.
At the bifurcation points, the similar phenomenon emerges, the entropy at
bifurcation points are nearly the same, see the horizontal line in Fig.7.. The
value of entropy at bifurcation point is about 4.32 for positive mass
parameter, and it is also the largest entropy for the scalarized black hole.
It is mentioned in the beginning of this section that every scalarized black
hole for fixed NUT parameter has a maximal value of mass parameter which is at
the bifurcation point. We conjecture that all the saclarized black holes with
the maximal mass parameter is greater than zero have a maximal entropy bound
of 4.32 which is located at the bifurcation point.

(a) $S/\lambda^{2}$ as a function of $m/\lambda$ (b) $S/\lambda^{2}$ as a
function of $r_{h}/\lambda$

Figure 7: The entropy $S$ of the black hole as functions of the mass parameter
$m$ is shown in the left figure, the entropy $S$ as functions of the horizon
radius $r_{h}$ is shown in the right figure. The entropy is nearly a constant
at the bifurcation point for positive mass parameter.

The temperature can be calculated through standard method, for our Taub-NUT
like metric, it is given by

$\displaystyle T=\frac{\sqrt{h^{\prime}(r_{h})f^{\prime}(r_{h})}}{4\pi}\,.$
(32)

The temperature is plotted as a function of mass parameter and event horizon
radius respectively in Fig.8. As is seen, the temperature of the hairy black
holes is always higher than that of scala-free black holes. In the right
panel, we can see that the temperature of one branch hairy black holes is
always higher than that of the other branch.

(a) $\lambda T$ as a function of $m/\lambda$ (b) $\lambda T$ as a function of
$r_{h}/\lambda$

Figure 8: The temperature $T$ of the black hole as functions of the mass
parameter $m$ and the horizon radius $r_{h}$ for different values of
$n/\lambda$.

Through the probe scalar analysis, we can find the bifurcation points of hairy
black holes and scalar-free black holes where the scalarized black holes
begin, connecting these bifurcation points gives the existence line. From all
the results we collect, the hairy black holes terminate at some points, these
points constitute the critical line. We plot these two line in Fig.9, between
these two lines is the ”habitat” of the new scalarized black hole solution.

Figure 9: The scalarized Taub-NUT blac holes start from the red line, the
value of mass parameter is positive for the dashed red line, while the solid
red line represents the part with negative mass parameter. The hairy Taub-NUT
black holes terminate at the blue line.

As we pointed out in the entropy part, the entropy and horizon area of black
hole in the bifurcation point are nealy constant for positive mass parameter,
which means $r_{h}^{2}+n^{2}$ remains constant. Since it happens in the
bifurcation points, we can trade $r_{h}$ with mass parameter $m$ using the
relation ${\frac{r_{h}^{2}-2mr_{h}-n^{2}}{r_{h}^{2}+n^{2}}}=0$ on the event
horizon, then we get the relation between $m$ and $n$

$\left(\frac{m}{\lambda}\right)^{2}+\left(\frac{n}{\lambda}\right)^{2}+\frac{m}{\lambda}\sqrt{\left(\frac{m}{\lambda}\right)^{2}+\left(\frac{n}{\lambda}\right)^{2}}\approx\frac{2.16}{\pi}\,.$
(33)

We plot the curve of this analytic relation together with the existence line
in the left of Fig.10. One can see that the two curves almost overlap for
positive mass, however, it deviates from each other when the mass parameter
becomes negative. We also plot the difference of NUT parameters of the the two
curves as a function of mass parameter in the right of Fig.10, the difference
is obvious for negative mass parameter, whilst it drops rapidly to zero as
mass parameter turns to positive value.

(a)  (b)

Figure 10: The left panel shows the bifurcation points where the scalarized
black hole begin and the curve for constant entroy, they lie on top of each
other for positive mass parameter. The right panel shows difference of NUT
parameter of existence line and constant entropy line, it is almost zero for
positive mass parameter.

##  5 conclusions

In this work we have studied the scalarization of Taub-NUT black holes in the
extended scalar-tensor-Gauss-Bonnet gravity theory, where the nontrivial
scalar is excited by the curvature of the spacetime through the Gauss-Bonnet
coupling term. We choose a coupling function such that the theory admits the
Ricci-flat Taub-NUT black hole as a solution. Guided by the probe scalar
analysis, we construct the hairy Taub-NUT black holes numerically. We observe
one novel phenomenon that there exists two branches of new hairy black hole
solutions which are smoothly connected to each other under certain parameter
range. As far as we know, there is no such phenomenon appears in the
scalarization of Schwarzschild black holes. The NUT parameter $n$ plays a
significant role in this phenomenon. As far as the NUT parameter is not zero,
the second branch of the hairy black hole emerges, whilst there is no trace of
the second branch as NUT parameter vanishes. When the NUT parameter is small,
the second branch hairy black hole is not obvious, and it becomes more obvious
as NUT parameter increases. However, when the NUT parameter is large enough,
the second branch solution disappears, again.

The obtained hairy Taub-NUT black holes bifurcate from the analytic Taub-NUT
black holes, but do not fall back to the scalar-free Taub-NUT black holes like
the scalarization of Schwarzschild. Instead, they terminate at some point for
all the cases we have studied. The bifurcation point is shifted to smaller
mass parameter as the NUT parameter increases.

We have also studied the entropy and temperature of the hairy black holes. The
temperature of scalarized black holes is always larger than scalar-free Taub-
NUT black holes. We found that the event horizon area of hairy black holes is
not always larger than that of scalar-free Taub-NUT black holes, but entropy
is, the reason is that the Gauss-Bonnet term has a additional contribution to
the entropy of hairy black holes. This implies that the new found black hole
solutions are thermodynamically more stable. We also observe a novel phenomena
in this system that the entropy of the black holes at the bifurcation point is
almost constant for positive mass parameter. We then conjecture a maximal
entropy bound for all the scalarized black hole whose mass parameter at the
bifurcation point is greater than zero. The physics under this phenomena is
worth exploring further.

The definition of mass and NUT charge of Taub-NUT black holes in Einstein’s
general relativity is still controversial, many progress has been made in this
direction[15, 20, 21, 22]. Turning to the higher curvature theories, such as
the theory we consider, the situation is even more complicated. The new
scalarized black holes supply one more platform to explore the definition of
mass and NUT charge in higher derivative gravity theories.

##  6 Acknowledgement

We are grateful to Yu-Qi Chen, Ze Li, Jun-Fei Liu and Hong Lu for useful
discussion. This work is supported in part by NSFC (National Natural Science
Foundation of China) Grant No. 12075166.

## References

  * [1] J. D. Bekenstein, Transcendence of the law of baryon-number conservation in black hole physics, Phys. Rev. Lett. 28, 452-455 (1972)
  * [2] J. D. Bekenstein, Novel “no-scalar-hair” theorem for black holes, Phys. Rev. D 51, no.12, R6608 (1995)
  * [3] D. D. Doneva and S. S. Yazadjiev, New Gauss-Bonnet Black Holes with Curvature-Induced Scalarization in Extended Scalar-Tensor Theories, Phys. Rev. Lett. 120, no.13, 131103 (2018) [arXiv:1711.01187 [gr-qc]].
  * [4] G. Antoniou, A. Bakopoulos and P. Kanti, Black-Hole Solutions with Scalar Hair in Einstein-Scalar-Gauss-Bonnet Theories, Phys. Rev. D 97, no.8, 084037 (2018) [arXiv:1711.07431 [hep-th]].
  * [5] H. O. Silva, J. Sakstein, L. Gualtieri, T. P. Sotiriou and E. Berti, Spontaneous scalarization of black holes and compact stars from a Gauss-Bonnet coupling, Phys. Rev. Lett. 120, no.13, 131104 (2018) [arXiv:1711.02080 [gr-qc]].
  * [6] M. Minamitsuji and T. Ikeda, Scalarized black holes in the presence of the coupling to Gauss-Bonnet gravity, Phys. Rev. D 99, no.4, 044017 (2019) [arXiv:1812.03551 [gr-qc]].
  * [7] P. Cunha, V.P., C. A. R. Herdeiro and E. Radu, Spontaneously Scalarized Kerr Black Holes in Extended Scalar-Tensor–Gauss-Bonnet Gravity, Phys. Rev. Lett. 123, no.1, 011101 (2019) [arXiv:1904.09997 [gr-qc]].
  * [8] Y. S. Myung and D. C. Zou, Instability of Reissner–Nordström black hole in Einstein-Maxwell-scalar theory, Eur. Phys. J. C 79, no.3, 273 (2019) [arXiv:1808.02609 [gr-qc]].
  * [9] T. Damour and G. Esposito-Farese, Nonperturbative strong field effects in tensor - scalar theories of gravitation, Phys. Rev. Lett. 70, 2220-2223 (1993)
  * [10] A. H. Taub, Empty space-times admitting a three parameter group of motions, Annals Math. 53, 472-490 (1951)
  * [11] E. Newman, L. Tamburino and T. Unti, Empty space generalization of the Schwarzschild metric, J. Math. Phys. 4, 915 (1963)
  * [12] R. A. Hennigar, D. Kubizňák and R. B. Mann,“Thermodynamics of Lorentzian Taub-NUT spacetimes,”Phys. Rev. D 100, no.6, 064055 (2019) doi:10.1103/PhysRevD.100.064055 [arXiv:1903.08668 [hep-th]].
  * [13] A. Ballon Bordo, F. Gray, R. A. Hennigar and D. Kubizňák,“The First Law for Rotating NUTs,”Phys. Lett. B 798, 134972 (2019) doi:10.1016/j.physletb.2019.134972 [arXiv:1905.06350 [hep-th]].
  * [14] A. Awad and S. Eissa, “Lorentzian Taub-NUT spacetimes: Misner string charges and the first law,” Phys. Rev. D 105, no.12, 124034 (2022) doi:10.1103/PhysRevD.105.124034 [arXiv:2206.09124 [hep-th]].
  * [15] H. S. Liu, H. Lu and L. Ma, Thermodynamics of Taub-NUT and Plebanski solutions, JHEP 10, 174 (2022) [arXiv:2208.05494 [gr-qc]].
  * [16] Y. Brihaye, C. Herdeiro and E. Radu, The scalarised Schwarzschild-NUT spacetime, Phys. Lett. B 788, 295-301 (2019) [arXiv:1810.09560 [gr-qc]].
  * [17] C. A. R. Herdeiro, E. Radu, N. Sanchis-Gual and J. A. Font, Spontaneous Scalarization of Charged Black Holes, Phys. Rev. Lett. 121, no.10, 101102 (2018) [arXiv:1806.05190 [gr-qc]].
  * [18] 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]].
  * [19] 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]].
  * [20] 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]].
  * [21] 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]].
  * [22] Z. Chen and J. Jiang, General Smarr relation and first law of a NUT dyonic black hole, Phys. Rev. D 100, no.10, 104016 (2019) [arXiv:1910.10107 [hep-th]].