# Faint Rapid Red Transients from
Neutron Star - CO White Dwarf Mergers

Yossef Zenati1, Alexey Bobrick2 and Hagai B. Perets1
1Physics Department, Technion - Israel Institute of Technology, Haifa 3200004,
Israel
2Lund University, Department of Astronomy and Theoretical physics, Box 43, SE
221-00 Lund, Sweden

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Mergers of neutron stars (NS) and white dwarfs (WD) may give rise to
observable explosive transient events. We use 3D hydrodynamical (SPH)
simulations, as well as 2D hydrodynamical-thermonuclear simulations (using the
FLASH AMR code) to model the disruption of CO-WDs by NSs, which produce faint
transient events. We post-process the simulations using a large nuclear
network and make use of the SuperNu radiation-transfer code to predict the
observational signatures and detailed properties of these transients. We
calculate the light-curves (LC) and spectra for five models of NS - CO-WD
mergers. The small yields of ${}^{56}{\rm Ni}$ (few$\times 10^{-3}\,{\rm
M_{\odot}}$) result in faint, rapidly-evolving reddened transients (RRTs) with
B (R) - peak magnitudes of at most $\sim-12$ ($-13$) to $\sim-13$ ($-15$),
much shorter and fainter than both regular and faint/peculiar type-Ia SNe.
These transients are likely to be accompanied by several months-long, $1$–$2$
mag dimmer red/IR afterglows. We show that the spectra of RRTs share some
similarities with rapidly-evolving transients such as SN2010x, although RRTs
are significantly fainter, especially in the I/R bands, and show far stronger
Si lines. We estimate that the upcoming Large Synoptic Survey Telescope could
detect RRTs at a rate of up to $\sim 10-70$ yr-1, through observations in the
R/I bands. The qualitative agreement between the SPH and FLASH approaches
supports the earlier hydrodynamical studies of these systems.

###### keywords:

stars: neutron – white dwarfs – supernovae: general

††pubyear: 2019††pagerange: Faint Rapid Red Transients from  Neutron Star - CO
White Dwarf Mergers–Appendix: Late-time fallback and fallback-powered
lightcurve

##  1 Introduction

# | ${M_{\rm WD}}[M_{\odot}]$ | ${\rho_{\rm max}[{\rm g}/{\rm cm}^{3}]}$ | ${R_{0}/r_{t}}$ | ${\rm\%He_{4}}$ | ${\rm\%C_{12}}$ | ${\rm\%O_{16}}$ | ${E_{\rm K}[{\rm erg}]}$ | ${\rm IGE_{U}[{\rm 10^{-3}M_{\odot}}]}$ | ${M_{\rm Tot-U}[{\rm 10^{-3}M_{\odot}}]}$
---|---|---|---|---|---|---|---|---|---
C | ${0.55}$ | ${8.5\times 10^{6}}$ | ${0.8}$ | ${-}$ | ${50}$ | ${50}$ | ${2.8\times 10^{48}}$ | ${9.0}$ | ${19}$
D | ${0.63}$ | ${6.4\times 10^{6}}$ | ${1}$ | ${9}$ | ${50}$ | ${41}$ | ${3.6\times 10^{49}}$ | ${9.6}$ | ${210}$
E | $0.63$ | $8.5\times 10^{6}$ | $0.8$ | $4$ | $49$ | $47$ | $1.4\times 10^{49}$ | $6.1$ | $90$
F * | $0.63$ | $2.3\times 10^{6}$ | $1.1$ | $-$ | $50$ | $50$ | $1.5\times 10^{48}$ | $3.3$ | $8.1$
J | $0.8$ | $4.4\times 10^{7}$ | $0.8$ | $-$ | $50$ | $50$ | $3.2\times 10^{49}$ | $45$ | $94$
Table 1:  The initial parameters of the simulated NS-WD merger models. The
columns show the names of the models, the masses of the WDs ($M_{\rm WD}$) in
solar masses, the maximum density in the torus in our FLASH simulations
($\rho_{\rm max}$), the torus radius in units of the tidal radius
($R_{0}/r_{t}$) and the initial mass fractions of ${\rm He}_{4}$, ${\rm
C}_{12}$ and ${\rm O}_{16}$. The last three columns show the kinetic energy at
the end of the FLASH simulations and the total mass of iron-group ejecta (IGE)
and all of the ejecta ($M_{\rm tot-U}$). The energy powering the outflows
comes almost entirely from accretion and not the nuclear processes involved,
e.g. Zenati et al. (2019).
*These models were run with ${M_{\rm NS}=2M_{\odot}}$

Figure 1: Snapshots for density distribution from the SPH simulation of model
D. The left column shows the density in the horizontal $x$-$y$ plane, and the
right column shows the density distribution in the vertical plane $x_{\rm
corot}$-$z$, corotating with the bulk of the white dwarf mass and centered on
the NS. The snapshots are taken at $34.4\,\textrm{s}$, $111.7\,\textrm{s}$ and
$303.6\,\textrm{s}$ from the beginning of the simulation. The simulation
starts with the donor relaxed into the state when the mass-transfer stream has
already developed. As the mass-transfer rate grows over time, the donor gets
increasingly deformed by the tidal field. At a certain point (here, at
$100\,{\rm s}$) the donor stops being gravitationally bound, gets stretched by
the tidal field and quickly relaxes into a debris-disk state (here, by about
$300\,{\rm s}$). The shock heating caused by self-crossings of the material
during this stage is responsible for most nuclear burning in the simulation.

Figure 2: The evolution of the WD debris for model E based on FLASH
simulations. The panels show the colour coded density distribution and
velocity fields (black arrows) throughout the simulation. The velocity scale
is $4500\,{\rm km}\cdot{\rm sec}^{-1}$. Magenta lines show the surfaces of
constant temperature.

Mergers of neutron stars (NSs) or black holes (BH) with white dwarfs (WDs)
have been explored through several hydrodynamical and nuclear-hydrodynamical
simulations, e.g. Fernández & Metzger (2013a); Bobrick et al. (2017); Zenati
et al. (2019); Fernández et al. (2019). All CO WD-NS mergers are expected to
lead to unstable mass transfer (Bobrick et al., 2017) where the WD is tidally
disrupted on dynamical timescales and forms an extended debris disk around the
NS (Papaloizou & Pringle, 1984; Fryer et al., 1999). The disc evolves mostly
viscously but is also affected by nuclear burning. Disc viscosity and nuclear
burning together power outflows throughout the disc (Metzger, 2012). Nuclear
burning in the disc proceeds steadily (with a possible weak detonation) and
produces small amounts of ${}^{56}{\rm Ni}$ (at most $10^{-2}$ ${\rm
M_{\odot}})$, e.g. Zenati et al. (2019). Such a small amount of ${}^{56}{\rm
Ni}$ is expected to lead to optical transients which are much fainter than
ordinary type Ia supernovae (SNe; which produce more than an order of
magnitude more ${}^{56}{\rm Ni}$). In this paper, we present the first
synthetic LCs and spectra for CO WD-NS mergers and show that they may give
rise to a novel class of faint, rapid red transients (RRTs), with no evidence
for hydrogen or helium.

To model the merger, we use two approaches. The first approach makes use of 2D
coupled hydrodynamical-thermonuclear simulations with the FLASH code, similar
to the method introduced by Fernández & Metzger (2013b). In these simulations,
the NS-WD merger is followed only from the advanced stage at which the debris
disk of the disrupted WD has already formed and given rise to an axisymmetric
structure which can be explored in 2D (Zenati et al., 2019, hereafter Paper
I). In the second approach, we applied a 3D SPH code to follow the full NS-WD
merger, including the early pre-disruption phase, but with nuclear burning
treated only in a post-processing step. As we discuss below, these two
approaches are complementary and generally produce very similar results.

As we and others have already found (Fernández & Metzger, 2013b; Zenati et
al., 2019), the evolution of the merger and the evolution of the WD debris
disk are dominated by the viscous evolution and the gravitational energy, with
nuclear burning playing only a minor role in the structural evolution of the
disk and the outflows, but playing an important role in driving the luminosity
evolution of the transient. These results motivated us to use 3D SPH
simulations which did not include any feedback from nuclear burning, which
were then post-processed using a 125-elements nuclear network to analyze
nuclear burning yields in such mergers. The nuclear yields and observational
predictions could then be compared with our results from the 2D FLASH
simulations.

To make detailed observational predictions for the transients arising from NS-
WD mergers, we have mapped the outputs from our hydrodynamic runs, post-
processed with a large nuclear network (in the SPH case), onto an input grid
for a 1D radiation transport code. We subsequently compared the obtained LCs
and spectra to LCs and spectra of various classes of peculiar transients. As
we show below, WD-NS mergers produce faint RRTs, fainter and redder than any
of the observed SNe thought to arise from thermonuclear explosions, and may,
therefore, give rise to a new class of potentially observable transients.
These transients may be accompanied by an even fainter red/infrared several
months-long afterglow.

Below we describe our numerical simulations and the various types of initial
conditions we explored in Section 2. We then present the main results in
Section 3 and discuss and summarize them in Section 4.

##  2 Methods

Our original FLASH simulations and their results are described in detail in
Zenati et al. (2019) (Paper I). Here we briefly summarize the main ingredients
of the model: the WD profiles produced with the MESA stellar evolution code,
the simulations of the debris disk evolution with the FLASH code, and the
nuclear post-processing - all of which we already described in Paper I. We
then focus on the new SPH simulations using the Water code and then describe
our radiation-transfer modelling with the SuperNu code and the resulting
observational signatures (light-curves and spectra), which are the main focus
of the current paper.

###  2.1 The neutron star and white dwarf models

We describe the properties of each of the NS-WD models we simulated in Table
1. We obtained the structural properties of the WDs through detailed stellar
evolution models of single and binary stars using the MESA code (Paxton et
al., 2011; Paxton et al., 2013, 2015). In all cases, we considered only
stellar progenitors of solar metallicity. Our models include both typical CO
WDs as well as hybrid HeCO WDs. The former are the outcomes of regular
evolution of single stars, which results in WDs composed of $\sim$50$\,\%$
carbon and $\sim$50$\,\%$ oxygen. The hybrid WDs, containing both CO and He,
are the outcomes of _binary_ stellar evolution, as described in Zenati et al.
(2018). Subsequently, we used these stellar profiles to set up our WD models
in the SPH code and their bulk compositions to set up the disrupted WD’s in
the FLASH code.

We do not resolve the NSs in our simulations and model them as point masses.
We considered two NS masses, $1.4\,M_{\odot}$ and $2\,M_{\odot}$, which
corresponds to the typical range of NS masses observed in binary radio
pulsars, e.g. Lattimer (2012).

###  2.2 FLASH simulations

During the early stages of disruption, the NS tidally shreds the WD into a
debris disc, e.g Papaloizou & Pringle (1984); Fryer et al. (1999). The
detailed initial conditions and the structure and evolution of the modelled
disk can be found in Paper I.

We simulate the evolution of such a disc using the publicly available FLASH
code v4.2 (Fryxell et al., 2000). We employ the unsplit ${\rm PPM}$ solver of
FLASH hydrodynamics code in ${\rm 2D}$ axisymmetric cylindrical coordinates on
a grid of size ${\rm 1\times 1\left[10^{10}cm\right]}$ using adaptive mesh
refinement. We follow the approach similar to the one adopted in other works
on thermonuclear SNe (e.g. Meakin et al. 2009). We handle detonations by the
reactive hydrodynamics solver in FLASH without the need for a front tracker,
which is possible since unresolved Chapman–Jouguet detonations retain the
correct jump conditions and propagation speeds.

In the cases of NS - HeCO-WD mergers (models D and E), we identify the
ignition and detonation of helium in the accretion disk. The helium ignition
occurs at sufficiently high accretion rates at densities of $\rho\gtrsim
8\times 10^{5}\,{\rm g}/{\rm cm}^{-3}$, which, through supersonic burning,
gives rise to a (weak) He-detonation in the HeCO mixed layer.

In the cases of NS - CO-WD mergers, the required pressure and density for a
weak detonation are $P\gtrsim 2\cdot 10^{25}\,{\rm erg}\cdot{\rm cm}^{-3}$ and
$\rho\gtrsim 2.5-3\times 10^{7}\,{\rm g}\cdot{\rm cm}^{-3}$, respectively,
which is higher than in the case of NS-He CO mergers. The timescale for ${\rm
C}_{12}$ burning is also much longer, in effect leading to a different and
delayed (by about one second) detonation compared to the HeCO WD cases.

In this work, we have also checked the resolution convergence following Fisher
et al. (2019) and have not observed any changes in the evolution of the disk.
We evolved our simulations for $(175$-$270)\,{\rm s}$. In one test case, we
ran a lower-resolution model for up to $400\,{\rm s}$. In comparison, the
orbital timescale at the circularisation radius, in case of model D, is about
$40$ seconds. The viscous timescale at the circularisation radius $R_{0}$,
given by:

$t_{\rm visc}\sim
260{\rm\,s}\left(\frac{0.1}{\alpha}\right)\left(\frac{R_{0}}{10^{9.3}\rm
cm}\right)^{3/2}\left(\frac{1.4M_{\sun}}{M_{\rm
c}}\right)^{1/2}\left(\frac{H_{0}}{0.5R_{0}}\right)^{-2},$ (1)

is about $260$ seconds, assuming the alpha-disc model parameter $\alpha=0.1$,
as adopted in our FLASH simulations. In the equation, $M_{\rm c}$ is the
central mass of the NS and $H_{0}$ is the scale height in the disc. In the
inner regions of the disc, where most nuclear burning takes place, due to a
smaller radius, the viscous timescale is only a few tens of seconds, which is
well covered by our runs.

While our FLASH simulations do not capture all the details of the initial
3-dimensional onset phase of mass transfer in the binaries, the orbital and
viscous timescales of the disc are well captured by our models. The simulation
timescales also agree with the timescales of our SPH simulations, which also
last for several hundreds of seconds. The overall similarity of the nuclear
abundance patterns and a reasonable, typically up to $1\,{\rm dex}$, agreement
of the nuclear yields supports the picture that both codes capture the phase
which is most relevant to nuclear burning.

As discussed in detail in Paper I, we use the Helmholtz EOS (Timmes & Swesty,
2000), and account for self-gravity of the disk. We employ a $19$-isotope
reaction network with correct treatment of the burning front (Fryxell et al.,
1989). To prevent the production of artificial unrealistic early detonation
that may arise from insufficient numerical resolution, we applied a limiter
approach following Kushnir et al. (2013). Further details of these simulations
can be found in Paper I. The initial conditions and the outcomes of the
simulations we use in the current paper are summarized in Table 1.

###  2.3 SPH simulations

We have complemented and verified our modelling with the FLASH code by using a
smoothed particle hydrodynamics (SPH) code Water (Bobrick et al., 2017). We
set up the SPH simulations and their post-processing aiming to reproduce the
same physical mergers like the ones modelled by the FLASH method. The SPH
approach allows us to evolve the models from early stages of mass transfer in
3D and follow the disruption of the WD and the formation of the debris disk,
while our 2D FLASH simulations can only be initiated at a late stage when a WD
debris disk is already assumed to have formed. On the other hand, our SPH
models only include nuclear burning as a post-processing step. The two
approaches, therefore, are different and complement each other.

The Water code is based on the most up-to-date SPH prescriptions available and
derives from the Oil-on-Water code used to model the onset stages of mass
transfer (Bobrick et al., 2017), see also Church et al. (2009). We refer to
Bobrick et al. (2017) for full implementation details and summarise the key
components of the code below.

We solve the Navier-Stokes equations by discretizing the Lagrangian as done in
Springel (2010). This way, we ensure exact energy and momentum conservation,
only limited by the accuracy of the gravity solver, which we implement as in
Benz et al. (1990). We limit the maximum radius of SPH particles in a
Lagrangian and continuous way, thus partially mitigating the so-called fall-
back problem when ejected particles become large and acquire too many
neighbours upon falling back on the donor. We use Wendland W6 kernels (Dehnen
& Aly, 2012), which prevent pairing instability by construction, and set the
number of neighbours to $400$. We use the artificial viscosity from Cullen &
Dehnen (2010), which switches on the viscosity only near shocks and does not
damp sound waves. We base the artificial conductivity on the prescription from
Bobrick et al. (2017); Hu et al. (2014), this way reducing discontinuities in
thermal energy near shocks. Finally, we use the KDK integrator as described in
Quinn et al. (1997); Springel (2005); Cullen & Dehnen (2010) with a time step
limiter by Saitoh & Makino (2009), which improves shock treatment.

We constructed the SPH models by using the stellar profiles prepared with the
MESA code, described in Section 2.1 (Table 1). We made the SPH models of the
WDs by mapping SPH particles onto appropriately spaced spherical shells (Saff
& Kuijlaars, 1997; Bobrick et al., 2017; Raskin & Owen, 2016). We used SPH
models with equal particle masses to avoid numerical artefacts reported in
earlier studies, e.g. Lorén-Aguilar et al. (2009); Dan et al. (2011). We
relaxed the single models in isolation over six dynamical times of the WD.
Then we placed the donors in a binary with the NS at $1.8\,a_{\rm RLOF}$ and
relaxed the binary in a co-rotating frame, continually spiralling in the donor
down to $0.975\,a_{\rm RLOF}$ over two orbital periods. Then we relaxed the
binary at the final separation over an additional half a period, to start the
main simulation with an existing mass-transfer stream, see Bobrick et al.
(2017). We removed the small number of particles accreted onto the NS during
the relaxation stage. Subsequently, we ran the simulations for up to $10$
orbital periods into the WD disruption until after the formation of a
developed disc. Since the SPH simulations allow us to model the merger from
the early stages to the advanced stages when a fully developed disc has
formed, we expect that they capture most of the phases relevant to
nucleosynthesis. Furthermore, as we show in Section 3.1, most of the nuclear
material, according to both FLASH and SPH, is produced in the advanced stages
of the merger when the disc has fully formed, which is captured by both codes.

In the SPH simulations, we used the same Helmholtz equation of state as we
used in FLASH, but did not include feedback from nuclear burning. However, as
mentioned above, the latter likely has only a little effect on the overall
dynamical evolution of the disk (Zenati et al., 2019). We softened the gravity
of the NS by a W6 Wendland kernel with a core radius of $0.1\,R_{\rm WD}$ and
have checked the convergence by simulating models in $20\,{\rm K}$, $50\,{\rm
K}$ and $100\,{\rm K}$ resolution.

###  2.4 Nuclear post-processing

We use a 19-isotope $\alpha$-chain network in the FLASH simulations and do not
evolve chemical compositions in the SPH simulations. The $\alpha$ network can
adequately capture the energy generated during nuclear burning (Timmes &
Swesty, 2000). To verify the abundances produced by the FLASH code and to
calculate the detailed nuclear yields for the SPH runs, we applied a
nucleosynthetic post-processing step with a large network. For the grid-based
FLASH simulations, we made use of $4000$-$10000$ tracer particles to follow
the evolution and to track the composition, velocity, density, and temperature
of the WD debris. These particles were evenly spaced throughout the WD-debris
disk with a step of ${2\times 10^{8}\,{\rm cm}}$. The detailed histories of
density and temperature from tracer particles (from FLASH) and SPH particles
(from the SPH Water code) were post-processed with MESA (version 10390) one-
zone burner (Paxton et al., 2015). We employed a 125-isotope network that
includes neutrons and composite reactions from JINA REACLIB (Cyburt et al.,
2010). We subsequently used the outputs of the original $\alpha$ network in
the FLASH case, to correctly weigh the nuclear abundances over the material,
and the outputs of the detailed MESA PPN network in the SPH case for further
analysis.

###  2.5 Radiation transfer modeling using SuperNu

Following the nuclear post-processing, we mapped the physical properties from
this step as inputs for the openly available radiation transfer code SuperNu
(Wollaeger et al., 2013; Wollaeger & van Rossum, 2014) in order to calculate
the light-curves (LC) and spectra expected from the mergers. SuperNu uses
Implicit Monte Carlo (IMC) and Discrete Diffusion Monte Carlo (DDMC) methods
to stochastically solve the special-relativistic radiative transport and
diffusion equations to first order in $v/c$ in three dimensions. The hybrid
IMC and DDMC scheme used in SuperNu makes it computationally efficient in
regions with high optical depth. This approach allows SuperNu to solve for
energy diffusion with very few approximations, which is very relevant for
supernova light curves.

SuperNu code follows the free-expansion phase of supernovae using a velocity
grid. We map the 3D velocities of the FLASH material and the SPH particles to
a 1D-spherical velocity grid of SuperNu, as they are, at the end of
simulations. In other words, we set the masses in each of the grid cells of
SuperNu equal to the sum of the masses of the corresponding elements (grid
cells for FLASH, particles for SPH). We set the chemical compositions and
electron fractions $Y_{e}$ as the mass-weighted average compositions in each
cell. Unlike the Eulerian-based FLASH code, SuperNu can handle true vacuum
within cells. Consequently, zero-mass grid cells containing no mapped
particles require no special treatment, although the velocity distributions
did not contain any significant gaps.

Since most of the material is gravitationally bound at the end of the
simulations, both in FLASH and SPH, mapping the 3D velocities directly to
SuperNu represents fast dynamical ejecta and is an approximation, as we
discuss in Section 4. To check for the effect of the additional ejecta
released as a slow wind, we additionally initialised SuperNu runs by randomly
assigning the element velocities typical values of $400\,{\rm km}/{\rm s}$
with a $1$-$\sigma$ spread of $200\,{\rm km}/{\rm s}$. As we discuss in
Section 4 the true fractions of material in fast and slow ejecta are
uncertain, and thus we obtain the characteristic signal expected from both
types of ejecta.

We initialised the SuperNu code at $0.5\,{\rm d}$ after the merger. The
homologously expanding material adiabatically cools down by at least six
orders of magnitude in temperature by this point. We set the temperature of
the material initialised in SuperNu to upper estimate value of
$10^{5}\,\textrm{K}$. Setting initial temperatures to values lower than
$10^{5}\,\textrm{K}$ led to essentially the same lightcurves or spectra. We
have also verified that the input distributions for SuperNu from SPH and FLASH
simulations broadly looked similar.

# | ${M_{\rm WD}[M_{\odot}}]$ | ${{}^{56}{\rm Ni}_{{\rm B}}\ [\times 10^{-3}{M_{\odot}}]}$ | ${\rm B_{peak}}$ | $\rm{\rm\varDelta B_{15}}$ | ${\rm R_{peak}}$ | ${\rm\varDelta R_{15}}$ | ${\rm I_{peak}}$ | ${\rm\varDelta I_{15}}$
---|---|---|---|---|---|---|---|---
C | ${0.55}$ | ${5.25}\ (6.51)$ | $-12.81\ (-12.68)$ | $4.50\ (2.85)$ | $-14.70\ (-14.66)$ | $3.29\ (1.81)$ | $-14.93\ (-14.89)$ | $2.34\ (1.78)$
D | ${0.62}$ | ${6.16\ (10.03)}$ | $-12.80\ (-12.66)$ | $3.41\ (2.85)$ | $-14.72\ (-15.11)$ | $2.75\ (1.54)$ | $-14.98\ (-15.26)$ | $2.01\ (1.57)$
E | $0.62$ | ${2.84\ (10.14)}$ | $-12.02\ (-12.92)$ | $3.48\ (1.93)$ | $-14.23\ (-14.90)$ | $2.52\ (1.44)$ | $-14.75\ (-15.01)$ | $2.10\ (1.40)$
F | $0.62$ | ${3.38\ (8.82)}$ | $-12.19\ (-12.40)$ | $3.64\ (2.47)$ | $-14.36\ (-14.79)$ | $2.77\ (1.51)$ | $-14.80\ (-14.91)$ | $2.02\ (1.47)$
J | $0.8$ | ${1.78\ (26.16)}$ | $-10.88\ (-12.28)$ | $3.35\ (2.67)$ | $-13.21\ (-14.35)$ | $1.91\ (1.52)$ | $-13.87\ (-14.85)$ | $1.82\ (1.42)$
Table 2:  The outcomes of the merger models explored in this study, based on
FLASH simulations and SPH results in close brackets. The columns correspond
to: model number (1); the white dwarf mass (2); the amount of ${}^{56}{\rm
Ni}$ produced (3); the B-peak luminosity (4); the decrease in B-magnitude over
15 days since the peak ($\Delta{\rm B}_{15}$) (5); the R-peak luminosity (6);
the decrease in R-magnitude over 15 days since the peak (7); the I-peak
luminosity (8) and the decrease in I-magnitude over 15 days since the peak
(9).

Figure 3:  The bolometric lightcurves. Top: The bolometric lightcurves for all
models based on the FLASH and SPH simulations. The subplot shows the light
curves for all the models based on the FLASH simulations. Bottom: Comparison
of the bolometric light curves from the SPH and the FLASH models for models D
and E.

Figure 4: The UBVRI light curve (LC) for model D (left) and model E (right).
The top panels show the results from the radiative transfer processing of the
SPH simulations. The bottom panels show for the same models from based 2D
FLASH simulations.

Figure 5: Spectra for models C, F and J from the FLASH simulations shown at
different times (top panel), spectra for model F superposed with spectra from
SN2010x and SN2005ek (middle panel) and a comparison of the FLASH- and SPH-
based spectra for model D (bottom panel).

Figure 6: Light curve (top panel), BVRI bands (middle panel) and spectra
(bottom panel) for the largest possible secondary transient for model D, in
which all the material is ejected as a slow wind $v=400\,{\rm km}/{\rm s}$,
and a $1$-$\sigma$ spread of $200\,{\rm km}/{\rm s}$. The upper panel: the
light curve zoom on the first 240 days, the inset figure shows the light curve
for up to 372 days. The middle panel: lightcurves in the B, V, R, I bands. U
band is not shown since its luminosity is below mag $-8$. The lower panel: the
spectra at $-1$, $50$, $100$, $150$ days after peak. The secondary red
transient is expected to have a wide peak with half-time of $\sim{50}\,{\rm
d}$, while its peak magnitude is expected to be no larger than in this
example.  Figure 7: Nuclear yields from the FLASH and the SPH simulations
(star and box symbols, correspondingly), with the bars showing the difference
between the SPH and the FLASH yields. The yields for 17 elements resulting
from nuclear post-processing of models D, E and F show an overall good
agreement: the abundances of most alpha-elements, including ${}^{56}{\rm Ni}$
are within $1\,{\rm dex}$ from the two codes, and the overall abundance
patterns are similar. Figure 8: Comparison of the R-band lightcurves of
SN2010x and SN2005ek with the modelled lightcurves based on the FLASH runs.
While sharing similar durations, the rapid red transients (RRTs) resulting
form NS-CO WD mergers show much fainter peak luminosities.

##  3 Results

###  3.1 Merger dynamics

We infer the three-dimensional picture of the merger and its long-term
behaviour from our SPH simulations (see Fig. 1). Initially, the mass-transfer
rates are low, and the evolution resembles that of stably transferring
binaries. However, as the Roche lobe gradually cuts off increasingly large
fractions of the donor, the mass transfer rate grows over time. At a certain
point, the mass transfer rate is so high, that the donor effectively becomes
gravitationally unbound and enters a free fall onto the accretor Margalit &
Metzger (2016). During this phase, the tidal forces stretch the white dwarf
(see figure 1), causing it to self-intersect, thus producing shocks and
circularising the material into a disc. The density and temperature rise in
the shocked regions, which leads to, initially mild, nuclear burning. Due to
the limitations of the 2D models, our 2D FLASH models (see Fig. 2) begin only
after the formation of the debris disk and subsequently follow its nuclear-
hydrodynamic evolution.

In the FLASH models, we see that as the density and temperature increase in
the inner regions, they attain critical conditions for thermonuclear burning,
which leads to a weak explosive detonation. We discuss this in more detail in
Paper I. However, only a small fraction of the accreting material participates
in thermonuclear burning. As a result, this produces only small amounts of
nuclear by-products. Both the FLASH and the post-processing of the SPH
modelling show that the mergers produce only up to $\sim
10^{-3}-10^{-2}\,{M_{\odot}}$ of intermediate or even iron group elements. In
particular, the amount of ${}^{56}{\rm Ni}$ is always at most
$10^{-2}\,{M_{\odot}}$. The bulk of the nuclear material is produced in the
developed disc stages, well represented both by SPH and FLASH. The total
released nuclear energy, as reported from the 19-isotope network in FLASH
simulation, reaches only up to $\sim{\rm few\,\times}10^{46}-10^{47}\,{\rm
erg}$. For comparison, the amount of the released gravitational energy,
calculated as the total increase in the kinetic and internal energies minus
the total generated nuclear energy, see, e.g. (Fernández & Metzger, 2013a;
Zenati et al., 2019; Fernández et al., 2019), is a few times $10^{47}$ to
$10^{49}\,{\rm erg}$. Thus, the released gravitational energy is $10$ to $100$
times larger than the released nuclear energy. Hence, gravitational dynamics
drives most of the outflows and heats most of the material in the systems.

Our models do not follow the long-term evolution of the NS-WD mergers and, in
particular, the evolution of the fallback material that could potentially
contribute to the transient luminosity. Nevertheless, we briefly discuss a
simplified analytical model for the possible fallback contribution in the
Appendix. We find that the fallback material is unlikely to contribute
significantly to the transient energetics and the observable lightcurves.

###  3.2 Light curve and spectra properties

Figures 3 and 4 show the resulting bolometric and UBVRI lightcurves from our
models, while Fig. 5 shows the predicted spectra from these models. Table 2
summarizes the basic properties of the transients. Figure 6 shows the
lightcurve and spectra for the possible secondary transient from the disk-wind
material.

As can be seen, our models predict that NS - CO-WD mergers give rise to vary
faint transients with absolute B magnitudes in the range between $-10.9$ and
$-12.8$. These transients are far brighter in the I/R bands, with peak I
magnitudes in the range between $-13.9$ and $-15$. These transients show
extremely fast decline in the B-band ($\Delta{\rm B}_{15}$ in the range
$3.35-4.5$) and rapid decline in the redder bands ($\Delta{\rm I}_{15}$ in the
range $1.8-2.34$). Such rapid red transients (RRTs), as predicted by our
models, are difficult to observe in current surveys, and have not been
observed yet to the best of our knowledge. Nevertheless, next-generation
surveys such as the upcoming Large Synoptic Sky Survey (LSST), might be able
to observe tens of such events every year (see the next section). An
observational feature specific to these transients may be a faint (up to mag
$-13$) ${}^{56}{\rm Ni}$-powered several-months long secondary red afterglow
transient. Due to lower ejection velocities compared to RRTs, the expansion
happens on roughly $10$ times longer timescales, resulting in larger diffusion
time, later nebular phase and several month-long peak half-times. The peak
luminosity of these secondary transients are an upper estimate and are likely
lower, depending on the fraction of the mass lost through the disc wind
ejecta.

Our results show comparable ranges of ejected masses, velocities and
${}^{56}{\rm Ni}$ production as found in the simplified 1D models by Margalit
& Metzger (2017a). These earlier models provided basic estimates for the
bolometric light curves which are broadly in agreement with our detailed
bolometric light curves.

###  3.3 Rate estimates

Despite the relatively small peak luminosities of such transients, upcoming
synoptic surveys such as LSST will observe them frequently. With the
$5$-$\sigma$ magnitude limits for single-visit detections of $24.0$ and $24.7$
in I and R bands, respectively (LSST Science Collaboration et al., 2009), LSST
will be able to detect such transients within $15$ days past the peak up to a
distance of $250\,\textrm{Mpc}$. In order to estimate the detection rates, we
assume that the merger rates of CO WD-NS binaries scale with the blue
luminosity (total luminosity, e.g. of galaxies, in the B-band). By using the
galactic merger rate of $1-2\cdot 10^{-4}\,\textrm{yr}^{-1}$ from our
population synthesis studies of NS-WD mergers (Toonen et al., 2018), the
galactic blue luminosity of $1.7\cdot 10^{10}L_{{\rm B},\odot}$ (Kalogera et
al., 2001) and the blue luminosity for the local universe $1.98\cdot
10^{8}\,{L_{{\rm B},\odot}}/{{\rm Mpc}^{3}}$ (Kopparapu et al., 2008), we find
that LSST will detect between $10$ and $70$ transients per year within its
field of view. The detection efficiency is similar for R and I bands. In
comparison, the ongoing ZTF survey is expected to observe between $0.05$ and
$0.3$ such transients per year, assuming R-band detection limit of $21$,
consistent with the current null-detection of such transients by ZTF.

###  3.4 Composition

The detailed composition yields from our FLASH models have already been
presented in detail in Paper I. The results from the post-processing of the
SPH simulations show overall agreement with the FLASH results. The nuclear
yields for most alpha-elements, as well as ${}^{56}{\rm Ni}$, are within
$1\,{\rm dex}$, as produced by the two codes, and the abundance patterns look
similar, although suggesting systematically somewhat lesser production of
burned material in FLASH. The detailed comparison for 17 elements can be seen
in Fig. 7.

##  4 Discussion and summary

In our study, we explored the observable manifestation of transients from NS-
WD mergers. We have made use of AMR and SPH hydrodynamical simulations to
model the mergers of NSs with CO WDs, followed their thermonuclear evolution
through concurrent and post-process analysis and then used radiative-transfer
models to predict the expected lightcurves and spectra of such transients. We
generally find that NS-WD mergers give rise to ultra-faint, rapidly evolving
reddened transients (RRTs), likely observable with next-generation surveys
such as the LSST.

In our study, we have qualitatively confirmed, for the first time with a 3D
code, the nuclear yields in NS-WD mergers obtained in earlier axisymmetric
simulations, e.g. (Fernández & Metzger, 2013a; Zenati et al., 2019; Fernández
et al., 2019). As discussed in Section 2.5, the fraction of the fast material
ejected dynamically, is the major uncertainty in our model, and our
synthesised lightcurves provide the maximal expected signal expected for RRT.
The FLASH simulations, dominated by $\alpha$-viscosity, may be representative
of the MRI-dominated disc stage following the merger and produce up to
$0.1$-$0.2\,M_{\odot}$ of fast ejecta, thus supporting the idea that the
amount of fast ejecta may be significant. The slow disc ejecta produced later
on may result in a faint red several months-long afterglow. The subsequent
evolution of the remnant object has been studied in detail by Margalit &
Metzger (2017b).

Our predicted RRTs are much fainter and faster evolving than both typical Ia
SNe as well as other peculiar faint and/or rapid transients observed to date.
Over the last two decades several classes of fast-evolving SNe had been
identified (de Vaucouleurs & Corwin, 1985; Poznanski et al., 2010; Tutukov et
al., 1992; Kasliwal et al., 2010; Perets et al., 2011; Drout et al., 2013,
2014; Ruiter, 2020). However, RRTs are too faint and typically redder than
these SNe. Even the faintest rapidly evolving SNe differ significantly from
our predicted RRTs, as can be seen in the lightcurve in figure 8 comparing SN
2010X and SN 2005ek with the RRTs. Moreover, their helium and aluminium
content may be too small to explain the He and Al lines identified in these
SNe (Kasliwal et al., 2010), and they show much stronger Si lines. Other faint
type Ia SNe such as Ca-rich 2005E-like SNe (Perets et al., 2010) or SN 2008ha-
like SNe (Foley et al., 2009) are still much brighter in the blue bands, are
not as red and decline much slower than our predicted RRTs.

Transients from NS-WD mergers may, therefore, represent a completely different
class of SNe that might be observed mostly in close-by galaxies using large
telescopes or possibly with next-generation surveys such as LSST, as discussed
above. As we show in Toonen et al. (2018) and discuss above, the rates of NS-
WD mergers could be sufficiently high as to be observed by LSST. The delay-
time distribution for such mergers peaks at early times of hundreds of Myr up
to one Gyr, and they are therefore expected to occur preferentially in late-
type host galaxies.

Finally, as we noted in Paper I, the total masses accreted onto the NS in our
models appear to be small, and we, therefore, expect that NS-WD mergers are
unlikely to produce regular GRBs. Nevertheless, one cannot exclude ultra-faint
long GRBs lasting over timescales comparable to accretion timescales; in this
case, an ultra-long GRB accompanied by an RRT would provide a very clear
smoking-gun signature for the origin of both these types of transients.

## Acknowledgements

We thank Dr. Daan Van Rossum, Prof. Silvia Toonen and Prof. Robert Fisher and
Prof. Ari Laor for stimulating discussions. The simulations with the Water
code and their post-processing were performed on the resources provided by the
Swedish National Infrastructure for Computing (SNIC) at the Lunarc cluster.

## References

  * Arnett (1979) Arnett W. D., 1979, ApJ, 230, L37
  * Benz et al. (1990) Benz W., Cameron A. G. W., Press W. H., Bowers R. L., 1990, ApJ, 348, 647
  * Bobrick et al. (2017) Bobrick A., Davies M. B., Church R. P., 2017, MNRAS, 467, 3556
  * Church et al. (2009) Church R. P., Dischler J., Davies M. B., Tout C. A., Adams T., Beer M. E., 2009, MNRAS, 395, 1127
  * Cullen & Dehnen (2010) Cullen L., Dehnen W., 2010, MNRAS, 408, 669
  * Cyburt et al. (2010) Cyburt R. H., et al., 2010, ApJS, 189, 240
  * Dan et al. (2011) Dan M., Rosswog S., Guillochon J., Ramirez-Ruiz E., 2011, ApJ, 737, 89
  * Dehnen & Aly (2012) Dehnen W., Aly H., 2012, MNRAS, 425, 1068
  * Dexter & Kasen (2013) Dexter J., Kasen D., 2013, The Astrophysical Journal, 772, 30
  * Drout et al. (2013) Drout M. R., et al., 2013, ApJ, 774, 58
  * Drout et al. (2014) Drout M. R., et al., 2014, ApJ, 794, 23
  * Fernández & Metzger (2013a) Fernández R., Metzger B. D., 2013a, ApJ, 763, 108
  * Fernández & Metzger (2013b) Fernández R., Metzger B. D., 2013b, ApJ, 763, 108
  * Fernández et al. (2019) Fernández R., Margalit B., Metzger B. D., 2019, MNRAS, p. 1661
  * Fisher et al. (2019) Fisher R., Mozumdar P., Casabona G., 2019, ApJ, 876, 64
  * Foley et al. (2009) Foley R. J., et al., 2009, The Astronomical Journal, 138, 376
  * Fryer et al. (1999) Fryer C. L., Woosley S. E., Herant M., Davies M. B., 1999, ApJ, 520, 650
  * Fryxell et al. (1989) Fryxell B. A., Arnett W. D., Müller E., 1989, in Bulletin of the American Astronomical Society. p. 1209
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Hu et al. (2014) Hu C.-Y., Naab T., Walch S., Moster B. P., Oser L., 2014, MNRAS, 443, 1173
  * Kalogera et al. (2001) Kalogera V., Narayan R., Spergel D. N., Taylor J. H., 2001, ApJ, 556, 340
  * Kasliwal et al. (2010) Kasliwal M. M., et al., 2010, ApJ, 723, L98
  * Kopparapu et al. (2008) Kopparapu R. K., Hanna C., Kalogera V., O’Shaughnessy R., González G., Brady P. R., Fairhurst S., 2008, ApJ, 675, 1459
  * Kushnir et al. (2013) Kushnir D., Katz B., Dong S., Livne E., Fernández R., 2013, ApJ, 778, L37
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201
  * Lattimer (2012) Lattimer J. M., 2012, Annual Review of Nuclear and Particle Science, 62, 485
  * Lorén-Aguilar et al. (2009) Lorén-Aguilar P., Isern J., García-Berro E., 2009, A&A, 500, 1193
  * Margalit & Metzger (2016) Margalit B., Metzger B. D., 2016, MNRAS, 461, 1154
  * Margalit & Metzger (2017a) Margalit B., Metzger B. D., 2017a, MNRAS, 465, 2790
  * Margalit & Metzger (2017b) Margalit B., Metzger B. D., 2017b, MNRAS, 465, 2790
  * Meakin et al. (2009) Meakin C. A., Seitenzahl I., Townsley D., Jordan IV G. C., Truran J., Lamb D., 2009, ApJ, 693, 1188
  * Metzger (2012) Metzger B. D., 2012, MNRAS, 419, 827
  * Papaloizou & Pringle (1984) Papaloizou J. C. B., Pringle J. E., 1984, Monthly Notices of the Royal Astronomical Society, 208, 721
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Perets et al. (2010) Perets H. B., et al., 2010, Nature, 465, 322
  * Perets et al. (2011) Perets H. B., Badenes C., Arcavi I., Simon J. D., Gal-yam A., 2011, ApJ, 730, 89
  * Poznanski et al. (2010) Poznanski D., et al., 2010, Science, 327, 58
  * Quinn et al. (1997) Quinn T., Katz N., Stadel J., Lake G., 1997, ArXiv Astrophysics e-prints,
  * Raskin & Owen (2016) Raskin C., Owen J. M., 2016, The Astrophysical Journal, 820, 102
  * Ruiter (2020) Ruiter A. J., 2020, arXiv e-prints, p. arXiv:2001.02947
  * Saff & Kuijlaars (1997) Saff E. B., Kuijlaars A. B., 1997, The Mathematical Intelligencer, 19, 5
  * Saitoh & Makino (2009) Saitoh T. R., Makino J., 2009, ApJ, 697, L99
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel (2010) Springel V., 2010, ARA&A, 48, 391
  * Timmes & Swesty (2000) Timmes F. X., Swesty F. D., 2000, ApJS, 126, 501
  * Toonen et al. (2018) Toonen S., Perets H. B., Igoshev A. P., Michaely E., Zenati Y., 2018, ArXiv:1804.01538,
  * Tutukov et al. (1992) Tutukov A. V., Yungelson L. R., Iben Icko J., 1992, ApJ, 386, 197
  * Wollaeger & van Rossum (2014) Wollaeger R. T., van Rossum D. R., 2014, ApJS, 214, 28
  * Wollaeger et al. (2013) Wollaeger R. T., van Rossum D. R., Graziani C., Couch S. M., Jordan IV G. C., Lamb D. Q., Moses G. A., 2013, ApJS, 209, 36
  * Zenati et al. (2018) Zenati Y., Toonen S., Perets H. B., 2018, Monthly Notices of the Royal Astronomical Society, 482, 1135
  * Zenati et al. (2019) Zenati Y., Perets H. B., Toonen S., 2019, MNRAS, 486, 1805
  * de Vaucouleurs & Corwin (1985) de Vaucouleurs G., Corwin Jr. H. G., 1985, ApJ, 295, 287

## Appendix: Late-time fallback and fallback-powered lightcurve

A significant fraction of the mass of the disrupted WD remains bound during
the entire merger. The fallback material may, therefore, accrete onto the NS
and potentially contribute to the emission from the merger over long
timescales after the initial disruption of the WD. However, our radiative
transfer models assume a homologous expansion of the material and therefore do
not account for any contribution from the fallback material.

To test how sensitive are our predictions to the way we initialize the ejecta
properties in SuperNu, we, therefore, considered a simplified analytical model
following Dexter & Kasen (2013). It describes how the longer-term fallback of
material onto the accretor powers the evolution of the light-curve, and
predicts

$L_{\rm
fallback}\left(t\right)=L_{0}\left(\frac{t_{0}}{2t_{d}}\right)^{n}e^{-t^{2}/2t_{d}^{2}}\left(-\frac{1}{2}\right)^{n/2}\times$
(2)
$\times\left[\gamma\left(1-\frac{n}{2},-\frac{t_{0}^{2}}{2t_{d}^{2}}\right)-\gamma\left(1-\frac{n}{2},-\frac{t^{2}}{2t_{d}^{2}}\right)\right],$

where $\gamma\left(s,x\right)$ is the incomplete gamma-function, $t_{0}=t_{\rm
dyn}$, $n={5/3}$, and $t_{d}$ follows from Arnett (1979):

${t_{d}=\sqrt{\frac{3}{4\pi}\frac{\left(M_{\rm ejecta}+M_{\rm
fallback}\right)\kappa}{v_{\rm final-ejecta}\ \cdot c}}}$ (3)

where $\rm c$ is speed of light; for the complete solution see Dexter & Kasen
(2013). Using this model we find that the potential contribution to the light-
curves from the fallback model are $3$ – $7$ orders of magnitude smaller in
luminosity compared to the ones predicted by our SuperNu modelling, or by our
$E_{\rm SN}$ model in Zenati et al. (2019). Given these results, we believe it
is safe to neglect the contribution from the emission from the fallback
accretion.

In the longer term, the remnants form a stable thick disk-like puffy structure
around the NS. In the yet longer-term, the leftover debris produce a more
spherical cloudy structure around the NS. Its detailed evolution is beyond the
scope of our models and is to be studied elsewhere.


Yan Liu111Email<EMAIL_ADDRESS>Hong-Da Lyu222Email:
<EMAIL_ADDRESS>and Avinash Raju333Email<EMAIL_ADDRESS>Center for Gravitational Physics, Department of Space Science,
and International Research Institute of Multidisciplinary Science,
Beihang University, Beijing 100191, China

###### Contents

  1. 1 Introduction
  2. 2 Set-up
    1. 2.1 Four dimensional Einstein-scalar theory
    2. 2.2 Thermodynamics of black holes
  3. 3 Behavior of the singularity across the phase transition
    1. 3.1 Second order phase transition
    2. 3.2 First order phase transition
    3. 3.3 Probes of the black hole singularity via geodesics
  4. 4 Conclusion and discussion
  5. A The boundary terms for the scalar field in holographic renormalization

##  1 Introduction

Black holes are robust predictions of Einstein’s gravitational theory, which
have attracted immense attention from both the theoretical community and
experimental observations. However, the physics of black hole interior is
still mysterious. It is believed that the appearance of singularities behind
the black hole horizon are unavoidable [1]. The description of the
mathematical structure of black hole singularities has been studied by
Belinski, Khalatnikov and Lifshitz (BKL) [2].

Among the families of spacelike BKL singularities, the most well-known
singularities are the ones in Schwarzschild black hole and FLRW cosmology.
More generally, in the simplest homogeneous cases, the BKL singularities are
known as Kasner singularities [3, 4] which take the form characterized by
Kasner exponents. Recently it was also shown that the geometry near the
singularity inside the black hole solutions of a class of Einstein-scalar
theories in AdS is of Kasner form [5].444Other work on the Kasner
singularities in gravitational systems can be found in e.g. [6, 7, 8, 9, 10,
11, 12]. Previous studies on singularities in asymptotic flat spacetime in
Einstein-scalar gravity can be found in [13].  More explicitly, in presence of
scalar field, close to the Kasner singularity, the $d+1$ dimensional geometry
and the scalar field are of the form

$\displaystyle
ds^{2}\sim-d\tau^{2}+\tau^{2p_{t}}dt^{2}+\tau^{2p_{i}}dx^{i}dx^{i}\,,~{}~{}~{}\phi\sim-\sqrt{2}p_{\phi}\log\tau$
(1.1)

where $\tau$ is a function of radial coordinate and the Kasner exponents
$p_{t},p_{i},p_{\phi}$ satisfy
$p_{t}+\sum_{i=1}^{d-1}p_{i}=p_{t}^{2}+\sum_{i=1}^{d-1}p_{i}^{2}+p_{\phi}^{2}=1$.
In general the Kasner metric has curvature singularity at $\tau=0$ except the
case where one of the Kasner exponents $p_{t},p_{i}$ is $1$ while the other
Kasner exponents vanishes.

Obviously the BKL picture is based on the classical gravitational theory which
would break down close to the singularity. It is expected that the singularity
will be resolved in the full quantum theory of gravity. Nevertheless, we
believe that analyzing the singularity within the framework of classical
gravity is still important and might shed light on the singularity resolution
mechanism. Here we shall discuss the singularity using classical gravity and
focus on the black hole systems with the singularities of Kasner form.555For
the generic case, for example we could consider inhomogenities, when we
approach the singularity, the Kasner exponent oscillates [4, 14] and we will
not consider this case.

The goal of this paper is to study the behavior of the black hole
singularities during the black hole phase transitions. It is well-known that
black hole phase transitions could be studied within framework of black hole
thermodynamics from the approach of Euclidean gravity which is related to the
physics outside of horizon. Therefore our study is expected to show a further
link between the physics inside the black hole horizon and that of the
outside. We shall study a four dimensional Einstein-scalar theories in AdS
with double trace deformations. By tuning the ratio between temperature and
the double trace deformation parameter (which is the only dimensionless
tunable parameter), the system undergoes different orders of phase transitions
by choosing different scalar potentials. Since the black holes solutions have
spacelike Kasner singularities inside the horizon, we shall study the
behaviors of the Kasner exponents during the phase transitions.

Although the black hole singularities are located behind the horizon of the
black holes, their information could can be extracted from the physical
quantities of the dual field theory in the context of AdS/CFT correspondence,
e.g. correlation functions [15, 16, 17], entanglement entropies [18] etc. The
fact that the information about singularity can be extracted from the dual
boundary field theory is consistent with the black hole complementarity [19].
We will discuss how to probe the Kasner exponents of the black holes from
spacelike and timelike geodesics.

Our paper is organized as follows. In Sec. 2, we set-up the holographic model
and present all the necessary ingredients for the numerical calculations to
study the phase transitions and black hole singularities. In Sec. 3 we discuss
the numerical results on the behaviors of singularities across the phase
transitions. We conclude and discuss some open questions in Sec. 4.

##  2 Set-up

We begin by collecting all the necessary ingredients for constructing black
hole solutions and studying the phase transitions in a four dimensional
Einstein scalar theory in AdS. We will first present the four dimensional
Einstein-scalar theory under consideration and show the numerical strategy to
solve the system from the boundary to the singularity. Then we will study the
thermodynamics of the black hole solutions in order to study the phase
diagrams of the system from which the singularity behavior could be uncovered.

###  2.1 Four dimensional Einstein-scalar theory

We consider a 3+1 dimensional theory of real scalar field coupled to
Einstein’s gravity666We have set the AdS radius $L=1$.

$\displaystyle S=\frac{1}{16\pi G}\int
d^{4}x\;\sqrt{-g}\,\Big{[}R+6-(\nabla\phi)^{2}-V(\phi)\Big{]}$ (2.1)

with a potential $V$ for the scalar field given by

$\displaystyle
V(\phi)=m^{2}\phi^{2}-\lambda_{3}\,\phi^{3}+\lambda_{4}\,\phi^{4}\,.$ (2.2)

To have a positive $V$ when $\phi\to\pm\infty$, we should have
$\lambda_{4}>0$. Note that when $\lambda_{3}\neq 0$, the $Z_{2}$ symmetry of
the system $\phi\to-\phi$ is broken, and the potential $V(\phi)$ has a W-shape
with two different local minima values located at
$\phi_{1,2}=\frac{3\lambda_{3}\pm\sqrt{9\lambda_{3}^{2}-32m^{2}\lambda_{4}}}{8\lambda_{4}}$.
We shall focus on the case $\lambda_{3}\geq 0$. We will be interested in a
class of solutions that asymptotically approach AdS4 at the boundary.
According to the AdS/CFT dictionary, the boundary CFT has a dual operator
whose scaling dimension is determined by the mass of the bulk scalar field.

We are interested in the finite temperature phases of the CFT, which are dual
to black holes in the bulk. We follow [5] to make following ansatz for the
asymptotically AdS4 metric and scalar field

$\displaystyle
ds^{2}=\frac{1}{r^{2}}\left(-f(r)e^{-\chi(r)}dt^{2}+\frac{dr^{2}}{f(r)}+dx^{2}+dy^{2}\right)\,,~{}~{}~{}\quad\phi=\phi(r)\,$
(2.3)

where $f$, $\chi$ are only functions of $r$. The AdS4 boundary is located at
$r\to 0$ while the singularity is located at $r\to\infty$. It is also useful
to work in the Eddington-Finkelstein coordinate where the above metric has the
following form

$\displaystyle
ds^{2}=\frac{1}{r^{2}}\left(-f(r)e^{-\chi(r)}dv^{2}+2e^{-\chi(r)/2}dvdr+dx^{2}+dy^{2}\right)\,,~{}~{}~{}\quad\phi=\phi(r)\,$
(2.4)

where $v$ is the infalling Eddington-Finkelstein coordinate. Obviously the
planar Schwarzschild black hole is a solution of the form (2.4) with
$f=1-\frac{r^{3}}{r_{h}^{3}},~{}\chi=\phi=0$. We will be interested in the
black hole solution with nontrivial scalar hair.

We focus on the case with $m^{2}=-2$. The equations of motion are given by

$\displaystyle\begin{split}\chi^{\prime}-r\phi^{\prime 2}&=0\,,\\\
f^{\prime}-\Big{(}\frac{r\phi^{\prime
2}}{2}+\frac{3}{r}\Big{)}\,f+\frac{3}{r}-\frac{V}{2r}&=0\,,\\\
\phi^{\prime\prime}+\Big{(}\frac{f^{\prime}}{f}-\frac{2}{r}-\frac{\chi^{\prime}}{2}\Big{)}\phi^{\prime}-\frac{\partial_{\phi}V}{2r^{2}f}&=0\,,\end{split}$
(2.5)

where the prime is the derivative with respect to the radial coordinate $r$.

From (2.5) we know that near the asymptotic AdS4 boundary the fields behave as
a power series in $r$ as

$\displaystyle\begin{split}f&=1+\frac{\alpha^{2}}{2}r^{2}+\left(m_{T}-2\alpha^{3}\lambda_{3}\log
r\right)r^{3}+\cdots\\\
\chi&=\chi_{0}+\frac{\alpha^{2}}{2}r^{2}+\frac{1}{3}\left(4\alpha\beta-\alpha^{3}\lambda_{3}-6\alpha^{3}\lambda_{3}\log
r\right)r^{3}+\cdots\\\ \phi&=\alpha
r+\left(\beta-\frac{3}{2}\alpha^{2}\lambda_{3}\log
r\right)r^{2}+\frac{1}{8}\left(-12\alpha\beta\lambda_{3}+\alpha^{3}(2-27\lambda_{3}^{2}+8\lambda_{4})+18\alpha^{3}\lambda_{3}^{2}\,\log
r\right)r^{3}+\cdots\end{split}$ (2.6)

The scalar field $\phi$ in the bulk is dual to a scalar operator
$\mathcal{O}$. We can interpret either the value $\alpha$ or $\beta$ as the
VEV of $\mathcal{O}$ for the current choice of mass parameter here, which are
known as alternative quantization and standard quantization respectively [20].

Here we use the alternative quantization for the scalar field, i.e. we set
$\alpha$ as the expectation value of the dual operator. Note that the
conformal dimension of the operator is one. In this case, we can consider a
relevant deformation of the boundary field theory by a double-trace operator
[21, 22]

$\displaystyle S\rightarrow S-\kappa\int d^{3}x\;\mathcal{O}^{2}\,.$ (2.7)

When $\kappa<0$, this deformation makes the dual system easier to condensate
and might induce a phase transition at finite temperature [23]. Therefore we
shall focus on the parameter regime with $\kappa<0$. With double trace
deformations, now the source of the operator in the boundary is
$\kappa\alpha-\beta$ while $\alpha$ is the expectation value of the operator
[21, 24].

Holography with double trace deformation has been widely studied in the
literature to explore the physics of phase transitions in the context of
AdS/CMT, see e.g. [23, 25, 26, 27]. In case of the symmetric potential in
(2.1), i.e. $\lambda_{3}=0$ and $\kappa<0$, there is a critical temperature
$T_{c}$ below which there is a nonzero value for $\langle\mathcal{O}\rangle$
[23, 27]. In the following we shall consider the general cases of the
potential and study the behavior of the phase transitions in the dual system
as well as the behaviors of the singularities during the phase transition in
the bulk.

At finite temperature, near horizon the fields can be expanded as

$\displaystyle\begin{split}f&=-\frac{V(\phi_{h})-6}{2\,r_{h}}(r_{h}-r)+\cdots\\\
\chi&=\chi_{h}+\frac{\phi_{h}^{2}(4+3\lambda_{3}\phi_{h}-4\lambda_{4}\phi_{h}^{2})^{2}}{r_{h}(V(\phi_{h})-6)}(r_{h}-r)+\cdots\\\
\phi&=\phi_{h}+\frac{\phi_{h}(4+3\lambda_{3}\phi_{h}-4\lambda_{4}\phi_{h}^{2})^{2}}{r_{h}(V(\phi_{h})-6)}(r_{h}-r)+\cdots\end{split}$
(2.8)

where
$V(\phi_{h})=-2\phi_{h}^{2}-\lambda_{3}\phi_{h}^{3}+\lambda_{4}\phi_{h}^{4}$.
The zero temperature near horizon condition will be discussed in the next
section.

There are three free parameters $\chi_{h},r_{h}$ and $\phi_{h}$ at the
horizon. The Hawking temperature of the black hole is determined by777Note
that to get the temperature, we should set $\chi_{0}$ in (2.6) to be zero.

$\displaystyle
T=\frac{|f^{\prime}(r_{h})|e^{-\chi_{h}/2}}{4\pi}=\frac{(6-V(\phi_{h}))e^{-\chi_{h}/2}}{8\pi\,r_{h}}\,,$
(2.9)

from which we have the relation between $\phi_{h}$ and temperature $T$

$\displaystyle 8\pi r_{h}Te^{\chi_{h}/2}+V(\phi_{h})-6=0\,.$ (2.10)

Given the fact that the temperature should be positive, the above relation
puts a constraint on the allowed $\phi_{h}$. For the geometry to flow to
asymptotic AdS4 with appropriate AdS radius, $\phi_{h}$ should satisfy
$\phi_{1}\leq\phi_{h}\leq\phi_{2}$, where $\phi_{1,2}$ are two minima of the
potential for the scalar field.

The following two scaling symmetries are useful for solving the system. The
first one is

$\displaystyle v\rightarrow bv\,,~{}~{}\chi\rightarrow\chi+2\log b\,,$ (2.11)

and the second one is

$\displaystyle(v,r,x,y)\rightarrow b(v,r,x,y)\,.$ (2.12)

We can set $\chi_{h}=0$ using (2.11) and set $r_{h}=1$ using (2.12), then near
the horizon we only have one free parameter $\phi_{h}$. We can take the
$\phi_{h}$ as a shooting parameter and integrate the equations from horizon to
AdS4 boundary numerically to obtain the black hole solution and therefore we
can get $\\{m_{T},\chi_{0},\alpha,\beta\\}$ near AdS4 boundary. From the
numerical solution we obtained, $\chi_{0}$ can be scaled to be zero using the
scaling symmetry (2.11) and then we have a nonzero $\chi_{h}$ near the
horizon. By imposing that the dual theory is sourceless with the double trace
deformation, i.e. $\kappa\alpha-\beta=0$, we obtain a nonzero $\kappa$.
Therefore, for a given $\phi_{h}$, we have the only dimensionless parameter
$T/(-\kappa)$ in dual field theory. Note that when $\lambda_{3}\neq 0$, the
logarithmic term in (2.6) indicates a scaling anomaly for the double trace
deformation parameter.888A similar scaling anomaly of the marginal double
trace deformation parameter was found in a different context [28, 29]. One
difference is that here the double trace deformation parameter is relevant.One
could introduce an explicit UV scale by replacing $\log r$ in (2.6) as
$\log(\Lambda r)$. With a specific scale $\Lambda$ which we set to be 1, and
we also fix the boundary value $\kappa=-1$ in this case. Thus for
$\lambda_{3}\neq 0$ the theories under consideration are defined at the same
renormalization scale with $\kappa=-1$.

It is necessary to know the interior region of the black hole solution to
study the behavior of Kasner exponents during black hole phase transition. So
we need to do the numerical integration from horizon to a large enough $r_{s}$
that can be thought as the singularity. From the equations of the system
(2.5), when $r\rightarrow\infty$, the solutions behave as

$\displaystyle\phi=\sqrt{2}c\log r+\cdots,~{}~{}~{}\chi=2c^{2}\log
r+\chi_{1}+\cdots,~{}~{}~{}f=-f_{1}r^{3+c^{2}}+\cdots\,,$ (2.13)

where $f_{1}>0$, and $c$ is a constant. This behavior is the same as [6]
because the higher power term of $\phi$ in $V$ is subleading near the
singularity. The values of $f_{1}$ and $c$ depend on the specific black hole
solution. For Schwarzschild black hole we have $c=0$. The near singularity
solution can be changed to Kasner form (1.1) by coordinate transformation
$r\rightarrow 1/(\tau^{\frac{2}{3+c^{2}}})$ and the Kasner exponents are
related to $c$ as

$\displaystyle
p_{x}=p_{y}=\frac{2}{3+c^{2}}\,,~{}~{}~{}p_{t}=\frac{c^{2}-1}{3+c^{2}}\,,~{}~{}~{}p_{\phi}=\frac{2\sqrt{2}c}{3+c^{2}}\,.$
(2.14)

Note that here we focus on the cases with isotropic and homogeneous spatial
directions in the dual theory and we do not see any oscillatory behavior of
the Kasner exponents in the examples considered.999One example where the
Kasner exponents oscillate is for holographic superconductor system [7, 8].
More generally, if we consider the spatial dependence of background fields,
usually the Kasner exponents show chaotic behavior [4].  Nevertheless it would
be extremely interesting to study the behavior of singularity during the
dynamical phase transition, e.g. [30].

####  2.1.1 Radially conserved charge

The system is invariant under the transformation (2.11). There is a radially
conserved charge associated to this symmetry which is given by

$\displaystyle Q$
$\displaystyle=\frac{1}{r^{2}}\,e^{-\frac{\chi}{2}}\,\left(f^{\prime}-f\chi^{\prime}\right)\,,$
(2.15)

satisfying $Q^{\prime}=0$. The conserved charge associated to the scaling
symmetry (2.12) can also be computed from the Noether theorem and using the
equations of motion (2.5) one can show that this conserved charge is exactly
the same as (2.15). The conserved charge (2.15) can be used to check the
numerical code.

Another utility of this conserved charge comes from the fact that it gives a
constraint relations between the parameters at the boundary, horizon and the
singularity. Explicitly evaluating (2.15) at the boundary, horizon and near
the singularity gives the relation

$\displaystyle 3m_{T}-4\alpha\beta+\alpha^{3}\lambda_{3}=-\frac{4\pi
Te^{-\chi_{h}/2}}{r_{h}^{2}}=f_{1}e^{-\chi_{1}/2}\,(c^{2}-3)\,,$ (2.16)

where we have used (2.14). Therefore we have $c^{2}\leq 3$ which gives a bound
of Kasner exponents

$\displaystyle\frac{1}{3}\leq p_{x}\leq\frac{2}{3}\,,~{}~{}-\frac{1}{3}\leq
p_{t}\leq\frac{1}{3}\,,~{}~{}-\sqrt{\frac{2}{3}}\leq
p_{\phi}\leq\sqrt{\frac{2}{3}}\,.$ (2.17)

Note that the singularity in Schwarzschild black hole saturates the upper
bound of $p_{x}$ and the lower bound of $p_{t}$. The plot of Kasner exponents
as a function of $c$ is shown in Fig.1. Although from the above relations,
both $p_{x}$ and $p_{t}$ could be larger than zero due to the existence of the
scalar field, it is interesting to note that for the parameters we considered
we only find the case with negative $p_{t}$.

Figure 1: The Kasner exponents $p_{t}$ (black line), $p_{x}$ (red line) and
$p_{\phi}$ (blue line) as functions of parameter $c$ in (2.14).

###  2.2 Thermodynamics of black holes

With the numerical solution constructed in the above subsection, we can
compute the free energy density of the dual system. When there are multiple
gravitational solutions, the system might experience phase transitions. In
this subsection, we compute the thermodynamical quantities in the system under
study.

The thermodynamics of the Einstein-Scalar system is encoded in the appropriate
renormalized Euclidean action. In our setup, the action must also take into
account the more generalized nature of the boundary conditions on the scalar
field due to alternative quantization and the double trace deformation. The
action takes the form [23, 24]

$\displaystyle\begin{split}S&=\int_{\mathcal{M}}d^{4}x\;\sqrt{-g}\left[R+6-g^{ab}\partial_{a}\phi\partial_{b}\phi-V(\phi)\right]+\int_{\partial\mathcal{M}}d^{3}x\;\sqrt{-\gamma}\,\left(2K\right)\\\
&\hskip
28.45274pt+\int_{\partial\mathcal{M}}d^{3}x\;\sqrt{-\gamma}\left(-4+2\phi\,n^{a}\partial_{a}\phi+(1+\kappa
r_{c})\phi^{2}+\sigma\phi^{3}-2\lambda_{3}\phi^{3}\log
r_{c}\right)\end{split}$ (2.18)

where we have set $16\pi G=1$, the boundary $\partial\mathcal{M}$ is at
$r=r_{c}$ with $r_{c}\to 0$, $\gamma_{\mu\nu}$ and $K$ are the induced metric
and exterior curvature of boundary respectively. The parameter $\sigma$ is
fixed by the fact that the cubic interaction term should not change the double
trace boundary conditions as shown in appendix A, and we have
$\sigma=-2\lambda_{3}$.

Using the ansatz (2.4) and the near-boundary expansion (2.6), it can be seen
that the Euclidean on-shell action is given by

$\displaystyle S^{E}_{\text{on-
shell}}=\frac{V}{T}\left(m_{T}-\alpha^{2}\kappa\right)$ (2.19)

where $V$ is spatial volume. The density of free energy is therefore given by

$\displaystyle f_{\kappa}\equiv\frac{F}{V}=\frac{TS^{E}_{\text{on-
shell}}}{V}=m_{T}-\alpha^{2}\kappa\,.$ (2.20)

Note that $m_{T}$ is the coefficient in the metric field near AdS4 boundary,
$\alpha$ is the VEV of the dual operator and $\kappa$ is the double trace
deformation parameter. Therefore, with numerical solution obtained above, one
can get the corresponding free energy density. For example, for the AdS planar
black hole solution, it can be verified that
$f_{\text{Sch}}=-\frac{1}{r^{3}_{h}}=-\left(\frac{4\pi}{3}\right)^{3}T^{3}$.

It can also be checked that the free energy $\eqref{free-energy}$ also
satisfies the thermodynamic relation

$\displaystyle f_{\kappa}=\epsilon-Ts$ (2.21)

where $\epsilon$ and $s$ are the energy and entropy densities. The expression
for the energy density can be read off from the Brown-York stress tensor at
the boundary

$\displaystyle
T_{\mu\nu}=2(K_{\mu\nu}-\gamma_{\mu\nu}K)+\frac{2}{\sqrt{-\gamma}}\frac{\delta
S_{\text{c.t.}}}{\delta\gamma^{\mu\nu}}$ (2.22)

where $S_{\text{c.t.}}$ is the last term in (2.18). We obtain the energy
density

$\displaystyle\epsilon=\lim_{r_{c}\to
0}\sqrt{-\gamma}T^{0}_{0}=-2m_{T}+4\alpha\beta-\sigma\alpha^{3}-3\lambda_{3}\alpha^{3}-\alpha^{2}\kappa\,.$
(2.23)

The expression for $Ts$ can be obtained in terms of the boundary parameters
using the relation (2.16) from the conserved charge, we have

$\displaystyle-Ts=3m_{T}-4\alpha\beta+\lambda_{3}\alpha^{3}\,.$ (2.24)

This concludes the proof of the thermodynamic relation (2.21).

##  3 Behavior of the singularity across the phase transition

In the previous section we have collected all the necessary ingredients to
solve the system, we shall show the numerical results of the system in this
section. We will first show the example with second order phase transition and
analyze the behavior of the singularity across the transition. Then we will
study the first order phase transitions and also analyze the singularity of
the black hole solutions. The probes of the black hole singularities in the
context of AdS/CFT will also be studied.

###  3.1 Second order phase transition

It has been shown in [23, 27] that when the potential $V(\phi)$ in (2.1) is
symmetric under $\phi\to-\phi$, i.e. $\lambda_{3}=0$, there is a second order
phase transition between Schwarzschild black hole and hairy black hole for
$\kappa<0$. Here we generalize these results with a different choice of the
potential for the scalar field and focus on the behaviors of the black hole
singularities.

We choose parameters $\lambda_{3}=0,\lambda_{4}=1/10$ in the scalar potential
(2.2) and solve the system with double trace deformation. With the boundary
condition that the source $\kappa\alpha-\beta$ should vanish, the hairy black
hole solution can only be found below a critical $T_{c}/(-\kappa)\approx
0.616$ which takes almost the same values as in [27]. This is due to the fact
that close to the phase transition, the value of the scalar field is very
small and the temperature can be analytically computed at the leading order in
the expansion of the scalar field. The free energy densities of the black hole
solutions and the expectation values of the scalar operator as a function of
$T/(-\kappa)$ are shown is Fig. 2. In the left plot of Fig. 2, the blue line
is the free energy of the hairy black hole while the grey line is the free
energy of the Schwarzschild black hole. It is seen that below
$T_{c}/(-\kappa)$ the hairy black hole is more stable. The right plot of Fig.
2 shows the behavior of the condensation, i.e. the order parameter of the
phase transition, as a function of $T_{c}/(-\kappa)$ in the hairy black hole.
When we lower the temperature, the condensate increases up to a constant value
at zero temperature. Close to the transition temperature $T\to T_{c}$, we have
$\langle{O}\rangle/(-\kappa)\propto(1-T/T_{c})^{1/2}$ and $\delta
f/(-\kappa)^{3}=(f_{\text{Sch}}-f_{\kappa})/(-\kappa)^{3}\propto(1-T/T_{c})^{2}$.
From these behaviors, we conclude that the phase transition is of second
order.

Figure 2: Left: The free energy across second order phase transition. The grey
and blue lines are the free energy densities of Schwarzschild black hole and
hairy black hole respectively. The red dot is the free energy of the hairy
black hole at zero temperature. Right: The order parameter $\langle
O\rangle/(-\kappa)$ as a function of $T/(-\kappa)$.

At zero temperature, the strategy for constructing the bulk solution should be
modified because now the near horizon conditions are different from (2.8). In
this case, the near horizon expansion as $r\to\infty$ are

$\displaystyle\begin{split}f&=1+\frac{1}{6\lambda_{4}}+\frac{a_{1}^{2}+96\lambda_{4}+6\lambda_{4}a_{1}^{2}}{48\lambda_{4}(a_{1}-3)}\,\phi_{0}^{2}\,r^{a_{1}}\,,\\\
\chi&=\frac{a_{1}}{4}\,\phi_{0}^{2}\,r^{a_{1}}\,,\\\
\phi&=\frac{1}{\sqrt{\lambda_{4}}}+\phi_{0}\,r^{\frac{a_{1}}{2}}\,,\end{split}$
(3.1)

where
$a_{1}=\frac{3+18\lambda_{4}-\sqrt{9+204\lambda_{4}+900\lambda_{4}^{2}}}{1+6\lambda_{4}}$
which is $3-\sqrt{15}$ for the choice of $\lambda_{4}=1/10$. At the leading
order the above near horizon geometry is AdS4, and the irrelevant deformation
flow the geometry to AdS4. Note that from the scaling symmetry (2.12),
$\phi_{0}$ can be rescaled to arbitrary value by scaling the radial coordinate
$r$. Therefore only the sign of $\phi_{0}$ is crucial to flow the near horizon
to AdS4 boundary and we have a unique solution at zero temperature, which is
an AdS4 to AdS4 domain wall with different AdS radius. The free energy of the
zero temperature solution is shown as a red dot in the left plot of Fig. 2. It
can be seen that the hairy black hole is always more stable than Schwarzschild
black hole. The trajectory of the near horizon value of the scalar field
$\phi$ when the temperature is lowered down can be seen from Fig. 3. At the
transition temperature $T_{c}/(-\kappa)$, we have $\phi_{h}=0$ and at zero
temperature we have $\phi_{h}=1/\sqrt{\lambda_{4}}$. Note that there is a
$Z_{2}$ symmetry $\phi\to-\phi$ for the system. We could have other branch of
solution with the same behavior as above.

Figure 3: When the temperature is lowered down, the flow of the near horizon
value of the scalar field is shown in the red line with arrow.

Next, we study the Kasner exponent across the second order phase transition
above. With the near horizon conditions for the hairy black hole solution, we
can integrate the system towards the singularity from the horizon, and obtain
the Kasner exponents numerically. One of the Kasner exponents, $p_{t}$ as a
function of $T/(-\kappa)$, is shown in Fig. 4, where the blue line is for
hairy black hole while the grey line is for Schwarzschild black hole. One
finds that during the second order phase transition at $T_{c}/(-\kappa)$, the
Kasner exponent $p_{t}$ is continuous while the first derivative of $p_{t}$
with respect to $T/(-\kappa)$ is not continuous. Since from (2.14) only one of
the three Kasner exponents is independent, we conclude that the Kanser
exponents are continuous and non-differentiable during the second order phase
transitions.

Figure 4: The Kasner exponent $p_{t}$ as a function of $T/(-\kappa)$ for
$\lambda_{3}=0,\lambda_{4}=1/10$. The grey line is $p_{t}$ of Schwarzschild
singularity which is $-1/3$. The blue line is for $p_{t}$ of the hairy black
hole solution. When $T/(-\kappa)$ goes to zero, $p_{t}$ goes to $-1/3$.

The above behavior of the Kasner exponents across the second order mean field
phase transition seems universal. Firstly, the same behavior can be found in
other holographic examples. The Kasner exponents are also studied in the
Einstein-Maxwell theory coupled to neutral scalar field in [6] and in Fig. 6
of [6] one example of phase transition was shown. From this we could see that
Kasner exponent is continuous while its first derivative with respect to
$T/\mu$ is not for a continuous phase transition if we identify $p_{t}=1$ for
the singularity of RN AdS black hole. Secondly, one naive analytical argument
could follow from the fact that the metric and matter fields are analytical
functions of radial coordinate $r$ at any temperature.101010We thank Run-Qiu
Yang for a helpful discussion on this point. Note that close to the transition
temperature the scalar field is very small. From the side of the hairy black
hole we have $\mathcal{O}\propto(1-T/T_{c})^{1/2}$, thus we expect that
$\phi(r)\propto(1-T/T_{c})^{1/2}$. From (2.13) we have
$c\propto(1-T/T_{c})^{1/2}$. Therefore, from (2.14) we have
$p_{t}+\frac{1}{3}\propto(1-T/T_{c})$. This means that in general for a second
order mean field phase transition, the Kasner exponent is continuous while the
first order derivative with respect to $T$ is not continuous.

###  3.2 First order phase transition

Now we consider the case with first order phase transition and study the
behavior of the singularity across the phase transitions. We choose a non-
symmetric scalar potential, and without loss of generality we set
$\lambda_{3}=1/8,\lambda_{4}=1/10$ in (2.2).

We find that the hairy black hole solution can exist at arbitrary value of
$T/(-\kappa)$. We plot the free energy of the hairy black hole solutions (in
blue) and Schwarzschild solution (in grey) in the left plot of Fig. 5. When we
lower the temperature, we find that there is a first order phase transition
$T_{c1}/(-\kappa)\approx 0.666$ where the first order derivative of free
energy density with respect to $T/(-\kappa)$ is discontinuous. The behavior of
the order parameter $\langle{O}\rangle/(-\kappa)$ as a function of
$T/(-\kappa)$ is shown in the right plot of Fig. 5 where the solid lines are
for stabler phases. At the transition temperature $T_{c1}/(-\kappa)$, the
condensation $\langle{O}\rangle/(-\kappa)$ has a jump to a negative value and
then decreasing until a constant value when we lower temperature. Therefore
the system goes through a first order phase transition.

Figure 5: Left: The free energy of the black hole solutions with
$\lambda_{3}=1/8$ and $\lambda_{4}=1/10$. The grey line is for Schwarzschild
solution and the other lines are for hairy black hole solutions. Right: The
VEV $\langle O\rangle$ as a function of $T/(-\kappa)$. In both plots, the
solid blue line is for stable hairy black hole solution while the dotted blue
lines are for unstable hairy black hole solutions.

At zero temperature, we have the near horizon boundary condition as
$r\to\infty$

$\displaystyle\begin{split}f&=c_{0}+c_{1}\,\phi_{0}^{2}\,r^{a_{2}}\,,\\\
\chi&=c_{2}\,\phi_{0}^{2}\,r^{a_{2}}\,,\\\
\phi&=\frac{3\lambda_{3}-\sqrt{9\lambda_{3}^{2}+64\lambda_{4}}}{8\lambda_{4}}+\phi_{0}\,r^{\frac{a_{2}}{2}}\,,\end{split}$
(3.2)

where $c_{0}$, $c_{1},c_{2}$ and $a_{2}$ are functions of $\lambda_{3}$ and
$\lambda_{4}$. For the case we considered, i.e. $\lambda_{3}=1/8$ and
$\lambda_{4}=1/10$, we have
$(c_{0},c_{1},c_{2},a_{2})=(2.135,-0.503,-0.235,-0.942)$. Note that $\phi_{0}$
can be rescaled to arbitrary value by the symmetry (2.12). We choose specific
$\phi_{0}$ to fix $\kappa=-1$. The free energy at zero temperature from the
above solution $f_{\kappa 0}/(-\kappa)^{3}=-3.125$ is smaller than the free
energy of Schwarzschild solution. Note that for boundary condition (3.2),
there are two $\phi_{0}^{\prime}$s corresponding to $\kappa=-1$. We have
chosen the more stabler one and the free energy of the other one is
$f_{\kappa}/(-\kappa)^{3}=-3.114$. Note that in (3.2) we chose the near
horizon value of $\phi$ to be the local minimal of the potential of the scalar
field. Since there exist two minima of the potential, there are two different
near horizon conditions and the other one has the form of
$\phi=\frac{3\lambda_{3}+\sqrt{9\lambda_{3}^{2}+64\lambda_{4}}}{8\lambda_{4}}+\tilde{\phi}_{0}\,r^{\frac{\tilde{a}_{2}}{2}}$.
However, the free energy of this one is $f_{\kappa}/(-\kappa)^{3}=-0.803$ and
it is bigger than $f_{\kappa 0}/(-\kappa)^{3}$. The trajectory of the near
horizon value of the scalar filed when we lower temperature is shown as the
solid red line with arrow in Fig. 6. Different from the previous symmetric
case, now the potential is not symmetric and the trajectory of the near
horizon value of the scalar filed is unique. When we lower down the
temperature, at $T_{c1}/(-\kappa)$, the near horizon value $\phi_{h}$ jumps
from 0 to a specific negative number and then flows to the left minimal of the
potential witch corresponding to the truly zero temperature solution. The
stabler phase can be seen from the horizon behavior of the scalar field. It is
interesting to compare with the study of holographic phase transitions in [31,
32] without double trace deformation, where a jump of the near horizon value
of scalar field was also found and is related to the existence of the
noncontinuous phase transition.

Figure 6: The red solid line with arrows denote the trajectory of the near
horizon value of the scalar field when we lower the temperature.

The Kasner exponents $p_{t}$ and $p_{\phi}$ of the singularities are shown in
Fig. 7 where the Kasner exponents of the hairy black hole and Schwarzschild
black hole are in blue and grey respectively. Below $T_{c1}/(-\kappa)$, the
Kasner exponents of the stable phases are in solid blue. Different from the
case of continuous phase transition, when the phase transition is of first
order, the Kasner exponents are discontinuous.

Figure 7: The Kasner exponents $p_{t}$ (left) and $p_{\phi}$ (right) of the
black hole solutions as a function of $T/(-\kappa)$. In both plots, the gray
line is for Schwartzschild black hole, the solid blue line is for stable hairy
black hole solution while the dotted blue lines are for unstable hairy black
hole solutions.

It is natural to expect that this discontinuity of the Kasner exponent is a
universal behavior for discontinuous phase transitions. Firstly, the exactly
same behavior can be found in other example of first order Hawking-Page
transition. For holographic CFTs on $S^{d-1}$ at finite temperature, there is
a well-known first order Hawking-Page transition [33] between the theory whose
dual geometry is described by a black hole in global AdSd+1 at high
temperature and the theory which is dual to a thermal AdS spacetime at low
temperature. The Schwarzschild black hole in AdS has Kasner exponents
$(p_{t},p_{\omega_{i}})=(\frac{2}{d}-1,\frac{2}{d})$ with $i=1,\dots,d-1$ and
$\omega_{i}$ the $i$-th spherical coordinate. It is easy to check that the
Kasner relations
$p_{t}+(d-1)p_{\omega_{i}}=p_{t}^{2}+(d-1)p_{\omega_{i}}^{2}=1$ are satisfied.
In the thermal AdS we do not have any singularity and the geometry is no
longer of Kasner form (1.1). In this sense the Kasner exponents of the
singularity are discontinuous during the first order phase transition. From
the holographic first order phase transition between the AdS Schwartzschild
black hole in planar coordinate and AdS-soliton one can obtain the same
conclusion. Secondly, in the Einstein-scalar theory with double trace
deformation, during the first order phase transition, the scalar field is
discontinuous, therefore we expect that the value $c$ in the near singularity
behavior of the scalar field (2.13) is discontinuous. From (2.14) we know that
the Kasner exponents are discontinuous during the first order phase
transitions.

###  3.3 Probes of the black hole singularity via geodesics

In this subsection, we will study some observable quantities of boundary field
theory from which we can read the information about Kasner exponents of the
black hole’s singularity. We shall focus on the geodesics since the bulk
geodesics can approach the singularity of the black hole in certain limit. In
the context of AdS/CFT duality, the length of bulk spacelike geodesics
corresponds to the correlation functions of the operators in the large
conformal dimensional limit of dual boundary field theory [15, 16]. Meanwhile,
the connection between the proper time from the black hole horizon to the
singularity and one-point functions were proposed in [17]. In the following we
will compute these two types of geodesics respectively.

We consider the radical geodesics for which
$g_{tt}\dot{t}^{2}+g_{rr}\dot{r}^{2}=\varepsilon$, where the dot denotes the
derivative with respect to the proper time $\tau$ and $\varepsilon=-1,0,1$
gives timelike, null and spacelike geodesics respectively. There is a
conserved charge $E=-g_{tt}\dot{t}$ along the geodesic. For timelike or null
geodesics $E$ characterize the conserved energy of the particle. The equation
of motion of the geodesic is

$\displaystyle\frac{E^{2}}{g_{tt}}+g_{rr}\dot{r}^{2}=\varepsilon$ (3.3)

from which we obtain

$\displaystyle\frac{dr}{d\tau}=\sqrt{\frac{\varepsilon
g_{tt}-E^{2}}{g_{tt}g_{rr}}}\,.$ (3.4)

In the following we will solve this equation for timelike geodesics and
spacelike geodesics respectively.

####  3.3.1 Radial timelike geodesics

Recently, it was proposed in [17] that the proper time $\tau_{s}$ from the
black hole horizon to the singularity is related to the expectation value of
certain operators in the large mass limit via
$\langle\mathcal{O}\rangle\propto e^{-im\tau_{s}-ml_{\text{hor}}}$, where $m$
is a complixified mass with $\text{Im}(m)<0$ and $l_{\text{hor}}$ is the
distance from horizon to the boundary. Therefore, from the one point function
of the operator one can extract the information of Kasner exponent once we
know the relation between $\tau_{s}$ and the Kasner exponents.

For timelike geodesics we have $\varepsilon=-1$ in (3.4). Plugging in the
ansatz for the metric we obtain

$\displaystyle\frac{d\tau}{dr}=\frac{e^{-\chi/2}}{r^{2}\sqrt{-E^{2}-(e^{-\chi}f/r^{2})}}\,.$
(3.5)

The proper time from horizon to singularity of a particle with $E=0$ is

$\displaystyle\tau_{s}=\int_{r_{h}}^{\infty}d\tau=\int_{r_{h}}^{\infty}\frac{dr}{r\sqrt{-f(r)}}\,.$
(3.6)

The proper time $\tau_{s}$ as a function of $T/(-\kappa)$ and $p_{t}$ are
shown in Fig. 8. The left two plots characterize that when we lower the
temperature, the behavior of the $\tau_{s}$, or equivalently the expectation
value of the operator $\mathcal{O}$. For this $\tau_{s}$, one can extract the
information of the Kasner exponent $p_{t}$ according to the right plots in
Fig. 8.

Figure 8: The proper time $\tau_{s}$ as a function of $T/(-\kappa)$ (left) and
$p_{t}$ (right) for the model in subsection 3.1 (top plots) and in subsection
3.2 (down plots). For these plots, the grey line is for Schwartzschild black
hole, the solid blue line is for stable hairy black hole solution while the
dotted blue lines are for unstable hairy black hole solutions.

####  3.3.2 Radial spacelike geodesics

The spacelike geodesics are well studied in the literature as they are the
relevant geodesics for the computation of correlators in AdS/CFT [15, 16]. A
significant difference between the spacelike geodesics from the timelike
geodesic is that they do not end in singularity. Instead the geodesics reach a
finite radius $r_{*}$, namely the turning point, which is controlled by the
energy $E$, after which it bounces back and goes to the other side of the
eternal black hole. Nevertheless, when $E\to\infty$, the turning point will
approach very close to the singularity and therefore in this limit one could
extract the information of the singularity with the spacelike geodesics.

For the spacelike geodesics, we set $\varepsilon=1$ in (3.4). The regularized
proper length is given by

$\displaystyle
L=2\int_{r_{c}}^{r_{*}}dr\frac{e^{-\chi/2}}{r^{2}\sqrt{E^{2}+\frac{1}{r^{2}}fe^{-\chi}}}+2\log
r_{c}$ (3.7)

where $r_{*}$ is the turning point which is given by $g_{tt}(r_{*})=E^{2}$ and
for real $E$, $r_{*}$ lies inside the horizon, and $r_{c}$ is the cutoff close
to the AdS boundary. As $E$ becomes large, $r_{*}$ tends towards the
singularity. At large $E$, most of the contribution comes from the boundary
and we may approximate the integral by an expansion in $1/E$.111111One
subtlety of the expansion is that the limit $r_{c}\to 0$ and $E\to\infty$ do
not commute with each other (see [5] for details). The resulting expansion is
of the form

$\displaystyle\begin{split}L&=2\log\frac{2}{E}+\frac{l_{1}}{E}+\frac{\alpha^{2}}{2}\frac{1}{E^{2}}+\frac{l_{3}}{E^{3}}+\Big{(}m_{T}-2\alpha\beta+\alpha^{3}\lambda_{3}-\alpha^{3}\lambda_{3}\log
2\Big{)}\,\frac{\log E}{E^{3}}\\\ &\hskip
113.81102pt-\frac{1}{2}\alpha^{3}\lambda_{3}\,\frac{\log^{2}E}{E^{3}}+l^{\prime}E^{\frac{1}{p_{t}}}\,,\end{split}$
(3.8)

where the coefficients $l_{1},l_{3}$ depend on the metric fields along the
radial direction and

$\displaystyle\begin{split}l^{\prime}&=\sqrt{\pi}(p_{t}-1)\frac{e^{\frac{\chi_{1}}{2p_{t}}}}{f_{1}^{\frac{p_{t}+1}{2p_{t}}}}\frac{\Gamma(\frac{p_{t}+1}{2p_{t}})}{\Gamma(\frac{1}{2p_{t}})}\,,\end{split}$
(3.9)

where $\chi_{1}$ and $f_{1}$ are defined in (2.13). The information of the
Kasner exponent of singularity is encoded in the non-analytic part of the
proper length $L$.

From the AdS/CFT correspondence, the proper length of the spacelike geodesic
is related to the correlators of the operators in the large conformal
dimension limit. Therefore, in principle, from the correlation functions of
the field theory one could obtain the Kasner exponent $p_{t}$.

##  4 Conclusion and discussion

We have studied the behaviors of the singularities in Einstein-scalar theory
with a double trace deformation in which the second order and the first order
phase transitions could be realized. We found that when the phase transition
is of second order, the Kasner exponents of the black hole singularity are
continous while the first derivative of the Kasner exponents with respect to
the temperature are not continuous. When the phase transition is of first
order, the Kasner exponents are not continuous. Our study confirms that the
physics of black hole interior is related to the physics outside the black
hole horizon. Finally, the Kasner exponents could be obtained from the
behavior of the field theoretical observable which are related to the
geodesics in the bulk.

There are lots of open questions to be explored. Firstly, since we only
considered a particular holographic model with a special choice of parameters,
it would be interesting to check how generic our results are. Especially, to
check if the behavior of Kasner exponents during the phase transition are
universal. Secondly, it would be interesting to study the evolution of black
hole singularities during the time-dependent dynamical phase transitions, e.g.
in [30, 34]. This might help us further understand the possible universal
physics of non-equilibrium physics. Thirdly, because the black hole interior
changes during the phase transition, it would be interesting to construct the
evaporating black hole solutions with phase transitions. In this way, when the
black evaporates, the changes of the black hole interior should have
nontrivial effect on the Page curve and checking the entropy of the Hawking
radiation in this case might help us further understand the physics of the
islands proposal [35].121212Studies on the Page curve for the eternal hairy
black hole can be found in e.g. [36]. Finally, the physics inside the black
hole interior might have some connections to black hole information paradox
[35, 37]. It would be interesting to understand further properties of dynamics
of quantum fields close to the black hole singularities and make connections
to the physics at the boundary.

### Acknowledgments

We would like to thank Li Li, Hong Lv, Hua-Jia Wang, Xin-Meng Wu, Run-Qiu
Yang, Xin-Yi Zhang for useful discussions. This work is supported by the
National Natural Science Foundation of China grant No.11875083. A.R is also
supported by Zhuoyue Postdoc Fellowship of Beihang University (ZYBH2018-01).

##  Appendix A The boundary terms for the scalar field in holographic
renormalization

In this appendix we describe a strategy to fix the free parameter $\sigma$ in
the counter-term action (2.18) for the scalar field in holographic
renormalization [24, 23].

Let us first consider the known case of symmetric potential, _i.e._
$\lambda_{3}=0$. In this case it is known that $\sigma=0$. The variation of
the (renormalized) action (2.18) with respect to the scalar field $\phi$ is,

$\displaystyle\delta
S\supset\int_{\mathcal{M}}d^{4}x\;\sqrt{-g}\,(\text{EOM})\,\delta\phi+\int_{\partial\mathcal{M}}d^{3}x\;\frac{e^{\chi/2}\sqrt{f}}{r^{3}}\left[-2r\sqrt{f}\;\phi\delta\phi^{\prime}+2(1+\kappa
r)\phi\delta\phi\right]\,,$ (A.1)

where $\partial\mathcal{M}$ is defined at the boundary $r=r_{c}$ with
$r_{c}\rightarrow 0$. Substituting the expansion (2.6) into (A.1), it reduces
to

$\displaystyle\delta
S\supset\int_{\partial\mathcal{M}}d^{3}x\;2\alpha\,\delta\left(\kappa\alpha-\beta\right)\,.$
(A.2)

The above expression shows that with double trace deformation the source for
the scalar operator is $\kappa\alpha-\beta$ while the expectation value is
$\alpha$.

In the non-symmetric potential case, i.e.$\lambda_{3}\neq 0$, we expect that
(A.2) continues to hold. Now the variation of the action with respect to the
scalar filed $\phi$ is given by

$\displaystyle\delta
S\supset\int_{\partial\mathcal{M}}d^{3}x\;\left(2\alpha\delta(\kappa\alpha-\beta)+3\sigma\alpha^{2}\delta\alpha+6\lambda_{3}\alpha^{2}\delta\alpha\right)\,.$
(A.3)

Comparing this to (A.2), we have $\sigma=-2\lambda_{3}$.

## References

  * [1] S. W. Hawking and R. Penrose, The Singularities of gravitational collapse and cosmology, Proc. Roy. Soc. Lond. A 314 (1970), 529-548.
  * [2] V. A. Belinsky, I. M. Khalatnikov, and E. M. Lifshitz, Oscillatory approach to a singular point in the relativistic cosmology, Adv. Phys. 19 (1970) 525; A general solution of the Einstein equations with a time singularity, Adv. Phys. 31, 639 (1982).
  * [3] E. Kasner, Geometrical theorems on Einstein’s cosmological equations, Am. J. Math. 43, 217-221 (1921).
  * [4] V. Belinski and M. Henneaux, The Cosmological Singularity, Cambridge Universe Press.
  * [5] A. Frenkel, S. A. Hartnoll, J. Kruthoff and Z. D. Shi, Holographic flows from CFT to the Kasner universe, JHEP 08, 003 (2020) [arXiv:2004.01192].
  * [6] S. A. Hartnoll, G. T. Horowitz, J. Kruthoff and J. E. Santos, Gravitational duals to the grand canonical ensemble abhor Cauchy horizons, JHEP 10 (2020), 102 [arXiv:2006.10056].
  * [7] S. A. Hartnoll, G. T. Horowitz, J. Kruthoff and J. E. Santos, Diving into a holographic superconductor, SciPost Phys. 10 (2021), 009 [arXiv:2008.12786].
  * [8] R. G. Cai, L. Li and R. Q. Yang, No Inner-Horizon Theorem for Black Holes with Charged Scalar Hair, [arXiv:2009.05520].
  * [9] Y. Q. Wang, Y. Song, Q. Xiang, S. W. Wei, T. Zhu and Y. X. Liu, Holographic flows with scalar self-interaction toward the Kasner universe, [arXiv:2009.06277].
  * [10] R. Q. Yang, R. G. Cai and L. Li, Constraining the number of horizons with energy conditions, [arXiv:2104.03012].
  * [11] S. A. H. Mansoori, L. Li, M. Rafiee and M. Baggioli, What’s inside a hairy black hole in massive gravity?, [arXiv:2108.01471].
  * [12] E. Shaghoulian and H. Wang, Timelike BKL singularities and chaos in AdS/CFT, Class. Quant. Grav. 33 (2016) no.12, 125020 [arXiv:1601.02599].
  * [13] V. A. Belinski and I. M. Khalatnikov, Effect of scalar and vector fields on the nature of the cosmological singularity, Sov. Phys. JETP, 36(4):591, 1973.
  * [14] T. Damour, M. Henneaux and H. Nicolai, Cosmological billiards, Class. Quant. Grav. 20 (2003), R145-R200 [arXiv:hep-th/0212256].
  * [15] L. Fidkowski, V. Hubeny, M. Kleban and S. Shenker, The Black hole singularity in AdS / CFT, JHEP 02, 014 (2004) [arXiv:hep-th/0306170].
  * [16] G. Festuccia and H. Liu, Excursions beyond the horizon: Black hole singularities in Yang-Mills theories. I., JHEP 04 (2006), 044 [arXiv:hep-th/0506202].
  * [17] M. Grinberg and J. Maldacena, Proper time to the black hole singularity from thermal one-point functions, JHEP 03 (2021), 131, [arXiv:2011.01004].
  * [18] T. Hartman and J. Maldacena, Time Evolution of Entanglement Entropy from Black Hole Interiors, JHEP 05 (2013), 014 [arXiv:1303.1080].
  * [19] L. Susskind, L. Thorlacius and J. Uglum, The Stretched horizon and black hole complementarity, Phys. Rev. D 48, 3743-3761 (1993) [arXiv:hep-th/9306069].
  * [20] I. R. Klebanov and E. Witten, AdS / CFT correspondence and symmetry breaking, Nucl. Phys. B 556 (1999), 89-114 [arXiv:hep-th/9905104].
  * [21] E. Witten, Multitrace operators, boundary conditions, and AdS / CFT correspondence, [arXiv:hep-th/0112258].
  * [22] M. Berkooz, A. Sever and A. Shomer, Double trace deformations, boundary conditions and space-time singularities, JHEP 0205 (2002) 034 [arXiv:hep-th/0112264].
  * [23] T. Faulkner, G. T. Horowitz and M. M. Roberts, Holographic quantum criticality from multi-trace deformations, JHEP 04 (2011), 051 [arXiv:1008.1581].
  * [24] I. Papadimitriou, Multi-Trace Deformations in AdS/CFT: Exploring the Vacuum Structure of the Deformed CFT, JHEP 05, 075 (2007) [arXiv:hep-th/0703152].
  * [25] J. H. She, B. J. Overbosch, Y. W. Sun, Y. Liu, K. Schalm, J. A. Mydosh and J. Zaanen, Observing the origin of superconductivity in quantum critical metals, Phys. Rev. B 84 (2011), 14452 [arXiv:1105.5377].
  * [26] N. Iqbal, H. Liu and M. Mezei, Quantum phase transitions in semilocal quantum liquids, Phys. Rev. D 91 (2015) no.2, 025024 [arXiv:1108.0425].
  * [27] E. Mefford and G. T. Horowitz, Simple holographic insulator, Phys. Rev. D 90, no.8, 084042 (2014) [arXiv:1406.4188].
  * [28] D. M. Hofman and N. Iqbal, Generalized global symmetries and holography, SciPost Phys. 4 (2018) no.1, 005 [arXiv:1707.08577].
  * [29] S. Grozdanov and N. Poovuttikul, Generalised global symmetries in holography: magnetohydrodynamic waves in a strongly interacting plasma, JHEP 04 (2019), 141 [arXiv:1707.04182].
  * [30] R. A. Janik, J. Jankowski and H. Soltanpanahi, Real-Time dynamics and phase separation in a holographic first order phase transition, Phys. Rev. Lett. 119 (2017) no.26, 261601 [arXiv:1704.05387].
  * [31] U. GÃ¼rsoy, E. Kiritsis, F. Nitti and L. Silva Pimenta, Exotic holographic RG flows at finite temperature, JHEP 10, 173 (2018) [arXiv:1805.01769].
  * [32] Y. Bea and D. Mateos, Heating up Exotic RG Flows with Holography, JHEP 08 (2018), 034 [arXiv:1805.01806].
  * [33] E. Witten, Anti-de Sitter space, thermal phase transition, and confinement in gauge theories, Adv. Theor. Math. Phys. 2 (1998), 505-532 [arXiv:hep-th/9803131].
  * [34] L. Bellantuono, R. A. Janik, J. Jankowski and H. Soltanpanahi, Dynamics near a first order phase transition, JHEP 10 (2019), 146 [arXiv:1906.00061].
  * [35] A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian and A. Tajdini, The entropy of Hawking radiation, [arXiv:2006.06872].
  * [36] E. Caceres, A. Kundu, A. K. Patra and S. Shashi, Page Curves and Bath Deformations, [arXiv:2107.00022].
  * [37] M. J. Perry, No Future in Black Holes, [arXiv:2106.03715].