% ****** Start of file apssamp.tex ******
%
%   This file is part of the APS files in the REVTeX 4.2 distribution.
%   Version 4.2a of REVTeX, December 2014
%
%   Copyright (c) 2014 The American Physical Society.
%
%   See the REVTeX 4 README file for restrictions and more information.
%
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
% as well as the rest of the prerequisites for REVTeX 4.2
%
% See the REVTeX 4 README file
% It also requires running BibTeX. The commands are as follows:
%
%  1)  latex apssamp.tex
%  2)  bibtex apssamp
%  3)  latex apssamp.tex
%  4)  latex apssamp.tex
%
\documentclass[%
 reprint,
%superscriptaddress,
%groupedaddress,
%unsortedaddress,
%runinaddress,
%frontmatterverbose, 
%preprint,
%preprintnumbers,
%nofootinbib,
%nobibnotes,
%bibnotes,
 amsmath,amssymb,
 aps,
%pra,
%prb,
%rmp,
%prstab,
%prstper,
%floatfix,
]{revtex4-2}

\usepackage{graphicx}
\usepackage{dcolumn}
\usepackage{bm}
\usepackage{amsmath,mathrsfs,amssymb}
\usepackage{subfigure}
\usepackage{float}
\usepackage{xcolor}
\usepackage{booktabs} 
\usepackage{adjustbox} 
\usepackage{array} 
\usepackage{siunitx} 
\begin{document}

\preprint{APS/123-QED}

\title{Black hole echos reflect the phase transition and fluctuations in Hawking radiation}


\author{
  Tianqi Yue$^{1}$ and 
  Jin Wang$^{2,}$
}
\email{Corresponding author:jin.wang.1@stonybrook.edu}
\affiliation{
  $^{1}$College of Physics, Jilin University, Changchun 130022, China \\
  $^{2}$Department of Chemistry, and Department of Physics and Astronomy, 
  Stony Brook University, Stony Brook, New York 11794, U.S.A.
}



%\date{\today}%

\begin{abstract}
Black holes are thermal objects. They can form thermodynamic phases and exhibit phase transitions. Furthermore, black holes can also radiate, termed as Hawking radiation. However, the signatures of these behaviors are challenging to observe. In this work,we consider Hawking radiation in black hole phase transitions. We uncovered that an echo can emerge from the correlations between individual single event and joint two events. This provides possible signature of black hole phase transition and fluctuations in Hawking radiation.
\end{abstract}




\maketitle
Black holes are established as thermodynamic systems with entropy and temperature via Hawking radiation \cite{hawking1975particle} and Bekenstein entropy \cite{bekenstein1973black}, leading to thermodynamics phases and phase transitions including the Hawking-Page transition in AdS spacetime \cite{hawking1983thermodynamics}. For RNAdS black holes, these transitions exhibit van der Waals-like behavior when the cosmological constant is treated as pressure \cite{Kastor2009,Dolan2011,Dolan2011Pressure}, with criticality analyzed in \cite{Kubiznak2012}. Recent kinetic approaches model transitions using free energy landscape (radius as order parameter) \cite{whiting1988action,Li2020,li2020thermal,li2022generalized}, where thermal fluctuations drive stochastic transitions between metastable states. To identify observable signatures of phase transitions and Hawking radiation, we investigate echo signals arising from event correlation differences within this free energy landscape framework incorporating evaporation effects \cite{yang2001two,cao2000event}. Our analysis of single and joint event probabilities reveal echo behavior through a correlation difference function—providing a potential probe of both phase transition dynamics and Hawking radiation.


\begin{figure*}[!ht]
    \centering
    \subfigure[]{
    \includegraphics[width=0.2\linewidth]{fg/Full_Reaction.jpg}\label{fig:Full_Reaction}
    }
    \subfigure[]{
    \includegraphics[width=0.15\linewidth]{fg/halfreactionA.jpg} \label{fig:halfreactionA}
    }
    \subfigure[]{
    \includegraphics[width=0.15\linewidth]{fg/halfreactionB.jpg } \label{fig:halfreactionB}
    }
     \subfigure[]{
    \includegraphics[width=0.16\linewidth]{fg/landscapewithreactionpath.jpg}\label{fig:landscapewithreactionpath}
    }
    \subfigure[]{
    \includegraphics[width=0.16\linewidth]{fg/locallanscape.jpg} \label{fig:locallanscape}
    }
    \caption{ Schematic of the reaction model and corresponding free energy landscape. (a) Full-reaction network. (b,c) Half-reactions A and B. (d,e) Landscapes with two basins. Stochastic trajectories within a basin represent local transitions, governed by the half-reaction dynamics and quantified by the Green's function.}
\end{figure*}


Consider a large system with two macrostates, each of which contains numerous substates. Assuming that these states are associated with certain probabilities, the dynamic equation that governs these probability evolutions is known as the master equation.
\begin{equation}
\begin{pmatrix}\label{master equation}
   \Dot{\rho}_a\\
   \Dot{\rho}_b 
\end{pmatrix}=
\begin{pmatrix}
   -K_a-K_{AB}&&&K_{BA}\\
   K_{AB} &&&-K_b-K_{BA}
\end{pmatrix}
\begin{pmatrix}
     \rho_a\\
   \rho_b
\end{pmatrix}
\end{equation}
This full-reaction system, comprising states A and B with internal substates, is decomposed into two half-reactions. The dynamics are governed by kinetic rate matrices (referring to Fig.\ref{fig:Full_Reaction},\ref{fig:halfreactionA},\ref{fig:halfreactionB}): \(K_{AB}\) and \(K_{BA}\) describe transitions between A and B, while \(K_a\) and \(K_b\) account for internal substate transitions. The master equations for each half-reaction are \(\dot{\rho}_a = (-K_a - K_{AB})\rho_a\) and \(\dot{\rho}_b = (-K_b - K_{BA})\rho_b\), with corresponding Green’s functions \(G_a(t) = e^{(-K_a - K_{AB})t}\) and \(G_b(t) = e^{(-K_b - K_{BA})t}\). At steady state, the probability flux is conserved, yielding normalized fluxes \(F_a = N_1^{-1} K_{BA} \rho_b\) and \(F_b = N_2^{-1} K_{AB} \rho_a\), where \(N_1 = N_2\). The single-event distribution for switching from A to B is \(f_a(t) = \sum K_{AB} G_a(t) F_a\), and similarly \(f_b(t) = \sum K_{BA} G_b(t) F_b\) for B to A. Joint distributions for consecutive transitions are given by \(f_{ab}(t_1, t_2) = \sum K_{BA} G_b(t_2) K_{AB} G_a(t_1) F_a\) and \(f_{ba}(t_1, t_2) = \sum K_{AB} G_a(t_2) K_{BA} G_b(t_1) F_b\). These expressions capture the probability fluxes and temporal evolution between and within states, with the Green’s functions playing a central role in governing the dynamical behavior. The probability distribution has clear physical meaning: for half-reaction A, \( G_a(t) \) governs the time evolution of probability density, while \( F_A \) represents substate fluxes. Evolution from \( A_0 \) under \( G_a \) gives \( A_1 \). Summing over substate transitions \( K_{AB,i} \) yields the total A→B probability. Other distributions follow analogously, highlighting the essential role of \( G_a(t) \) in dynamics.




To quantify the degree of correlation between individual single switching event and joint two switching event, we define a function known as the difference function.
\begin{equation}
    \delta(t_1,t_2) = \left|f_{ab}(t_1,t_2) - f_a(t_1) f_b(t_2)\right|
    \label{delta}
\end{equation}
\(\delta\) quantifies the correlation between two switching events. A larger \(\delta\) indicates stronger correlation in the sequence \(A \xrightarrow{G_a(t_1)} B \xrightarrow{G_b(t_2)} A\), analogous to a matter wave echoing from A to B and back. The echo time corresponds to the extremum of \(\delta\), marking the point of highest correlation where the second event is most influenced by the first.




In the case of an RNAdS black hole, the metric has the following form,
\(ds^2 = -f(r)dt^2 + f(r)^{-1}dr^2 + r^2 d\Omega^2\), where \(f(r)=1 - \frac{2M}{r} + \frac{Q^2}{r^2} + \frac{r^2}{L^2}\), \(M\) is the mass, \( Q \) is the charge, and \( L =\sqrt{-3/\Lambda} \) is the AdS curvature radius (\( \Lambda\) being the cosmological constant). RNAds black hole can have thermodynamic phase transitions \cite{Chamblin1999, Chamblin1999Holography, Wu2000}. To explore this phase transition process, we can define the generalized free energy landscape for the black hole state in order parameter \(r_+\), which is given by \cite{whiting1988action,Li2020, li2020thermal, li2022generalized}
\begin{equation}
    G=M-T S=\frac{r_+}{2}\left(1+\frac{r_+^2}{L^2}+\frac{Q^2}{r_+^2}\right)- \pi T r_+^2
    \label{G}
\end{equation}
The thermodynamic pressure is related to the cosmological constant or AdS curvature radius \cite{Kastor2009, Dolan2011, Dolan2011Pressure}: \(P = \frac{3}{8 \pi} \frac{1}{L^2}\), with \(r_+\) the black hole horizon radius, simplified as \(r\). The ensemble or environmental temperature, ranging from the minimum to maximum Hawking temperature \(T_H\), results in two stable states on the free energy landscape. The kinetic rate can be approximated analytically \cite{Zwanzig2001}, and a Taylor expansion of the free energy landscape can be performed near the stable and transition states (local barrier site):\(G(r) \approx G(r_A) + \frac{1}{2}\omega_{A}^2(r-r_{A})^2\) and \(G(r) \approx G(r_m)- \frac{1}{2}\omega_{max}^2(r-r_{m})^2\)
    
The kinetic time for the phase transition can be estimated using the transition state theory, where \( \langle t_{\text{mfp}} \rangle \) is expressed as \( \frac{2\pi \eta}{\omega_A \omega_{\text{max}}} e^{\beta \Delta G_{mA}} \), with \( \beta = \frac{1}{k_B T} \) and \( \Delta G_{mA} \) representing the barrier height between the initial A and transition state. Here, \( \omega_A \) is the fluctuation frequency around basin A, \( \omega_{\text{max}} \) is the frequency at the top of barrier, and \( \eta \) is the friction coefficient derived from the diffusion coefficient \( D = \frac{k_B T}{\eta} \), where we set \( k_B = 1 \). The kinetic rate for the transition from A to B is given by \( K_{AB} = \frac{\omega_A \omega_{\text{max}}}{2\pi \eta} e^{-\beta \Delta G_{mA}} \), with a similar expression for \( K_{BA} \) in the reverse transition.

Hawking radiation in this model is represented as a reaction event characterized by the kinetic rate. According to the Stefan-Boltzmann law, the mass reduction of a black hole follows the differential equation \( \frac{dM}{dt} = -\sigma A T_H^4 \), where \( T_H \) is the Hawking temperature. This temperature depends on the black hole’s horizon radius \( r_+ \), the AdS curvature \( L \), and the charge \( Q \). Using the relationship \( \text{rate} = \text{flux}/\text{density} \), the Hawking radiation rate for RNAdS black holes can then be derived \cite{li2021kinetics}.

\begin{equation}
    K_{HR} = \left|\frac{dM}{Mdt}\right|=\frac{ \left(1+8 \pi P r^2 - \frac{Q^2}{r^2}\right)^4}{7680 \pi r^3 \left(1 + \frac{8 \pi Pr^2}{3 } + \frac{Q^2}{r^2}\right)}
    \label{H-R rate}
\end{equation}
where we have set \(c = k_B = \hbar = 1\).
The Hawking evaporation process is relatively slow, with usually longer timescale compared to the kinetic rates of phase transitions at normal temperatures. Thus, energy, entropy, and Hawking temperature can be approximated as quasi-steady states over time. In the following sections, we demonstrate that the Hawking radiation rate can produce echoes. This mechanism can be summarized as follows: the Hawking radiation rate, coupled with the fluctuations of the black hole and associated phase transitions, ultimately determines the echo behavior.

Black hole phase transition process can be described by diffusion on the free energy landscape while black hole radiation can be treated as reactions. In this reaction-diffusion framework, the system's probability evolution depends on continuous order parameters and time, with kinetic rates being coordinate-dependent. The Smoluchowski operator \(\hat{L}_D = \lambda \theta \frac{\partial}{\partial x} \left( \frac{\partial}{\partial x} + \frac{x}{\theta} \right)\) governs the diffusion under a linearized force (harmonic potential approximation), where \(\theta\) gives the equilibrium variance and \(\lambda \theta = D\) is the diffusion constant. This steady-state treatment linearizes the free energy landscape near stable states, expanding to quadratic order. Consequently, the complete probability distribution separates into independent local distributions for large and small black holes - analogous to matrix models where off-diagonal terms are negligible, permitting decomposition into independent half-reactions.


At this point, the operator \(\hat{L}_D\) describes only the local behavior of the black hole in steady state, excluding the effects of probability outflows due to phase transitions or Hawking radiation. When considering phase transitions and Hawking radiation, characterized by the kinetic rate \(K(x)\), the probability of the half-reaction should take the following form:
\begin{equation}
   \frac{\partial \rho(x,t)}{\partial t} = -K(x) \rho(x,t) + \hat{L}_D \rho(x,t)
   \label{Diffusion raction equation}
\end{equation}
This model in Fig.\ref{fig:locallanscape} corresponds to matrix model described in Fig.\ref{fig:halfreactionA} and Fig.\ref{fig:halfreactionB}. From certain perspective, the reaction-diffusion model can be regarded as an infinite-dimensional limit of the matrix model, introducing appropriate structures to characterize the properties of reactions. When the local order parameter of a steady-state black hole changes, this scenario resembles a single state comprising multiple substates (infinitely many substates defined by continuous parameters). Similarly, the reaction-diffusion model describes the diffusion driven by the underlying free energy landscape, as well as kinetic rates characterizing the reaction behavior. In addition, in this case, the local stationary state follows a Gaussian distribution, \(F = \rho = e^{-\frac{x^2}{2\theta}} / \sqrt{2 \pi \theta}\) under linear force. We will not delve into further details here, as these will be elaborated upon in the application of this model to black hole kinetics.



At this point, we have accounted for the majority of behaviors in a system with two stable-state black holes. The specific steps for applying the kinetic rate of black hole phase transitions and the kinetic rate of Hawking radiation are outlined as follows. \(K_{AB},K_{BA}\) is the interaction term (switching speeds between A and B) respectively. Here we make a transformation for the order parameter from the radius of the black hole r to near the stable state, where it becomes the difference from the stable state, for example, the small black hole A, \(x=| r-r_A |\). We immediately get that on the free energy landscape, the relaxation rate is \(\lambda_A=\omega_A^2/\eta\),and the variance is \(\theta_A=k_B T/\omega_A^2\) (note that it is actually the same form for A or B). 
Since the Hawking radiation rate is local, we can also perform the Taylor expansion for the Hawking radiation rate. The order parameter is also expanded around \(r_A\),  we use Gaussian functions to simulate \(\delta\) functions at a specific radius to smooth out Hawking radiation rates in the stable state for ease of analytic treatment.
\begin{equation}
    \begin{aligned}
    \frac{\partial \rho_a(x,t)}{\partial t} &= -K_{HR_a}(x)\delta_a(x) \rho_a(x,t) \\&+ \hat{L}_{D_{a}} \rho_a(x,t)-K_{AB}\rho_a(x,t)
      \label{Diffusion raction equation RNAdsa1}
\end{aligned}
\end{equation}
\(\delta_a(x) \approx e^{-\frac{x^2}{b_a^2}} / \sqrt{\pi}b_a\)  function are approximated by Gaussian functions and the higher order terms are ignored by doing Taylor expansion
\begin{equation}
    \begin{aligned}
    \frac{\partial \rho_a(x,t)}{\partial t} &= -\frac{K_{HR}(r_A)}{\sqrt{\pi}b_a}(1-\frac{x^2}{b_a^2}) \rho_a(x,t) \\&+ \hat{L}_{D_{a}} \rho_a(x,t)-K_{AB}\rho_a(x,t)
   \label{Diffusion raction equation RNAdsGsa1}
\end{aligned}
\end{equation}
The Green's function of half-reaction A from equation \eqref{Diffusion raction equation RNAdsa1}(refer to \cite{cao2000event,risken1984solutions} and End Matter \ref{Green's Function}) is given as
\begin{equation}
    \begin{aligned}
        G_a(x,y,t) &= e^{-K_{eff_a} t}\left[ \frac{s_a}{2 \pi \theta_a (1 - e^{-2 \lambda_a s_a t})}  \right]^{1/2} \\
    &\exp{\left[-B_a(x - y e^{-\lambda_a s_a t})^2+\alpha_a(x^2 - y^2)\right]}\label{Ga of RNAdS} 
    \end{aligned}  
\end{equation}    
where \(B_a\) is \(\frac{s_a}{2 \theta_a (1 - e^{-2 \lambda_a s_a t})} \), effective rate is \(K_{eff_a} = K_{AB} + K_{a_1} + \frac{\lambda_a (s_a-1)}{2}\),
\(s_a = \sqrt{1-\frac{4K_{a_2}\theta_a}{\lambda_a}}\),\(\alpha_a=\frac{s_a - 1}{4 \theta_a}\),\(K_{a_1}=\frac{K_{HR}(r_A)}{\sqrt{\pi}b_a}\),\(K_{a_2}=\frac{K_{HR}(r_A)}{\sqrt{\pi}b_a^3}\). One can see that \(\alpha_a\) is an important parameter, which is obtained by performing Taylor expansion \(\alpha_a=-K_{a_2}/2\lambda_a\). This parameter quantifies the fluctuation of Hawking radiation relative to the relaxation rate.In fact, it can also be seen from the effective rate \(K_{eff_a} = K_{AB} + K_{a_1} - K_{a_2}\theta_a\). It is clear that Hawking radiation promotes the black hole to deviate from the stable state. The fluctuation of Hawking radiation rate is closely related to the echo as seen later.

In the same way, one can obtain the evolution law of the half-reaction of large black hole B. The order parameter as the radius of the black hole switches to the deviation of \(r_B\), \(x=|r-r_B|\) .
The Green's function for half-reaction B can be obtained using a similar method.  
\begin{equation}
    \begin{aligned}
    G_b(x,y,t) &= e^{-K_{eff_b} t}\left[ \frac{s_b}{2 \pi \theta_b (1 - e^{-2 \lambda_b s_b t})} \right]^{1/2}\\
    &\exp{\left[-B_b(x - y e^{-\lambda_b s_b t})^2 + \alpha_b(x^2 - y^2)\right]}
    \end{aligned}
\end{equation} 
where \(B_b\) is \(\frac{s_b}{2 \theta_b (1 - e^{-2 \lambda_b s_b t})}\)  
effective rate is \(K_{eff_b} = K_{BA} + K_{b_1} +\frac{\lambda_b(s_b-1)}{2}\),relative rate is \(\alpha_a=\frac{s_a - 1}{4 \theta_a}\),\(K_{b_1}=\frac{K_{HR}(r_B)}{\sqrt{\pi}b_b}\),\(K_{b_2}=\frac{K_{HR}(r_B)}{\sqrt{\pi}b_b^3}\).


One can obtain the distribution of the switching event \( f_a(t) =K_{AB}\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} dx \, dy \,  G_a(x, y, t) \rho_b(y)\). The kinetic event of transition \(f_b(t) =K_{BA}\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} dx \, dy \,  G_b(x, y, t) \rho_a(y)\). The joint distribution of switching events \(f_{ba}\) can be obtained
\begin{equation}
    \begin{aligned}
     f_{ba}(t_1,t_2)&=
     K_{AB}K_{BA}\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}dx \, dy \, dz \,\\
      &G_a(x, y, t_2)G_b(y, z, t_1) \rho_a(z) 
    \label{RNAdS fba}
    \end{aligned}
\end{equation}


The kinetic rates \(K_{AB}\) and \(K_{BA}\), which characterize the phase transition driven by environmental thermal fluctuations and depend on \(T\), \(Q\), and \(P\), form the essential framework of the echo phenomenon. While Hawking radiation contributes a correction to the effective rates, it is the underlying phase transition kinetics that are indispensable: if \(K_{AB} = K_{BA} = 0\), the phenomenon vanishes entirely. The difference function, computed via multi-dimensional Gaussian integration or numerical methods using the distributions \(f_a(t) = \Delta_a K_{AB}e^{-K_{eff_a}t}\), \(f_b(t) = \Delta_b K_{BA}e^{-K_{eff_b}t}\), and \(f_{ba}(t_1,t_2) = \Delta_{ba} K_{AB}K_{BA}e^{-K_{eff_b}t_1}e^{-K_{eff_a}t_2}\) (see \ref{distribution}), thus probes the combined effect of the black hole phase transition kinetics and the Hawking radiation. Analyzing how \(T\), \(Q\), and \(P\) influence the echo through this framework allows the echo itself to be used as a probe of black hole characteristics such as phase transition and Hawking radiation.






\begin{figure}[!ht]
    \centering
    \includegraphics[width=0.95\linewidth]{fg/Twoevent3MFPT.pdf}
    \caption{ Mean first-passage time (MFPT), calculated from the distribution \(f_{ba}(t,t)\), as a function of temperature. Curves correspond to dissipation coefficients \(\eta = 100\) (red), \(10^4\) (blue), and \(10^5\) (purple), at fixed \(Q=0.1\), \(P=0.003/(8\pi)\), and \(b_a=b_b=50\). }
    \label{fig:ScalingBehavior}
\end{figure}

\begin{figure*}[!ht]
    \centering
 
    \subfigure[]{
    \includegraphics[width=0.39\linewidth]{fg/Logtvsecho.pdf}\label{fig:3RNecho}
    }
    \subfigure[]{
    \includegraphics[width=0.38\linewidth]{fg/relativerate.pdf}\label{fig:relativerate}
    }
    \caption{
 (a) Same-time difference function \(\delta(t)\) versus \(\log(t)\) for temperatures \(T = 0.0312\) (blue), \(0.0313\) (red), and \(0.0314\) (green), with fixed \(Q=0.1\), \(P=0.003/(8\pi)\), \(\eta=100\), and \(b_a=b_b=50\). The echo time is defined as the location of the maximum following an initial rapid decrease. (b) Relative Hawking radiation rate \( |\alpha_b|\) for the large black hole (state B) as a function of temperature.
}
    \label{fig:3}
\end{figure*}

\begin{figure*}[!ht]
    \centering
 
    \subfigure[]{
    \includegraphics[width=0.4585\linewidth]{fg/echopeakvsTandQ.jpg}\label{fig:echopeakvsTandQ}
    }
    \subfigure[]{
    \includegraphics[width=0.5\linewidth]{fg/echopeakvsTandP.jpg}\label{fig:echopeakvsTandP}
    }
    \caption{Variation of the echo peak height with parameters. (a,b) The peak increases with \(Q \sim 0.1\text{--}0.15\) and \( P \sim \frac{3}{8 \pi}0.01\text{--}\frac{3}{8 \pi}0.0101 \), while showing a rise-and-fall behavior with \(T\), for fixed \(\eta=100\) and \(b_a=b_b=50\).
 }
    \label{fig:4}
\end{figure*}


\begin{table*}[!ht]
\caption{Echo peak amplitudes and echo times for different parameters with \(T=0.03,Q=0.1,P=\frac{3}{8 \pi} 0.01\).}
\label{tab:echo_data}
\begin{ruledtabular}
\begin{tabular}{@{}*{9}{c}@{}}
\toprule
$\eta$ & $b_a$ & $b_b$ & $K_{BA}$ & $K_{AB}$ & $K_{a_2}$ & $K_{b_2}$ & Echo Peak & Echo Time \\
\midrule
50    & 50    & 50  & $8.97394 \times 10^{-9}$  & $1.15006 \times 10^{-13}$ & $1.9637 \times 10^{-13}$ & $1.8728 \times 10^{-11}$ & $2.09332 \times 10^{-34}$ & 13.5656      \\
100   & 50    & 50  & $4.48697 \times 10^{-9}$  & $5.7503 \times 10^{-14}$ & $1.9637 \times 10^{-13}$ & $1.8728 \times 10^{-11}$ & $1.04431 \times 10^{-34}$ & 26.3525      \\
200   & 50    & 50  & $2.24349 \times 10^{-9}$  & $2.87515\times 10^{-14}$ & $1.9637 \times 10^{-13}$ & $1.8728 \times 10^{-11}$ & $5.22625 \times 10^{-35}$ & 52.0387     \\
100   & 100   & 50  & $4.48697 \times 10^{-9}$  & $5.7503 \times 10^{-14}$ & $2.45463 \times 10^{-14} $ & $1.8728 \times 10^{-11}$ & $1.35417 \times 10^{-35}$ & 30.3996      \\
100   & 50    & 100 & $4.48697 \times 10^{-9}$  & $5.7503 \times 10^{-14}$ & $1.9637 \times 10^{-13}$ & $2.341 \times 10^{-12}$ & $1.04384 \times 10^{-34}$ & 23.1872      \\
\bottomrule
\end{tabular}
\end{ruledtabular}
\end{table*}

Figure \ref{fig:ScalingBehavior} shows kinetic turnover emerges at small friction coefficients as temperature increases. For RNAdS large black holes, stronger Hawking radiation at higher temperatures drives this turnover from kinetics-dominated to radiation-dominated phase transitions. With large friction coefficients, this turnover vanishes, revealing distinct dynamical regimes. These timescales reflect competition between phase transitions and Hawking radiation, characterized by the effective rate \(K_{eff}\) governing exponential decay. The joint switching distribution \(f_{ba}(t,t)\) captures these dynamics, while the difference function provides additional radiation characteristics, as we demonstrate below.

The same time difference function, defined as \(t_1 = t_2\), can characterize most of the properties of the echo. For instance, this function can characterize the echo peak and the timescale of the echo.  








In Fig. \ref{fig:3RNecho}, the echo peak shows strong temperature dependence, linked to Hawking radiation fluctuations. The difference function takes the form:
\(
\delta = \left| \Delta_{ba} - \Delta_b \Delta_a \right| \exp\left[-K_{eff_b} t_1\right] \exp\left[-K_{eff_a} t_2\right],
\)
comprising an exponential decay envelope controlled by effective rates and a fluctuation term involving relaxation rate \(\lambda\), variance \(\theta\), and relative rate \(\alpha\)—or a Hawking radiation fluctuation function of \(t_1\) and \(t_2\). Without Hawking fluctuations (i.e., \(K_{a_2} = K_{b_2} = 0\)), the distributions simplify to exponential forms, yielding \(\delta = 0\), indicating event independence. These results suggest that fluctuations in Hawking radiation significantly influence echo behavior, where the difference function is proportional to the variance of the stochastic rate \cite{cao2000event}. The echo vanishes when \(\alpha = 0\), showing that relative fluctuation rates—not absolute values—govern the echo amplitude. The effective rates set the exponential decay scale, while relative fluctuations drive the signal.

Fig. \ref{fig:relativerate} illustrates the relative rate \(\alpha_b\) for large black holes, decreasing at low temperatures (due to rising relaxation rates) and increasing at high temperatures (dominated by evaporation), reflecting the interplay between Hawking radiation and phase transition kinetics \cite{li2021kinetics}.






The echo peak height is governed by (\(K_{AB}\) and \(K_{BA}\)) the phase transition kinetics between states A and B, exhibiting strong dependence on parameters such as the relative fluctuation rate \(\alpha\) (Fig. \ref{fig:echopeakvsTandQ}). It reaches a maximum near the critical temperature, where the two states become energetically degenerate. This behavior stems from two competing factors: (1) the comparable transition rates \(K_{AB}\) and \(K_{BA}\) near degeneracy enhance the probability of correlated \(B \rightarrow A \rightarrow B\) sequences, while (2) larger energy barriers away from criticality suppress such transitions. The peak magnitude thus reflects the combined effect of phase transition dynamics and Hawking radiation, both being shaped by the underlying free energy landscape.



As shown in Fig.\ref{fig:4}, the echo peak increases with charge \(Q\) or pressure \(P\), while its temperature dependence is non-monotonic. These behaviors originate from parameter-induced changes in the free energy landscape, which alter the black hole radius, curvature, Hawking radiation rate, effective kinetic rates, and relative fluctuation rate \(\alpha\). The echo phenomenon is thus structurally determined by the landscape.

In Table \ref{tab:echo_data},further mechanistic insight comes from selectively varying parameters: changing the Gaussian fluctuation factor \(b\) only affects Hawking radiation contributions and modifies \(\alpha\), thereby altering the echo amplitude and position, although varying parameter \(b\) alters the contribution of Hawking radiation to the effective rates, we can at least observe distinct changes in both the echo peak height and its temporal position solely due to modifications in Hawking radiation. In contrast, adjusting the dissipation coefficient \(\eta\) only influences kinetic rates \(K_{AB}\) and \(K_{BA}\), affecting both the amplitude and temporal position of the echo peak. This decoupling confirms that Hawking fluctuations and kinetic rates distinctly shape echo behavior. When the difference function is analyzed for a fixed ratio between the two time arguments, \( k = t_1/t_2 \), its behavior closely mirrors the single-time scenario, demonstrating a similar strong dependence on phase transitions and Hawking radiation fluctuations. (refer to End Matter \ref{two dimensional difference function}).




Echo signals originate from the interplay between phase transition kinetics and Hawking radiation fluctuations, providing a dynamical probe of black hole thermodynamics. Their amplitude reflects fluctuation strength while their timescale encodes transition rates. These correlation dynamics may extend to other contexts such as gravitational waves from cosmological phase transitions\cite{hogan1986gravitational} and analogue gravity systems\cite{unruh1981experimental,weinfurtner2011measurement,kolobov2021observation}.Using experimental parameters from the acoustic black hole analogy, as realized in Bose-Einstein condensate analogues \cite{steinhauer2016observation}, the effective Hawking evaporation rate is estimated as \(K_{HR} \sim \alpha^4 c_{\text{out}}/L\). Here, \(c_{\text{out}}\) is the speed of sound just outside the analogue event horizon, \(L\) represents the characteristic size of the analogue black hole, and \(\alpha\) is a dimensionless coupling parameter characterizing the interaction strength in the condensate. This rate reaches \(\mathcal{O}(0.1)\,\mathrm{s}^{-1}\) for typical configurations where the Hawking temperature satisfies \(k_{B}T_{H} \sim \alpha\, m c_{\text{out}}^2\) (with \(m\) representing the mass of an atom in the condensate), suggesting potential observability of echo-like correlations. This underscores the need to unify phase transition dynamics with radiation fluctuations within the free energy landscape for a complete microscopic description.


\begin{acknowledgments}
T. Y. acknowledges support from the National Natural Science Foundation of China (Grant No. 12234019).
\end{acknowledgments}







\nocite{*}
\bibliography{apssamp}


\appendix 
\section{End matter}
\subsection{Green's function }\label{Green's Function}

From the partial differential equation, one can get the equation satisfied by its Green function
\begin{equation}
\begin{aligned}
       \frac{\partial}{\partial t} G(x, y, t) &= - (K_{a_1}  + K_{AB})G(x, y, t)+K_{a_2}x^2 G(x, y, t)\\&+ \lambda_a \theta_a \frac{\partial}{\partial x} \left( \frac{\partial}{\partial x} +\frac{x}{\theta_a}   \right) G(x, y, t)
    \label{A G}
\end{aligned}
\end{equation}
For simplicity, some of the lower corner symbols a are omitted from the above symbols.
The initial condition is
\begin{equation}
    G(x, y, 0) = \delta(x - y) \label{initial condition a}
\end{equation}
Applying the transformation
\begin{equation}
   G(x, y, t) = g(x, y, t) e^{\alpha_a(x^2 - y^2)}
   \label{transformation}
\end{equation}
where \( \alpha_a=\frac{s_a - 1}{4 \theta_a},
    s_a= \sqrt{1-\frac{4K_{a_2}\theta_a}{\lambda_a}}\)
The \(g(x,y,t)\) satisfies the Fokker-Planck equation for the Ornstein-Uhlenbeck process
\begin{equation}
\begin{aligned}
    \frac{\partial}{\partial t} g(x, y, t) &= \left[ \lambda_a s_a \frac{\partial}{\partial x} x+ \lambda_a \theta_a \frac{\partial^2}{\partial x^2} \right] g(x, y, t) \\&- K_{eff_a} g(x, y, t)
\end{aligned}
\label{Ornstein-Uhlenbeck process a}
\end{equation}

where \(K_{eff_a}=K_{AB}+K_{a_1}+\frac{\lambda_a (s_a-1)}{2}\)
The initial condition is 
\begin{equation}
    g(x, y, t) = \delta(x - y)
    \label{initial condition a1}
\end{equation}
Rewriting \(g(x,y,t)\) as
\begin{equation}
    g(x, y, t) = g_1(x, y, t)e^{-K_{eff_a}t} 
    \label{Rewriting a}
\end{equation}
where \(g1(x,y,t)\) is the Green function for the standard OrnsteinUlenbeck process
\begin{equation}
    \frac{\partial P}{\partial t} = \gamma \frac{\partial}{\partial x} (xP) + D \frac{\partial^2 P}{\partial x^2}
    \label{standard OrnsteinUlenbeck process a}
\end{equation}
with \(\gamma= \lambda s\) and \(D = \lambda \theta\). The standard solution is given  as
\begin{equation}
\begin{aligned}
    g(x, y, t) &= e^{-K_{eff_a}t}\sqrt{\frac{s_a}{ 2 \pi \theta_a (1 - e^{-2\lambda_a s_a t})}}\\&\exp\left[-\frac{s_a(x - ye^{-\lambda_a s_a t})^2}{2\theta_a (1 - e^{-2\lambda_a s_a t})}\right]
    \label{ga}
\end{aligned}
\end{equation}
Thus, one have a complete Green function in equation.
\begin{equation}
    \begin{aligned}
        G_a &= e^{-K_{eff_a} t}\left[ \frac{s_a}{2 \pi \theta_a (1 - e^{-2 \lambda_a s_a t})}  \right]^{1/2} \\
    &\exp{\left[-B_a(x - y e^{-\lambda_a s_a t})^2+\alpha_a(x^2 - y^2)\right]}\label{Ga of RNAdS} 
    \end{aligned}  
\end{equation}    
where \(B_a\) is \(\frac{s_a}{2 \theta_a (1 - e^{-2 \lambda_a s_a t})} \), effective rate is \(K_{eff_a} = K_{AB} + K_{a_1} + \frac{\lambda_a (s_a-1)}{2}\),
\(s_a = \sqrt{1-\frac{4K_{a_2}\theta_a}{\lambda_a}}\),\(\alpha_a=\frac{s_a - 1}{4 \theta_a}\). 



















\subsection{Distributions of single event and joint two events}\label{distribution}
\begin{figure}[!ht]
    \centering
    \subfigure[]{
    \includegraphics[width=0.99\linewidth]{fg/deltavsk.jpg} 
    } 
    \caption{ The difference function \(\delta(t_1, t_2/k)\) at fixed  parameters: \(T = 0.0312\), \(Q = 0.1\), \(P = 3/(8\pi) \times 0.01\),\(\eta=100\) and \(b_a=b_b=50\).}\label{fig:deltavsk}
\end{figure}
One obtains the distribution of the kinetic event of transition \(f_a(t)\)
\begin{equation}
\begin{aligned}
    f_a(t) &= \frac{\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} dx \, dy \, K_{AB} G_a(x, y, t) K_{BA}\rho_b(y) }{\int_{-\infty}^{\infty}dx K_{BA}\rho_b} \\& =K_{AB}\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} dx \, dy \,  G_a(x, y, t) \rho_b(y) 
    \label{RNAdS fa}
\end{aligned}   
\end{equation}
where \(\rho_b(x) = e^{-\frac{x^2}{2\theta_b}} / \sqrt{2 \pi \theta_b}\). In fact, the second equal sign in the equation \eqref{RNAdS fa} is normalization process of flux:\(F_a=K_{BA} \rho_b/\sum K_{BA} \rho_b\).
However, since the kinetic rate is independent of \(x\), we have \(F_b = \rho_b(x)\). When combined with the Green’s function, equation \eqref{RNAdS fa} represents a Gaussian integral with complex coefficients, yielding the result:
\begin{subequations}
    \begin{align}
        f_a(t)&=\Delta_aK_{AB}e^{-Keff_at} \\
        \Delta_a&=\sqrt{\frac{s_a}{4\theta_a\theta_b(1-e^{-2 \lambda_as_at})}}\sqrt{\frac{1}{A_{a_1}B_{a_1}}}\\
        A_{a_1}&=B_ae^{-2\lambda_a s_at} +\frac{1}{2\theta_b}+\alpha_a\\
        B_{a_1}&=B_a-\frac{B_a^2e^{-2\lambda_a s_a t}}{A_{a_1}}-\alpha_a
    \end{align}
\end{subequations}
In summary,\(f_a\) is proportional to \(\Delta_a\) and decaying exponentially with the effective kinetic rate \(K_{eff_a}\),\(\Delta_a\) is also naturally a function of time \(t\),\(\theta\),\(\lambda\) and \(s\).


\begin{equation}
\begin{aligned}
      f_{ba}(t_1,t_2)  &=K_{AB}K_{BA}\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}dx \, dy \, dz \, \\&G_a(x, y, t_2)G_b(y, z, t_1) \rho_a(z) 
\end{aligned}
    \label{RNAdS fba}
\end{equation}
The result of integration is
\begin{subequations}
    \begin{align}
        f_{ba}(t_1,t_2)&=\Delta_{ba}K_{AB}K_{BA}e^{-Keff_bt_1}e^{-Keff_at_2} \\
        \Delta_{ba}=&
        \sqrt{\frac{s_a}{8\theta_a^2\theta_b(1-e^{-2 \lambda_as_at_2})(1-e^{-2 \lambda_bs_bt_1})}}\\&*\sqrt{\frac{1}{A_{ba}B_{ba}C_{ba}}}\\
        A_{ba}&=B_be^{-2\lambda_b s_b t_1} +\frac{1}{2\theta_a}+\alpha_b\\
        B_{ba}&=B_b-\frac{B_b^2e^{-2\lambda_b s_b t_1}}{A_{ba}}+B_ae^{-2\lambda_a s_a t_2}+\alpha_a-\alpha_b\\
        C_{ba}&=B_a-\frac{B_a^2e^{-2\lambda_a s_a t_2}}{B_{ba}}-\alpha_a
    \end{align}
\end{subequations}


\subsection{Difference function in time under different parameters }\label{two dimensional difference function}

It is important to note that these parameters are also functions of time. For example, \(B_b\) depends on \(t_1\). In summary, the specific dependence on either \(t_1\) or \(t_2\) is not explicitly indicated here; however, it can be inferred from the subscript of \(\lambda\), where \(a\) corresponds to 2 and \(b\) corresponds to 1. A complete characterization of black hole echoes is achieved by investigating the two-dimensional difference function for trajectories defined by \( t_1 = k t_2 \) across different values of \( k \). Figure \ref{fig:deltavsk} illustrates that the echo features evolve with \( k \), in a manner similar to the one-dimensional scenario. The parameter \( k \) specifically influences the decay rate and the amplitude of the prefactor \( |\Delta_{ba} - \Delta_b\Delta_a| \). Therefore, this method enables a full investigation of the echo behaviors.

\end{document}

