# Three-dimensional de Sitter holography and bulk correlators at late time

Heng-Yu Chena,b and Yasuaki Hikidac aDepartment of Physics, National Taiwan
University, Taipei 10617, Taiwan  bPhysics Division, National Center for
Theoretical Sciences, Taipei 10617, Taiwan  cCenter for Gravitational Physics
and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto 606-8502, Japan

###### Abstract

We propose an explicitly calculable example of holography on 3-dimensional de
Sitter space by providing a prescription to analytic continue a higher-spin
holography on 3-dimensional anti-de Sitter space. Applying the de Sitter
holography, we explicitly compute bulk correlation functions on 3-dimensional
de Sitter space at late time in a higher-spin gravity. These expressions are
consistent with recent analysis based on bulk Feynman diagrams. Our explicit
computations reveal how holographic computations could provide fruitful
information.

††preprint: YITP-22-37

##  I Introduction

Quantum gravity on de Sitter (dS) background is important to understand our
early universe. In particular, the late time correlation functions on dS have
been argued to give a clue to understand what happened at inflation era
Maldacena (2003); Maldacena and Pimentel (2011); Arkani-Hamed and Maldacena
(2015); Arkani-Hamed  _et al._ (2020). It is expected that dS/CFT
correspondence Strominger (2001); Witten (2001a); Maldacena (2003) is useful
to formulate quantum gravity on dS as the better explored AdS/CFT
correspondence Maldacena (1998) does for quantum gravity on anti-de Sitter
(AdS) space. Until recently, however, there has been only one concrete
explicit prototype of dS/CFT correspondence, which uses higher-spin gravity on
4-dimensional dS Anninos  _et al._ (2017). The duality can be understood as an
analytic continuation of higher-spin AdS4 holography by Klebanov and Polyakov
(2002), see Sundborg (2001); Mikhailov (2002); Witten (2001b); Sezgin and
Sundell (2002) for previous works. Recently, another explicit proposal was
made on 3-dimensional de Sitter holography Hikida  _et al._ (2022a, b), see
Ouyang (2011) for a previous attempt. The purpose of this letter is to compute
bulk correlation functions on 3-dimensional dS at late time by elaborating the
proposal furthermore. To our best knowledge, this is the first computation of
bulk dS correlators at late time based on an explicit concrete holographic
setup. We compare our results with the recent analysis based on direct bulk
Feynman diagram computations Sleight and Taronna (2021a, b). See also Shukla
_et al._ (2016); Di Pietro  _et al._ (2022) as well. Since lower dimensional
examples are much tractable, we expect to learn a lot about mysterious dS/CFT
correspondence through them.

For our purpose, we develop holographic method for bulk computations and
higher-spin dS3 holography. We first identify the phases relating the dS
boundary operators with their AdS counterparts due to the analytic
continuation from AdS. We then elaborate the dS3/CFT2 correspondence in Hikida
_et al._ (2022a, b). In order to analyze bulk correlators, we consider an
analytic continuation of Gaberdiel-Gopakumar duality Gaberdiel and Gopakumar
(2011), which is between 3d Prokushkin-Vasiliev theory Prokushkin and Vasiliev
(1999) and 2d coset model

$\displaystyle\frac{\text{SU}(N)_{k}\times\text{SU}(N)_{1}}{\text{SU}(N)_{k+1}}\,.$
(1)

This coset was proven to describe WN-minimal model Arakawa  _et al._ (2019).
We provide a prescription to perform an analytic continuation, which is
different from the previous one 111In Hikida  _et al._ (2022a, b) a limit of
$k\to-N$ is considered, while in this letter an ’t Hooft limit is taken with
large $|N|,|k|$ but finite $\lambda=N/(N+k)$. We then apply the Maldacena’s
holographic prescription to compute bulk dS correlators at late time given in
Maldacena (2003) and find explicit expressions for 2- and 3-pt. correlators.
We also examine a simple 4-pt. correlator. We finally compare our results to
generic arguments on bulk Feynman diagrams in Sleight and Taronna (2021a, b)
and comment on the advantage of holographic approach.

##  II Preliminary

We start by reviewing the AdS/CFT correspondence in order to explain how
dS/CFT correspondence is obtained from an analytic continuation. We consider a
gravity theory on $d+1$-dimensional AdS and $d$-dimensional CFT on the AdS
boundary. In the bulk theory, we assume there exist symmetric tensor fields
$\sigma_{i_{1}\cdots i_{s}}^{\text{AdS}}$ with mass $m$ and spin $s$, and the
scalar fields correspond to the ones with $s=0$. The CFT operators dual to the
bulk fields are denoted by $J^{i_{1}\cdots i_{s}}_{\text{AdS}}$. For the
Euclidean AdS metric, we use the Poincaré coordinates,

$\displaystyle\mathrm{d}s^{2}=\ell_{\text{AdS}}^{2}\,\frac{\mathrm{d}y^{2}+\mathrm{d}\vec{x}^{2}}{y^{2}}$
(2)

with the AdS radius $\ell_{\text{AdS}}$. We consider the region with $y\geq 0$
and the boundary is located at $y=0$. The bulk field behaves near $y=0$ as

$\displaystyle\sigma_{i_{1}\cdots
i_{s}}^{\text{AdS}}(y,\vec{x})\sim\sigma_{i_{1}\cdots
i_{s}}^{+,\text{AdS}}(\vec{x})\,y^{\Delta_{+}-s}+\sigma_{i_{1}\cdots
i_{s}}^{-,\text{AdS}}(\vec{x})\,y^{\Delta_{-}-s}$ (3)

with

$\displaystyle\ell_{\text{AdS}}^{2}\,m^{2}=-(\Delta_{+}\Delta_{-}+s)\,,\qquad\Delta_{-}=d-\Delta_{+}\
.$ (4)

We use the flat boundary metric as $\mathrm{d}s^{2}=\mathrm{d}\vec{x}^{2}$,
and the coupling between the bulk field and its dual operator is

$\displaystyle\ell^{d}_{\text{AdS}}\int\mathrm{d}^{d}\vec{x}\,\sigma^{\pm,\text{AdS}}_{i_{1}\cdots
i_{s}}\,J^{i_{1}\cdots i_{s}}_{\pm,\text{AdS}}\,,$ (5)

where $\ell_{\text{AdS}}$-dependent factor is explicitly written.

We then consider dS/CFT correspondence. The Poincaré patch for Lorentzian
dSd+1 is described by

$\displaystyle\mathrm{d}s^{2}=\ell^{2}\,\frac{-\mathrm{d}\eta^{2}+\mathrm{d}\vec{x}^{2}}{\eta^{2}}$
(6)

with dS radius $\ell$. We consider the region $\eta\leq 0$ and a boundary is
located at the future infinity $\eta\to-0$. The metric can be related to (2)
by

$\displaystyle y=-i\eta\ ,\qquad\ell_{\text{AdS}}=-i\ell\,.$ (7)

As in the case of AdSd+1, we consider bulk tensor fields on dSd+1 denoted by
$\sigma_{i_{1}\cdots i_{s}}$, which are dual to boundary operators
$J^{i_{1}\cdots i_{s}}$. In particular, the conformal dimension of the dual
operator is

$\displaystyle\ell^{2}\,m^{2}=\Delta_{+}\Delta_{-}+s\,,\qquad\Delta_{-}=d-\Delta_{+}\
.$ (8)

In this letter, we consider the case with small $\ell^{2}m^{2}$ in order to
avoid subtlety associated with complex $\Delta_{\pm}$. Near $\eta\to-0$, the
bulk field behaves as

$\displaystyle\begin{aligned} &\sigma_{i_{1}\cdots i_{s}}(\eta,\vec{x})\\\
&\sim\sigma^{+}_{i_{1}\cdots
i_{s}}(\vec{x})\,(-\eta)^{\Delta_{+}-s}+\sigma^{-}_{i_{1}\cdots
i_{s}}(\vec{x})\,(-\eta)^{\Delta_{-}-s}\,.\end{aligned}$ (9)

Following Maldacena (2003), we may identify $\sigma_{i_{1}\cdots
i_{s}}^{\text{AdS}}$ with $\sigma_{i_{1}\cdots i_{s}}$ by generalizing the map
(7) as follows. Notice that bulk indices are lowered and raised by the bulk
metric $g_{\mu\nu}\sim\ell^{2}$ but the boundary metric is independent of
$\ell$, see, e.g., Anninos  _et al._ (2017). It is thus natural to assign

$\displaystyle\ell^{-s}\,\sigma^{\pm}_{i_{1}\cdots
i_{s}}(\vec{x})=i^{\Delta_{\pm}-s}\,\ell_{\text{AdS}}^{-s}\,\sigma^{\pm,\text{AdS}}_{i_{1}\cdots
i_{s}}(\vec{x})\,,$ (10)

and we have

$\displaystyle\begin{aligned} &\sigma^{\pm}_{i_{1}\cdots
i_{s}}(\vec{x})=e^{i\frac{\pi}{2}\Delta_{\pm}}\,\sigma^{\pm,\text{AdS}}_{i_{1}\cdots
i_{s}}(\vec{x})\,,\\\ &J^{\pm}_{i_{1}\cdots
i_{s}}(\vec{x})=e^{i\frac{\pi}{2}(d-\Delta_{\pm})}\,J^{\pm,\text{AdS}}_{i_{1}\cdots
i_{s}}(\vec{x})\,.\end{aligned}$ (11)

For conserved currents with $\Delta_{+}=s+d-2$, symmetric traceless tensors
are

$\displaystyle J^{+}_{i_{1}\cdots
i_{s}}(\vec{x})=e^{i\frac{\pi}{2}(2-s)}\,J^{+,\text{AdS}}_{i_{1}\cdots
i_{s}}(\vec{x})\,.$ (12)

In particular the energy momentum tensor with $s=2$ does not receive any phase
factor and the standard operator product expansions are preserved.

As mentioned above, we would like to compute bulk correlators on dS at late
time by following the prescription of Maldacena (2003), which may be
summarized as

$\displaystyle\Psi[\sigma^{\pm}_{i_{1}\cdots
i_{s}}]=\left\langle\exp\left(\ell^{d}\int\mathrm{d}^{d}\vec{x}\,\sigma^{\pm}_{i_{1}\cdots
i_{s}}\,J^{i_{1}\cdots i_{s}}_{\pm}\right)\right\rangle\,.$ (13)

The left hand side is the wave functional of universe for a fixed
$d$-dimensional metric $g_{\mu\nu}=h_{\mu\nu}$ at late time $\eta\to-0$. The
right hand side is computed by a certain CFT with operators $J^{i_{1}\cdots
i_{s}}_{\pm}$ coupled with their dual bulk fields $\sigma^{\pm}_{i_{1}\dots
i_{s}}$. The bulk correlation functions are then computed as expectation
values as

$\displaystyle\left\langle\,\prod_{j=1}^{m}\psi_{j}(\vec{x}_{j})\,\right\rangle=\int[{\cal
D}\psi_{i}]\,\left|\Psi\left[\psi_{l}\right]\right|^{2}\,\prod_{j=1}^{m}\psi_{j}(\vec{x}_{j})\,,$
(14)

where we set $\psi_{j}=\sigma^{\pm}_{i_{1}\cdots i_{s_{j}}}$. See appendix A
and, e.g., Arkani-Hamed  _et al._ (2020) for more details.

##  III Higher-spin AdS3 holography

In this letter, we construct an explicit example of dS3/CFT2 correspondence by
suitable analytic continuation of Gaberdiel-Gopakumar duality Gaberdiel and
Gopakumar (2011) as in Hikida  _et al._ (2022a, b). The gravity side of this
duality is given by the Prokushkin-Vasiliev theory Prokushkin and Vasiliev
(1999) on AdS3, which consists of higher-spin gauge fields and two complex
scalar fields. The higher-spin sector can be described by Chern-Simons gauge
fields $A=A_{\mu}\mathrm{d}x^{\mu},\bar{A}=\bar{A}_{\mu}\mathrm{d}x^{\mu}$.
The gauge fields take values in higher-spin algebra hs$[\hat{\lambda}]$. The
algebra can be defined such as to be truncated to sl$(N^{\prime})$ at
$\hat{\lambda}=N^{\prime}$ with positive integer $N^{\prime}$ even though we
use $0<\hat{\lambda}<1$. The Chern-Simons action is given by 222We use the
definition of Tr for the infinite dimensional algebra given in Vasiliev
(1991); Gaberdiel and Hartman (2011).

$\displaystyle\begin{aligned} &S=S_{\text{CS}}[A]-S_{\text{CS}}[\bar{A}]\,,\\\
&S_{\text{CS}}=\frac{k_{\text{CS}}}{4\pi}\int\text{Tr}\left(A\wedge\mathrm{d}A+\frac{2}{3}A\wedge
A\wedge A\right)\,.\end{aligned}$ (15)

Here $k_{\text{CS}}$ is the level of the Chern-Simons theory, and it relates
to gravity parameters as $k_{\text{CS}}=\ell_{\text{AdS}}/(4G_{N})$, where
$G_{N}$ is the Newton constant. The asymptotic symmetry near the AdS boundary
is found to be a W-algebra Henneaux and Rey (2010); Campoleoni  _et al._
(2010). The central charge is the same as the Brown-Henneaux one Brown and
Henneaux (1986) as

$\displaystyle c=6k_{\text{CS}}=\frac{3\ell_{\text{AdS}}}{2G_{N}}\,.$ (16)

Note that the classical gravity limit with $G_{N}\to 0$ corresponds to the
large central charge limit $c\to\infty$. More precisely speaking, the symmetry
algebra is W${}_{\infty}[\hat{\lambda}]$, which can be truncated to
W${}_{N^{\prime}}$-algebra with spin-$s\,(=2,3,\ldots,N^{\prime})$ currents at
$\hat{\lambda}=N^{\prime}$ Gaberdiel and Gopakumar (2012). The two complex
scalar fields $\phi_{\pm}$ have masses given by
$\ell^{2}_{\text{AdS}}\,m^{2}=\hat{\lambda}^{2}-1$.

The CFT dual to the higher-spin theory on AdS3 was proposed to be the WN-
minimal model at the ’t Hooft limit Gaberdiel and Gopakumar (2011). The WN-
minimal model can be described by a coset CFT (1) with the central charge

$\displaystyle c=(N-1)\left(1-\frac{N(N+1)}{(N+k)(N+k+1)}\right)\,.$ (17)

The classical limit of higher-spin gravity is proposed to be dual to the ’t
Hooft limit, where $N,k\to\infty$ but ’t Hooft parameter $\lambda=N/(N+k)$ is
kept finite. Note that $0<\lambda<1$ for real positive integer $N,k$. The
proposal of Gaberdiel and Gopakumar (2011) is that $\lambda$ is identified
with $\hat{\lambda}$ appearing in the higher-spin theory. In the following, we
only use $\lambda$ to express the parameter.

The dual operators are given by two complex scalar operators
$\mathcal{O}_{\pm}^{\text{AdS}}(z,\bar{z})\,(\equiv\mathcal{O}^{\text{AdS}}_{\pm}(z))$
and holomorphic conserved spin-$s$ currents $J_{(s)}^{\text{AdS}}(z)$ with
$\Delta_{+}=s=2,3,\ldots$ (and anti-holomorphic currents
$\bar{J}_{(s)}^{\text{AdS}}(\bar{z})$). The conformal dimensions of
$\mathcal{O}_{\pm}^{\text{AdS}}(z)$ are $2h_{\pm}=1\pm\lambda$ at the ’t Hooft
limit. The 3-pt. functions of scalar-scalar-higher-spin current are computed
from the bulk Vasiliev equations in Ammon  _et al._ (2012) as

$\displaystyle\begin{aligned}
&\langle\,\mathcal{O}_{\pm}^{\text{AdS}}(z_{1})\,\bar{\mathcal{O}}_{\pm}^{\text{AdS}}(z_{2})\,J_{(s)}^{\text{AdS}}(z_{3})\,\rangle\\\
&=C_{\pm}^{(s)}\left(\frac{z_{12}}{z_{13}z_{23}}\right)^{s}\,\langle\,\mathcal{O}_{\pm}^{\text{AdS}}(z_{1})\,\bar{\mathcal{O}}_{\pm}^{\text{AdS}}(z_{2})\,\rangle\end{aligned}$
(18)

with

$\displaystyle
C_{\pm}^{(s)}=\frac{\eta^{s}_{\pm}}{2\pi}\,\frac{\Gamma(s)^{2}}{\Gamma(2s-1)}\frac{\Gamma(s\pm\lambda)}{\Gamma(1\pm\lambda)}$
(19)

at the leading order in $1/c$. The phase factors can be chosen arbitrary but
here we set $\eta^{s}_{+}=1$ and $\eta^{s}_{-}=(-1)^{s}$ as in Hikida and
Uetoko (2017). The holomorphic higher-spin currents are normalized as

$\displaystyle\begin{aligned}
&\langle\,J_{(s)}^{\text{AdS}}(z)\,J_{(s)}^{\text{AdS}}(0)\,\rangle=\frac{cB^{(s)}}{z^{2s}}\,,\\\
&B^{(s)}=\frac{1}{2^{2s}\pi^{\frac{5}{2}}}\,\frac{\sin(\pi\lambda)}{\lambda(1-\lambda^{2})}\,\frac{\Gamma(s)\,\Gamma(s-\lambda)\,\Gamma(s+\lambda)}{\Gamma\left(s-\frac{1}{2}\right)}\,.\end{aligned}$
(20)

Since the higher-spin theory has one dimension-less parameter
$k_{\text{CS}}\,(=c/6)$, the coupling to scalar fields is also organized by
the same parameter. Thus the canonical normalization of 2-pt. functions of
scalar operators is

$\displaystyle\langle\,\mathcal{O}_{\pm}^{\text{AdS}}(z_{1})\,\bar{\mathcal{O}}_{\pm}^{\text{AdS}}(z_{2})\,\rangle\sim\frac{c}{|z_{12}|^{4h_{\pm}}}$
(21)

up to an overall real factor as shown in Anninos  _et al._ (2017).

Here we remark that the 3-pt. correlators (18) can also be obtained from the
bulk 3-pt. Witten diagram:

$\displaystyle\int\mathrm{d}^{3}pK_{2h_{\pm},0}(z_{1};p)K_{2h_{\pm},0}(z_{2};p)K_{s,s}(z_{3};p)\,,$
(22)

see fig. 1. Here we represent a bulk point by $p=(y,z)$ and a bulk-to-boundary
spin-$s$ propagator with dual dimension $\Delta$ by
$K_{\Delta,s}(z^{\prime};p)$. While the $z_{i}$-dependence was fixed
kinematically, this integral was evaluated in Costa  _et al._ (2014) to
produce the overall dynamical factor, and we can show that it matches with
(19) up to an $s$-independent factor.

Figure 1: Witten diagrams for 3- and 4-pt. functions. The exchange diagrams
may be expressed as products of 3-pt. functions via split representation of
bulk-to-bulk propagators.

##  IV Analytic continuation to dS3/CFT2

We would like to map the bulk fields on AdS3 to those on dS3 applying (7). Due
to the map, the level of Chern-Simons theory in (15) has to be changed as
Witten (1988)

$\displaystyle
k_{\text{CS}}=\frac{\ell_{\text{AdS}}}{4G_{N}}=-i\frac{\ell}{4G_{N}}=-i\kappa$
(23)

with $\kappa\in\mathbb{R}$. The symmetry at late time is supposed to be
$W_{\infty}[\lambda]$-algebra with pure imaginary $c=-i\,6\kappa$ Strominger
(2001); Ouyang (2011), see (16). The masses of two complex scalars
$\phi_{\pm}$ become $\ell^{2}m^{2}=1-\lambda^{2}$. Here we again assign
$0<\lambda<1$ but now for $\ell^{2}m^{2}$ to be positive.

The dual CFT is given by (1) but the central charge (17) is now pure
imaginary. We would like to provide a concrete prescription to perform an
analytic continuation in the coset model. The symmetry underlying the
holography is $W_{\infty}[\lambda]$-algebra, which was shown to be uniquely
determined with two finite parameters $\lambda,c$ Gaberdiel and Gopakumar
(2012); Linshaw (2021). The symmetry of the coset is realized by a special
choice of $\lambda,c$ in terms of positive integer $k,N$. These facts imply
that the correct analytic continuation is keeping $\lambda$ fixed but setting
$c=-ic^{(g)}$ with $c^{(g)}>0$ 333We use the notation $c=-ic^{(g)}$ by
following the previous works Hikida  _et al._ (2022a, b).. Here we compute
coefficient functions in the wave functional of universe via the relation
(13), i.e., in terms of coset CFT (1), which can be also described by Toda
CFT, see Creutzig and Hikida (2021). We then change the parameters from $N,k$
to $c,\lambda$, which can be easily analytically continued.

The CFT has two complex scalar operators
$\mathcal{O}_{\pm}(z,\bar{z})(\equiv\mathcal{O}_{\pm}(z))$ with conformal
dimensions $2h_{\pm}=1\pm\lambda$ at the ’t Hooft limit and conserved higher-
spin currents $J_{(s)}(z)$ (and its anti-holomorphic partner
$\bar{J}_{(s)}(\bar{z})$) with $s=2,3,\ldots$. The relation to bulk fields on
AdS3 can be read off from (9). For the scalar operator, we assign

$\displaystyle\mathcal{O}_{\pm}=-e^{-i\pi
h_{\pm}}\mathcal{O}_{\pm}^{\text{AdS}}\,,\quad\tilde{\mathcal{O}}_{\pm}=-e^{-i\pi
h_{\pm}}\bar{\mathcal{O}}_{\pm}^{\text{AdS}}\,.$ (24)

While $\bar{\mathcal{O}}_{\pm}^{\text{AdS}}$ is the complex conjugate of
$\mathcal{O}_{\pm}^{\text{AdS}}$, the equation (24) imply that
$\mathcal{O}_{\pm}$ and $\tilde{\mathcal{O}}_{\pm}$ are no longer complex
conjugate to each other. For higher-spin current, we have

$\displaystyle J_{(s)}=-e^{-\frac{i}{2}\pi
s}J_{(s)}^{\text{AdS}}\,,\quad\tilde{J}_{(s)}=-e^{-\frac{i}{2}\pi
s}\bar{J}_{(s)}^{\text{AdS}}\,.$ (25)

We then compute bulk correlators at late time from (14). The 2-pt. correlators
can be computed as

$\displaystyle\langle\,\phi_{\pm}(z_{1})\,\tilde{\phi}_{\pm}(z_{2})\,\rangle=-\frac{1}{2\text{Re}\langle\mathcal{O}_{\pm}(z_{1})\tilde{\mathcal{O}}_{\pm}(z_{2})\rangle}\,.$
(26)

Note that only real part appears due to the square of wave functional as in
(14) and $\tilde{\phi}_{\pm}$ are complex conjugate of $\phi_{\pm}$. Here the
functional inverse is defined via (see, e.g., Dolan and Osborn (2011); Osborn
(2012))

$\displaystyle
k_{h,\bar{h}}\frac{1}{\pi}\int\mathrm{d}^{2}y\,\frac{(x-y)^{2h-2}(\bar{x}-\bar{y})^{2\bar{h}-2}}{(y-z)^{2h}(\bar{y}-\bar{z})^{2\bar{h}}}=\delta^{(2)}(x-z)\,,$
$\displaystyle
k_{h,\bar{h}}=\frac{\Gamma(2-2\bar{h})}{\Gamma(2h-1)}=(-1)^{2(h-\bar{h})}\frac{\Gamma(2-2h)}{\Gamma(2\bar{h}-1)}\,.$
(27)

Noticing (21) and (24), we have

$\displaystyle\begin{aligned}
&\langle\,\mathcal{O}_{\pm}(z_{1})\,\tilde{\mathcal{O}}_{\pm}(z_{2})\,\rangle=-i\,e^{-2i\pi
h_{\pm}}\,\frac{1}{|z_{12}|^{4h_{\pm}}}\end{aligned}$ (28)

by changing the overall factors of operators. From the formula (27), we have

$\displaystyle\langle\,\phi_{\pm}(z_{1})\,\tilde{\phi}_{\pm}(z_{2})\,\rangle$
$\displaystyle=-a_{h_{\pm}}\frac{1}{\pi}\frac{\Gamma(2-2h_{\pm})}{\Gamma(2h_{\pm}-1)}\frac{1}{|z_{12}|^{4-4h_{\pm}}\,}$
(29)

with

$\displaystyle a_{h_{\pm}}=\frac{1}{2\sin(2\pi h_{\pm})}\,.$ (30)

The factor $a_{h_{\pm}}$ arises due to the analytic continuation and indeed it
reproduces the corresponding factor $c_{\Delta}^{\text{dS-AdS}}$ with
$\Delta=2h_{\pm}$ given in (2.15) of Sleight and Taronna (2021b). Note that
the factor $c_{\Delta}^{\text{dS-AdS}}$ can be used also for spinning fields
with conformal dimension $\Delta$ and spin $s$.

We then consider the higher-spin gauge field and describe the deviation from
the background values by $\mu^{(s)}$ with $s=2,3,\ldots$. Using the relation
(25), the 2-pt. coefficient functions are written as

$\displaystyle\begin{aligned}
&\langle\,J_{(s)}(z_{1})\,J_{(s)}(z_{2})\,\rangle=(-1)^{s+1}\frac{ic^{(g)}}{c}\,\,\langle\,J_{(s)}^{\text{AdS}}(z_{1})\,J_{(s)}^{\text{AdS}}(z_{2})\,\rangle\,.\end{aligned}$
(31)

Thus the bulk correlators of higher spin fields at late time can be written as

$\displaystyle\begin{aligned}
&\langle\mu_{(s)}(z_{1})\mu_{(s)}(z_{2})\rangle\\\
&\quad=a_{(s)}\frac{1}{\pi\Gamma(2s-1)c^{(g)}B^{(s)}}\frac{1}{(z_{12})^{2-2s}(\bar{z}_{12})^{2}}\,.\end{aligned}$
(32)

The factor

$\displaystyle a_{(s)}=(-1)^{s+1}\frac{i}{2}$ (33)

arises due to the analytic continuation as in the scalar case. However, the
phase factor (33) does not match with $c_{\Delta}^{\text{dS-AdS}}$ for a
massive spin-$s$ field in Sleight and Taronna (2021b) as it diverges for
higher-spin gauge fields with $\Delta=s$. This implies that the massless limit
of higher-spin field is quite subtle. Let us see an example of massive spin-1
field on Lorentzian dSd+1 space-time. Its bulk two-point function was computed
in (3.18)-(3.20) of Allen and Jacobson (1986), but it diverges in the massless
limit $m\to 0$. In the limit, gauge symmetry appears and analysis has to be
redone from the beginning as in section 4 of Allen and Jacobson (1986). In
Sleight and Taronna (2021b), a factor was inserted such that bulk massive
spin-$s$ propagator on dSd+1 has correct behavior near short distance. To
obtain (33) by their method, we have to work directly with massless higher-
spin gauge field or take a massless limit with a special care.

We then examine 3-pt. correlators at late time. Using (24) and (25), the 3-pt.
coefficient functions are related to (18) as

$\displaystyle\begin{aligned}
&\langle\,\mathcal{O}_{\pm}(z_{1})\,\tilde{\mathcal{O}}_{\pm}(z_{2})\,J_{(s)}(z_{3})\,\rangle\\\
&\quad=i\,e^{-i\frac{\pi}{2}(4h_{\pm}+s)}C_{\pm}^{(s)}\left(\frac{z_{12}}{z_{13}z_{23}}\right)^{s}\frac{1}{|z_{12}|^{4h_{\pm}}}\end{aligned}$
(34)

and similarly for
$\langle\mathcal{O}_{\pm}\tilde{\mathcal{O}}_{\pm}\tilde{J}_{(s)}\rangle$. The
sum of (34) and the conjugate of
$\langle\mathcal{O}_{\pm}\tilde{\mathcal{O}}_{\pm}\tilde{J}_{(s)}\rangle$
becomes

$\displaystyle-2\sin\left((4h_{\pm}+s)\frac{\pi}{2}\right)C_{\pm}^{(s)}\left(\frac{z_{12}}{z_{13}z_{23}}\right)^{s}\frac{1}{|z_{12}|^{4h_{\pm}}}\,.$
(35)

In order to read off the bulk 3-pt. correlators, we need to multiply two of
(29) and one of (32) and integrate over the positions, see (50). We then
obtain

$\displaystyle\begin{aligned}
&\langle\,\phi_{\pm}(z_{1})\,\tilde{\phi}_{\pm}(z_{2})\,\mu_{(s)}(z_{3})\rangle\\\
&\quad=-\lambda_{h_{\pm},h_{\pm},s}\frac{\Gamma(\mp\lambda)\Gamma(1\mp\lambda)}{\Gamma(s)^{2}\Gamma(1-s\pm\lambda)\Gamma(s\pm\lambda)}\\\
&\quad\times\frac{C_{\pm}^{(s)}}{c^{(g)}B^{(s)}}\left(\frac{z_{12}}{z_{13}z_{23}}\right)^{1-s}\left(\frac{\bar{z}_{12}}{\bar{z}_{13}\bar{z}_{23}}\right)\frac{1}{|z_{12}|^{4h_{\mp}}}\end{aligned}$
(36)

with

$\displaystyle\lambda_{h_{\pm},h_{\pm},s}=2a_{h_{\pm}}^{2}a_{(s)}\sin\left((4h_{\pm}+s)\frac{\pi}{2}\right)\,.$
(37)

The same factor appears in (3.24) of Sleight and Taronna (2021b) up to the
subtlety associated to the massless limit mentioned above.

##  V Four-point correlators

In the coset model (1), the 4-pt. functions of scalar operators were computed
in Papadodimas and Raju (2012) exactly in all orders of $N,k$. In this letter,
we analyze a simple one given by

$\displaystyle\begin{aligned}
G_{-+}(z)&=\langle\,\mathcal{O}^{\text{AdS}}_{-}(\infty)\,\mathcal{O}^{\text{AdS}}_{+}(1)\,\bar{\mathcal{O}}^{\text{AdS}}_{+}(z)\,\bar{\mathcal{O}}^{\text{AdS}}_{-}(0)\,\rangle\\\
&=|1-z|^{-2\Delta_{+}}|z|^{\frac{2}{N}}\left|1+\frac{1-z}{Nz}\right|^{2}\,.\end{aligned}$
(38)

Following the analytic continuation procedures explained above, we should be
able to obtain 4-pt. bulk correlator

$\displaystyle\langle\,\phi_{-}(\infty)\,\phi_{+}(1)\,\tilde{\phi}_{+}(z)\,\tilde{\phi}_{-}(0)\,\rangle\,$
(39)

in dS3 at late time.

The expansion of $G_{-+}(z)$ in (38) in terms of global conformal blocks was
given in Hikida and Uetoko (2017) as

$\displaystyle
G_{-+}(z)=|1-z|^{-2{\left(1+\lambda-\frac{1-\lambda^{2}}{c}\right)}}+|1-z|^{-2{(1+\lambda)}}$
(40)
$\displaystyle\times\left[\sum_{s=2}^{\infty}(-1)^{s}\frac{C^{(s)}_{-}C^{(s)}_{+}}{cB^{(s)}}\,(1-z)^{s}{}_{2}F_{1}\left(s,s;2s;1-z\right)+\text{c.c.}\right]$

up to the order $1/c^{2}$. In particular, the conformal dimension of
$\mathcal{O}_{+}$ can be expanded as
$\Delta_{+}=1+\lambda-(1-\lambda^{2})/c+\mathcal{O}(c^{-2})$. The expansion
now can be written as a sum over holomorphic and anti-holomorphic conserved
spin-$s$ exchange with $s=2,3,\ldots,\infty$ and the coefficients are given as
products of 3-pt. functions. In terms of Witten diagram, the 4-pt. function
should be computed from

$\displaystyle\int\mathrm{d}^{3}p_{1}\mathrm{d}^{3}p_{2}K_{2h_{-},0}(\infty;p_{1})K_{2h_{+},0}(0;p_{1})$
(41) $\displaystyle\quad\times
G_{s,s}(p_{1};p_{2})K_{2h_{+},0}(z;p_{2})K_{2h_{-},0}(1;p_{2})\,,$

see fig. 1. Here $G_{\Delta,s}(p_{1};p_{2})$ represents the bulk-to-bulk
spin-$s$ propagator with dual dimension $\Delta$. In the split representation,
it can be expressed as

$\displaystyle
G_{\Delta,s}(p_{1};p_{2})=\int_{-\infty}^{\infty}\frac{\mathrm{d}\nu}{\nu^{2}+(\Delta-1)^{2}}\Omega_{\nu,s}(p_{1};p_{2})\,,$
(42)
$\displaystyle\Omega_{\nu,s}(p_{1};p_{2})=\frac{\nu^{2}}{\pi}\int\mathrm{d}^{2}zK_{1+i\nu,s}(p_{1},z)K_{1-i\nu,s}(p_{2},z)$

up to contact term contributions Costa  _et al._ (2014). Picking up poles in
(42), the 4-pt. function becomes a product of 3-pt. functions, and the
integration over the boundary coordinates leads to a conformal partial wave

$\displaystyle\mathcal{I}_{s,0}(z)=z^{s}{}_{2}F_{1}(s,s;2s;z)+\frac{\Gamma(2s-1)\Gamma(2s)}{\Gamma(s)^{4}}$
$\displaystyle\times
z^{1-s}{}_{2}F_{1}(1-s,1-s;2-2s;z)\bar{z}{}_{2}F_{1}(1,1;2;\bar{z})$ (43)

up to an unimportant overall constant. See (57) and Simmons-Duffin (2014) for
details. In terms of conformal partial waves, the 4-pt. function can be
expanded as:

$\displaystyle
G_{-+}(z)=|1-z|^{-2{\left(1+\lambda-\frac{1-\lambda^{2}}{c}\right)}}+|1-z|^{-2(1+\lambda)}$
(44)
$\displaystyle\times\left[\sum_{s=2}^{\infty}(-1)^{s}\frac{C^{(s)}_{-}C^{(s)}_{+}}{cB^{(s)}}\,\mathcal{I}_{s,0}(1-z)-\frac{1-\lambda^{2}}{c}\ln\bar{z}+\text{c.c.}\right]\,.$

As in (44), there are three kinds of terms in the expansions of 4-pt. function
(38), each has a natural interpretation in the AdS and dS bulks. The first
term is just the product of 2-pt. functions, corresponding to disconnected
Witten diagrams. The second term is the sum of conformal partial waves, which
can be expressed in lower-pt. coefficient functions as mentioned above. When
we perform the path integral over $\phi_{\pm},\bar{\phi}_{\pm}$, they will be
over counted if we include them in the 4-pt. coefficient function. Thus, it is
convenient to remove it beforehand. The third term corresponds to the contact
4-pt. interaction in the bulk. Analytically continuing from AdS to dS spaces,
the connected parts of their respective 4-pt. CFT correlators can be related
via:

$\displaystyle\begin{aligned}
&\langle\,\mathcal{O}_{-}(\infty)\,\mathcal{O}_{+}(1)\,\tilde{\mathcal{O}}_{+}(z)\,\tilde{\mathcal{O}}_{-}(0)\,\rangle_{\text{c}}\\\
&=-i\,\langle\,\mathcal{O}^{\text{AdS}}_{-}(\infty)\,\mathcal{O}^{\text{AdS}}_{+}(1)\,\bar{\mathcal{O}}^{\text{AdS}}_{+}(z)\,\bar{\mathcal{O}}^{\text{AdS}}_{-}(0)\,\rangle_{\text{
c }}\,,\end{aligned}$ (45)

where the subscript “c” indicates the connected part. Notice that while using
(24) yields a phase factor $e^{-i\pi(2h_{+}+2h_{-})}=1$, the overall $-i$
comes from the fact that the connected part is proportional to $c$ if we use
the normalization (21) 444We can show that any connected $n$-pt. function is
proportional to $c$ with the normalization (21) by repeating operator product
expansions to reduce it to products of 3-pt. functions. This is consistent
with the bulk perturbation by $k_{\text{CS}}$.. The contribution from the
third term vanishes in this case after taking its real part.

Thus, the only non-trivial contribution comes from that corresponding to the
dS bulk exchange diagrams, which can be written as

$\displaystyle\langle\,\phi_{-}(\infty)\,\phi_{+}(1)\,\tilde{\phi}_{+}(z)\,\tilde{\phi}_{-}(0)\,\rangle_{\text{c}}=\sum_{s=2}^{\infty}\frac{\lambda_{h_{+},h_{+},s}\lambda_{h_{-},h_{-},s}}{a_{(s)}c^{(g)}}$
$\displaystyle\times\frac{\Gamma(-\lambda)\Gamma(1-\lambda)\Gamma(\lambda)\Gamma(1+\lambda)}{\Gamma(1-s+\lambda)\Gamma(s+\lambda)\Gamma(1-s-\lambda)\Gamma(s-\lambda)}$
(46)
$\displaystyle\times\frac{1}{|1-z|^{2(1-\lambda)}}\left[(-1)^{s}\frac{C^{(s)}_{-}C^{(s)}_{+}}{B^{(s)}}\,\mathcal{I}_{s,0}(1-z)+\text{c.c.}\right]\,.$

The factor $\lambda_{h_{\pm},h_{\pm},s}$ comes from that in (LABEL:bulk3pt)
and the division by $a_{(s)}$ eliminates the over counting of (32). The result
is consistent with the generic expression, say, in (4.39) of Sleight and
Taronna (2021b).

##  VI Discussion

In this letter, we computed bulk dS3 correlators at late time by developing
holographic method and dS3/CFT2 correspondence. The expressions are consistent
with the previous analysis of Sleight and Taronna (2021a, b) based on bulk
Feynman diagrams in the in-in formulation, whose review may be found in
Weinberg (2005). Let us end by briefly commenting these two complementary
approaches exploring the dS/CFT holography and where our results fit in here.

Notice that, in the in-in formulation, there are two types of interaction
vertices with time-ordering and anti-time ordering and the corresponding
propagators connecting them in the dS bulk. This is related to consider both
signs for analytic continuation instead of our prescription (7), see Sleight
and Taronna (2021a, b). Moreover, as in the case of AdS computations, bulk
4-pt. dS correlators can reduce to the evaluation of the product of two 3-pt.
functions by applying a formula analogous to (42), see, e.g., Sleight and
Taronna (2021a, b). However in such computations, further integrals over
spectral parameter $\nu$ are needed to make full comparisons with CFT
correlators, which are usually very difficult.

In holographic method, wave functional and its complex conjugation are
involved as in (14) and the integration over boundary fields provides
connections between them. We can see how the two methods complement each other
from our explicit calculations. As an advantage of holographic method, there
are no integrals over spectral parameter $\nu$ in (14), and the difficulty is
avoided by directly evaluating the coefficient functions via dual CFT. In our
explicit example, the non-trivial information is included in the the third
term of (44). The contribution vanishes at the end of computation, but the
corresponding term survives for different 4-pt. correlators obtained from
$\langle\mathcal{O}_{\pm}^{\text{AdS}}\bar{\mathcal{O}}_{\pm}^{\text{AdS}}\mathcal{O}_{\pm}^{\text{AdS}}\bar{\mathcal{O}}_{\pm}^{\text{AdS}}\rangle$.
See Papadodimas and Raju (2012); Hikida and Uetoko (2017) for their exact
forms and conformal block expansions. In addition, coefficient functions can
be obtained even with finite $\lambda,c$ from the explicit dual CFT. For these
reasons, holographic method would enable us to evaluate more complicated bulk
correlators at loop levels, on dS3 cosmological backgrounds Krishnan  _et al._
(2014), and so on. We are planing to report on more details on the relation
between two methods and further analysis of bulk correlators in near future.

Acknowledgements We are grateful to Thomas Creutzig, Kenta Suzuki, Tadashi
Takayanagi, Yusuke Taki, Seiji Terashima and Takahiro Uetoko for useful
comments. We would like to particularly thank Tatsuma Nishioka for the
collaboration at the early stage and other related work. The work of H. Y. C.
is supported in part by Ministry of Science and Technology (MOST) through the
grant 110-2112-M-002-006-. This work of Y. H. was supported by JSPS Grant-in-
Aid for Scientific Research (B) No. 19H01896, Grant-in-Aid for Scientific
Research (A) No. 21H04469, and Grant-in-Aid for Transformative Research Areas
(A) “Extreme Universe” No. 21H05187.

## References

  * Maldacena (2003) Juan Martin Maldacena, “Non-Gaussian features of primordial fluctuations in single field inflationary models,” JHEP 05, 013 (2003), arXiv:astro-ph/0210603 .
  * Maldacena and Pimentel (2011) Juan M. Maldacena and Guilherme L. Pimentel, “On graviton non-Gaussianities during inflation,” JHEP 09, 045 (2011), arXiv:1104.2846 [hep-th] .
  * Arkani-Hamed and Maldacena (2015) Nima Arkani-Hamed and Juan Maldacena, “Cosmological collider physics,” (2015), arXiv:1503.08043 [hep-th] .
  * Arkani-Hamed  _et al._ (2020) Nima Arkani-Hamed, Daniel Baumann, Hayden Lee, and Guilherme L. Pimentel, “The cosmological bootstrap: Inflationary correlators from symmetries and singularities,” JHEP 04, 105 (2020), arXiv:1811.00024 [hep-th] .
  * Strominger (2001) Andrew Strominger, “The dS/CFT correspondence,” JHEP 10, 034 (2001), arXiv:hep-th/0106113 .
  * Witten (2001a) Edward Witten, “Quantum gravity in de Sitter space,” in  _Strings 2001: International Conference_ (2001) arXiv:hep-th/0106109 .
  * Maldacena (1998) Juan Martin Maldacena, “The Large $N$ limit of superconformal field theories and supergravity,” Adv. Theor. Math. Phys. 2, 231–252 (1998), arXiv:hep-th/9711200 .
  * Anninos  _et al._ (2017) Dionysios Anninos, Thomas Hartman, and Andrew Strominger, “Higher spin realization of the dS/CFT correspondence,” Class. Quant. Grav. 34, 015009 (2017), arXiv:1108.5735 [hep-th] .
  * Klebanov and Polyakov (2002) I. R. Klebanov and A. M. Polyakov, “AdS dual of the critical O$(N)$ vector model,” Phys. Lett. B 550, 213–219 (2002), arXiv:hep-th/0210114 .
  * Sundborg (2001) Bo Sundborg, “Stringy gravity, interacting tensionless strings and massless higher spins,” Nucl. Phys. B Proc. Suppl. 102, 113–119 (2001), arXiv:hep-th/0103247 .
  * Mikhailov (2002) Andrei Mikhailov, “Notes on higher spin symmetries,” (2002), arXiv:hep-th/0201019 .
  * Witten (2001b) Edward Witten, “Spacetime reconstruction,”  (2001b), the John Schwarz 60-th Birthday Symposium, California Institute of Technology, U.S.A., 3-4 Nov 2001\.
  * Sezgin and Sundell (2002) E. Sezgin and P. Sundell, “Massless higher spins and holography,” Nucl. Phys. B 644, 303–370 (2002), [Erratum: Nucl.Phys.B 660, 403–403 (2003)], arXiv:hep-th/0205131 .
  * Hikida  _et al._ (2022a) Yasuaki Hikida, Tatsuma Nishioka, Tadashi Takayanagi, and Yusuke Taki, “Holography in de Sitter Space via Chern-Simons Gauge Theory,” Phys. Rev. Lett. 129, 041601 (2022a), arXiv:2110.03197 [hep-th] .
  * Hikida  _et al._ (2022b) Yasuaki Hikida, Tatsuma Nishioka, Tadashi Takayanagi, and Yusuke Taki, “CFT duals of three-dimensional de Sitter gravity,” JHEP 05, 129 (2022b), arXiv:2203.02852 [hep-th] .
  * Ouyang (2011) Peter Ouyang, “Toward higher spin dS3/CFT2,” (2011), arXiv:1111.0276 [hep-th] .
  * Sleight and Taronna (2021a) Charlotte Sleight and Massimo Taronna, “From AdS to dS exchanges: Spectral representation, Mellin amplitudes and crossing,” Phys. Rev. D 104, L081902 (2021a), arXiv:2007.09993 [hep-th] .
  * Sleight and Taronna (2021b) Charlotte Sleight and Massimo Taronna, “From dS to AdS and back,” JHEP 12, 074 (2021b), arXiv:2109.02725 [hep-th] .
  * Shukla  _et al._ (2016) Ashish Shukla, Sandip P. Trivedi, and V. Vishal, “Symmetry constraints in inflation, $\alpha$-vacua, and the three point function,” JHEP 12, 102 (2016), arXiv:1607.08636 [hep-th] .
  * Di Pietro  _et al._ (2022) Lorenzo Di Pietro, Victor Gorbenko, and Shota Komatsu, “Analyticity and unitarity for cosmological correlators,” JHEP 03, 023 (2022), arXiv:2108.01695 [hep-th] .
  * Gaberdiel and Gopakumar (2011) Matthias R. Gaberdiel and Rajesh Gopakumar, “An AdS3 dual for minimal model CFTs,” Phys. Rev. D 83, 066007 (2011), arXiv:1011.2986 [hep-th] .
  * Prokushkin and Vasiliev (1999) S. F. Prokushkin and Mikhail A. Vasiliev, “Higher spin gauge interactions for massive matter fields in 3-D AdS space-time,” Nucl. Phys. B 545, 385 (1999), arXiv:hep-th/9806236 .
  * Arakawa  _et al._ (2019) Tomoyuki Arakawa, Thomas Creutzig, and Andrew R. Linshaw, “W-algebras as coset vertex algebras,” Invent. Math. 218, 145–195 (2019), arXiv:1801.03822 [math.QA] .
  * Note (1) In Hikida  _et al._ (2022a, b) a limit of $k\to-N$ is considered, while in this letter an ’t Hooft limit is taken with large $|N|,|k|$ but finite $\lambda=N/(N+k)$.
  * Note (2) We use the definition of Tr for the infinite dimensional algebra given in Vasiliev (1991); Gaberdiel and Hartman (2011).
  * Vasiliev (1991) Mikhail A. Vasiliev, “Higher spin algebras and quantization on the sphere and hyperboloid,” Int. J. Mod. Phys. A 06, 1115–1135 (1991).
  * Gaberdiel and Hartman (2011) Matthias R. Gaberdiel and Thomas Hartman, “Symmetries of holographic minimal models,” JHEP 05, 031 (2011), arXiv:1101.2910 [hep-th] .
  * Henneaux and Rey (2010) Marc Henneaux and Soo-Jong Rey, “Nonlinear $W_{\infty}$ as asymptotic symmetry of three-dimensional higher spin anti-de Sitter gravity,” JHEP 12, 007 (2010), arXiv:1008.4579 [hep-th] .
  * Campoleoni  _et al._ (2010) Andrea Campoleoni, Stefan Fredenhagen, Stefan Pfenninger, and Stefan Theisen, “Asymptotic symmetries of three-dimensional gravity coupled to higher-spin fields,” JHEP 11, 007 (2010), arXiv:1008.4744 [hep-th] .
  * Brown and Henneaux (1986) J. David Brown and M. Henneaux, “Central charges in the canonical realization of asymptotic symmetries: An example from three-dimensional gravity,” Commun. Math. Phys. 104, 207–226 (1986).
  * Gaberdiel and Gopakumar (2012) Matthias R. Gaberdiel and Rajesh Gopakumar, “Triality in minimal model holography,” JHEP 07, 127 (2012), arXiv:1205.2472 [hep-th] .
  * Ammon  _et al._ (2012) Martin Ammon, Per Kraus, and Eric Perlmutter, “Scalar fields and three-point functions in $D=3$ higher spin gravity,” JHEP 07, 113 (2012), arXiv:1111.3926 [hep-th] .
  * Hikida and Uetoko (2017) Yasuaki Hikida and Takahiro Uetoko, “Three point functions in higher spin AdS3 holography with $1/N$ corrections,” Universe 3, 70 (2017), arXiv:1708.02017 [hep-th] .
  * Costa  _et al._ (2014) Miguel S. Costa, Vasco Gonçalves, and João Penedones, “Spinning AdS propagators,” JHEP 09, 064 (2014), arXiv:1404.5625 [hep-th] .
  * Witten (1988) Edward Witten, “(2+1)-dimensional gravity as an exactly soluble system,” Nucl. Phys. B 311, 46 (1988).
  * Linshaw (2021) Andrew R. Linshaw, “Universal two-parameter $\mathcal{W}_{\infty}$-algebra and vertex algebras of type $\mathcal{W}(2,3,\dots,N)$,” Compos. Math. 157, 12–82 (2021), arXiv:1710.02275 [math.RT] .
  * Note (3) We use the notation $c=-ic^{(g)}$ by following the previous works Hikida  _et al._ (2022a, b).
  * Creutzig and Hikida (2021) Thomas Creutzig and Yasuaki Hikida, “Correlator correspondences for Gaiotto-Rapčák dualities and first order formulation of coset models,” JHEP 12, 144 (2021), arXiv:2109.03403 [hep-th] .
  * Dolan and Osborn (2011) F. A. Dolan and H. Osborn, “Conformal partial waves: Further mathematical results,” (2011), arXiv:1108.6194 [hep-th] .
  * Osborn (2012) H. Osborn, “Conformal blocks for arbitrary spins in two dimensions,” Phys. Lett. B 718, 169–172 (2012), arXiv:1205.1941 [hep-th] .
  * Allen and Jacobson (1986) Bruce Allen and Theodore Jacobson, “Vector two point functions in maximally symmetric spaces,” Commun. Math. Phys. 103, 669 (1986).
  * Papadodimas and Raju (2012) Kyriakos Papadodimas and Suvrat Raju, “Correlation functions in holographic minimal models,” Nucl. Phys. B 856, 607–646 (2012), arXiv:1108.3077 [hep-th] .
  * Simmons-Duffin (2014) David Simmons-Duffin, “Projectors, shadows, and conformal blocks,” JHEP 04, 146 (2014), arXiv:1204.3894 [hep-th] .
  * Note (4) We can show that any connected $n$-pt. function is proportional to $c$ with the normalization (21) by repeating operator product expansions to reduce it to products of 3-pt. functions. This is consistent with the bulk perturbation by $k_{\text{CS}}$.
  * Weinberg (2005) Steven Weinberg, “Quantum contributions to cosmological correlations,” Phys. Rev. D 72, 043514 (2005), arXiv:hep-th/0506236 .
  * Krishnan  _et al._ (2014) Chethan Krishnan, Avinash Raju, Shubho Roy, and Somyadip Thakur, “Higher spin cosmology,” Phys. Rev. D 89, 045007 (2014), arXiv:1308.6741 [hep-th] .

##  Appendix A Appendix A: Bulk correlators from wave functional

In this appendix, we collect formulas to compute bulk correlators from wave
functional of universe. Expanding the wave functional (13) by the fields
$\psi_{j}$ at the late time, we may write

$\displaystyle\Psi[\psi_{l}]$ (47) $\displaystyle=\exp\left[\sum_{m\geq
2}\int\mathrm{d}^{d}\vec{x}_{1}\cdots\mathrm{d}^{d}\vec{x}_{m}C_{m}(\\{\vec{x}_{l}\\})\prod_{j=1}^{m}\psi_{j}(\vec{x}_{j})\right]\,.$

The proposal here is that the coefficient functions $C_{m}(\\{\vec{x}_{l}\\})$
can be computed by certain dual conformal field theory.

The explicit forms of correlation functions in terms of CFT correlators are
given as follows (see, e.g., Maldacena (2003); Maldacena and Pimentel (2011);
Arkani-Hamed  _et al._ (2020) for their momentum space expressions). The 2-pt.
functions are written as

$\displaystyle\begin{aligned} \langle\,\sigma_{i_{1}\cdots
i_{s}}(\vec{x}_{1})\,\sigma^{j_{1}\cdots
j_{s}}(\vec{x}_{2})\,\rangle=-\frac{\Pi^{j_{1}\cdots j_{s}}_{i_{1}\cdots
i_{s}}(\hat{x}_{12})}{2\,\text{Re}\,\langle\,J_{s}(\vec{x}_{1})\,J_{s}(\vec{x}_{2})\,\rangle}\,,\end{aligned}$
(48)

where we have used the expression

$\displaystyle\langle\,J_{i_{1}\cdots i_{s}}(\vec{x}_{1})\,J^{j_{1}\cdots
j_{s}}(\vec{x}_{2})\,\rangle=\Pi^{j_{1}\cdots j_{s}}_{i_{1}\cdots
i_{s}}(\hat{x}_{12})\,\langle\,J_{s}(\vec{x}_{1})\,J_{s}(\vec{x}_{2})\,\rangle$
(49)

with the unit vector $\hat{x}=\vec{x}/|\vec{x}|$. The 3-pt. functions used in
the main context are

$\displaystyle\langle\,\phi_{1}(\vec{x}_{1})\,\phi_{2}(\vec{x}_{2})\,\sigma_{i_{1}\cdots
i_{s}}(\vec{x}_{3})\,\rangle$ (50) $\displaystyle\
=-\int\prod_{l=1}^{3}\mathrm{d}^{d}\vec{x}_{l}^{\prime}\text{Re}\,\langle\,\mathcal{O}_{1}(\vec{x}_{1}^{\prime})\,\mathcal{O}_{2}(\vec{x}_{2}^{\prime})\,J_{j_{1}\cdots
j_{s}}(\vec{x}_{3}^{\prime})\,\rangle$
$\displaystyle\times\frac{\Pi^{j_{1}\cdots j_{s}}_{i_{1}\cdots
i_{s}}(\vec{x}_{33^{\prime}})}{4\left[\prod_{i=1}^{2}\text{Re}\,\langle\,\mathcal{O}_{i}(\vec{x}_{i})\,\mathcal{O}_{i}(\vec{x}_{i}^{\prime})\,\rangle\right]\text{Re}\,\langle\,J_{s}(\vec{x}_{3})\,J_{s}(\vec{x}_{3}^{\prime})\,\rangle}\
.$

Here we set for $s=0$ as $\phi=\sigma_{i_{1}\cdots i_{s}}$ and
$\mathcal{O}=J^{i_{1}\cdots i_{s}}$. The scalar 4-pt. function can be
decomposed as

$\displaystyle\begin{aligned}
&\langle\,\phi_{1}(\vec{x}_{1})\,\phi_{2}(\vec{x}_{2})\,\phi_{3}(\vec{x}_{3})\,\phi_{4}(\vec{x}_{4})\,\rangle\\\
&=\int\prod_{l=1}^{4}\mathrm{d}^{d}\vec{x}_{l}^{\prime}\frac{\langle\mathcal{O}^{4}\rangle_{\text{c}}+\langle\mathcal{O}^{4}\rangle_{\text{d}}}{\prod_{i=1}^{4}(-2\,\text{Re}\,\langle\,\mathcal{O}_{i}(\vec{x}_{i})\,\mathcal{O}_{i}(\vec{x}_{i}^{\prime})\,\rangle)}\end{aligned}$
(51)

with

$\displaystyle\langle\,\mathcal{O}^{4}\,\rangle_{\text{c}}=2\,\text{Re}\,\langle\,\mathcal{O}_{1}(\vec{x}_{1}^{\prime})\,\mathcal{O}_{2}(\vec{x}_{2}^{\prime})\,\mathcal{O}_{3}(\vec{x}_{3}^{\prime})\,\mathcal{O}_{4}(\vec{x}_{4}^{\prime})\,\rangle\,,$
(52) $\displaystyle\langle\,\mathcal{O}^{4}\,\rangle_{\text{d}}$
$\displaystyle=\sum_{s}\int\mathrm{d}^{d}\vec{y}_{1}\mathrm{d}^{d}\vec{y}_{2}2\,\text{Re}\,\langle\,\mathcal{O}_{1}(\vec{x}_{1}^{\prime})\,\mathcal{O}_{2}(\vec{x}_{2}^{\prime}),J_{j_{1}\cdots
j_{s}}(\vec{y}_{1})\,\rangle$ $\displaystyle\times\frac{\Pi^{j_{1}\cdots
j_{s}}_{i_{1}\cdots i_{s}}(\vec{y}_{12})\text{Re}\,\langle\,J^{i_{1}\cdots
i_{s}}(\vec{y}_{2})\,\mathcal{O}_{3}(\vec{x}_{3}^{\prime})\,\mathcal{O}_{4}(\vec{x}_{4}^{\prime})\,\rangle}{\text{Re}\,\langle\,J_{s}(\vec{y}_{1})\,J_{s}(\vec{y}_{2})\,\rangle}\,.$
(53)

For $\langle\,\mathcal{O}^{4}\,\rangle_{d}$, we have to sum up all
intermediate operators.

##  Appendix B Appendix B: Conformal integrals

In order to extract information of bulk correlators from the coefficient
functions of wave functional, we need to perform some conformal integrals.
Here we summarize some formulas, see, e.g., Dolan and Osborn (2011); Osborn
(2012).

For 2-pt. correlators, we have used the formula (27). For 3-pt. correlators,
it is convenient to apply the formula

$\displaystyle\begin{aligned}
&k_{h,\bar{h}}\frac{1}{\pi}\int\mathrm{d}^{2}y\frac{1}{(z-y)^{2-2h}(\bar{z}-\bar{y})^{2-2\bar{h}}}\mathcal{F}^{h,\bar{h}}_{12}(y,\bar{y})\\\
&=\frac{\Gamma(1-h-h_{12})\Gamma(1-\bar{h}+\bar{h}_{12})}{\Gamma(\bar{h}+\bar{h}_{12})\Gamma(h-h_{12})}\mathcal{F}^{1-h,1-\bar{h}}_{12}(y,\bar{y})\end{aligned}$
(54)

with

$\displaystyle\begin{aligned}
&\mathcal{F}^{h,\bar{h}}_{12}(z,\bar{z})=\frac{1}{z_{12}^{h_{1}+h_{2}-h}(z_{1}-z)^{h+h_{12}}(z_{2}-z)^{h-h_{12}}}\\\
&\quad\times\frac{1}{\bar{z}_{12}^{\bar{h}_{1}+\bar{h}_{2}-\bar{h}}(\bar{z}_{1}-\bar{z})^{\bar{h}+\bar{h}_{12}}(\bar{z}_{2}-\bar{z})^{\bar{h}-\bar{h}_{12}}}\,.\end{aligned}$
(55)

Here we have used the notation such as $z_{12}=z_{1}-z_{2}$,
$h_{12}=h_{1}-h_{2}$ and so on.

For 4-pt. correlator, we expand it in terms of conformal partial waves as in
(44). We define a shadow operator $\hat{\Phi}_{1-h,1-\bar{h}}$ by

$\displaystyle\begin{aligned} &\hat{\Phi}_{1-h,1-\bar{h}}(x)\\\
&=k_{h,\bar{h}}\frac{1}{\pi}\int\mathrm{d}^{2}y\,\frac{1}{(x-y)^{2-2h}(\bar{x}-\bar{y})^{2-2\bar{h}}}\,\Phi_{h,\bar{h}}(y)\,,\end{aligned}$
(56)

where $k_{h,\bar{h}}$ was defined in (27). The conformal partial wave of
spin-$s$ exchange can be obtained from integrals

$\displaystyle\frac{k_{1-s,1}k_{s,0}}{\pi^{2}cB^{(s)}}\int\mathrm{d}^{2}x\,\mathrm{d}^{2}y\,\langle\,\mathcal{O}_{-}^{\text{AdS}}(\infty)\,\bar{\mathcal{O}}_{-}^{\text{AdS}}(0)\,J_{(s)}^{\text{AdS}}(x)\,\rangle\,$
$\displaystyle\times\frac{1}{(x-y)^{2-2s}(\bar{x}-\bar{y})^{2}}\,\langle\,J_{(s)}^{\text{AdS}}(y)\,\mathcal{O}_{+}^{\text{AdS}}(1)\,\bar{\mathcal{O}}_{+}^{\text{AdS}}(z)\,\rangle$
$\displaystyle=\frac{\Gamma(2s)}{\pi
cB^{(s)}}\int\mathrm{d}^{2}x\,\langle\,\mathcal{O}_{-}^{\text{AdS}}(\infty)\,\bar{\mathcal{O}}_{-}^{\text{AdS}}(0)\,J_{(s)}^{\text{AdS}}(x)\,\rangle\,$
$\displaystyle\times\,\langle\,\hat{J}_{(s)}^{\text{AdS}}(x)\,\mathcal{O}_{+}^{\text{AdS}}(1)\,\bar{\mathcal{O}}_{+}^{\text{AdS}}(z)\,\rangle$
(57)
$\displaystyle=(-1)^{s}\frac{1}{|1-z|^{2\Delta_{+}}}\frac{C_{-}^{(s)}C_{+}^{(s)}}{cB^{(s)}}\mathcal{I}_{s,0}(1-z)\,,$

where $\mathcal{I}_{s,0}(z)$ is given in (43), see, e.g., Simmons-Duffin
(2014).


# Dark matter halo mass density profiles around an accreting supermassive
blackhole in the Schrödinger-Newton approach

Chon-Fai Kam  Email<EMAIL_ADDRESS>Department of Physics and The Hong
Kong Institute of Quantum Information Science and Technology, The Chinese
University of Hong Kong, Shatin, New Territories, Hong Kong, China  Iat-Neng
Chan  Email<EMAIL_ADDRESS>Department of Physics and Chemistry, Faculty of
Science and Technology, University of Macau, Avenida da Universidade, Taipa,
Macau, People’s Republic of China

###### Abstract

We study how supermassive blackhole accretion affects the mass density profile
of a fuzzy dark matter soliton core at the centre of a dark matter halo. The
supermassive blackhole at the center of a galaxy was regarded as a point mass.
The Schrödinger-Newton equation for the scalar field was solved numerically.
We find that the time-dependent perturbation has a significant squeezing
effect on the soliton density profile, which both reduces the size of the core
and increases the central density.

##  I Introduction

Despite a strong circumstantial evidence for its existence, the exact particle
nature of dark matter, which makes up 27% of the Universe’s energy density, is
still yet unknown and remains one of the biggest mysteries in physics
bertone2018new . As dark matter is approximately five times more abundant than
baryonic matter such as neutrons and protons, which adds up to less than 5% of
the energy density in the universe, it dominates the formation of large-scale
structure of the Universe frenk2012dark , and governs the formation and
evolution of galaxies within merging dark matter haloes
vogelsberger2020cosmological ; de2020dark . Although the nature of dark matter
is still unknown, there is overwhelming evidence that dark matter cannot be
composed of any particles described by the standard model of particle physics,
based on a variety of observations bertone2018new . Among the various possible
solutions to the dark matter problem such as neutrinos, light bosons, or
modified gravity, the cold dark matter (CDM) model with weakly interacting
massive particles (WIMPs) seems the most promising one, since it can
successfully simulates the large-scale structure of the Universe
davis1985evolution ; springel2005simulations .

Despite the CDM model being successful in describing the formation of large-
scale structure of the Universe, it presents several difficulties to overcome
on small scales ($\lesssim$ hundreds of kpc) weinberg2015cold ;
bullock2017small . The biggest challenge is the so called “cusp-core” (CC)
problem moore1994evidence , where the mass density profile of cold dark matter
halos is expected to has a cusp from collisionless $N$-body simulations of
low-mass galaxies, i.e., $\rho\propto r^{-x}$ with $1\lesssim x\lesssim 1.5$,
which deeply contradicts the observed core-like mass density profile, i.e.,
$\rho\propto r^{0}$ of dwarf galaxies, irregulars, and low surface brightness
galaxies navarro1996structure ; navarro2010diversity ; saburova2014surface .
Another challenge is the “missing satellite” (MSP) problem, which refers to a
discrepancy between the number of predicted number of halo substructures in
collisionless $N$-body simulations and the observed number of satellite
galaxies, i.e., the CDM model predicts that there are around 1000 satellite
galaxies in the Local Group, much higher than the observed number ($\sim$ 50)
of dwarf galaxies moore1999dark ; klypin1999missing . Due to the discrepancies
between the predictions of the CDM model and the actual observations on small
scales, and also the lack of direct evidence supporting cold dark matter
particles such as the WIMPs with masses of order a few GeV to a few hundred
GeV agnes2018low , alternative solutions to the dark-matter problem need to be
fully exploited.

One possibility is that the galactic dark matter is an ultralight scalar field
with negligible self-interactions, which consists of ultralight light scalar
bosons with masses of order a few $10^{-22}$eV hu2000fuzzy ;
goodman2000repulsive . In this scenario, the ultralight scalar bosons form a
Bose-Einstein condensate, in which all particles share one wave function. Such
an alternative theoretical model, so called fuzzy cold dark matter (FCDM)
model, provides intriguing solutions to the cusp-core problem: the ultralight
scalar particles have Compton wavelengths of order a few light years, so that
wave-like behaviors are manifest at the galactic scale, and the cusps in the
center of the FCDM halos could be avoided due to the Heisenberg uncertainty
principle mocz2018schrodinger ; mocz2019first . In this regard, the FCDM model
yields a smoother structure compared to the classical collisionless CDM model
on scales below a few kpc due to quantum pressure mocz2019first . For example,
ultralight scalar bosons with mass of order $10^{-22}$eV transform dark matter
halo cusps of galaxies into soliton cores of order a few kpc
mocz2018schrodinger ; mocz2019first . Beside this the existence of dwarf
galaxies in dark matter halos with masses of $\sim 10^{9}$ times the mass of
the Sun in the Local Universe constrain the FCDM model. For smaller boson
masses, i.e., $m\lesssim 10^{-24}$eV, the observed large-scale cosmic
structure will be erased; while for larger boson masses, i.e., $m\gtrsim
10^{-20}$eV, the classical collisionless CDM behavior on all cosmological
scales recovers, which makes it difficult to restrict the mass of FCDM
particles mocz2018schrodinger .

Observational evidence indicates that at the center of nearly every large
galaxy there is a supermassive black hole with a mass which is million to
billion times the mass of the Sun ferrarese2000fundamental ;
gebhardt2000relationship . Hence the study of how dark matter halos evolve
under the influence of its central supermassive black hole is a necessary step
toward understanding the connection between the galaxies and their surrounding
dark-matter environments reddick2013connection . Although a fully relativistic
analysis of scalar field dark matter evolution around a supermassive blackhole
is possible barranco2012schwarzschild ; barranco2017self ;
avilez2018possibility ; hui2019black , we are here mainly interested in the
dark matter halo mass density profile far from the Schwarzschild radius, so
that the supermassive blackhole can be represented by a point mass
davies2020fuzzy . Such a simplified approach offers the simplest possible
theoretical explanation for existing data. As most galactic cores are known to
host a massive black hole that is accreting matter and displaying a
sufficiently strong luminosity, the impact of blackhole accretion on the dark
matter halo mass density profile also merits study. This is the main goal of
this work.

##  II The Model

To describe the dark matter halo density profiles around an accreting
supermassive black hole, we consider a scalar field dark matter model, in
which the dark matter are comprised of ultralight fuzzy dark matter (FDM)
particles. As the occupation numbers in galactic halos are high, the
ultralight fuzzy dark matter can be described by a scalar field
$\Psi(\mathbf{r},t)$ minimally coupled to the gravity, obeying the
Schrödinger-Newton equation tod1999analytical

$i\hbar\partial_{t}\Psi(\mathbf{r},t)=\left(-\frac{\hbar^{2}\nabla^{2}}{2m}+m\Phi(\mathbf{r},t)+mV(r,t)\right)\Psi(\mathbf{r},t),$
(1)

where $m\approx 10^{-22}$eV is the mass of the fuzzy dark matter particle,
$\Phi(\mathbf{r},t)$ is the self-potential of the wave function which obeys
Poisson’s equation describing the evolution of the gravitational field in the
Newtonian limit

$\nabla^{2}\Phi(\mathbf{r},t)=4\pi G|\Psi(\mathbf{r},t)|^{2},$ (2)

and $V(r,t)$ is a time-dependent black hole potential given by

$V(r,t)=-\frac{GM(t)}{r}.$ (3)

Here, $M(t)$ is a time-dependent mass of the black hole, and $G$ is the
gravitational constant. According to the standard Shakura-Sunyaev disk (SSD)
model for accreting black holes, we assume that the mass of the black hole
varies in time as smole2015smbh

$M(t)\equiv
M_{0}\exp\left(\lambda\frac{1-\epsilon}{\epsilon}\frac{t}{t_{E}}\right)\equiv
M_{0}e^{\alpha t},$ (4)

where $\lambda\equiv L/L_{E}$ is the Eddington ratio, which relates the
luminosity $L$ of an accreting super-massive black hole with the Eddington
luminosity $L_{E}\equiv 4\pi cGMm_{p}/\sigma_{T}$, with $\sigma_{T}$ being the
Thomson scattering cross-section for the electron, and $m_{p}$ being the mass
of the proton. Here, $\epsilon\approx 0.1$ is the radiative efficiency of
accretion, and $t_{E}\approx 4.5\times 10^{7}$yrs is the Salpeter time scale —
an e-folding time-scale for super-massive black hole growth. As the Poisson
equation for the self-potential $\Phi(\mathbf{r},t)$ is solved by the integral

$\Phi(\mathbf{r},t)=-G\int\frac{|\Psi(\mathbf{r}^{\prime},t)|^{2}}{|\mathbf{r}-\mathbf{r}^{\prime}|}d^{3}r^{\prime},$
(5)

the scalar field $\Psi(\mathbf{r},t)$ obeys the following differential-
integral equation

$\displaystyle i\hbar\partial_{t}\Psi(\mathbf{r},t)$
$\displaystyle=\left(-\frac{\hbar^{2}\nabla^{2}}{2m}-Gm\int\frac{|\Psi(\mathbf{r}^{\prime},t)|^{2}}{|\mathbf{r}-\mathbf{r}^{\prime}|}d^{3}r^{\prime}\right.$
$\displaystyle\left.-Gm\frac{M(t)}{r}\right)\Psi(\mathbf{r},t).$ (6)

For spherically symmetric solutions, one obtains a 1+1 dimensional
differential integral equation which governs the spherically symmetric
solutions of the scalar field

$\displaystyle i\hbar\partial_{t}\Psi(r,t)$
$\displaystyle=-\frac{\hbar^{2}}{2mr}\frac{\partial^{2}}{\partial
r^{2}}(r\Psi(r,t))-Gm\left(\frac{4\pi}{r}\int_{0}^{r}|\Psi(r^{\prime},t)|^{2}r^{\prime
2}dr^{\prime}\right.$
$\displaystyle\left.+4\pi\int_{r}^{\infty}|\Psi(r^{\prime},t)|^{2}r^{\prime}dr^{\prime}+\frac{M_{0}e^{\alpha
t}}{r}\right)\Psi(r,t)$ $\displaystyle\equiv\hat{H}(r,t)\Psi(r,t).$ (7)

##  III The Numerical Method

We now use the Crank-Nicolson scheme to solve the above 1+1 dimensional
differential integral equation. To do so, we define the spatial and temporal
grid sizes as $\Delta r$ and $\Delta t$ respectively, and use the index
notation $\Psi_{k}^{n}(r,t)\equiv\Psi(k\Delta r,n\Delta t)$. Then, we can
write Eq. (II) in the Cayley form as salzman2005investigation ;
franklin2016dynamics

$\exp\left(\frac{i\Delta
t\hat{H}}{2\hbar}\right)\Psi_{k}^{n+1}=\exp\left(-\frac{i\Delta
t\hat{H}}{2\hbar}\right)\Psi_{k}^{n}.$ (8)

After linearization, it can be written as

$\Psi_{k}^{n+1}=(\hat{Q}^{-1}-1)\Psi_{k}^{n},\>\hat{Q}\equiv\frac{1}{2}\left(1+\frac{i\Delta
t}{2\hbar}\hat{H}\right).$ (9)

Or equivalently, a system of linear equations

$\Psi^{n+1}=\chi^{n}-\Psi^{n},\>\chi^{n}\equiv\hat{Q}^{-1}\Psi^{n}.$ (10)

Here, the radial part of the Laplacian has the form

$\displaystyle\nabla^{2}$ $\displaystyle\equiv\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial r}\>\>\mbox{for}\>\>r>0,$ (11a)
$\displaystyle\nabla^{2}$ $\displaystyle=3\frac{\partial^{2}}{\partial
r^{2}}\>\>\>\>\>\>\>\>\>\>\>\>\>\mbox{for}\>\>r=0,$ (11b)

where the finite difference discretization scheme for the radial part of the
Laplacian is (see Appendix B)

$\displaystyle\nabla^{2}\Psi_{k}^{n}$ $\displaystyle=\frac{1}{(\Delta
r)^{2}}\left(\frac{k+1}{k}\Psi_{k+1}^{n}-2\Psi_{k}^{n}+\frac{k-1}{k}\Psi_{k-1}^{n}\right)\>\>\mbox{for}\>\>k>0,$
(12a) $\displaystyle\nabla^{2}\Psi_{k}^{n}$ $\displaystyle=\frac{6}{(\Delta
r)^{2}}\left(\Psi_{1}^{n}-\Psi_{0}^{n}\right)\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\>\mbox{for}\>\>k=0.$
(12b)

Besides, the discretization of the self-potential and the black hole potential
have the form

$\displaystyle\Phi_{k}^{n}$ $\displaystyle\approx-4\pi G(\Delta
r)^{2}\left(\frac{1}{k}\sum_{i=0}^{k-1}|\Psi_{i}^{n}|^{2}i^{2}+\sum_{i=k}^{N-1}|\Psi_{i}^{n}|^{2}i\right)\equiv-4\pi
G(\Delta r)^{2}f_{k}^{n},$ (13a) $\displaystyle V_{k}^{n}$
$\displaystyle=-GM_{0}\frac{\exp(\alpha n\Delta t)}{k(\Delta r)}\equiv-
GM_{0}g_{k}^{n}.$ (13b)

Away from the origin and the point before spatial infinity of the spatial
grid, i.e., $k\neq 0,N-1$, one obtains the finite difference equations

$\displaystyle\hat{Q}\chi_{k}^{n}\equiv
a_{k}\chi_{k-1}^{n}+b_{k}\chi_{k}^{n}+c_{k}\chi_{k+1}^{n},$ (14a)
$\displaystyle
a_{k}\equiv\beta\frac{k-1}{k},\>b_{k}\equiv\frac{1}{2}-2\beta+\gamma
f_{k}^{n}+\delta g_{k}^{n},\>c_{k}\equiv\beta\frac{k+1}{k},$ (14b)
$\displaystyle\beta\equiv-\frac{i\hbar}{8m}\frac{\Delta t}{(\Delta
r)^{2}},\gamma\equiv\frac{-i\pi Gm(\Delta r)^{2}\Delta
t}{\hbar},\delta=\frac{-iGM_{0}m\Delta t}{4\hbar}.$ (14c)

Similarly, at the point directly before spatial infinity of the spatial grid,
i.e., $k=N-1$, one obtains the finite difference equation
$\hat{Q}\chi_{N-1}^{n}=a_{N-1}\chi_{N-2}^{n}+b_{N-1}\chi_{N-1}^{n}$, where the
coefficients $a_{N-1}$ and $b_{N-1}$ are given by Eqs. (14b) - (14c) by
substituting $k$ for $N-1$. Finally, at the origin, i.e., $k=0$, one obtains
the finite difference equation $\hat{Q}\chi_{0}^{n}\equiv
b_{0}\chi_{0}^{n}+c_{0}\chi_{0}^{n}$, where
$b_{0}\equiv\frac{1}{2}-6\beta+\gamma f_{0}^{n}+\delta g_{0}^{n}$, and
$c_{0}\equiv 6\beta$. Denoting $\tau\equiv t/t_{E}$, one obtains the
dimensionless finite difference equation at the $(n+1)$th time step,
$\tilde{\Psi}_{k}^{n+1}=(\hat{Q}^{-1}-1)\tilde{\Psi}_{k}^{n}$, where

$\displaystyle\hat{Q}\tilde{\Psi}_{k}^{n}=\tilde{a}_{k}\tilde{\Psi}_{k-1}^{n}+\tilde{b}_{k}\tilde{\Psi}_{k}^{n}+\tilde{c}_{k}\tilde{\Psi}_{k+1}^{n},\tilde{\Psi}_{k}^{n}\equiv\sqrt{4\pi
G}t_{E}\Psi_{k}^{n},$ (15a)
$\displaystyle\tilde{a}_{k}\equiv\tilde{\beta}\frac{k-1}{k},\tilde{b}_{k}\equiv\frac{1}{2}-2\tilde{\beta}+\tilde{\gamma}\tilde{f}_{k}^{n}+\delta\tilde{g}_{k}^{n},\tilde{c}_{k}\equiv\tilde{\beta}\frac{k+1}{k},$
(15b) $\displaystyle\tilde{\alpha}\equiv\alpha
t_{E},\tilde{\beta}=\frac{-i\Delta\tau}{8(\Delta
u)^{2}},\tilde{\gamma}\equiv\frac{-i(\Delta
u)^{2}\Delta\tau}{4},\tilde{\delta}\equiv\frac{-i\lambda\Delta\tau}{4}.$ (15c)

Hence, $\hat{Q}$ can be represented as a sparse $N$ by $N$ matrix in a
tridiagonal form

$\hat{Q}\equiv\begin{pmatrix}b_{0}&c_{0}&0&0&\cdots&0\\\
a_{1}&b_{1}&c_{1}&0&\cdots&0\\\ 0&a_{2}&b_{2}&c_{2}&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\cdots&a_{N-2}&b_{N-2}&c_{N-2}\\\
0&0&\cdots&0&a_{N-1}&b_{N-1}\end{pmatrix}.$ (16)

Notice that only the coefficients $b_{k}$ are time-dependent, since they
contain a sum over the self-potential and the time-dependent black hole
potential. Here, the initial condition for the scalar field $\Psi(r,t)$ is
fixed by setting $\Psi(r,0)=\phi(r)$, where $\phi(r)$ is the steady-state
spherical solution for $\alpha=0$, which corresponds to a time independent
black hole potential with a mass $M(t)\equiv M_{0}$.

##  IV The steady-state spherical solution

The steady-state spherical solution $\phi(r)$ can be solved by the ansatz
$\Psi(r,t)=e^{-iEt/\hbar}\phi(r)$, so that the Schrodinger-Poisson equation
for a time-independent black hole potential becomes

$\displaystyle E\phi(r)$
$\displaystyle=\left(-\frac{\hbar^{2}}{2m}\nabla^{2}+m\Phi(r)+mV(r)\right)\phi(r),$
(17a) $\displaystyle\nabla^{2}\Phi(r)$ $\displaystyle=4\pi G|\phi(r)|^{2}.$
(17b)

Or equivalently

(a) $\lambda=0$ (b) $\lambda=0.1$

(c) $\lambda=0.5$ (d) $\lambda=1.0$

Figure 1: Schematic of the asymptotic normalized densities
$\tilde{\phi}(u)^{2}$ for $u\rightarrow\infty$ as functions of the parameter
$\epsilon$ for different $\lambda$ values, where $\lambda\equiv GM_{0}m/(\hbar
c)$ is the dimensionless mass of the black hole. Here, we fixed the
dimensionless distance $u$ to $15$ when evaluating the asymptotic normalized
densities $\tilde{\phi}(u)^{2}$.

$\displaystyle\frac{\hbar^{2}}{2m}\frac{d^{2}}{dr^{2}}\left(r\phi(r)\right)$
$\displaystyle=\left(rm\Phi(r)-rE-GmM_{0}\right)\phi(r),$ (18a)
$\displaystyle\frac{d^{2}}{dr^{2}}\left(r\Phi(r)\right)$ $\displaystyle=4\pi
Gr\phi(r)^{2},$ (18b)

which are subjected to the initial conditions $\phi(0)=1$,
$\phi^{\prime}(0)=\Phi(0)=\Phi^{\prime}(0)=0$ and the boundary condition
$\phi(\infty)=0$. We use the shooting method to solve the above set of
differential equations, which starts by a trial solution at $r=0$, and tries
to reach the asymptotic boundary condition at $r=\infty$. To simplify
numerical computation, one may introduce a set of dimension-less variables via

$u\equiv\frac{r}{L_{0}},\tilde{\Phi}\equiv\frac{\Phi}{c^{2}},\epsilon\equiv\frac{E}{mc^{2}},$
(19)

so that the renormalized time-independent coupled SP equations becomes

(a) $\lambda=0$ (b) $\lambda=0.1$

(c) $\lambda=0.5$ (d) $\lambda=1.0$

Figure 2: Schematic of the normalized densities $\tilde{\phi}^{2}_{n}\equiv
4\pi G\hbar^{2}\phi_{n}^{2}/(m^{2}c^{4})$ as functions of the normalized
distance $u\equiv mcr/\hbar$ for different $\lambda$ values, where $n$ is an
integer which labels the discrete eigenvalues $\epsilon_{n}$, and
$\lambda\equiv GM_{0}m/(\hbar c)$ is the dimensionless mass of the black hole.
Here, we only plot the normalized densities for the first four eigenvalues.

$\displaystyle\frac{d^{2}}{du^{2}}(u\tilde{\phi})$
$\displaystyle=2u\left(\tilde{\Phi}-\epsilon-\frac{\lambda}{u}\right)\tilde{\phi},$
(20a) $\displaystyle\frac{d^{2}}{du^{2}}(u\tilde{\Phi})$
$\displaystyle=u\tilde{\phi}^{2},$ (20b)

where $L_{0}\equiv\hbar/(mc)$, $\lambda\equiv GM_{0}m/(\hbar c)$, and
$\tilde{\phi}\equiv\sqrt{4\pi G}L_{0}\phi/c$. As an example, for $m\approx
10^{-22}$eV, one obtains $L_{0}\approx 0.1380$ m, $\lambda\approx 5.3832\times
10^{-27}M_{0}$, and $\tilde{\phi}\approx 1.3326\times 10^{-14}\phi$. The above
coupled SP equations can be written as a system of first-order ordinary
differential equation as

$\lambda$ | $\epsilon_{0}$ | $\epsilon_{1}$ | $\epsilon_{2}$ | $\epsilon_{3}$
---|---|---|---|---
0 | 0.649599900660242 | 0.855382257832184 | 0.950376 | 1.027724
0.1 | 0.537515694614713 | 0.742953148283380 | 0.836921 | 0.918315
0.5 | 0.116385322639802 | 0.360002326932728 | 0.459582 | 0.560518
1.0 | -0.411778217857106 | 0.043028169337089 | 0.171107 | 0.281542
Table 1: Table of the four smallest eigenvalues $\epsilon_{n}$ which obey the
asymptotic boundary condition $\tilde{\phi}(\infty)=0$ for different values of
$\lambda$.

(a) $\lambda=0$ (b) $\lambda=0.1$

(c) $\lambda=0.5$ (d) $\lambda=1.0$

Figure 3: Schematic of the normalized densities
$|\tilde{\Psi}(u,\tau)|^{2}\equiv 4\pi G\hbar^{2}/(m^{2}c^{4})|\Psi(r,t)|^{2}$
as functions of radical distance and time in the dimensionless units for an
initial 0-node solution (ground state) for $\lambda\equiv GM_{0}m/(\hbar
c)=0,0.1,0.5$ and $1.0$ respectively. Here, we fix the value of $\alpha=0.1$
for different values of $\lambda$.

$\displaystyle y_{1}^{\prime}$ $\displaystyle=y_{2},$ (21a) $\displaystyle
y_{2}^{\prime}$
$\displaystyle=-\frac{2}{u}y_{2}+2\left(y_{3}-\epsilon-\frac{\lambda}{u}\right)y_{1},$
(21b) $\displaystyle y_{3}^{\prime}$ $\displaystyle=y_{4},$ (21c)
$\displaystyle y_{4}^{\prime}$ $\displaystyle=-\frac{2}{u}y_{4}+y_{1}^{2},$
(21d)

where $y_{1}\equiv\tilde{\phi}$, $y_{3}\equiv\tilde{\Phi}$. One may also
simplify the the coupled SP equations by introducing $y_{1}\equiv
u\tilde{\phi}$ and $y_{3}\equiv u\tilde{\Phi}-\lambda$, which obey the
following set of coupled differential equations

$\displaystyle y_{1}^{\prime\prime}$
$\displaystyle=2\left(\frac{y_{3}}{u}-\epsilon\right)y_{1},$ (22a)
$\displaystyle y_{3}^{\prime\prime}$ $\displaystyle=\frac{1}{u}y_{1}^{2},$
(22b)

where the initial conditions are $y_{1}(0)=0$, $y_{1}^{\prime}(0)=1$,
$y_{3}(0)=-\lambda$, and $y_{3}^{\prime}(0)=0$, and the asymptotic boundary
condition at infinity is $\lim_{u\rightarrow\infty}u^{-1}y_{1}=0$. Eqs. (22a)
- (22b) can be written as a system of first-order ordinary differential
equations as

(a) $\lambda=0$ (b) $\lambda=0.1$

(c) $\lambda=0.5$ (d) $\lambda=1.0$

Figure 4: Schematic of the normalized densities
$|\tilde{\Psi}(u,\tau)|^{2}\equiv 4\pi G\hbar^{2}/(m^{2}c^{4})|\Psi(r,t)|^{2}$
as functions of radical distance and time in the dimensionless units for an
initial 1-node solution (lowest excited state) for $\lambda\equiv
GM_{0}m/(\hbar c)=0,0.1,0.5$ and $1.0$ respectively. Here, we fix the value of
$\alpha=0.1$ for different values of $\lambda$.

$\displaystyle y_{1}^{\prime}$ $\displaystyle=y_{2},$ (23a) $\displaystyle
y_{2}^{\prime}$ $\displaystyle=2\left(\frac{y_{3}}{u}-\epsilon\right)y_{1},$
(23b) $\displaystyle y_{3}^{\prime}$ $\displaystyle=y_{4},$ (23c)
$\displaystyle y_{4}^{\prime}$ $\displaystyle=\frac{1}{u}y_{1}^{2}.$ (23d)

where the initial conditions are $y_{1}(0)=0$, $y_{2}(0)=1$,
$y_{3}(0)=-\lambda$, and $y_{4}(0)=0$.

##  V Restrictions on fuzzy dark matter particle mass

###  V.1 The mass-radius relation for the soliton/ boson star

For an intuiting understanding of nonlinear, bound objects like solitons or
boson stars, one can start with the analysis of different energy scales. In
particular, solitons are those compact objects in which quantum pressure
balances gravitational potential energy. i.e.,

$\frac{GM_{s}}{R_{s}}\propto\frac{1}{2m^{2}}\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}\propto\frac{1}{m^{2}R_{s}^{2}},$
(24)

where $\rho\equiv m|\phi|^{2}$ is the mass density, $R_{s}$ is the radius of
the soliton, and $M_{s}$ is the total mass of the soliton given by
davies2020fuzzy

$M_{s}\equiv\int_{0}^{\infty}|\phi|^{2}4\pi r^{2}dr=\frac{\hbar
c}{Gm}\int_{0}^{\infty}\tilde{\phi}^{2}u^{2}du\equiv\frac{\hbar
c\tilde{M}_{s}}{Gm},$ (25)

where $u\equiv r/L_{0}$, $L_{0}\equiv\hbar/(mc)$,
$\tilde{\phi}\equiv\sqrt{4\pi G}L_{0}\phi/c$, and $\tilde{M}_{s}$ is the
normalized total mass of the soliton. The constant $\hbar c/(Gm)$ sets the
mass of the soliton and is given by

$\frac{\hbar c}{Gm}\approx 1.336\times
10^{12}M_{\odot}\left(\frac{10^{-22}\mbox{eV}}{m}\right).$ (26)

One may introduce the dimensionless ratio of the black hole mass to the
soliton mass, $\eta\equiv\lambda/\tilde{M}_{s}=M_{0}/M_{s}$, where $M_{0}$ is
the mass of the supermassive black hole, and $\lambda$ sets the dimensionless
mass of the black hole. Clearly, the mass ratio $\eta$ is invariant under
scaling, and is the main factor which determines the density profile of the
soliton.

From Eq. (24), one deduces that the total mass of a soliton is inversely
proportional to its core size, i.e., hui2021wave

$M_{s}R_{s}\propto\frac{1}{Gm^{2}}\propto 2.2\times 100\mbox{pc}\times
10^{9}M_{\odot}\left(\frac{10^{-22}\mbox{eV}}{m}\right)^{2}.$ (27)

Using the data $1\mbox{pc}=3.085677581\times 10^{16}\mbox{m}$, one obtains the
mass-radius relation for a soliton in the dimensionless unit

(a) $\lambda=0.1$ (b) $\lambda=0.1$

Figure 5: Schematic of the normalized densities
$|\tilde{\Psi}(u,\tau)|^{2}\equiv 4\pi G\hbar^{2}/(m^{2}c^{4})|\Psi(r,t)|^{2}$
as functions of radical distance and time in the dimensionless units for
$\alpha=0,0.1$ and $0.2$ respectively. In the left panel we plot the
normalized densities for an initial 0-node solution (ground state), while in
the right panel we plot the normalized densities for an initial 1-node
solution (lowest excited state). Here, we fix the value of $\lambda=0.1$ for
different values of $\alpha$.

$\tilde{M}_{s}U_{s}=2.57436836743.$ (28)

###  V.2 Soliton condensation: the relationship between the soliton mass and
the halo mass

On the large scale of the Universe, it was first pointed out by Schive et al.
schive2014cosmic ; schive2014understanding  from the cosmological simulations
that, virialized halos tend to have a core with a soliton mass that scales
with the halo mass as:

$M_{s}\propto a^{-1/2}M_{h}^{1/3},$ (29)

where $a$ is the cosmic scale factor, $M_{s}$ is the soliton mass, and $M_{h}$
is the dark matter halo mass. Employing the scaling relation at the present
age of the Universe ($a=1$), one obtains

$M_{s}\approx 1.25\times
10^{9}M_{\odot}\left(\frac{M_{h}}{10^{12}M_{\odot}}\right)^{1/3}\left(\frac{10^{-22}\mbox{eV}}{m}\right).$
(30)

As a remark, a dark matter halo is said to be virialized when it has reached
an equilibrium between the kinetic and potential energy, and and obeys the
virial theorem hui2021wave . The virial radius of a dark matter halo is the
radius that which the density of the halo equals $200$ times the critical
density of the Universe. The mass contained inside the virial radius is
regarded as the total mass of the dark matter halo.

###  V.3 The relationship between the soliton mass and the FDM particle mass

The restriction of the range of dark matter soliton mass mainly comes from the
cosmological setting. In order that the soliton does not collapse directly to
a black hole, there is an upper limit of the soliton mass hui2021wave :
$GM/R\lessapprox 1$. Substitution in the mass-radius relation for the soliton
yields the maximum soliton mass $M_{2}$ hui2017ultralight

Figure 6: The log-log plot of the soliton mass $M_{s}$ against the fuzzy dark
matter mass $m$. The total mass of the soliton at the center of the Milky Way
and M87 are shown with dark and blue solid curves respectively. The soliton
mass is obtained from the relationship between the soliton mass and the mass
of halos where the galaxies reside. The upper and the lower bounds of the dark
matter soliton mass are shown by red and magenta broken lines respectively.

$M_{2}\propto\frac{1}{Gm}\propto 8.46\times
10^{11}M_{\odot}\left(\frac{m}{10^{-22}\mbox{eV}}\right).$ (31)

Similarly, there is a lower limit of the soliton mass. The soliton ceases to
exist if the central density is lower than $200$ times the critical density of
the Universe, as we have remarked above. Numerically, the minimum soliton mass
$M_{1}$ is hui2017ultralight

$M_{1}\propto 1.40\times
10^{7}M_{\odot}\left(\frac{m}{10^{-22}\mbox{eV}}\right)^{-3/2}.$ (32)

Finally, using the following data

$\displaystyle\frac{GM_{\odot}}{c^{2}}$ $\displaystyle=14.7662563825\times
10^{8}\mu\mbox{m},$ (33a) $\displaystyle\hbar c$
$\displaystyle=0.1973269804\mbox{eV}\times\mu\mbox{m},$ (33b)

one readily obtains from Eq. (25) the numerical relationship between the
soliton mass and the FDM particle mass

$M_{s}/M_{\odot}=1.336337765\times 10^{-x-10}\times\tilde{M}_{s},$ (34)

where $x\equiv\log_{10}(m/\mbox{eV})$, and $\tilde{M}_{s}$ is the normalized
total mass of the soliton.

##  VI The Result

From Figs. 1 and 2, one can see that, for a given dimensionless black hole
mass $\lambda$, there are a unique set of discrete eigenvalues $\epsilon_{n}$
which obey the asymptotic boundary condition $\phi(\infty)=0$. In the
dimensionless unit, the eigenvalues $\epsilon_{n}$ and the associated
eigenfunctions $\tilde{\phi}_{n}$ are labeled by the number of nodes of
$\tilde{\phi}(u)$ as functions of the normalized distance $u\equiv mcr/\hbar$.
The $0$-node solution in the fuzzy dark matter theory is regarded as the
soliton, or boson star. In other words, the dimensionless black hole mass
determines the density profile of the soliton.

In Table. 1, we list the first few lowest discrete eigenvalues of the
nonlinear Schrödinger-Newton equation. Unlike its linear counterpart, the
asymptotic boundary condition $\phi(\infty)=0$ cannot be truly satisfied as it
requires an infinite accuracy of the nonlinear eigenvalues. To ensure
sufficient accuracy, we retain 16 digits in the computation of the lowest
eigenvalues.

Figure 7: The log-log plot of the enclosed soliton mass $M$ within a radius of
$r\approx 0.01\mbox{pc}$ near the central supermassive blackhole against the
fuzzy dark matter mass $m$. Here, the dimensionless mass of the supermassive
blackhole is chosen to be $0$, $0.1$, $0.5$ and $1.0$ respectively.

In Fig. 3, we plot the FDM soliton mass density profile with respect to the
radical distance in the dimensionless unit for an initial 0-node solution. The
results show that when the black hole is accreting at a fixed accretion rate
(e.g., $\alpha=0.1$), the dark matter soliton density profile for an initial
0-node solution is squeezed near the central supermassive black hole. For a
fixed accretion rate, the squeezing effects for larger black holes are more
significant than those for smaller black holes. In Fig. 4, the schematic of
the FDM soliton mass density profile for an initial 1-node solution also shows
significant squeezing effects surrounding the central black hole.

In Fig. 5, we compare the FDM soliton-like mass density profiles with a fixed
dimensionless black hole mass $\lambda=0.1$ for different accretion rates,
i.e., $\alpha=0,0.1$ and $0.2$. The result shows that the squeezing effects
are more significant when the accretion rate for the central supermassicve
increases, and is weaken for an initial 1-node solution compared to an initial
0-node FDM soliton.

In Fig. 6, we see that for a fuzzy dark matter particle with a mass
$m\lessapprox 10^{-26}\mbox{eV}$, the total mass of the FDM soliton at the
center of Milky Way already fails below the lower bound. Hence, the mass of
FDM particle has to be $m\gtrapprox 10^{-26}\mbox{eV}$.

In Fig. 7, we compare the mass profile of different solitons near the galactic
central supermassive black hole ($M_{0}\propto 10^{5}\sim 10^{10}M_{\odot}$)
with dimensionless masses $\lambda\equiv GM_{0}m/(\hbar c)$ ranging from $0$
to $1$, where $M_{0}$ is the mass of the supermassive black hole. The mass
profile is mainly unaffected by the FDM particle mass when $m\lessapprox
10^{-22.5}\mbox{eV}$, and is significantly modulated by the FDM particle mass
when $m\gtrapprox 10^{-22.5}\mbox{eV}$. Surrounding a supermassive black hole,
the mass of enclosed soliton can be lowered by two orders of magnitude, when
the black hole mass is comparable to the total mass of the soliton.

##  VII Conclusion

We numerically computed the mass density profiles of the soliton cores in the
center of fuzzy dark matter halos by solving the Schrödinger-Newton equation
with the presence of a time-dependent point mass perturber which represents
the supermassive blackhole accretion at the center of galaxy. We showed that
there is a significant squeezing effects on the density profile in both the
core size and the central density.

Here our computation for the density profiles is non-relativistic. In future
works, we may perform a fully relativistic analysis of scalar field dark
matter evolution around an accreting supermassive blackhole, as well as self-
consistent time-dependent simulations of the blackhole-soliton-halo system.

##  Appendix A Derivation of the 1+1 dimensional differential integral
equation for spherically symmetric solutions

For spherically symmetric solutions, one may use the Laplace expansion of the
inverse distance between the two points $\mathbf{r}$ and $\mathbf{r}^{\prime}$

$\frac{1}{|\mathbf{r}-\mathbf{r}^{\prime}|}=\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}\frac{r_{>}^{l}}{r_{<}^{l+1}}Y_{l}^{-m}(\theta,\phi)Y_{l}^{m}(\theta^{\prime},\phi^{\prime}),$
(35)

where $r_{<}=\mbox{min}(r,r^{\prime})$, $r_{>}=\mbox{max}(r,r^{\prime})$, and
$Y_{l}^{m}(\theta,\phi)$ is a normalized spherical harmonic function. Then, an
integration of the inverse distance between the two points $\mathbf{r}$ and
$\mathbf{r}^{\prime}$ weighted by the density $|\Psi(r^{\prime},t)|^{2}$
becomes

$\displaystyle\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}\int_{0}^{\infty}|\Psi(r^{\prime},t)|^{2}r^{\prime
2}dr^{\prime}\frac{r_{<}^{l}}{r_{>}^{l+1}}Y_{l}^{-m}(\theta,\phi)\int
Y_{l}^{m}(\theta^{\prime},\phi^{\prime})d\Omega^{\prime}$
$\displaystyle=\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}\int_{0}^{\infty}|\Psi(r^{\prime},t)|^{2}r^{\prime
2}dr^{\prime}\frac{r_{<}^{l}}{r_{>}^{l+1}}Y_{l}^{-m}(\theta,\phi)\sqrt{4\pi}\delta_{l,0}\delta_{m,0}$
$\displaystyle=\frac{4\pi}{r}\int_{0}^{r}|\Psi(r^{\prime},t)|^{2}r^{\prime
2}dr^{\prime}+4\pi\int_{r}^{\infty}|\Psi(r^{\prime},t)|^{2}r^{\prime}dr^{\prime},$
(36)

where
$d\Omega^{\prime}\equiv\sin\theta^{\prime}d\theta^{\prime}d\phi^{\prime}$ is
the differential solid angle in the direction
$(\theta^{\prime},\phi^{\prime})$. Using the composite trapezoidal rule, one
may discretize the first integral along the radial spatial direction as

$\displaystyle\frac{4\pi(\Delta
r)^{2}}{k}\left(\sum_{j=1}^{k-1}|\Psi_{j}(t)|^{2}j^{2}+\frac{k^{2}}{2}|\Psi_{k}(t)|^{2}\right)\>\>$
$\displaystyle\mbox{for}\>\>k>1,$ (37a) $\displaystyle 4\pi(\Delta
r)^{2}\cdot\frac{1}{2}|\Psi_{1}(t)|^{2}\>\>\>\>\>\>\>\>\>\>\>\>\>\>$
$\displaystyle\mbox{for}\>\>k=1,$ (37b)

where $r\equiv k\Delta r$ and $\Psi_{j}(t)\equiv\Psi(j\Delta r,t)$. Similarly,
one may discretize the second integral along the radial spatial direction as

$\displaystyle 4\pi(\Delta
r)^{2}\left(\frac{k}{2}|\Psi_{k}(t)|^{2}+\sum_{j=k+1}^{N-1}|\Psi_{j}(t)|^{2}j+\frac{N}{2}|\Psi_{N}(t)|^{2}\right)\>\>$
$\displaystyle\mbox{for}\>\>k<N-1,$ (38) $\displaystyle 4\pi(\Delta
r)^{2}\left(\frac{N-1}{2}|\Psi_{N-1}(t)|^{2}+\frac{N}{2}|\Psi_{N}(t)|^{2}\right)\>\>\>\>\>\>\>\>$
$\displaystyle\mbox{for}\>\>k=N-1,$ (39)

where $r_{max}\equiv N\Delta r$ for a large positive integer $N$ represents
the spatial infinity. Hence, for $0<k<N-1$, the discretization of the self-
potential along the radial spatial direction has the form

$\displaystyle\Phi_{k}(t)$ $\displaystyle=-4\pi G(\Delta
r)^{2}\left(\frac{1}{k}\sum_{j=1}^{k}|\Psi_{j}(t)|^{2}j^{2}\right.$
$\displaystyle\left.+\sum_{j=k+1}^{N-1}|\Psi_{j}(t)|^{2}j+\frac{N}{2}|\Psi_{N}(t)|^{2}\right),$
(40)

where $\Phi_{k}(t)\equiv\Phi(k\Delta r,t)$ and we used the boundary condition
$\Psi_{N}(t)=0$ at the spatial infinity. In particular, for $k=N-1$, the
discretization of the self-potential in the radial direction has the form

$\Phi_{1}(t)=-4\pi G(\Delta
r)^{2}\left(\frac{1}{k}\sum_{j=1}^{N-1}|\Psi_{j}(t)|^{2}j+\frac{N}{2}|\Psi_{N}(t)|^{2}\right).$
(41)

For a non-uniform partition of $[0,r_{k}]$ such that
$0=r_{0}<r_{1}<\cdots<r_{k-1}<r_{k}$, the trapezoidal rule yields the
discretization of the first integral along the spatial direction

$\frac{4\pi}{r_{k}}\left[\frac{1}{2}\sum_{j=1}^{k}\left(|\Psi_{j-1}(t)|^{2}r_{j-1}^{2}+|\Psi_{j}(t)|^{2}r_{j}^{2}\right)\Delta
r_{j}\right],$ (42)

where $\Delta r_{j}\equiv r_{j}-r_{j-1}$. Similarly, for a non-uniform
partition of $[r_{k},r_{max}]$ such that
$r_{k}<r_{k+1}<\cdots<r_{N-1}<r_{N}=r_{max}$, the trapezoidal rule yields the
discretization of the second integral along the spatial direction

$4\pi\left[\frac{1}{2}\sum_{j=k+1}^{N}\left(|\Psi_{j-1}(t)|^{2}r_{j-1}+|\Psi_{j}(t)|^{2}r_{j}\right)\Delta
r_{j}\right].$ (43)

Using the central difference formula, the discretization of the second order
derivative of $\Psi(r,t)$ with respect to $r$ at $r_{k}$ reads

$\displaystyle\left.\frac{\partial^{2}\Psi(r,t)}{\partial
r^{2}}\right|_{r=r_{k}}=\frac{\Psi^{\prime}(r_{k+1/2},t)-\Psi^{\prime}(r_{k-1/2},t)}{r_{k+1/2}-r_{k-1/2}}$
$\displaystyle=\frac{2}{\Delta r_{k+1}+\Delta
r_{k}}\left(\frac{\Psi(r_{k+1},t)-\Psi(r_{k},t)}{\Delta
r_{k+1}}-\frac{\Psi(r_{k},t)-\Psi(r_{k-1},t)}{\Delta r_{k}}\right)$
$\displaystyle=\frac{2\Psi_{k+1}(t)}{(\Delta r_{k+1}+\Delta r_{k})\Delta
r_{k+1}}-\frac{2\Psi_{k}(t)}{\Delta r_{k+1}\Delta
r_{k}}+\frac{2\Psi_{k-1}(t)}{(\Delta r_{k+1}+\Delta r_{k})\Delta r_{k}}.$ (44)

Similarly, the discretization of the first order derivative of $\Psi(r,t)$
with respect to r at $r_{k}$ reads

$\left.\frac{\partial\Psi(r,t)}{\partial
r}\right|_{r=r_{k}}=\frac{\Psi_{k+1}(t)-\Psi_{k-1}(t)}{\Delta r_{k+1}+\Delta
r_{k}}.$ (45)

Hence, the finite difference formula for the radial part of the Laplacian
reads

$\displaystyle\left.\nabla^{2}\Psi(r,t)\right|_{r=r_{k}}\equiv\left.\left(\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial
r}\right)\right|_{r=r_{k}}=\frac{2\Psi_{k+1}(t)}{(\Delta r_{k+1}+\Delta
r_{k})\Delta r_{k+1}}\frac{r_{k+1}}{r_{k}}$
$\displaystyle-\frac{2\Psi_{k}(t)}{\Delta r_{k+1}\Delta
r_{k}}+\frac{2\Psi_{k-1}(t)}{(\Delta r_{k+1}+\Delta r_{k})\Delta
r_{k}}\frac{r_{k-1}}{r_{k}}.$ (46)

##  Appendix B Boundary conditions for the scalar field at $r\rightarrow 0$

In this work, the boundary conditions for the scalar field $\Psi(r,t)$ are
given by

$\partial_{r}\Psi(0,t)=0\>\>\mbox{and}\>\>\Psi(\infty,t)=0.$ (47)

For the radical part of the Laplacian, one obtains

$\lim_{r\rightarrow 0}\nabla^{2}\Psi(r,t)\equiv\lim_{r\rightarrow
0}\left(\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial
r}\right)\Psi(r,t)=3\frac{\partial^{2}\Psi}{\partial r^{2}}(r,t),$ (48)

where we have used the L’Hospital’s Rule to determine the indeterminate form
of $0/0$

$\lim_{r\rightarrow 0}\frac{1}{r}\frac{\partial\Psi(r,t)}{\partial
r}=\frac{\partial^{2}\Psi(r,t)}{\partial r^{2}}.$ (49)

The finite difference discretization of the boundary condition at
$r\rightarrow 0$ has the form

$\lim_{r\rightarrow 0}\frac{\partial\Psi}{\partial
r}(r,t)\approx\frac{\Psi_{1}^{n}-\Psi_{-1}^{n}}{2\Delta r}=0,$ (50)

which yields the boundary condition $\Psi_{1}^{n}=\Psi_{-1}^{n}$ at arbitrary
time steps $n$. Hence, the finite difference discretization of the radical
part of the Laplacian has the form

$\lim_{r\rightarrow 0}\nabla^{2}\Psi(r,t)=\frac{3}{(\Delta
r)^{2}}\left(\Psi_{1}^{n}-2\Psi_{0}^{n}+\Psi_{-1}^{n}\right)=\frac{6}{(\Delta
r)^{2}}\left(\Psi_{1}^{n}-\Psi_{0}^{n}\right).$ (51)

## References

  * (1) G. Bertone and T. M. P. Tait, Nature 562, 51 (2018).
  * (2) C. S. Frenk and S. D. M. White, Ann. Phys. (Berl.) 524, 507 (2012).
  * (3) M. Vogelsberger, F. Marinacci, P. Torrey, and E. Puchwein, Nat. Rev. Phys. 2, 42 (2020).
  * (4) I. de Martino, S. S. Chakrabarty, V. Cesare, A. Gallo, L. Ostorero, and A. Diaferio, Universe 6, 107 (2020).
  * (5) M. Davis, G. Efstathiou, C. S. Frenk, and S. D. M. White, Astrophys. J. 292, 371 (1985).
  * (6) V. Springel, S. D. M. White, A. Jenkins, C. S. Frenk, N. Yoshida, L. Gao, J. Navarro, R. Thacker, D. Croton, J. Helly, J. A. Peacock, S. Cole, P. Thomas, H. Couchman, A. Evrard, J. Colberg, and F. Pearce, Nature 435, 629 (2005).
  * (7) D. H. Weinberg, J. S. Bullock, F. Governato, R. K. de Naray, and A. H. Peter, Proc. Natl. Acad. Sci. 112, 12249 (2015).
  * (8) J. S. Bullock and M. Boylan-Kolchin, Annu. Rev. Astron. Astrophys. 55, 343 (2017).
  * (9) B. Moore, Nature 370, 629 (1994).
  * (10) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996).
  * (11) J. F. Navarro, A. Ludlow, V. Springel, J. Wang, M. Vogelsberger, S. D. M. White, A. Jenkins, C. S. Frenk, and A. Helmi, MNRAS 402, 21 (2010).
  * (12) A. Saburova and A. D Popoloć, MNRAS 445, 3512 (2014).
  * (13) B. Moore, S. Ghigna, F. Governato, G. Lake, T. Quinn, J. Stadel, and P. Tozzi, Astrophys. J. Lett. 524, L19 (1999).
  * (14) A. Klypin, A. V. Kravtsov, O. Valenzuela, and F. Prada, Astrophys. J. 522, 82 (1999).
  * (15) P. Agnes, I. F. M. Albuquerque, T. Alexander, A. K. Alton, G. R. Araujo, D. M. Asner, M. Ave, H. O. Back, B. Baldin, G. Batignani, et al., Phys. Rev. Lett. 121, 081307 (2018).
  * (16) W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000).
  * (17) J. Goodman, New Astron. 5, 103 (2000).
  * (18) P. Mocz, L. Lancaster, A. Fialkov, F. Becerra, and P. H. Chavanis, Phys. Rev. D 97, 083519 (2018).
  * (19) P. Mocz, A. Fialkov, M. Vogelsberger, F. Becerra, A. A. Amin, S. Bose, M. Boylan-Kolchin, P. H. Chavanis, L. Hernquist, L. Lancaster, F. Marinacci, V. H. Robles, and J. Zavala, Phys. Rev. Lett. 123, 141301 (2019).
  * (20) L. Ferrarese and D. Merritt, Astrophys. J. Lett. 539, L9 (2000).
  * (21) K. Gebhardt, R. Bender, G. Bower, A. Dressler, S. M. Faber, A. V. Filippenko, R. Green, C. Grillmair, L. C. Ho, J. Kormendy, T. D. Lauer, J. Magorrian, J. Pinkney, D. Richstone, and S. Tremaine, Astrophys. J. Lett. 539, L13 (2000).
  * (22) R. M. Reddick, R. H. Wechsler, J. L. Tinker, and P. S. Behroozi, Astrophys. J. 771, 30 (2013).
  * (23) J. Barranco, A. Bernal, J. C. Degollado, A. Diez-Tejedor, M. Megevand, M. Alcubierre, D. Núñez, and O. Sarbach, Phys. Rev. Lett. 109, 081102 (2012).
  * (24) J. Barranco, A. Bernal, J. C. Degollado, A. Diez-Tejedor, M. Megevand, D. Núñez, and O. Sarbach, Phys. Rev. D 96, 024049 (2017).
  * (25) A. A. Avilez, T. Bernal, L. E. Padilla, and T. Matos, MNRAS 477, 3257 (2018).
  * (26) L. Hui, D. Kabat, X. Li, L. Santoni, and S. S. Wong, MNRAS 2019, 2019 (2019).
  * (27) E. Y. Davies and P. Mocz, MNRAS 492, 5721 (2020).
  * (28) P. Tod and I. M. Moroz, Nonlinearity 12, 12 (1999).
  * (29) M. Smole, M. Micic, and N. Martinović, MNRAS 451, 1964 (2015).
  * (30) J. Franklin, Y. Guo, K. C. Newton, and M. Schlosshauer, Classical Quant Grav 33, 075002 (2016).
  * (31) P. J. Salzman, Investigation of the time-dependent Schroedinger-Newton equation (University of California, Davis, 2005).
  * (32) L. Hui, Annu. Rev. Astron. Astrophys. 59, 247 (2021).
  * (33) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D 95, 043541 (2017).
  * (34) H. Y. Schive, T. Chiueh, and T. Broadhurst, Nature Physics 10, 496 (2014).
  * (35) H. Y. Schive, M. H. Liao, T. P. Woo, S. K. Wong, T. Chiueh, T. Broadhurst, and W. P. Hwang, Phys. Rev. Lett. 113, 261302 (2014).