\section{Introduction}\label{sec:Intro}

Increasingly precise constraints on $\Lambda$CDM parameters from different probes have revealed tensions between measurements from early and late time Universe, including the well-known $H_0$ tension. The value of $H_0$ inferred from Cosmic Microwave Background (CMB) temperature, polarization and lensing anistropies by the Planck Collaboration is 67.36 $\pm$ 0.54 km/s/Mpc~\cite{Planck2018:cosmo}. High resolution CMB observations from ground based experiments such as ACTPol and SPT3G have also independently measured $H_0$ and found 67.9 $\pm$ 1.5 km/s/Mpc~\cite{ACT:A20} and 68.8 $\pm$ 1.5 km/s/Mpc~\cite{dutcher2021measurements} respectively.  All of these measurements are inconsistent with the latest measurement from the Cepheids-calibrated cosmic distance ladder $H_0 = 73.2 \pm 1.3$ km/s/Mpc~\cite{Riess2021}.\\

This discrepancy could indicate the need for new physics beyond $\Lambda$CDM. Many models have already been proposed to solve the Hubble tension. Modifications to early time physics, such as new physics that change the physical size of the sound horizon at recombination, are investigated and the different ways of solving the current tension have been reviewed in~\cite{Knox2020, Divalentino2021}. Upcoming cosmic variance limited measurements of CMB polarization for a wide multipole range from Simons Observatory~\cite{Ade2019} and CMB-S4~\cite{Abazajian2016} will put severe constraints on these models.\\

Another hypothesis which could explain these tensions is the presence of remaining systematic effects not accounted for in the data analysis. The presence of bias in the cosmic distance ladder data analysis is actively studied, in particular biases coming from supernovae light curves standardization~\cite{Rigault2018, Saunders2018}. An alternative calibration of cosmic distances, using the Tip of the Red Giant Branch (TRGB) is also explored, giving another measurement of the Hubble parameter~\cite{Freedman2019, Freedman2020}.\\

In this paper we focus on possible instrumental systematic effects that can affect the measurement of CMB anisotropies. A way of testing this possibility is to use observables that are less sensitive to the exact instrument model. One of these observables is the Pearson's correlation coefficient of T and E modes $\mathcal{R}_\ell^{\rm TE}=C_\ell^\mathrm{TE} / \sqrt{C_\ell^\mathrm{TT}C_\ell^\mathrm{EE}}$. It was mentioned in~\cite{Planck2015:ps}, in which the authors give geometrical interpretations of the correlation coefficient as a cosine or as a decorrelation angle. Correlation coefficients have also been used to look at the correlations between galaxy surveys and CMB-derived lensing power spectrum in an unbiased way~\cite{Peacock2018}.  $\mathcal{R}_\ell^{\rm TE}$ statistical properties have been discussed in details in~\cite{Louis2019}. This observable has been shown to be robust against multiplicative instrumental systematics, such as beam error, calibration, polarization efficiency uncertainties or unmodelled transfer functions.\\

This work presents the first estimation of the cosmological parameters from $\mathcal{R}_\ell^{\rm TE}$. The nearly Gaussian nature of the correlation coefficient, at scales where we can ensure a high $EE$ signal-to-noise ratio, makes possible the use of a Gaussian likelihood function to constrain cosmological parameters and to compare with constraints from the combination of $C_\ell^{\rm TT}$, $C_\ell^{\rm TE}$ and $C_\ell^{\rm EE}$ power spectra.\\

The paper is organized as follows. In section~\ref{sec:syst} we illustrate the effect of different systematic effects on cosmological parameters constrained from the CMB power spectra. In section~\ref{sec:like} we discuss the construction of a $\mathcal{R}_\ell^{\rm TE}$-based likelihood. In section~\ref{sec:data} we show the results of a $\mathcal{R}_\ell^{\rm TE}$-based analysis for the Planck PR4 dataset. We conclude in section~\ref{sec:Conclusion}.

\section{Effect of systematics in power spectra based likelihood}\label{sec:syst}

In this section, we introduce the dataset and likelihood used in the rest of the paper, we then study the impact of systematics on cosmological parameter estimation using a set of biased simulations.

\subsection{Dataset and $C_\ell$-based likelihood}\label{subsec:dataset}

We use the latest Planck data release (PR4). The maps were produced using the \texttt{NPIPE} processing pipeline that jointly analyzes data from Planck High Frequency Instrument (HFI) and Planck Low Frequency Instrument (LFI)~\cite{PlanckNPIPE}.\\

Cosmological parameters are estimated using the Planck \texttt{HiLLiPoP} (High-L Likelihood Polarized for Planck) likelihood. The code is publicly available on github\footnote{\url{https://github.com/planck-npipe}}. This is a multifrequency likelihood for the Planck cosmology channels : 100, 143 and 217 GHz. The dataset consists of two split map for each frequency. From these maps, we obtain 15 cross power spectra from which we compute the six cross-frequency spectra used in the likelihood. The foreground model includes galactic dust, CIB, tSZ \& kSZ contributions, tSZ-CIB correlation and Poisson-like distributed point source foregrounds. A more detailled description is given in~\cite{Couchot2017}. The likelihood assumes that the $TT$, $TE$ and $EE$ CMB power spectra follow a Gaussian distribution, which is a good assumption at high multipoles ($\ell>30$). The likelihood has been tested extensively against the Plik likelihood~\cite{Planck2015:ps}.\\

In section~\ref{subsec:syst}, we study the impact of systematics on simulated Planck data. Cosmological constraints from the Planck PR4 dataset are presented in section~\ref{sec:data}.
\subsection{Bias on the cosmological parameters}\label{subsec:syst}
The correlation coefficient is designed to be insensitive to any multiplicative bias. To illustrate the effect of this kind of bias on cosmological parameters estimated from the combination of the $TT$, $TE$ and $EE$ power spectra, we generate a set of biased simulations and fit for the cosmological parameters. The simulations are generated at the power spectra level from the best fit \texttt{HiLLiPoP} cosmology and foreground model, we then use the Planck PR4 covariance matrix to add scatter to the spectra. Fiducial cosmology has been set to $100\theta_{\mathrm{MC}} = 1.04065$, $\Omega_bh^2 = 0.02231$, $\Omega_ch^2 = 0.1193$, $\mathrm{ln}(10^{10}A_\mathrm{s}) = 3.045$, $n_\mathrm{s} = 0.9619$ and $\tau = 0.0566$.\\

We include systematics into our simulated dataset such that the observed temperature $\Tilde{a}_{\ell m}^\mathrm{T}$ and polarization E-modes $\Tilde{a}_{\ell m}^\mathrm{E}$ are given by $\Tilde{a}_{\ell m}^\mathrm{T} = \epsilon_\ell^\mathrm{T}a_{\ell m}^\mathrm{T}$ and $\Tilde{a}_{\ell m}^\mathrm{E} = \epsilon_\ell^\mathrm{P}a_{\ell m}^\mathrm{E}$. The measured power spectra are
\begin{align}\label{eq:bias}
    \Tilde{C}_\ell^{TT} &= ({\epsilon_\ell^{\mathrm{T}}})^2C_\ell^{TT}\nonumber\\
    \Tilde{C}_\ell^{TE} &= {\epsilon_\ell^{\mathrm{T}}}\epsilon_\ell^{\mathrm{P}}C_\ell^{TE}\nonumber\\
    \Tilde{C}_\ell^{EE} &= ({\epsilon_\ell^{\mathrm{P}}})^2C_\ell^{EE}
\end{align}
Systematic effects can be constant over multipoles (e.g. polarization efficiency) or scale-dependant (e.g. transfer functions).\\

We obtain the posterior distributions of cosmological parameters using the Markov Chain Monte Carlo (MCMC) algorithm implemented in \texttt{cobaya}~\cite{Torrado2020} with the $C_\ell$-based \texttt{HiLLiPoP} likelihood. Figure~\ref{fig:bias_syst} displays the distributions of the six standard $\Lambda$CDM parameters for a set of different transfer functions. We consider three different kinds of systematics :
\begin{enumerate}[a.]
    \item $\epsilon^{\rm P}_\ell = \mathrm{cte}$, $\epsilon^{\rm T}_\ell = 1$ (Fig.~\ref{fig1:poleff})
    \item $\epsilon^{\rm P}_\ell = \epsilon^{\rm P}(\ell)$, $\epsilon^{\rm T}_\ell = 1$ (Fig.~\ref{fig1:polTF})
    \item $\epsilon^{\rm P}_\ell = 1$, $\epsilon^{\rm T}_\ell = \epsilon^{\rm T}(\ell)$ (Fig.~\ref{fig1:tempTF})
\end{enumerate}

For illustration purpose we choose a transfer function such as
\begin{equation}\label{eq:tf_model}
    \epsilon(\ell) =
    \begin{cases}
        \epsilon_{\mathrm{min}} & \text{if $\ell < \ell_{\mathrm{min}}$} \\
        \epsilon_{\mathrm{min}} + \Delta\epsilon\cdot\mathrm{sin}^2\left(\frac{\pi}{2}\frac{\ell - \ell_{\mathrm{min}}}{\Delta\ell}\right) & \text{if $\ell_{\mathrm{min}}\le\ell\le\ell_{\mathrm{max}}$} \\
        \epsilon_{\mathrm{max}} & \text{if $\ell > \ell_{\mathrm{max}}$}
    \end{cases}
\end{equation}
where $\Delta\epsilon = \epsilon_{\rm max} - \epsilon_{\rm min}$ and $\Delta\ell = \ell_{\rm max} - \ell_{\rm min}$ with $\ell_{\rm min} = 100$. The transfer function will smoothly increase from $\epsilon_\mathrm{min} = 0.95$ to $\epsilon_\mathrm{max} = 1$ between $\ell_\mathrm{min}$ and $\ell_\mathrm{max}$.\\

Adding an inconsistency between temperature and polarization leads to significant shifts for all the cosmological parameters for the three different models used in this analysis. The effect is clearly noticeable for baryon density and particularly for the polarization efficiency case~(\ref{fig1:poleff}). Using the biased simulations, we obtain constraints on $\Omega_bh^2$ that are more than $3\sigma$ discrepant with the unbiased $C_\ell$s-derived constraint for $\epsilon_{\rm P}\le0.97$. $H_0$ is particularly affected by a temperature transfer function : the $H_0$ constraint is shifted towards lower values~(\ref{fig1:tempTF}). The $\mathcal{R}_\ell^{TE}$ correlation coefficient could therefore be an interesting consistency test, it could help avoiding large bias in our cosmological parameter measurements arising from multiplicative systematic effects.\\

\begin{figure*}
\centering
\subfloat[Polarization efficiency]{\label{fig1:poleff}\includegraphics{bias_pol.pdf}}

\subfloat[Polarization transfer function]{\label{fig1:polTF}\includegraphics{bias_tf.pdf}}

\subfloat[Temperature transfer function]{\label{fig1:tempTF}\includegraphics{bias_tt.pdf}}
\caption{\label{fig:bias_syst} 1D posterior distributions of the six $\Lambda$CDM parameters using a simulated dataset ($TT$, $TE$, $EE$) biased with a constant polarization efficiency~(\ref{fig1:poleff}), with a polarization transfer function~(\ref{fig1:polTF}) and with a temperature transfer function~(\ref{fig1:tempTF}). The transfer function model is described in Eq.~(\ref{eq:tf_model}). We vary the~$\ell_{\rm max}$ parameter from $200$ to $1000$. The dashed black lines correspond to the fiducial model. Gray bands correspond to $\pm 1\sigma$ $\mathcal{R}_\ell^{\rm TE}$ constraints.}
\end{figure*}

\section{A likelihood for the Pearson's correlation coefficient of T and E modes}\label{sec:like}

In this section we recall some of the statistical properties of $\mathcal{R}_\ell^{\rm TE}$ and we propose and validate a multifrequency likelihood for estimating cosmological parameters based on the correlation coefficients.
\begin{figure*}
\includegraphics[scale = 0.9]{sim_check.pdf}
\caption{\label{fig:check_like} $\Lambda$CDM parameters constraints for a set of $N_{\rm sim}=100$ simulations (gray). The average posterior distribution is shown in red. The distance from the average to the input value (in sigma units) is displayed in the top right corner of each panel. Mean posterior distributions are consistent with input values of the simulations}
\end{figure*}
\subsection{Statistic of the multifrequency correlation coefficient}\label{subsec:stat}
An estimator for the correlation coefficient can be obtained using the estimated $TT$, $EE$ and $TE$ power spectra : $\hat{\mathcal{R}}_b^{TE} = {\hat{C}_b^{TE}}/{\sqrt{\hat{C}_b^{TT}\hat{C}_b^{EE}}}$. By definition, the measured $\Tilde{\mathcal{R}}_\ell^{\rm TE}$ correlation coefficient is unaffected by any of the systematic effects we discuss in section~\ref{sec:syst}. (i.e $\Tilde{\mathcal{R}}_\ell^{\rm TE} = \mathcal{R}_\ell^{\rm TE}$).\\

This estimator is affected by a subdominant bias that can be corrected at first order as : $\hat{\mathcal{R}}_b^{\mathrm{TE},c} = \hat{\mathcal{R}}_b^{\mathrm{TE}}(1-\alpha_b)$. An analytical expression for $\alpha_b$ is given in appendix~\ref{app:statistics}. The expression for the covariance matrix of the $\mathcal{R}_b^{\rm TE}$ estimator was also derived in~\cite{Louis2019}, the generalization of the expression to a multifrequency case is given by
\begin{align}\label{eq:covmat}
    \mathbf{\Gamma}(\mathcal{R}_b^{\rm TE}, \mathcal{R}_{b'}^{\rm TE}) &= \mathbf{\Gamma}(C_{b}^{\rm{TE}}, C_{b'}^{\rm{TE}})\nonumber\\
    &+ \frac{1}{4}\big[\mathbf{\Gamma}(C_{b}^{\rm{TT}}, C_{b'}^{\rm{TT}}) + \mathbf{\Gamma}(C_{b}^{\rm{EE}}, C_{b'}^{\rm{EE}})\big]\nonumber\\
    &- \frac{1}{2}\big[\mathbf{\Gamma}(C_{b}^{\rm{TE}}, C_{b'}^{\rm{TT}}) + \mathbf{\Gamma}(C_{b}^{\rm{TT}}, C_{b'}^{\rm{TE}})\nonumber\\
    &\quad+ \mathbf{\Gamma}(C_{b}^{\rm{TE}}, C_{b'}^{\rm{EE}}) + \mathbf{\Gamma}(C_{b}^{\rm{EE}}, C_{b'}^{\rm{TE}})\big]\nonumber\\
    &+ \frac{1}{4}\big[\mathbf{\Gamma}(C_{b}^{\rm{TT}}, C_{b'}^{\rm{EE}}) + \mathbf{\Gamma}(C_{b}^{\rm{EE}}, C_{b'}^{\rm{TT}})\big]
\end{align}
where $\mathbf{\Gamma}(\rm{X}, \rm{Y}) = {\text{cov}(X^{\nu_1\times\nu_2}, Y^{\nu_3\times\nu_4})}/{(X^{\nu_1\times\nu_2}\cdot Y^{\nu_3\times\nu_4})}$. Derivation of Eq.~(\ref{eq:covmat}) is developed in appendix~\ref{app:statistics}.

\subsection{$\mathcal{R}_\ell^{\rm TE}$-based likelihood construction}\label{subsec:Rlike}
We adapt the \texttt{HiLLiPoP} $C_\ell$-based likelihood introduced in section~\ref{sec:syst} in order to fit cosmological parameters from the correlations coefficient. We use the multi-frequency covariance matrix given in Eq.~(\ref{eq:covmat}) and assume that the correlation coefficients follow a multivariate Gaussian distribution.\\

\noindent The likelihood is defined as follows :
\begin{align}\label{eq:like}
    \mathrm{ln}\mathcal{L} \simeq &-\frac{1}{2}\left(\Delta\mathcal{R}^{\mathrm{vec}}\right)^{\mathrm{T}}\mathbf{\Xi}^{-1}\left(\Delta\mathcal{R}^{\mathrm{vec}}\right)
\end{align}
where $\mathbf{\Xi}$ is the $\mathcal{R}_\ell^{TE}$ multi-frequency covariance matrix and $\Delta\mathcal{R}^{\mathrm{vec}} = \mathcal{R}^{\mathrm{vec},\mathrm{data}} - \mathcal{R}^{\mathrm{vec},\mathrm{th}}$ is the residual correlation coefficient vector including the 6 cross-frequency spectra. The $\mathcal{R}_\ell^{\rm TE}$ covariance matrix is computed from the Planck PR4 power spectra covariance matrices using Eq.~(\ref{eq:covmat}).\\
\begin{figure*}
\centering
\includegraphics{RTE_cosmo.pdf}
\caption{\label{fig:cosmo} 2D posterior distributions for $H_0$, $\Omega_bh^2$ and $\Omega_ch^2$ constrained using $C_\ell$-based likelihood (blue) or the $\mathcal{R}_\ell^{TE}$-based likelihood (red). The expected 1$\sigma$ scatter between parameters estimated from the two different likelihoods is displayed with gray dotted lines.}
\end{figure*}

The data vector is computed from the $C_\ell$s data vector using the unbiased estimator $\hat{\mathcal{R}}_b^{\mathrm{TE},c}$. The model for the correlation coefficient is constructed using theoretical and foregrounds $TT$, $EE$, and $TE$ power spectra. For a cross frequency $\nu_1\times\nu_2$, we model the $XY$ power spectra as
\begin{align}\label{eq:model_like}
    C_\ell^{XY, \mathrm{th},\nu_1\times\nu_2} = A_{\rm Pl}A^{XY, \nu_1\times\nu_2}[&C_\ell^{XY, \mathrm{CMB}}(\theta)\nonumber\\ &+ C_\ell^{XY, \mathrm{fg}, \nu_1\times\nu_2}(\theta_{\rm fg})]
\end{align}
where $(X, Y)\in\{T, E\}^2$, $A^{XY, \nu_1\times\nu_2}$ depends on the map calibration parameters and $A_{\rm Pl}$ is a global calibration parameter. $\theta_{\rm fg}$ is the set of parameters describing the foregrounds power spectra amplitudes and $\theta$ are the $\Lambda$CDM parameters. We give more details about Eq.~(\ref{eq:model_like}) in appendix~\ref{app:ps_model}.\\

\begin{table*}[t]
\begin{ruledtabular}
\begin{tabular}{lp{1.5cm}p{1.5cm}p{1.5cm}p{1.5cm}p{1.5cm}p{1.5cm}||p{1.5cm}p{1.5cm}}
 &\multicolumn{2}{c}{$C_\ell^{\rm TT}$}&\multicolumn{2}{c}{$C_\ell^{\rm EE}$}&\multicolumn{2}{c||}{$C_\ell^{\rm TE}$}&\multicolumn{2}{c}{$\mathcal{R}_\ell^{\rm TE}$}\\
Frequencies (GHz) & \hfil$\ell_{\mathrm{min}}$ & \hfil$\ell_{\mathrm{max}}$ & \hfil$\ell_{\mathrm{min}}$ & \hfil$\ell_{\mathrm{max}}$ & \hfil$\ell_{\mathrm{min}}$ & \hfil$\ell_{\mathrm{max}}$ & \hfil$\ell_{\mathrm{min}}$ & \hfil$\ell_{\mathrm{max}}$\\
\hline
    100$\times$100 &  \hfil30 & \hfil1200 & \hfil30 & \hfil1000 & \hfil30 & \hfil1200 & \hfil50 & \hfil1000 \\
    100$\times$143 & \hfil30 & \hfil1500 & \hfil30 & \hfil1250 & \hfil30 & \hfil1500 & \hfil50 & \hfil1250 \\
    100$\times$217 & \hfil100 & \hfil1500 & \hfil400 & \hfil1250 & \hfil300 & \hfil1500 & \hfil400 & \hfil1250 \\
    143$\times$143 & \hfil30 & \hfil2000 & \hfil30 & \hfil1750 & \hfil30 & \hfil1750 & \hfil50 & \hfil1500 \\
    143$\times$217 & \hfil100 & \hfil2500 & \hfil400 & \hfil1750 & \hfil300 & \hfil1750 & \hfil400 & \hfil1500 \\
    217$\times$217 & \hfil100 & \hfil2500 & \hfil400 & \hfil2000 & \hfil300 & \hfil2000 & \hfil400 & \hfil1500\\
\end{tabular}
\end{ruledtabular}
\caption{\label{tab:multipoles}Multipole ranges used for each power spectrum and for the $\mathcal{R}_\ell^{\rm{TE}}$ correlation coefficient.}
\end{table*}

While the approximation of Gaussianity is excellent for CMB power spectra at multipoles $\ell \ge 30$, it is not straightforward that it holds for a non linear combination of power spectra. However, the Gaussianity of the correlation coefficient has been shown to be a good assumption in the high $EE$ signal-to-noise regime~\cite{Louis2019}. We check the robustness of this approximation running MCMC chains for the set of simulations described in section~\ref{sec:syst}. Figure~\ref{fig:check_like} displays the distributions and the mean posterior for each $\Lambda$CDM parameters for $N_{\rm sim} = 100$ simulations. Recovered posterior distributions are statistically consistent with expectations. The use of a $\tau$-prior explains the small dispersion of the posterior distributions for this parameter.

\section{Application to Planck data}\label{sec:data}

In this section, we estimate cosmological parameters using the $\mathcal{R}_\ell^{\rm TE}$-likelihood on Planck PR4 data. We also compute a CMB-only correlation coefficient, by marginalizing over the foreground parameters.

\subsection{Cosmological results}\label{subsec:cosmo}
We use \texttt{cobaya} to sample the likelihoods and derive the posterior of the cosmological and foregrounds parameters. We run MCMC chains on Planck PR4 data using the $C_\ell$-likelihood and the $\mathcal{R}_\ell^{\rm TE}$-likelihood. We use wide flat priors for cosmological parameters, and a Gaussian prior on the optical depth of reionization $\tau = 0.054\pm 0.007$~\cite{Planck2018:cosmo}. The derivation of the statistical properties of the correlation coefficient (cf. appendix~\ref{app:statistics}) relies on a second order development, which is valid only on scales with a high $EE$ signal-to-noise ratio. We apply cuts in multipole to avoid the foreground-dominated scales (low-$\ell$) and noise dominated scales (high-$\ell$). The multipole ranges used in the two likelihoods for the different cross-frequency spectra are displayed in table~\ref{tab:multipoles}.\\

The posterior distributions for $H_0$, and the density parameters $\Omega_bh^2$ and $\Omega_ch^2$ are shown in Fig.~\ref{fig:cosmo}. We display the expected 1 $\sigma$ fluctuations between $C_\ell$-derived parameters and $\mathcal{R}_\ell^{TE}$-derived ones in gray. This was derived from the set of simulations introduced in section~\ref{sec:syst}. Cosmological parameters estimated from the correlation coefficient $\mathcal{R}_\ell^{TE}$ are consistent with the one estimated from power spectra. We do not detect the effect of any multiplicative systematic in the data.\\

\begin{table}
\begin{ruledtabular}
\begin{tabular}{lcc}
Parameter & $C_\ell^{\rm TT}$, $C_\ell^{\rm TE}$, $C_\ell^{\rm EE}$ & $\mathcal{R}_\ell^{\rm TE}$ \\
\hline
$H_0$ [km/s/Mpc] & $67.3 \pm 0.5$ & $67.5 \pm 1.3$\\
$\Omega_bh^2$ & $0.02233 \pm 0.00014$ & $0.02235 \pm 0.00037$\\
$\Omega_ch^2$ & $0.1194 \pm 0.0012$ & $0.1192 \pm 0.0028$\\
$\mathrm{ln}(10^{10}A_s)$ & $3.040 \pm 0.015$ & $3.098 \pm 0.152$
\end{tabular}
\end{ruledtabular}
\caption{\label{tab:results} Cosmological parameter constraints (mean values and $1\sigma$ errors) derived from the $C_\ell$-likelihood and the $\mathcal{R}_\ell^{\rm TE}$-likelihood.}
\end{table}

Table~\ref{tab:results} presents mean values and 1$\sigma$ errors for the cosmological parameters discussed in this paper. As expected, using $\mathcal{R}_\ell^{\rm TE}$ worsens the constraints on cosmological parameters with respect to the $C_\ell$-based likelihood. The $\mathcal{R}_\ell^{\rm TE}$-derived errors on $H_0$ and $\Omega_bh^2$ are $2.6$ times wider than the $C_\ell$-derived ones and the error on $\Omega_ch^2$ is $2.3$ times larger. The correlation coefficient is poorly sensitive to the parameters describing the shape of the initial matter power spectrum, $A_s$ and $n_s$. Interestingly, while the amplitude parameter $A_s$ appears to cancel in the ratio of power spectra, it can be measured through the effect of lensing on the power spectra~\cite{Lewis2006}. The posterior distribution of $\mathrm{log}(10^{10}A_s)$ is displayed in Fig.~\ref{fig:logA}. The $\mathcal{R}_\ell^{\rm TE}$-derived error on $\mathrm{ln}(10^{10}A_s)$ is 10 times larger than the $C_\ell$-derived one. \\

We note that our value of the Hubble parameter $H_0 = 67.5 \pm 1.3 \;\mathrm{km/s/Mpc}$, derived from the correlation coefficient is consistent with other CMB based measurements of $H_0$. The Hubble parameter we obtain using $\mathcal{R}_\ell^{\rm TE}$ is still $3.1\sigma$ discrepant with the latest measurement from Cepheids-calibrated cosmic distance ladder~\cite{Riess2021}.

\begin{figure}[t]
    \centering
    \includegraphics{logA_R.pdf}
    \caption{1D posterior distribution for $\mathrm{log}(10^{10}A_s)$ using the $C_\ell$-based likelihood (blue) or the $\mathcal{R}_\ell^{\rm TE}$-based likelihood (red). The correlation coefficient constraint is sensitive to $A_s$ only through the effect of lensing on the power spectra.}
    \label{fig:logA}
\end{figure}
\subsection{CMB-only $\mathcal{R}_\ell^{TE}$}\label{subsec:cmbonly}
\begin{figure*}
\includegraphics[scale = 0.9]{R_act_planck.pdf}
\caption{\label{fig:cmb_only} (Top) : CMB-only $\mathcal{R}_\ell^{\rm{TE}}$ correlation coefficient for Planck PR4 data (red) and ACT data (green) after foreground marginalization. The best-fit correlation coefficient (gray dashed line) is computed with \texttt{CAMB} from the \texttt{HiLLiPoP} $C_\ell$-likelihood best fit cosmology --- (Bottom) : Residual plot with respect to the binned correlation coefficient best-fit.}
\end{figure*}

We produce CMB-only power spectra, marginalizing over foreground parameters. We follow the method used in~\cite{Dunkley2013}. The multi-frequency vector is modelled as
\begin{equation}\label{eq:cmb_model}
    C_\mathrm{b}^{\mathrm{model}} = \mathbf{A}\cdot C_\mathrm{b}^{\mathrm{CMB}} + C_\mathrm{b}^{\mathrm{fg}}(\theta_{\mathrm{fg}})
\end{equation}
with $\mathbf{A}$ a matrix that projects the CMB bandpowers to their corresponding elements in each cross-frequency spectrum, and $C_\mathrm{b}^{\mathrm{fg}}(\theta_{\mathrm{fg}})$ the foregrounds power spectra corresponding to a set of foreground parameters $\theta_\mathrm{fg}$.\\

Instead of sampling simultaneously CMB bandpowers and nuisance parameters, we split the joint probability distribution $p(C_\mathrm{b}^{\mathrm{CMB}}, \theta_\mathrm{fg}|\mathrm{data})$ into two conditional probability densities $p(C_\mathrm{b}^{\mathrm{CMB}}|\theta_\mathrm{fg},\mathrm{data})$ and $p(\theta_\mathrm{fg}|C_\mathrm{b}^{\mathrm{CMB}}, \mathrm{data})$. The main advantage here is that the CMB bandpowers sampling is simple because $C_\mathrm{b}^\mathrm{CMB}$ follows a multivariate normal distribution with mean $\hat{C}_\mathrm{b}$ and covariance $\mathbf{Q}$ given by

\begin{align}
    &\hat{C}_\mathrm{b} = [\mathbf{A}^\mathrm{T}\mathbf{\Sigma}^{-1}\mathbf{A}]^{-1}[\mathbf{A}^\mathrm{T}\mathbf{\Sigma}^{-1}(C_\mathrm{b}^\mathrm{data} - C_\mathrm{b}^\mathrm{fg}(\theta_\mathrm{fg}))]\label{eq:gibbs}\\
    &\mathbf{Q}^{-1} = \mathbf{A}^\mathrm{T}\mathbf{\Sigma}^{-1}\mathbf{A}\label{eq:covgibbs}
\end{align}
with $\mathbf{\Sigma}$ the $C_\ell$s multifrequency covariance matrix. We alternate CMB bandpowers sampling (using the known probability distribution described in Eqs.~[\ref{eq:gibbs},\ref{eq:covgibbs}]) and Metropolis-Hastings sampling for the nuisance parameters using a Gaussian likelihood. We then compute the correlation coefficient and its covariance matrix using the CMB-only power spectra and the $\mathbf{Q}$ covariance matrix.\\

In Fig.~\ref{fig:cmb_only}, we compare the CMB-only $\mathcal{R}_\ell^{\rm TE}$ with the best-fit model obtained using the $C_\ell$-based likelihood. The correlation coefficient is compatible with the model with $\chi^2/\mathrm{dof} = 52.2/52$ (PTE = 0.47).\\

Higher resolution experiments such as Atacama Cosmology Telescope (ACT) bring more information about the correlation coefficient at small scales. We display on Fig~\ref{fig:cmb_only} the ACT CMB-only $\mathcal{R}_\ell^{\rm TE}$, computed from the power spectra provided by the ACT collaboration~\cite{Choi2020} and publicly available.\footnote{\url{https://github.com/ACTCollaboration/}} Before computing the correlation coefficient, we remove points for which the $EE$ SNR is lower than $3$. We find that the ACT correlation coefficient is in good agreement with the Planck PR4 best-fit, with $\chi^2/\mathrm{dof} = 39.0/36$ (PTE = 0.33).

\section{Conclusion}\label{sec:Conclusion}
Given the current context of tensions between early and late Universe measurements, it is important to assess the robustness of the cosmological parameter constraints. In this work, we have studied the impact of remaining multiplicative systematic effects that would have been neglected in the data treatment, resulting in a bias at the power spectra level. We have shown that this kind of bias can significantly shift the cosmological parameter constraints determined from the combination of $C_\ell^{\rm TT}$, $C_\ell^{\rm TE}$ and $C_\ell^{\rm EE}$.\\

To prevent the cosmological parameter measurements from being biased, we have proposed to use a likelihood based on the Pearson's correlation coefficient $\mathcal{R}_\ell^{\rm TE}$, an observable that is insensitive to multiplicative biases. We have obtained the first constraints on cosmology only from the correlation between $T$ and $E$ modes. We have shown that it gives cosmological parameters that are consistent with other CMB-based measurements. Using $\mathcal{R}_\ell^{\rm TE}$ increases the error on the cosmological parameters with respect to the errors derived from the combination of $TT$, $TE$, and $EE$ power spectra. However, we have measured the Hubble parameter $H_0 = 67.5 \pm 1.3$ km/s/Mpc with an associated error which is similar to the error obtained by ground-based CMB experiments such as ACTPol~\cite{ACT:A20}. The correlation coefficient provides a good consistency check on Planck data. We have not observed the effect of any multiplicative bias in Planck PR4 data and we have obtained a value for the Hubble parameter which is still discrepant with late Universe measurements.\\

Next generation of ground based CMB experiments such as Simons Observatory~\cite{Ade2019} or CMB-S4~\cite{Abazajian2016} will produce data with a high EE signal-to-noise ratio for a wide range of multipoles. This will allow the computation of $\mathcal{R}_\ell^{\rm TE}$ up to smaller scales, increasing its constraining power and relevance. Ground-based experiments are also typically more affected by transfer functions due to atmospheric and ground pick up filtering. We expect that observables such as $\mathcal{R}_\ell^{\rm TE}$ will provide a strong consistency check on these data.\\

\subsection*{Acknowledgements}
We thank Gilles Weymann for useful comments and discussion. The theoretical power spectra used in this paper were computed using \texttt{CAMB} Boltzmann solver~\cite{Lewis2000, Howlett2012}. We gratefully acknowledge the IN2P3 Computer Center (\url{http://cc.in2p3.fr}) for providing the computing resources and services needed to this work.

\onecolumngrid
\begin{appendix}

\section{Statistical properties of $\mathcal{R}_\ell^{\rm TE}$ estimator}\label{app:statistics}

In the high signal-to-noise regime, the correlation coefficient estimator can be developed as :
\begin{align}\label{eq:lin_dev}
    \hat{\mathcal{R}}_b^{\rm TE} &= \mathcal{R}_b^{\rm TE}\frac{1+ \frac{\Delta C_b^{\rm TE}}{C_b^{\rm TE}}}{\sqrt{\left(1+\frac{\Delta C_b^{\rm TT}}{C_b^{\rm TT}}\right)\left(1+\frac{\Delta C_b^{\rm EE}}{C_b^{\rm EE}}\right)}}\nonumber\\
    &= \mathcal{R}_b^{\rm TE}\left(1 + \frac{\Delta C_b^{\rm TE}}{C_b^{\rm TE}}\right)\cdot\left(1 - \frac{1}{2}\frac{\Delta C_b^{\rm TT}}{C_b^{\rm TT}} + \frac{3}{8}\left(\frac{\Delta C_b^{\rm TT}}{C_b^{\rm TT}}\right)^2\right)\cdot\left(1 - \frac{1}{2}\frac{\Delta C_b^{\rm EE}}{C_b^{\rm EE}} + \frac{3}{8}\left(\frac{\Delta C_b^{\rm EE}}{C_b^{\rm EE}}\right)^2\right)
\end{align}
where $\mathcal{R}_b^{\rm TE}$, $C_b^{\rm TT}$, $C_b^{\rm EE}$ and $C_b^{\rm TE}$ are the spectra we want to estimate.\\

Such an estimator is not an unbiased estimator (i.e. $\langle\hat{\mathcal{R}}_b^{\rm TE}\rangle \neq \mathcal{R}_b^{\rm TE}$), but is affected by a bias term such that : $\langle\hat{\mathcal{R}}_b^{\rm TE}\rangle = \mathcal{R}_b^{\rm TE}(1+\alpha_b)$. The bias term $\alpha_b$ can be easily computed taking the mean value of Eq.~(\ref{eq:lin_dev}) :
\begin{align}\label{eq:bias_cc}
    \alpha_b = \frac{3}{8}\left(\frac{\mathrm{cov}(C_b^{\rm TT}, C_b^{\rm TT})}{{C_b^{\rm TT}}^2} + \frac{\mathrm{cov}(C_b^{\rm EE}, C_b^{\rm EE})}{{C_b^{\rm EE}}^2}\right) - \frac{1}{2}\left(\frac{\mathrm{cov}(C_b^{\rm TT}, C_b^{\rm TE})}{C_b^{\rm TT}C_b^{\rm TE}} + \frac{\mathrm{cov}(C_b^{\rm EE}, C_b^{\rm TE})}{C_b^{\rm EE}C_b^{\rm TE}}\right) + \frac{1}{4}\frac{\mathrm{cov}(C_b^{\rm TT}, C_b^{\rm EE})}{C_b^{\rm TT}C_b^{\rm EE}}.
\end{align}

From Eq.~(\ref{eq:lin_dev}) we can compute the $\mathcal{R}_\ell^{TE}$ covariance matrices defined by
\begin{equation}\label{eq:r_covmat}
    \mathrm{cov}(\mathcal{R}_b^{\mathrm{TE}, \nu_1\times\nu_2},
                 \mathcal{R}_{b'}^{\mathrm{TE}, \nu_3\times\nu_4}) =
        \langle\Delta\mathcal{R}_b^{\mathrm{TE}, \nu_1\times\nu_2}
                \Delta\mathcal{R}_b^{\mathrm{TE}, \nu_3\times\nu_4}
        \rangle,
\end{equation}
with $\nu_1\times\nu_2$ and $\nu_3\times\nu_4$ two cross frequencies and $\Delta\mathcal{R}_b^{\rm TE} = \hat{\mathcal{R}}_b^{\rm TE} - \mathcal{R}_b^{\rm TE}$ the deviation to the mean value. $\Delta\mathcal{R}_b^{\rm TE}$ can be expressed analytically (at second order) using Eq.~(\ref{eq:lin_dev})

\section{\texttt{HiLLiPoP} power spectra model}\label{app:ps_model}
In section~\ref{sec:syst} we present the \texttt{HiLLiPoP} $C_\ell$-based likelihood. In section~\ref{sec:like}, we construct a $\mathcal{R}_\ell^{\rm TE}$-based likelihood. In this appendix, we explain how the model in Eq.~(\ref{eq:model_like}) is computed.\\

We work with three frequencies : 100, 143 and 217 GHz (indexed by $\nu_i$, $\nu_j$ for $(i, j)\in\{0,1,2\}^2$) and two split maps (indexed by $R$, $S$ with $(R,S)\in\{A,B\}^2$. We also define $X$ and $Y$ such as $(X,Y)\in\{T, E\}^2$. For two frequencies ($\nu_i$, $\nu_j$) with $j\ge i$ and two maps ($R$, $S$), we model the power spectra as

\begin{equation}\label{eq:ps_map}
    C_\ell^{XY, \mathrm{th}, \nu_i^R\times\nu_j^S} = A_{\rm Pl}c^{\nu_i^R}c^{\nu_j^S}\left[C_\ell^{XY, \mathrm{CMB}} + C_\ell^{XY, \mathrm{fg}, \nu_i\times\nu_j}\right],
\end{equation}
where $A_{\rm Pl}$ is a global amplitude calibration parameter and $c^{\nu_i^R}$, $c^{\nu_j^S}$ are calibration parameters at map level. This set of parameter is sampled in the likelihood as nuisance parameters.\\

We take the weighted average of cross-map power spectra to compute the cross-frequency power spectra, ignoring the auto power spectra.

\begin{equation}\label{eq:ps_freq}
    C_\ell^{XY, \mathrm{th}, \nu_i\times\nu_j} = A_{\rm Pl}\left[\sum_{\substack{(R,S)\in\\\{A,B\}^2}} w_\ell^{XY, \nu_i^R\times\nu_j^S}c^{\nu_i^R}c^{\nu_j^S}\left(C_\ell^{XY, \mathrm{CMB}} + C_\ell^{XY, \mathrm{fg},\nu_i\times\nu_j}\right)\left(1 - \delta_{\nu_i\nu_j}(1 - \delta_{RA}\delta_{SB})\right)\right],
\end{equation}
where $w_\ell^{XY,\nu_i^R\times\nu_j^S}$ are the weights associated to the $\nu_i^R\times\nu_j^S$ $XY$ power spectrum and $\delta$ is the Kronecker delta. We can express Eq.~(\ref{eq:ps_freq}) as Eq.~(\ref{eq:model_like}) using the following definition

\begin{align}\label{eq:fg_freq}
    A^{XY, \nu_i\times\nu_j} = \sum_{\substack{(R,S)\in\\\{A,B\}^2}} w_\ell^{XY, \nu_i^R\times\nu_j^S}c^{\nu_i^R}c^{\nu_j^S}\left(1 - \delta_{\nu_i\nu_j}(1 - \delta_{RA}\delta_{SB})\right).
\end{align}
We obtain the equation describing the model used to compute the power spectra in the likelihood
\begin{equation}
    C_\ell^{XY, \mathrm{th},\nu_i\times\nu_j} = A_\mathrm{Pl}A^{XY, \nu_i\times\nu_j}\left[C_\ell^{XY, \rm CMB} + C_\ell^{XY, \mathrm{fg}, \nu_i\times\nu_j}\right].
\end{equation}

\end{appendix}
\twocolumngrid


\section{introduction}
\label{sec:introduction}
Since the discovery of Hawking radiation \cite{Hawking:1974sw},
quantum black holes have been one of the most intriguing physical objects in quantum theory of gravity.
In the semiclassical approximations,
the total number of particles emitted from black holes
turned out to be the Planckian distribution of thermal radiation
\cite{Hawking:1974sw,Hartle:1976tp,Gibbons:1977mu,Akhmedov:2015xwa}.

On the other hand, there has been much attention to the collapse of the dust shell for
the Schwarzschild black hole \cite{Boulware:1975fe,Gerlach:1976ji,10.1143/PTP.81.826,Alberghi:1998xe,Alberghi:2001cm}
and lower-dimensional black holes \cite{Mann:1992my,Ross:1992ba,Peleg:1994wx,Crisostomo:2003xz,Mann:2006yu,Hyun:2006xt}.
In particular, the functional Schr\"odinger formalism  was employed
to investigate quantum radiation from the collapsing shell for the Schwarzschild metric
\cite{Vachaspati:2006ki,Vachaspati:2007hr,Greenwood:2008zg,Greenwood:2010sx,Kolopanis:2013sty,Saini:2015dea,Saini:2016rrt}.
In fact,  the advantage of the functional Schr\"odinger
equation would be
that wave functionals
living on a space of field configurations explicitly depend on time \cite{Kiefer:1991xy}.
This functional Schr\"odinger formalism
was also applied to the formation of various dynamical geometries such as
the Reissner-Nordstr\"om black hole \cite{Greenwood:2009pd,Das:2019iru},
the black string in anti-de Sitter(AdS) space \cite{Greenwood:2009gp},
the three-dimensional AdS black hole \cite{Saini:2017mur},
and the regular black hole \cite{Um:2019qgc}.
In the formation of the black holes,
the spectrum of quantum radiation is approximately thermal
as the shell approaches its own horizon,
and the temperature of the radiation is in agreement with
the Hawking temperature of the corresponding black holes.
In asymptotically flat black holes formed by collapsing shells,
spacetime for the interior region of the shell
was naturally prescribed by Minkowski spacetime because the interior region is assumed to be empty \cite{Vachaspati:2006ki,Vachaspati:2007hr,Greenwood:2008zg,Greenwood:2009pd,Greenwood:2010sx,Kolopanis:2013sty,Saini:2015dea,Saini:2016rrt,Das:2019iru,Um:2019qgc}.

In the asymptotically non-flat spacetime such as
the nonrotating
Ba\~nados-Teitelboim-Zanelli (BTZ) black hole \cite{Banados:1992wn}
described by the line element
\begin{equation}
\label{eq:metric of BTZ bh}
ds^2=-f(r) dt^2 +\frac{1}{f(r)}dr^2 +r^2 d\phi^2,
\end{equation}
where $f(r) =-M +r^2/\ell^2$
with $M$ and $\ell$ being the mass parameter and the radius of curvature,
respectively, the BTZ black hole is defined for $M>0$ and the massless black hole is for $M=0$.
In addition, the global AdS space of a thermal soliton emerges as a bound state
for $M=-1$ which has a lower energy than that of the massless black hole.
In the formation of the BTZ black hole, the interior region of the shell could be chosen
as the global AdS space ($M=-1$) \cite{Peleg:1994wx,Crisostomo:2003xz} or the massless black hole geometry ($M=0$)
\cite{Saini:2017mur}. For $M=0$,
the radiation temperature was shown to be close to the Hawking temperature by explicit
numerical analysis \cite{Saini:2017mur}.

In this paper,
we will revisit the quantum radiation by solving the functional Schr\"odinger equation
analytically in the formation of the BTZ black hole.
For this purpose,
the exterior and interior regions of the dust shell will be described by the AdS geometry \eqref{eq:metric of BTZ bh}
characterized by different exterior and interior mass parameters of $M_{\rm out}$ and $M_{\rm in}$.
The mass parameters are chosen as
$M_{\rm out} >0$ and $M_{\rm in} =0,-1$.
Namely, the vacuum state for the interior region
is described by  the massless black hole or the global AdS space.
Importantly, in the incipient limit of $R\to R_H$
where $R$ is a shell radius, the occupation number will be expressed by
explicit Bessel functions
and it turns out to be independent of $M_{\rm in}$.
After all, the radiation temperature is found to be close to the Hawking temperature.

The organization of this paper is as follows.
In Sec.~\ref{sec:general},
we will find a mass relation from the junction condition on the dust shell
for arbitrary mass parameters of $M_{\rm out}$ and $M_{\rm in}$.
Then, the functional Schr\"odinger equation from the action for a massless scalar field
will be derived for the formation of the BTZ black hole.
In Sec.~\ref{sec:cases},
the wave functions from the functional Schr\"odinger equations
will be obtained analytically and the corresponding spectra of quantum radiation will be investigated.
Finally, conclusion and discussion will be given in Sec.~\ref{sec:conclusion}.

\section{Schr\"odinger equations for collapsing dust shell}
\label{sec:general}
\subsection{Classical collapse of the shell}
\label{sec:general classical}

In this section, we provide a general formalism for the classical collapse of the shell in AdS space.
Let us start with  spacetime divided by an infinitely thin shell with a radius $R$ in three-dimensional AdS space.
For the exterior and interior coordinates denoted by $\left(t,r,\phi\right)$
and $\left(T,r,\phi\right)$,
the radius would be given as
$R(t)$ which is a decreasing function of the time $t$.

For the exterior region of the shell, the spacetime is given by the line element
    \begin{equation}
    \label{eq:exterior}
    ds_{\rm (out)}^2 = -f_{\rm out}(r) dt^2+\frac{1}{f_{\rm out}(r)} dr^2 +r^2 d\phi^2~~(r>R(t)),
    \end{equation}
where $f_{\rm out}(r) = -M_{\rm out}+r^2/\ell^2$.
For the interior region of the shell, the line element is also described by
	\begin{equation}
    \label{eq:interior}
	ds_{\rm (in)}^2  = -f_{\rm in}(r) dT^2 + \frac{1}{f_{\rm in}(r)} dr^2 +r^2 d\phi^2~~(r<R(t)),
	\end{equation}
where $f_{\rm in} (r) =-M_{\rm in}+r^2/\ell^2$.
Thus, the line element at the shell is assumed to be
    \begin{equation}
    \label{eq:shell}
    ds_{\rm (shell)}^2 =-d\tau^2+ R^2(t) d\phi^2~~(r=R(t))
    \end{equation}
with $\tau$ being the proper time of the collapsing shell.
Since the line element should be continuous at the shell,
Eqs.~\eqref{eq:exterior} and \eqref{eq:interior} are coincident with Eq.~\eqref{eq:shell} at $r=R(t)$,
which gives the relations between the time coordinates $t$ and $T$ as
    \begin{equation}
    \label{eq:time rel}
    \lambda(t)=\frac{dT}{dt}  \bigg|_{r=R(t)}= \sqrt{ \frac{f_{\rm out} (R)}{f_{\rm in} (R)}-
    \frac{1}{f_{\rm in}(R)} \left\{ \frac{1}{f_{\rm out}(R)} -\frac{1}{f_{\rm in}(R)} \right\} \left(\frac{dR}{dt}\right)^2
    },
    \end{equation}
where we did not fix $f_{\rm out}(R)$ and $f_{\rm in}(R)$ for our general argument.

Next, we consider a combined metric tensor defined as
$g_{\mu\nu} =\Theta(r-R(t)) g^{\rm (out)}_{\mu\nu} + \Theta(-r+R(t)) g^{\rm (in)}_{\mu\nu}$
where $\Theta(x)=1 (x>0)$ and $\Theta(x)=0 (x<0)$.
Employing the first junction condition that the combined metric tensor is
continuous on the shell,
one can derive the projection of the Einstein equation onto the shell as \cite{Israel:1966rt}
    \begin{equation}
    \label{eq:Einstein eq}
    -\left[K_{\mu\nu}\right]+h_{\mu\nu} \left[K\right] = 8\pi S_{\mu\nu},
    \end{equation}
where $K_{\mu\nu}=h_\mu^\gamma \nabla_\gamma n_\nu$, $K=g_\mu^\gamma K_{\gamma\nu}$, and $S_{\mu\nu}=h_\mu^\gamma T_{\gamma \nu}$
are the projections of the extrinsic curvature and the energy-momentum tensor of the shell, respectively.
$[...]$ denotes
$[A] = \lim_{r\to R+} A(r) - \lim_{r\to R-} A(r)$.
And $h_{\mu\nu}=(g_{\mu\nu} -n_\mu n_\nu) \big|_{r=R(t)}$
is a projection operator onto the shell with a normal vector of the shell $n_\mu$,
where $h_{\mu\nu}$ automatically satisfies $h_\mu^\gamma h_{\gamma\nu} = h_{\mu\nu}$.
We consider the dust shell and so its energy-momentum tensor is given as $S^{\mu\nu}=\sigma u^\mu u^\nu$, where $\sigma$ is the surface energy density of
the shell and $u^\mu$ is the three-velocity at the shell.
Taking the inner product to Eq.~\eqref{eq:Einstein eq} with $u^\mu$ and $u^\nu$,
one can finally get the so-called mass relation,
    \begin{equation}
    \label{eq:mass rel}
    \sqrt{\dot R^2 + f_{\rm out}(R)}-\sqrt{\dot R^2 + f_{\rm in}(R)} =-8\pi\sigma R,
    \end{equation}
where $\dot R=dR/d\tau$. Note that Eq.~\eqref{eq:mass rel}
certainly shows that
$\sigma>0$ for $M_{\rm out} > M_{\rm in}$
while $\sigma<0$ for $M_{\rm out} < M_{\rm in}$.

Plugging $f_{\rm out}(r)$ and $f_{\rm in}(r)$ in Eqs.~\eqref{eq:exterior} and \eqref{eq:interior} into Eq.~\eqref{eq:mass rel}, one can find
    \begin{equation}
    \label{eq:R dot square}
    \dot R^2 =\frac 1 2 \left(M_{\rm out}+ M_{\rm in} \right)
    +\left(\frac{f_{\rm out}(R)- f_{\rm in} (R)}{16\pi\sigma R} \right)^2
    +\left(\left(4\pi\sigma\right)^2-\frac{1}{\ell^2} \right) R^2,
    \end{equation}
where we supposed that $(4\pi\sigma)^2 > \ell^{-2}$ to make the shell collapse continuously even
for an arbitrary large $R$.
Using the relation of $\left(dR/dt\right)^2 = \left(dt/d\tau\right)^{-2} \dot R^2$
for the collapsing shell, we get
    \begin{equation}
    \label{eq:dR over dt}
    \frac{dR}{dt}=-\frac{f_{\rm out}(R) \dot R}{\sqrt{f_{\rm out}(R)+\dot R^2}}.
    \end{equation}

\subsection{The functional Schr\"odinger equation}
\label{sec:general quantum}
To study the quantum radiation
from the collapsing shell, we consider a minimally coupled scalar field $\Phi$.
On the spacetime background described by Eqs.~\eqref{eq:exterior}, \eqref{eq:interior}, and \eqref{eq:shell},
we consider the total action which consists of the sum of the exterior and interior actions as
    \begin{align}
    \label{eq:total action}
    & S_{\rm (tot)} = S_{\rm (out)} +S_{\rm (in)}, \\ \label{eq:exterior action}
    & S_{\rm (out)} =\pi \int dt \int_{R(t)}^\infty  \!\!\!\!dr~\left(\frac{1}{f_{\rm out}(r)} \left(\partial_t \Phi\right)^2 -f_{\rm out}(r)\left(\partial_r \Phi\right)^2 \right),\\ \label{eq:interior action}
    & S_{\rm (in)} =\pi \int dT \int^{R(t)}_0 \!\!\!\!\!\!dr~ \left(\frac{1}{f_{\rm in}(r)} \left(\partial_T \Phi\right)^2 -f_{\rm in}(r)\left(\partial_r \Phi\right)^2 \right).
    \end{align}
For convenience, let us rewrite the interior action by using Eq.~\eqref{eq:time rel} as
    \begin{equation}
    \label{eq:interior action}
    S_{\rm (in)} =\pi \int dt \int^{R(t)}_0 \!\!\!\!\!\!dr~\left(\frac{1}{\lambda(t) f_{\rm in}(r)} \left(\partial_t \Phi\right)^2 - \lambda(t) f_{\rm in}(r)\left(\partial_r \Phi\right)^2 \right).
    \end{equation}
From Eq.~\eqref{eq:total action}, the functional Schr\"odinger equation will be derived
in the limit of $\lambda(t) \ll 1$,
which will be shown to be of relevance to the incipient limit of the BTZ black hole of $R\to R_H$.

In order to study the formation of the BTZ black hole, we
consider the configuration of $M_{\rm out}> 0$ and $M_{\rm in}=0,-1$.
In particular, for $M=-1$,
nonsingular vortex solutions in AdS$_3$ background were numerically found \cite{Edery:2020kof}.
However, we will consider the global AdS space as a thermal soliton
without any classical fields.
As the radius of the shell approaches its horizon, $i.e.$, $R(t) \to R_H =\sqrt{M_{\rm out}}\ell$ called the incipient limit,
Eq.~\eqref{eq:R dot square} reduces to a finite value as
$\dot R^2 \approx \left.\left(f_{\rm in}(R)/(8\mu)-2\mu\right)^2\right|_{R=R_H}$,
where $\mu$ is the constant of motion defined by $\mu =2 \pi R \sigma$ \cite{Saini:2017mur}.
Thus, Eq.~\eqref{eq:dR over dt} in this limit also reduces to a simpler form as
    \begin{equation}
    \label{eq:dR over dt sim}
    \frac{dR}{dt} \approx -f_{\rm out} (R),
    \end{equation}
which yields the radius function as
    \begin{equation}
    \label{eq:radius function}
    R(t) =
        \begin{cases}
        R_0 & (t<0), \\
        R_H +\left(R_0-R_H \right) e^{-f_{\rm out}'(R_H) t} & ( 0\le t\le t_f),\\
        R_f & (t>t_f)~,
        \end{cases}
    \end{equation}
where $t=0$ and $R_0$ are the initial time and radius, respectively.
For a nice mathematical manipulation of infinite time,
we artificially take the collapse to stop at some time
$t=t_f$ with the radius $R_f$
defined as $R(t_f)=R_H +(R_0 -R_H) e^{-f'_{\rm out} (R_H) t_f}$,
where Eq.~\eqref{eq:radius function} is continuous at $t=t_f$, $i.e.$, $R_f=R(t_f)$.
Then we take the limit of $t_f \to\infty$ to investigate the eternal collapse case
which implies that the shell collapses to form the black hole ($i.e.$, $R(t)\to R_H$),
in other words, $R_f \to R_H$ in Eq.~\eqref{eq:radius function}.
Note that,  if one considered up to the order of $(R-R_H)^2$ in Eq.~\eqref{eq:mass rel},
$M_{\rm in}$ would appear in $(dR/dt)$ in such a way that
$(dR/dt)\approx -f_{\rm out}(R)+f_{\rm out}^2(R)/(2 \dot R^2)$.
However, as long as we are concerned with the lowest order of $(R-R_H)$,
$i.e.$, in the incipient limit, $(dR/dt)$ turns out to be independent of $M_{\rm in}$.
Therefore, in the formation of the BTZ black hole,
the choice of the interior region is not crucial as long as $M_{\rm out} > M_{\rm in}$.

From Eq.~\eqref{eq:dR over dt},
Eq.~\eqref{eq:time rel} can be cast in the incipient limit as
    \begin{equation}

    \label{eq:approx lambda}
    \lambda(t)

    \approx \frac{16\mu^2+f_{\rm in}(R_H)}{\left(16\mu^2-f_{\rm in}(R_H)\right)f_{\rm in}(R_H)} f_{\rm out}(R(t))
    \end{equation}
so that $\lambda(t)$ can be made very small
because $f_{\rm out}\to 0$ near the horizon.
For $\lambda(t)\ll 1$, the kinetic term in the exterior action \eqref{eq:exterior action} and the potential term in the interior action \eqref{eq:interior action}
can be ignored, and thus, the total action \eqref{eq:total action} is simplified as
    \begin{equation}
    \label{eq:tot action sim}
    S_{\rm (tot)}=\pi \int dt \left(\frac{1}{\lambda(t)} \int_0^{R(t)} \!\!\!\!\!\!dr~ \frac{1}{f_{\rm in} (r)} \left(\partial_t \Phi\right)^2-
    \int_{R(t)}^\infty \!\!\!\!dr~ f_{\rm out}(r) \left(\partial_r \Phi\right)^2 \right).
    \end{equation}
After decomposition of $\Phi$ in terms of $\Phi =\sum_{k} a_k (t) d_k (r)$
and the simultaneous diagonalization of the following Hermitian operators \cite{goldstein:mechanics,Vachaspati:2006ki},
$A_{k,k'} =2\pi\int_0^{R(t)} dr~ d_k(r) d_{k'}(r)/f_{\rm in} (r)$ and
$B_{k,k'} =2\pi\int_{R(t)}^\infty dr~f_{\rm out}(r)(dd_k (r)/dr)(dd_{k'} (r)/dr)$,
we obtain the total action \eqref{eq:tot action sim} as
    \begin{eqnarray}
    \notag
    S_{\rm (tot)} & =&\int dt L_\Phi \\
    \label{eq:tot action diag}
    &=& \sum_k \int dt \left(\frac{\alpha_k}{2\lambda(t)} \left(\frac{d b_k(t)}{dt}\right)^2 -\frac 1 2 \beta_k b_k^2(t) \right),
    \end{eqnarray}
where $\alpha_k$ and $\beta_k$ are eigenvalues of $A_{k,k'}$ and $B_{k,k'}$
for a set of eigenvectors $\{b_k\}$, respectively.
For a mode amplitude $b_k(t)$,
the conjugate momentum is defined as $\Pi_k=\partial L_\Phi/\partial \left(d b_k(t)/dt\right)$
$=\alpha_k \lambda^{-1}(t) \left(db_k (t)/dt\right)$.
Then, from the Legendre transformation of Eq.~\eqref{eq:tot action diag}, the Hamiltonian can easily be obtained
as $H_\Phi=\sum_k \left\{\lambda(t)/\left(2\alpha_k\right) \Pi_k^2+(1/2) \beta_k b_k^2\right\}$.
Imposing the quantization rule $\left[b_k(t),\Pi_{k'}\right]=i\delta_{k,k'}$
where $\Pi_k \to -i\partial/\partial b_k$, we get the Hamiltonian for quantum radiation as
    \begin{equation}
    H_\Phi = \sum_k H_{\Phi,k}
    =\sum_k \left\{-\frac{\lambda(t)}{2\alpha_k} \frac{\partial^2}{\partial b_k^2}+\frac 1 2 \beta_k b_k^2\right\}.
    \end{equation}
It is worthy to note that the Hamiltonian is separable for each mode $k$
so that the functional Schr\"odinger equation of $H_\Phi \Psi= i \partial \Psi /\partial t$ \cite{Vachaspati:2006ki,Vachaspati:2007hr,Greenwood:2008zg,Greenwood:2009gp,Greenwood:2009pd,Greenwood:2010sx,Kolopanis:2013sty,Saini:2015dea,Saini:2016rrt,Saini:2017mur,Das:2019iru,Um:2019qgc}
can be rewritten as the family of the equation of $H_{\Phi,k} \psi_k= i \partial \psi_k /\partial t$
for a certain mode $k$,
where $\Psi = \prod_k \psi_k (t,b_k)$.
Therefore, for one eigenvector $b \in\{b_k\}$ for which the wave function $\psi (t,b) \in \{\psi_k(t,b_k)\}$,
the functional Schr\"odinger equation is written as \cite{Vachaspati:2006ki}
    \begin{equation}
    \label{eq:Sch eq eta}
    \left\{-\frac{1}{2\alpha} \frac{\partial^2}{\partial b^2}
    +\frac 1 2 \alpha \frac{\omega^2_0}{\lambda(\eta(t))} b^2 \right\}\psi(\eta,b)= i\frac{\partial \psi}{\partial \eta} (\eta, b)
    \end{equation}
in terms of a time parameter $\eta = \int_0^t d t'~ \lambda(t')$,
where $\omega_0 = \sqrt{\beta / \alpha }$
with $\alpha\in\{\alpha_k \}$ and $\beta\in\{\beta_k \}$.
Using Eqs.~\eqref{eq:dR over dt sim} and \eqref{eq:approx lambda},
$R(\eta)$ can also be written in the form of
    \begin{equation}
    \label{eq:Radius eta}
    R(\eta) =\left\{
    \begin{matrix}
        R_0 & (\eta<0), \\
        R_0 -\dfrac{\left(16\mu^2-f_{\rm in}(R_H)\right)f_{\rm in}(R_H)}{16\mu^2+f_{\rm in}(R_H)}\eta & ( 0\le\eta\le \eta_f),\\
        R_f & (\eta> \eta_f),~
    \end{matrix}
    \right.
    \end{equation}
which will be used in Sec.~\ref{sec:cases}
in order to solve the functional Schr\"odinger equation \eqref{eq:Sch eq eta} analytically.

\section{Quantum radiation from the collapsing shell}
\label{sec:cases}
\subsection{Exact solution of the functional Schr\"odinger equation}
\label{sec:exact sol BH}

The functional Schr\"odinger equation \eqref{eq:Sch eq eta}
can be used to study the quantum radiation of the collapsing BTZ black hole.

For $\eta<0$,
the radius \eqref{eq:Radius eta} is constant, and thus, the wave function $\psi(\eta,b)$ can be given by
the linear combination of the eigenstates of the simple harmonic oscillator as
    \begin{equation}
    \label{eq:eigenstate SHO}
    \varphi_n (b) =\left( \frac{\alpha \bar \omega}{\pi}\right)^{\frac 1 4} \frac{1}{\sqrt{2^n n!}} H_n (\sqrt{\alpha \bar \omega} b) e^{-\frac 1 2 \alpha \bar \omega b^2},
    \end{equation}
where $H_n (x)$ is a Hermite polynomial and
$\bar \omega=\omega_0$.
Starting from a vacuum state represented by a quantum state with the lowest energy,
we set $\psi(\eta,b)=e^{-i \omega_0 \eta} \varphi_0(b)\big|_{\bar \omega =\omega_0}$.

For $0\leq\eta \leq\eta_f$,
let us obtain the exact solution to satisfy Eq.~\eqref{eq:Sch eq eta}.
By plugging $R(\eta)$ in Eq.~\eqref{eq:Radius eta} into Eq.~\eqref{eq:approx lambda}, $\lambda(t)$ in Eq.~\eqref{eq:approx lambda} can be rewritten as
    \begin{equation}
    \label{eq:lambda eta}
    \lambda(\eta)=\frac{d\eta}{dt}
    =\frac{f'_{\rm out}(R_H)}{H} \left(1-H\eta\right)
    \end{equation}
where
$H=\{\left(16\mu^2+f_{\rm in}(R_H)\right)(R_0-R_H)\} / \{\left(16\mu^2-f_{\rm in}(R_H)\right)f_{\rm in}(R_H)\}$.
Using Eq.~\eqref{eq:lambda eta} with the following redefined relations of
$\hat \alpha = H\alpha$, $\hat \omega_0 = \omega_0 /\sqrt{f'_{\rm out}(R_H) H}$,
and $\hat \eta=H \eta$,
one can simplify the functional Schr\"odinger equation in the incipient limit as \cite{Kolopanis:2013sty}
    \begin{equation}
    \label{eq:Sch eq eta hat}
    \left\{-\frac{1}{2\hat \alpha} \frac{\partial^2}{\partial b^2} +
    \frac 1 2 \hat \alpha \frac{\hat \omega_0^2}{1-\hat \eta}\right\} \psi=
    i\frac{\partial\psi}{\partial \hat\eta}.
    \end{equation}
The exact solution to Eq.~\eqref{eq:Sch eq eta hat}
 is known as \cite{Dantas:1990rk}
    \begin{equation}
    \label{eq:solution}
    \psi(\hat \eta,b)=e^{i \theta(\hat \eta)} \left(\frac{\hat \alpha}{\pi \zeta^2}\right)^{\frac 1 4} e^{i\left(\frac{1}{\zeta} \zeta_{\hat \eta}+\frac{i}{\zeta^2}\right) \frac{\alpha b^2}{2}},
    \end{equation}
where $\theta(\hat \eta)=(-1/2)\int_0^{\hat \eta} d \hat \eta ' \zeta^{-2} (\hat \eta ')$.
Note that $\zeta(\hat \eta)$ is the solution to the inhomogeneous nonlinear equation of
    \begin{equation}
    \label{eq:soluble diff eq}
    \frac{d ^2 \zeta}{d \hat \eta^2}+\frac{\hat \omega_0 ^2}{1-\hat \eta} \zeta=\frac{1}{\zeta^3}
    \end{equation}
with the initial conditions: $\zeta(0)=1/\sqrt{\hat \omega_0}$ and
$\zeta_{\hat \eta}(0)=0$,
where $\zeta_{\hat \eta}=d\zeta/d\hat \eta$.
In the previous work for the formation of the BTZ black hole with $M_{\rm in}=0$
\cite{Saini:2017mur},
Eq.~\eqref{eq:soluble diff eq} was solved numerically.
Instead, we will find the analytic expression for Eq.~\eqref{eq:soluble diff eq}.

Let us now suppose that the solution satisfying Eq.~\eqref{eq:soluble diff eq} is written as \cite{Kolopanis:2013sty}
    \begin{equation}
    \label{eq:solution1}
    \zeta(\hat \eta)=\frac{1}{\sqrt{\hat \omega_0}} \sqrt{\xi^2 (\hat \eta)+\chi^2 (\hat \eta) }
    \end{equation}
with the appropriate initial conditions:
$\xi(0)=1,~\xi_{\hat \eta}(0)=0,~\chi(0)=0,~\chi_{\hat \eta} (0)=\hat \omega_0$,
and its derivative with respect to $\hat \eta$ is also written as
$\zeta_{\hat \eta} =\left(\xi \xi_{\hat \eta}+\chi \chi_{\hat \eta}\right)/(\hat \omega_0 \zeta)$.
In Eq.~\eqref{eq:soluble diff eq}, the two linearly independent homogeneous solutions
denoted by $y_1(\hat \omega_0,\hat \eta)$ and $y_2(\hat \omega_0,\hat \eta)$ are obtained as
the Bessel function of the first kind $J_1 (2\hat \omega_0 \sqrt{1-\hat \eta})$ and
the Bessel function of the second kind $Y_1 (2\hat \omega_0 \sqrt{1-\hat \eta})$, respectively \cite{10.5555/1098650}.
Then, we can construct $\xi(\hat \eta)$ and $\chi(\hat \eta)$ in Eq.~\eqref{eq:solution1} as
    \begin{eqnarray}
    \label{eq:wrons1}
    \xi(\hat \omega_0,\hat \eta) =& W^{-1}
    \left( y_{2\hat \eta}(\hat \omega_0,0) y_1 (\hat \omega_0,\hat \eta)\right.
    \left. -y_{1\hat \eta}(\hat \omega_0,0) y_2 (\hat \omega_0,\hat \eta)\right),\\
    \chi(\hat \omega_0,\hat \eta) =& -\hat \omega_0 W^{-1}
    \left( y_{2}(\hat \omega_0,0) y_1 (\hat \omega_0,\hat \eta)\right.
    \left. -y_{1}(\hat \omega_0,0) y_2 (\hat \omega_0,\hat \eta)\right),\\
    \xi_{\hat \eta} (\hat \omega_0,\hat \eta) = &W^{-1}
    \left( y_{2\hat \eta}(\hat \omega_0,0) y_{1\hat \eta} (\hat \omega_0,\hat \eta)\right.
    \left. -y_{1\hat \eta}(\hat \omega_0,0) y_{2 \hat \eta} (\hat \omega_0,\hat \eta)\right),\\
    \label{eq:wrons4}
    \chi_{\hat \eta} (\hat \omega_0,\hat \eta) = &-\hat \omega_0 W^{-1}
    \left( y_2(\hat \omega_0,0) y_{1\hat \eta} (\hat \omega_0,\hat \eta)\right.
    \left. -y_1(\hat \omega_0,0) y_{2 \hat \eta} (\hat \omega_0,\hat \eta)\right),
    \end{eqnarray}
where $W$ is the Wronskian at the initial time calculated as $W=y_{2\hat \eta}(\hat\omega_0,0) y_1(\hat\omega_0,0)-
y_{1\hat \eta}(\hat\omega_0,0) y_2(\hat\omega_0,0)=-\pi^{-1}$.
From Eqs.~\eqref{eq:wrons1}-\eqref{eq:wrons4},
we explicitly get
    \begin{eqnarray}
    \label{eq:sol1}
    \xi (x,z) =& \phantom{-}\dfrac{\pi z}{2}&\left(Y_0 (x) J_1(z)-J_0(x) Y_1(z) \right),\\ \label{eq:sol2}
    \chi (x,z) =& \phantom{-}\dfrac{\pi z}{2}&\left(Y_1 (x) J_1(z)-J_1(x) Y_1(z) \right),\\ \label{eq:sol3}
    \xi_{\hat \eta} (x,z) =& -\dfrac{ \pi x^2}{4} &\left(Y_0 (x) J_0(z) -J_0 (x) Y_0 (z)\right),\\ \label{eq:sol4}
    \chi_{\hat \eta} (x,z)  =& -\dfrac{ \pi x^2}{4} &\left(Y_1 (x) J_0(z) -J_1 (x) Y_0 (z)\right),
    \end{eqnarray}
where the new variables are defined by
    \begin{equation}
    \label{eq:x and z BH}
    x=2\hat \omega_0,\qquad z=2\hat \omega_0 \sqrt{1-\hat \eta}.
    \end{equation}
Consequently, we obtain the analytic solution $\zeta(\hat \eta)$ and its derivative $\zeta_{\hat \eta}$,
    \begin{equation}
    \label{eq:sol tot}
    \zeta(x,z)=\sqrt{\frac 2 x \big(\xi^2 +\chi^2\big)},\quad
    \zeta_{\hat \eta}(x,z)=\sqrt{\frac 2 x}~ \frac{\xi \xi_{\hat \eta}+\chi \chi_{\hat \eta}}{\sqrt{\xi^2 +\chi^2 }}
        \end{equation}
in terms of the combinations of the Bessel functions.

Finally, if $\eta>\eta_f$, then $R(\eta)= R_f$,
and thus, Eq.~\eqref{eq:Sch eq eta} becomes the Schr\"odinger equation of the time-independent simple harmonic oscillator
so that the wave function is simply given as
$\psi(\eta,b)=\left(\alpha \omega_f /\pi\right)^{\frac 1 4} e^{-i \omega_f \eta_f}  e^{-\frac 1 2 \alpha  \omega_f b^2}$,
where $\omega_f=\omega_0/\sqrt{\lambda(\eta_f)}$.

\subsection{Spectrum of the quantum radiation}
\label{sec:spectrum BH}
In order to study the quantum radiation from the shell,
we should calculate the occupation number defined by
    \begin{equation}
    \label{eq:occ num def}
    N=\sum_{n=0}^\infty n |c_n (\eta) |^2,
    \end{equation}
where the probability amplitude in the $n$th eigenstate \eqref{eq:eigenstate SHO} for the frequency
$\bar \omega =\omega_f$
can be calculated as
    \begin{align}
    c_{n} (\eta)
    &= \int_{-\infty}^{\infty} \textrm d b~ \varphi_n (b)  \psi(\eta,b)\\
    &=\left\{
        \begin{array}{ll}
        e^{i \theta(\hat \eta)} \frac{(n-1)!!}{\sqrt{n!}} \frac{1}{\sqrt[4]{\hat \omega_f \zeta^2 } }
    \sqrt{(-1)^n \frac 2 P \left(1-\frac 2 P\right)^n}  & \qquad (n=0,2,4,...), \\
        0 & \qquad (n=1,3,5,...),\label{eq:prob amp}
        \end{array}
    \right.
    \end{align}
where $P=1-i(\zeta_{\hat \eta}/\hat \omega_f \zeta) + 1/(\hat \omega_f \zeta^2)$ with $\hat \omega_f = \hat \omega_0 /\sqrt{1-H \eta_f}$.
After some calculations, we obtain
    \begin{equation}
    \label{eq:occ num formal}
    N(\hat \eta_f,\hat \omega_0)=\frac{\hat \omega_f \zeta^2}{4}
    \left\{\left(1-\frac{1}{\hat \omega_f \zeta^2}\right)^2
    +\left(\frac{\zeta_{\hat \eta}}{\hat \omega_f \zeta}\right)^2\right\}.
    \end{equation}
Using the exact
expression \eqref{eq:sol tot},
we can rewrite Eq.~\eqref{eq:occ num formal} as
    \begin{equation}
    \label{eq:occ num BH finite}
    N(x,z) = \left\{\frac x z (\xi^2 +\chi^2)\right\}^{-1}
    \left[~\frac 1 4 \!\left(\frac{x}{z} (\xi^2 +\chi^2)-1\right)^2
    +\left(\frac{\xi\xi_{\hat \eta}+\chi\chi_{\hat \eta}}{x}\right)^2\right].
    \end{equation}
By integrating Eq.~\eqref{eq:lambda eta},
the relation between $t$ and $\eta$ in the incipient limit can easily be obtained as
$1-H\eta = e^{-f'_{\rm out}(R_H)\phantom{`}t}$
and $\Omega_f$ is expressed by
$\Omega_f = \lambda(t_f) \omega_f = (f'_{\rm out}(R_H)/H) e^{-f'_{\rm out} (R_H)}\omega_f$.
In this limit, Eq.~\eqref{eq:x and z BH} can also be rewritten in terms of $t_f$ and $\Omega_f$ as
    \begin{equation}
    \label{eq:x and z}
    x=\frac{2\Omega_f}{f'_{\rm out} (R_H)} e^{\frac 1 2 f'_{\rm out}(R_H) t_f},\quad
    z=\frac{2\Omega_f}{f'_{\rm out} (R_H)}.
    \end{equation}
Thus, the occupation number becomes the functions of $t_f$ and
$\Omega_f$, {\it i.e.}, $N(x,z)=N(t_f,\Omega_f)$.
\begin{figure}[b]
\centering
\includegraphics[width=0.6\textwidth]{N_BH.jpg}
\caption{The spectrum of the occupation number is plotted
with respect to the dimensionless frequency $\Omega_f \ell$
for two different times.
The black and grey curves are for a finite time \eqref{eq:occ num BH finite} ($t_f/\ell=1.5$) and
for infinite time \eqref{eq:occ num BH infinite}, respectively.
Although the spectrum for the finite time is damped oscillating, the spectrum turns out to be
approximately thermal
for $t_f \to \infty$ when the shell approaches the horizon.
}\label{fig:black hole}
\end{figure}

Note that
the functional Schr\"odinger equation \eqref{eq:Sch eq eta}
tells us that
the wave function is apparently responsible for the exterior and interior geometries through Eq.~\eqref{eq:approx lambda};
however, the occupation number $N(t_f,\Omega_f)$ in the incipient limit is of no relevance to the interior geometry
because $x$ and $z$ are independent of $f_{\rm in}$ as seen from Eq.~\eqref{eq:x and z}.
For the asymptotic observer, the shell does not collapse to $R_H$ in a finite time, but the collapsing shell takes infinite time to reach $R_H$.
So, we should take the limit of $t_f \to \infty$ to study the quantum radiation
when the black hole forms.
For $t_f \to\infty$, we finally obtain the occupation number \eqref{eq:occ num BH finite} as
    \begin{equation}
    \label{eq:occ num BH infinite}
    N_{\infty} (z)= \frac \pi 8 \frac{z}{J_1^2(z)+Y_1^2(z)}
    \left[
    \left(J_1^2(z)+Y_1^2(z)-\frac{2}{\pi z}
    \right)^2
    +\left(J_0(z) J_1(z)+Y_0(z) Y_1(z)\right)^2
    \right],
    \end{equation}
where we used the following asymptotic forms of Bessel functions
    \begin{eqnarray}
    & \dfrac x z \left(\xi^2 +\chi^2\right) \to \dfrac{\pi z}{2}\left(J_1^2(z)+Y_1^2(z)\right), &\\
    & \dfrac{\xi\xi_{\hat \eta}+\chi\chi_{\hat \eta}}{x} \to -\dfrac{\pi z}{4}\left(J_0(z) J_1(z)+Y_0(z) Y_1(z)\right).&
    \end{eqnarray}

For convenience, the spectrum of the occupation number is plotted in Fig.~\ref{fig:black hole}.
In a finite time, the radiation is nonthermal
since the spectrum is quite different from the
Planckian distribution.
When the time goes to infinity,
the spectrum approaches the Planckian distribution.
In particular, in a low frequency region,  Eq.~\eqref{eq:occ num BH infinite} is written as
    \begin{equation}
    \label{eq:Hawking temp}
    N_\infty (\Omega_f)\approx \frac{f'_{\rm out}(R_H)}{4\pi \Omega_f}
    =\frac{\sqrt{M_{\rm out}}}{2\pi\ell~\Omega_f}.
    \end{equation}
Hence, the temperature of the quantum radiation can be read off
from Eq.~\eqref{eq:Hawking temp} by comparing it to
the Planckian distribution in a low frequency as
\begin{equation}
\label{eq:temperature result}
T=\frac{\sqrt{M_{\rm out}}}{2\pi\ell}.
\end{equation}
For $M_{\rm out}> 0$ and $M_{\rm in}=0$, the present analytic account
is compatible with the numerical result~\cite{Saini:2017mur}
that the temperature of the emitted radiation
is very close to the Hawking temperature
as the shell
approaches its own horizon.
Of course, Eq.~\eqref{eq:temperature result} is also valid for $M_{\rm out}> 0$ and $M_{\rm in}=-1$.

\section{conclusion and discussion}
\label{sec:conclusion}
In conclusion, we considered the collapsing dust shell where
the exterior and interior regions of the shell were
described by the AdS geometry
characterized by the different exterior and interior mass parameters of $M_{\rm out}$ and $M_{\rm in}$,
and studied the quantum radiation from the collapsing shell
by using the functional Schr\"odinger formalism.
In the formation of the BTZ black hole where $M_{\rm out} >0$ and $M_{\rm in}=0,-1$,
the wave function in the incipient BTZ black hole limit was fortunately expressed by the analytic form and
so the occupation number of the excited states was calculated exactly.
The spectrum of the occupation number is nonthermal for a finite time, but it is getting close to
the thermal Hawking radiation as time goes to infinity.

Until now, we have studied the collapse of the BTZ black hole for which $M_{\rm out} >0$.
Let us speculate what happens for the collapse where $M_{\rm out} <0$.
As a matter of fact, for $-1< M_{\rm out}< 0$,
the end points of collapse depend on the initial velocity and position, and
can lead to a naked
singularity \cite{Peleg:1994wx,Crisostomo:2003xz,Mann:2006yu}.
To evade the naked singularity attributed to the exterior geometry,
we can instead consider the collapse of the global AdS space as a thermal soliton of $M_{\rm out} =-1$ without the horizon.
Firstly, assuming that the surface energy density $\sigma$ of the shell is positive
in the formation of the global AdS space, then the interior mass parameter should be chosen as $M_{\rm in} <-1$
from Eq.~\eqref{eq:mass rel} so that the interior geometry leads necessarily to the naked singularity.
Of course, if the surface energy density is negative, then
the interior mass parameter can be $M_{\rm in}=0$ for a particular value of $\mu$
where $\mu$ is the constant of motion defined by $\mu =2 \pi R \sigma$; however,
the negative energy density would be unnatural in the classical collapse.
Secondly, the spectrum of the occupation number will appear to be nonthermal
at very late times
when the radius of the shell shrinks to $R \to 0$.
This fact might be expected from the previous work for the collapse of a regular black hole in the absence of the horizon \cite{Um:2019qgc}.
Eventually, the occupation number would not be well-defined when encountering the naked singularity,
which deserves further study.

\acknowledgments
This research was supported by Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by
the
Korea government(MSIP) (NRF-2017R1A2B2006159) and
the Ministry of Education through the Center for Quantum Spacetime (CQUeST) of Sogang University (NRF-2020R1A6A1A03047877).

\bibliographystyle{JHEP}