I. INTRODUCTION

The (nonlinear) gravitational wave (GW) memory effect [1][2][3][4][5][6][7][8][9], representing a physical manifestation of the nonlinearity intrinsic to general relativity (GR), is a particularly interesting aspect of the theory and has yet to be confirmed experimentally. This memory reveals itself by a permanent displacement of the freely floating test masses within an idealized GW detector after it has been traversed by gravitational radiation. It has been subject to many investigations and plays a key role, for instance, in the Bondi-Metzner-Sachs formalism of asymptotically flat spacetimes [10][11][12], Weinbergs'ss o f tt h e o r e m s [13],a n d the infrared triangle [14,15]. The memory carries a manifold of analytically and observationally interesting phenomenology, ranging from connections to the fundamental symmetries of spacetime (as established by Bondi and others) to intrinsic (kinematic) properties of binary black hole (BBH) mergers [16]. Beyond that, however, it can also pick up more subtle features such as corrections due to modified gravity theories [17][18][19][20][21][22][23] and extra compact dimensions [24]. One aspect that has only sparsely been addressed in literature is the remnant of potential black hole (BH) quantum features in the GW memory. While there is a small body of work aiming to treat gravitational radiation from a quantum mechanical point of view, for instance [25,26], the quantumness in GW physics comes up majorly in the context of the GW echo effect [27][28][29][30][31] (see also [32][33][34][35][36]; we point out [37] for a more generalized treatment of quantum modifications such as echos and [38][39][40][41][42] for earlier works on echoes). As will be demonstrated in this work, the echo induces a plentitude of features within gravitational waveforms, including a contribution to the memory. Generally, the echo effect refers to residual GW signals from a BBH merger arriving at the detector with a constant (but potentially frequency-dependent) time delay after the primary waveform has passed through. This echolike signature is produced when the ingoing radiation from the ringdown phase of a binary merger is reflected by a (quantum) surface close to the classical event horizon. A delayed GW echo is predicted by a range of theoretical models, including those involving deviations from GR [43,44]. In this work, the cases of interest are echoes arising due to the presence of near-horizon (quantum) structure surrounding the BH [30,[45][46][47][48][49][50][51] as well as the existence of exotic compact objects (ECO) replacing the traditional BH by substituting the would-be horizon with a physical boundary [52][53][54][55][56][57]. In particular in the context of near-horizon quantum structure, the echo manifests inadmissible evidence of an interplay between quantum physics and gravity [36,47,51]. Recent studies also suggest the possibility of echoes arising due to simple quantization arguments, containing imprints of the specific underlying theory of quantum gravity [31]. ECOs on the other hand challenge the very foundations of our understanding of compact stellar objects and provide new avenues to resolve long-standing conceptual issues that have been present in BH physics. Despite the existence of stable ECOs and their production of echoes being under debate (see for instance [58][59][60][61][62][63]), searching for GW echoes in upcoming measurement data has the potential to enlighten the fundamental dark spots in our current understanding of ECOs and, in particular, (quantum) BHs.

The variety of phenomenological models sourcing a GW echo generally complicates the computation of contribution in an agnostic manner. Nonetheless, the prospect of new data collected by future space-and ground-based instruments has facilitated the efforts toward constructing accurate gravitational waveforms including echoes [64][65][66]. Attaching an echo to the existing waveform model templates further increases the already large parameter space for numerical waveform simulations. Even for echoless waveforms, more efficient waveform modeling is already an active research topic as the demands for the precision of such models will increase for future GW data [67,68]. Adding another layer of complexity to the template construction process by attaching a model-dependent echo to the waveform seems counterproductive. One can therefore ask, whether simulating the exact shape of the echo time series can be avoided without losing important astrophysical and BH related information, or, in other words, does the echo also leave traces in the waveform which are less (or not at all) sensitive toward the particular assumptions that determine its the echo's morphology?

A common feature of all strain signals passing through a GW detector is that they carry gravitational energy. This implies that for each signal, there exists a finite interval in the strain time series hðtÞ during which j ḣðtÞj 2 > 0.E v e n though the precise form of the strain time series may be uncertain due to precision limitations in the data analysis or regarding the waveform templates, the features appearing in detector data due to the passing energy flux are well defined: The flux through the detector will cause a contribution to the nonlinear memory, leading to a permanent displacement of the (ideals) detector's test masses. Consequently, each echo contributes to the accumulated GW memory, simplifying the complex time series into a steplike increase in the measured strain. In this work, we explore this claim and estimate the detectability of echo-induced memory, examine the information that can be extracted from such observations, and investigate how these findings are connected to the dynamics that generate GW echoes.

To date, neither GW memory nor the echo have been detected. Regarding the memory, the measurements of current ground-based detectors have not achieved sufficient signal-to-noise ratio (SNR) for a claimed detection, as they typically respond to GWs with timescales much shorter than those of typical memory signals. However, future groundbased detectors [69][70][71][72][73] and space-based missions [74,75] hold great potential for detecting the GW memory. In particular, the Laser Interferometer Space Antenna (LISA) [76], with its freely floating test masses and excellent sensitivity in the low-frequency band, where memory signals tend to be stronger, stands out as a promising candidate to probe the memory effect. Indeed, early estimates indicate a promising SNR for massive binary systems at redshifts up to z ≈ 2 [77]. Similarly, while there has been no direct detection of GW echoes despite numerous efforts (see [78,79] for a review and also [80]), LISA's prospects for detecting echoes have been demonstrated to be promising in a companion paper [81]. These prospects remain strong even when echoes are treated as separate signals, detached from their "sourcing" waveform-a scenario that has been suggested as plausible [82].

Building up on the tentative findings for echoes in [81], we begin our investigation into the echo-induced memory by examining the reflective properties of ECOs and BHs in Sec. II. For the BH reflectivity, we rely on the model developed in [81], which is based on area quantization arguments and referred to as quantum black holes (QBH) throughout this work. For ECOs, we use the model adopted in [66]. Section II B details the methodology for constructing echo templates from numerical waveforms, comparing the results for both reflectivity models. In Sec. III,w e explore the properties of echo-induced memory and derive semianalytical expressions for its dominant components. Additionally, we investigate the associated flux balances. To complement the analytical findings, we conduct numerical simulations to explore the echo-induced memory and its related flux laws, utilizing a set of waveforms from the Simulating eXtreme Spacetimes (SXS) dataset (see [83] for the latest update). Distinctive features of memory contributions resulting from different theoretical origin are explored. Furthermore, we estimate the SNR for the memory contributions due to echoes using a recently developed LISA data analysis pipeline [77]. A comprehensive discussion of the numerical results and their implications is provided in Sec. V.

II. GRAVITATIONAL WAVES FROM QUANTUM BLACK HOLES: THEORY

The mechanism leading to the creation of GW echoes is based on a perturbative treatment of the spacetime metric outside the compact object formed during a binary merging event. After the merger, the remnant body adopts a perturbed Kerr (or Schwarzschild) metric, which gradually settles into an unperturbed state. The decay of these perturbations is governed by the Teukolsky equations [84,85]. The latter can be solved analytically, revealing both an ingoing and outgoing solution in terms of quasinormal modes (QNMs). The outgoing solution corresponds to the ringdown phase of the gravitational waveform observed by the GW detector. The ingoing solution, in the classical framework, crosses the event horizon of the BH and is lost forever inside the horizon. This part of the solution is typically considered irrelevant in classical treatments of BH physics. However, this picture changes dramatically if the horizon is no longer treated as a fully transmissive, "unspectacular" surface.

A. Reflective properties of quantum black holes and exotic compact objects

In the quantum regime, the BH horizon plays a significant role for physical processes, as it has first been established by Hawking when he demonstrated that the BH in fact emits thermal radiation. To date, the exact quantum structure of a BH's horizon remains unknown. Nevertheless, numerous theories have been proposed to model its quantum properties, including the membrane paradigm [86][87][88] (see also [36,89,90]), effective field theories [91], phenomenological approaches [51] and quantum mechanics-inspired models [92]. In its original form, the membrane paradigm postulates a replacement of the BH horizon by a fictitious (physical) membrane. This framework has been particularly useful for modeling the outer layers of ECOs. In this article, we focus on two recent proposals in which (i) the membrane paradigm is used to impose boundary conditions near an ECO's surface [36,50,51] and a Boltzmann reflectivity applies (see also [57] for similar approaches avoiding membrane paradigm), and (ii) a BH reflectivity is constructed based purely on area quantization arguments [31,81]. We outline the relevant results from both approaches in the following sections.

Boltzmann reflectivity and membrane paradigm

Replacing the BH horizon with a (quantum) membrane fundamentally alters the behavior of linear perturbations of the BH spacetime. In [50,51] it was argued that quantum effects near the horizon, which lead to thermal emission at the Hawking temperature T H , modify the equations governing the evolution of these linear perturbations. Concretely, the membrane paradigm and other effects in convolution with the fluctuation-dissipation theorem lead to

as the evolution equation for QNM functions, ψ ω ,ofGWs under a stochastic fluctuation field ξ ω . In the latter, ΩðxÞ ≔ ω= ffiffiffiffiffiffi g 00 p is the blueshifted frequency, E Pl the Planck energy, γ the dimensionless dissipation parameter, and VðxÞ the BH potential. The dissipation term resembles the viscous dissipation for sound waves but is suppressed by the gravitational coupling Ω=E Pl . The latter nicely matches with phenomenology resulting from quantum gravity theories where constraints on spacetime viscosity can be implemented via viscous terms [93].

Following (1), the fluctuation-dissipation theorem predicts a thermal spectrum for the mode functions ψ ω .Near the horizon, one can thus compute the analytical solution of the QNM functions assuming a constant surface gravity of κ ¼ 2πT H . Physically, this can be interpreted as the energy flux of ingoing GWs being unable to penetrate the BH horizon and getting either absorbed or reflected. This result is consistent with the models for which the BH interior cannot be entered, such as the membrane paradigm or the string theory inspired fuzzball complementary conjecture [94,95].

In the near horizon limit, logðE Pl =γjωjÞ ≪ κx ≪-1, one thus finds the Boltzmann (flux) reflectivity [51] jRj¼

We note this equation is independent of the dissipation parameter γ appearing in the dissipation term of Eq. ( 1). This is due to taking the absolute value. The parameter γ determines the time-separation of the reflected echoes [51] and thus enters as a complex phase. Concretely, including dissipation effects, the position of reflection is marked by the distance x 0 for which γΩðxÞ ∼ E Pl . Thus, in this model,

For a more exhaustive derivation of the latter, we refer to [50]. The time delay can be rewritten as a complex phase and incorporated into R such that

Note here that while the arguments leading to the Boltzmann reflectivity are less model dependent, the phase of reflectivity results purely from Eq. ( 1), in particular, the structure of the dissipation term. Thus, the time delay between echoes is highly sensitive to the underlying (quantum) theory. 1 Similar results to [36,50,51] have also been obtained following approaches with very different physical motivation, e.g., [91]. In this work, the reflectivity (4) will serve as a toy model for the reflective properties of ECOs. Thus, we henceforth label Eq. ( 4) with R ECO .W e generalize the ECO's reflectivity slightly by introducing an effective horizon temperature T QH , in place of the standard Hawking temperature T H ¼ 1=8π. This allows us to explore how varying the horizon temperature influences the properties of GW echoes, providing insights into the impact of different quantum horizon effects on the observed signals.

Physically, the first factor in Eq. ( 4) represents the position of the membrane located near the horizon from which the GW reflects. In radial coordinates, this location can be expressed as r echo ∼ ln γ 2πT QH . The second exponential factor in (4) governs the frequency range over which the reflectivity is significantly large. Thus, while the dissipation factor γ characterizes the time separation of the echoes, the effective horizon temperature T QH primarily determines the width of the reflected frequency band. Generally, the exponential suppression of larger frequencies is motivated by a quantum mechanical perspective, treating an isolated BH as a multilevel system [51]. In this "giant atom" model, Hawking radiation serves as the mechanism through which the BH can (spontaneously) de-excite, while the BH's reflectivity can be interpreted as a form of stimulated emission. Consequently, the frequency scale ω ≲ T H naturally emerges as the regime where stimulated emission is possible, while for frequencies ω ≫ T H the BH behaves as effectively all absorbing.

Area quantization and black hole energy states

Even when ECOs are not permissable, 2 simple quantization arguments can still support the existence of GW echoes. A compelling proposition was made in [31], which is based on Bekenstein and Mukhanov's area quantization formula A N ¼ αl Pl N, where l Pl is the Planck length and N a positive integer. This approach leads to a discrete mass spectrum for a given BH, which in turn quantizes its emission and absorption processes. Considering the spin of the BH-though we generally focus on nonspinning remnant BH with dimensionless angular momentum a ≪ 1-and the fact that, in astrophysical context, we encounter macroscopic BHs only, i.e., N ≫ 1, one finds for the characteristic frequencies being absorbed by the BH [31] that

The coefficients κ and Ω H denote the surface gravity and the angular momentum, respectively. With the frequency ω N scaling as 1=M, Planck-scale effects are magnified and elevated into the frequency regime relevant for present-day GW interferometers.

In principle, the characteristic frequencies serve as the sole, narrow pathway for GWs to enter a BH by crossing its horizon. However, if the BH is spinning rapidly enough, the width of the quantized energy levels, denoted as Γ, becomes significant. This width is inversely proportional to the decay rate associated with the spontaneous emission of Hawking radiation, which leads to the de-excitation of the BH [31].If the energy states were to overlap, BHs would behave as true absorbers, allowing any frequency ω to cross the horizon. This scenario, however, does not hold, even for highly spinning remnant BH, as long as the phenomenological constant α exceeds a critical value α > α crit . Notably, this critical value α crit is much smaller than the lowest phenomenological constant typically considered in the literature [31],i . e . ,i sα ¼ 4 log 2. Consequently, the overlap of energy levels is generally not anticipated and even for α ¼ 4 log 2 they remain quite narrow.

The latter implies that a considerable amount of the GW ringdown's modes content cannot be absorbed by a QBH. As argued in [30], the remaining modes could be reflected, producing a late-time echo in the gravitational waveform. Modeling the reflectivity of QBH is highly nontrivial and requires careful consideration of various (quantum) effects. In this work, we utilize the simplified, phenomenologically motivated toy model described in [81]. The latter includes a broad range of potential echo morphologies and is based on the key assumption that only the characteristic frequencies ω N can be absorbed. Including phenomena like line broadening and exponential decay, the characteristic frequencies' sharp lines in the QBH absorption spectrum turn into broader cusplike features with exponential damping. The final reflectivity coefficient for the QBH is found to be [81]

Just as for the reflectivity of the ECO, R ECO , Eq. ( 6) contains an exponential decay term governed by T QH .I n addition, R QBH acquires the sine term with exponent δ which effectively determines the sharpness of the absorption lines corresponding to the characteristic frequencies 1 In [36], for instance, the membrane paradigm is treated more rigorously and via a replacement of the horizon with a fictitious quantum membrane. This membrane constitutes an ensemble of microscopic degrees of freedom in the ground state following a Gaussian wave function. As a result, a frequency-dependent reflectivity can be defined (see Fig. 1 in [36]) and similarly, a time separation between echoes is defined as

Both the reflectivity and Δt are largely dependent on the variance of the quantum state describing the microscopic degrees of freedom, σ 2 .

2 An detailed discussion on this matter is presented in [78].

ω N ðαÞ. The sine-term itself assures the equidistant distribution of absorption lines. 3 The complex exponential, similar to R ECO , encapsulates the time shift of the echoes compared to the QNM of the initial GW's ringdown. For the QBH, the exponent is constructed based on the following considerations: The time delay between the onset of the ringdown and the arrival of the first echo at the GW detector is determined by the time duration it takes for the ingoing gravitational radiation-resulting from solving the Teukolsky equation [84,85] and morally emitted from the BH potential barrier 4 to travel to the reflective surface and return to the BH's potential barrier, i.e.,

where r à ¼ r þ 2 ln ð r 2 -1Þ is the tortoise coordinate and r Shell >r H for r H ¼ 2M. Thereby, the distance between the reflective surface and the BH horizon is crucial for an outside observer's perception due to time dilation GW modes experience close to the horizon. The latter distance is quantified by the parameter β. For a detailed calculation of Eq. ( 7), we refer to [32,96].

In summary, the reflectivity R QBH is determined by four parameters: the phenomenological constant α (which in turn unequivocally determines ω N ), β, δ, and T QH .F o r R QBH , we reparametrize T QH in R QBH by T QH ¼ ϵT H such that ϵ steers the exponential decay of the reflectivity coefficient. To better distinguish between the models in the analysis of this work, we keep T QH as a free parameter for R ECO such that it is completely determined by γ and ϵ.It is important to note that the parameters γ and β only impact the time separation of the echo and the GW in the strain time series picked up by the GW detector. Additionally, the parameter α does not constitute a freely chosen model parameter of reflectivity; rather, it is derived from fundamental properties of the underlying quantum theory of gravity. In Bekenstein's original proposal, the constant takes the value α ¼ 8π [97] and subsequent computations relying on different approaches have confirmed his findings [98]. Nevertheless, alternative values are present in the literature, and we therefore treat α as a free parameter. The actual significance of these model parameters concerning the echo-induced features in the gravitational wave memory will be explored below.

B. Echo generation via the hybrid method and reconstruction of spacetime

The GW echo is primarily driven by the QNMs associated with the ringdown phase of the GW signal. Despite presenting a significant challenge, various techniques of extracting (analytical) echo templates have been established [99,100]. In this work, we employ the methodology outlined in [66], where the waveforms near the horizon of an ECO or QBH are reconstructed using asymptotic information contained in the Newman-Penrose scalars Ψ °0 and Ψ °4 at future null infinity, I þ . Below, we summarize the key concepts underlying the reconstruction of GW echoes.

The construction of [66] largely follows the hybrid approach [101,102] in which spacetime is divided into two regions as depicted in Fig. 1. For the larger, blueshaded region, BH perturbation theory is valid, while below the three-dimensional timelike tube Σ Shell , the system enters the strong-field regime. The exterior region is treated as linearly perturbed Schwarzschild, and for the interior, post-Newtonian (PN) theory is applied, matching the perturbed Schwarzschild metric at Σ Shell . 5 For the blue-shaded region, the no-incoming wave condition at I -is prescribed, while outgoing waves are imposed at I þ . The dynamics in this region are described by the Teukolsky equation. Based on the outgoing gravitational waveform, the perturbative field close to H þ is determined and subsequently, the echo waveform is derived. In general, the homogeneous solutions of the Teukolsky equation for Ψ 0 and Ψ 4 can be decomposed using spinweighted spheroidal harmonics. By setting the dimensionless spin parameter a ¼ 0, i.e., considering only mergers with nonspinning remnants, the solutions can be written using spin-weighted spherical harmonics s Y lm ,a s

The radial functions s R lmω ðrÞ can be determined by matching the solutions of the radial Teukolsky equations with solutions determined in BBH coalescence spacetimes, asymptotically given by [66]

where r à , again, is the tortoise coordinate. In Eq. ( 9), the top solutions approximate the behavior close to I þ , the bottom solutions yield the dynamics close to the BH horizon/ECO surface. A relation connecting the waves escaping to infinity, Z ∞ lmω and Y ∞ lmω (corresponding to Ψ °4 and Ψ °0, respectively), with the ingoing waves at the future horizon, H þ , Z in lmω and Y in lmω ,i sg i v e nb y[66]

where the coefficients D in lmω ;D out lmω ;C in lmω ;C out lmω are computed using the Black-Hole Perturbation Toolkit [103]. They encapsulate the physics of the BH potential barrier: the transitivity for radiation traveling from H -to I þ is given by 1=D out lmω (1=C out lmω ), the reflectivity of the potential barrier displayed in Fig. 1 by D in lmω =D out lmω (C in lmω =C out lmω ). Together, the reflective surface of the QBH or ECO and the BH potential barrier form a cavity. The amplitudes of inand outgoing radiation at null infinity, Z ∞ lmω and Y ∞ lmω , can be obtained from numerical relativity (NR) simulations using the Cauchy-Characteristic-Evolution (CCE) simulation pipeline [104,105] incorporated in the numerical relativity code SpECTRE [106,107].

The strain detected by a GW detector, morally located at I þ , can be decomposed into two components: the primary waveform of the BBH merger, where we denote its amplitude as Z ∞ lmω , and subsequent echoes, here denoted as Z echo lmω . The conversion from radial amplitude solution of the Teukolsky equation to GW strain at I þ is easily written down in Fourier space as

For the transformation of ingoing to outgoing solution, we rely on the Teukolsky-Starobinsky (TS) relations

where

An overview of the relations between the asymptotic information at I þ and at the horizon H þ for both ingoing and outgoing radiation is provided in Fig. 3 of [66].

For the computation of the echo, the relevant information is encapsulated in the ingoing radiation at the future horizon H þ , here represented by Y in lmω (see Fig. 1). This ingoing radiation travels toward the (partially) reflective surface, where part of it trapped in the cavity. To determine how each frequency mode interacts with the BH or ECO, both a reflectivity parameter and an appropriate boundary condition are required. With the reflectivity coefficients already established in the previous subsection, we now turn to the selection of a suitable boundary condition: For the ECO, we apply the same approach as in [66] relying on an earlier work considering the tidal forces applied on zero-angular-momentum fiducial observers [108].B y computing the transverse component of the tidal tensor field from the point of view of these fiducial observers, the ECO's boundary condition can be derived as

where the upper index on the left-hand side signifies that Eq. ( 14) computes the first echo only (one for each cycle the radiation undergoes within the cavity). Y in ECO lmω is determined by the ingoing Ψ 0 which takes the form of curvature perturbations of fiducial observers near the horizon and includes only the fraction of Y in lmω associated with the QNM of the ringdown phase. 6 In the observer's frame, the ECO gives rise to a local reflectivity R ECO and an additional conversion factor of ð-1Þ lþmþ1 =4. For the outgoing amplitude defined in (14) to reach I þ , this first echo has to be transmitted through the BH barrier and, thus,

A similar procedure can be applied to the QBH reflectivity model, for which the boundary condition has not yet been specified. From a phenomenological standpoint, there is no compelling reason to assume significant alterations to the ingoing wave beyond the possibility of a phase shift. All other effects are encapsulated in the model-dependent reflectivity coefficient R QBH . Therefore, inspired by the arguments that lead to Eq. ( 14), we postulate the boundary condition for the QBH to read 7

As mentioned above, the boundary conditions ( 14) and ( 16) represent only the first echo of the GW. However, as the potential barrier is not perfectly transmitting, a portion of the radiation Z out ECO lmω (Z out QBH lmω ) is reflected back (with the reflectivity given by D in lmω =D out lmω ) and the cycle starts over (see Fig. 5). The reflected quantity can again be transformed into ingoing radiation via the relations above. For simplicity, we include all necessary transformations as well as the numerical factor of ð-1Þ lþmþ1 into

representing the reflectivity of the potential barrier. Using the latter, we can compute the full echo as

In the latter equation, the total echo is rewritten as a sum over individual contributions, each associated with n cycles in the cavity formed by the barrier and the reflective surface.

To avoid instabilities of the ECO, one must ensure that jR ECO R BH j < 1 at all times. An exception to the inequality is the QNM, ω n , of the ECO (QBH) for which R ECO ðω n ÞR BH ðω n Þ¼1 (R QBH ðω n ÞR BH ðω n Þ¼1). The QNMs therefore appear as poles of the corresponding transfer function

The transfer functions for R ECO and R QBH are displayed in Fig. 2. We clearly see the difference in the QNMs for ECO and QBH. The global maxima, however, align as they both represent the fundamental QNM frequency of the Schwarzschild BH. The roots of the QBH's reflectivity are indicated by dashed green lines. For large frequencies, both transfer functions are suppressed exponentially. For small frequencies, the 1=D out lmω factor in Eq. ( 18) dominates the transfer function. For the ECO, both T QH and γ FIG. 2. Transfer functions for R ECO and R QBH with

The green dashed lines mark the characteristic absorption frequencies for the QBH. The dotted orange lines mark the QNM of the ECO. Similarly, the peaks of K QBH mark the QNMs of the QBH. For better readability, the transfer function of the QBH, K QBH , is multiplied with an overall factor of 5. 6 As Y in lmω is inferred from the full waveform at I þ , Y ∞ lmω , via (10), it naturally carries information about all regimes of the waveform. However, physically, only the ringdown part exists within the potential barrier and can fall down the horizon H þ . The procedure rectifying this issue is explained below. 7 We acknowledge that the boundary condition in ( 14) is fundamentally tied to the tidal deformability of ECOs. As classical BHs do not exhibit tidal deformability, the analogy is lacking here. Consequently, it is possible that BH could respond in a completely different manner. However, for the purposes of this study, we adopt a slightly generalized ECO-like boundary condition and leave the intriguing question of the appropriate QBH boundary condition to future investigations. determine its QNMs. Additionally, T QH determines the exponential suppression. For the QBH, β, δ impact the location of the QBH'sQ N M s ,ω QBH QNM ,w h i l eα shifts the roots and ϵ governs the exponential decay. The overall amplitude is suppressed by ζ. All parameters and their correspondence between reflectivity models are summarized in Table I.

The factors jR ECO j 2 and jR QBH j 2 represent the corresponding energy reflectivity of the ECO and QBH, respectively. Similarly, jR BH j 2 represents the energy reflectivity of the potential barrier [109]. Based on this interpretation, we can equally define a coefficient of transmission for the amplitude of ingoing radiation at H þ . We define, for ECO and QBH analogously,

such that, for the initial amplitude penetrating the ECO's surface (the QBH's reflective shell) and propagating toward the horizon H þ , we find

As for the reflected radiation sourcing the echo in Eq. ( 18), the transmitted portion receives a contribution for every cycle. Adding all subsequently transmitted radiation, we obtain

Equations ( 18) and ( 22) will be used in subsequent sections to determine the energy and angular momentum flux across both horizons I þ and H þ . We stress at this point that Eq. ( 20) does not fully determine T ECO as a complex phase remains undetermined. Technically, this issue can only be resolved by adding further assumptions to the model. However, in this work, we are solely interested in the energy and angular momentum fluxes. As Eq. ( 20) precisely describes the energy transmitivity, its definition is sufficient for the subsequent analysis.

Both the amplitude falling into the BH, Y horizon lmω , and the one escaping to I þ , Z ∞ lmω , are determined by Y in ECO lmω , the ingoing radial component of the Ψ 0 wave falling down the future horizon in the classical picture. Following [66], Y in ECO lmω can be determined by computing the late-time portion of Y in lmω that corresponds to the QNM, i.e., the ringdown phase of the full waveform. The relevant portion is determined by the retarded time parameter, v H Σ , associated to the shell Σ shell separating the two spacetime regions in Fig. 1. With v H

Σ at hand, one can define [66]

with F ðvÞ being the Planck-Taper filter as defined in [66].

Physically speaking, the Planck-Taper filter guarantees the desired inclusion of only the ringdown-affiliated part of the waveform. The shell time v H Σ is determined via minimizing the NR waveforms mismatch with a sample waveform purely constructed from QNM overtones with different initial times (see [66] for details).

With the algorithm sketched above, determining the final echo-corrected waveform numerically requires solely the NR data for Ψ °0 and Ψ °4 defined at I þ . As mentioned above, we obtain the NR data through CCE simulations. The latter adopts the conventional normalization, that is,

Here, M f is the mass of the remnant compact object. Equation ( 18) provides an expression for the strain of the GW echo, determined by a given reflectivity model and the properties of the original GW signal at I þ ,u pt oa frequency-dependent prefactor. At this stage, it is worth mentioning that this contribution adds on top of the original waveform. As these signals and the process of origin are disentangled, one can simply compute the complete waveform as

where h ∞ results from Z ∞ lmω , the raw asymptotic waveform at I þ . A depiction of h eff is provided in Fig. 3, showing the combination of echo and main strain in the top panel and

TABLE I. Model parameters for the reflectivity of ECOs and QBHs. An implementation of the transfer functions KðωÞ is displayed in Fig. 2. As the ECO's transfer function does not obtain roots, there is no correspondence between α, δ for QBHs. KðωÞ attribute Parameter ECO Parameter QBH Time separation echo γ;T QH β Exponential supp. for large ω T QH ϵ Location QNMs γ;T QH β; ϵ Separation of roots ÁÁÁ α "Sharpness" of roots ÁÁÁ δ Boundary suppression 1=41 =ζ

the isolated echo, i.e., without h ∞ in the bottom one. Thereby, h echo can be explicitly computed as a convolution of all aforementioned substeps, using Eqs. ( 10), ( 11), ( 18), (23), and (24) as

where

and we define Ψ 0;lm such that

where FðΨ 0;lm Þ ≕ Y ∞ lmω and Fð•Þ denotes the Fourier transform. Note also that C in lm represents the Fourier transform of C in lmω . An exemplary echo is displayed in the bottom plot, Fig. 3. The parameter dependence of QBH echoes is displayed in Fig. 10 of Appendix B (compare with Table I). A similar plot for ECOs is found in [66].As expected, the two types of echoes exhibit only slight variations in their morphology.

The statement in Eq. ( 25) is independent of the time separation that is determined by the corresponding parameters of a given reflectivity model, i.e., independent of the complex phase factors in ( 4) and (6). Therefore, the echo can easily be added to (numerical) waveforms by adding it on top of the raw strain time series in a postprocessing step similar to [110]. Crucially, the echo, as well as its memory, can be treated as individual "events." 8

III. MEMORY AND FLUX FOR (QUANTUM) GRAVITATIONAL WAVEFORMS

In the literature, numerous studies have explored the GW memory based on distinct approaches [4][5][6][7][8][9]. The equivalence of the resulting expressions for memory has been well established. For applications involving numerical waveforms, a specific formulation of the gravitational memory based on the so-called energy-momentum or Bondi-Metzner-Sachs (BMS) balance laws [67] (see also [111,112] for more technical details) has become widely adopted [110,113,114]. These balance laws provide a series of constraint equations that relate the GW strain at future null infinity to physical properties of the remnant compact object, such as its velocity and mass. Notably, they include expressions for both the linear and nonlinear GW memory. These flux balance laws have been applied in prior studies for waveform analysis and NR applications, among others [68,115,116].

A. Computing the memory

For the numerical analysis in this work, we utilize the memory components derived from the BMS balance laws as implemented in [110]. This choice is driven by the applicability of these equations to NR waveforms and the convenient format in which they are presented. Specifically, the memory components, as outlined in [110], enable direct computation of time-dependent memory corrections for each individual harmonic strain mode. Since the BMS laws are derived from full, nonlinear GR, we expect these FIG. 3. Comparison of the full waveform h eff ¼ h echo þ h ∞ (top) for ECO and QBH, and isolated echo (bottom) for QBH simulated for event SXS:BBH:0207. For the ECO, we choose T QH ¼ 2T H and γ ¼ 10 -15 . For the QBH, the parameter choice is displayed on top of the bottom plot. The latter displays the isolated h echo from the upper plot. 8 This approach aligns with recent proposals suggesting that echoes should be treated as independent signals rather than as extensions of the BBH merger waveform due to the potentially substantial time delays involved [82].

equations to apply equally well to any strain (quantum) corrections that arise within the regime of validity of GR as an effective theory.

The memory equation in [110] is directly derived from the balance flux laws as described in [67]. The time series containing only the memory-induced contributions to the strain, after a few simplifications, can be written as

where D -1 is defined via

and ð defines the spin-weighted derivative operator whose precise definition can be found in [67,110,112]. The "Δ" for a tensor or scalar thereby denotes the change in time, i.e., Δf ≡ fðu 2 Þfðu 1 Þ. The h appearing in Eq. ( 29) is the full GW strain. It is fundamentally dependent on the line of sight, spanning from the detector to the binary frame. The latter can be characterized by two angles θ; ϕ, marking a point on the celestial sphere. The initial time u i in Eq. ( 29) marks the beginning of the gravitational waveform. As, in this work, we consider numerical waveforms, by construction, we do not use the precise definition of Eq. ( 29) in which one would take the limit u i → -∞ [110] but instead identify u i with the starting time of the provided NR waveforms. This truncation must be compensated by an angular dependent constant κ. 9The full memory given by Eq. ( 29) can be separated into a dominating term associated with the energy flux and an oscillatory term,

While the latter equations both describe the displacement memory that can be directly related to the supertranslation part of the asymptotic symmetry group at I þ [e.g., [117,118], see [119] for a full derivation of (29) and following equations], one may also consider spin memory contributions [117,[120][121][122] for completeness. The spin memory in turn, is connected to the superrotations subgroup. Both superrotations and supertranslations each generalize the known translation and rotation symmetries included in the Poincaré group respectively by an infinite tower of "new" generators. For details, we refer to [117]. 10The spin memory can be brought to a similar form as Eq. ( 29). Then, rearranging the individual contributions of spin and displacement memory, one can write them as in [110,113],

Here, m and N correspond to the Bondi mass aspect and the angular momentum aspect, respectively [110], i.e., in terms of strain and Newman Penrose scalar they read

While, in principle, we can compute the deviations due to echoes present in the strain time series for all components given in (33), we are mainly interested in the components associated with the energy and angular momentum flux across the horizon I þ ,i . e . ,h E and h J . Both components can be related to the null memory (also called nonlinear memory) but have different parity, i.e., h E is electric and h J is magnetic. For practical purposes, it is even sufficient to consider solely h E as it largely dominates the nonlinear memory. Being interested in the kinematics of the considered binary, we nevertheless include the computation of the memory associated with the angular momentum flux. Applying the notation of Sec. II B, we can identify and separate contributions in h E and h J . We start by considering h E , where Eq. ( 25) leads to mixed terms of the form

Both h echo and h ∞ being time series waveforms, we assume that the echo, and thus also its time derivative, is time-wise clearly separated from the main waveform h ∞ as, for instance, depicted in Fig. 3. In this case, we find ḣ∞ ḣecho ¼ ḣ∞ ḣecho ¼ 0 and the memory h E separates into a waveform and an echo contribution. The assumption of time-wise separability requires a sufficiently large phase shift incorporated in the reflectivities R ECO , R QBH . The relevant parameter bounds can be quantified as γ ≲ 10 -4 (for T QH ¼ T H , and T QH ≲ 5T H for γ ¼ 10 -15 ) for the echo of an ECO, and for echoes induced by QBHs one finds β ≲ 10 -7 . For the remainder of this work, we restrict ourselves to a range in parameter space where this nooverlap assumption holds valid, such that

where we omit the constant α.

Next, we consider h J . Structurally, the main difference with respect to h E is the angular derivative ðh. We argue that the derivative operator only acts on the angular dependent part of the strain, i.e., we can decompose hðu; θ; ϕÞ¼

such that the strain modes, and thus Z in lmω (Y in lmω ), as functions of retarded (advanced) time, are not affected. Consequently, even with the angular derivative operator ð applied on the strain, the separation argument still holds, and terms like ðh ∞ ðh echo vanish due to the absence of nontrivial time series overlap. Thus, inserting (25) in h J we find that h echo does not mix with h ∞ and thus

where h echo J is simply

A comparison between h echo J , h echo E and h ∞ J and h ∞ E is displayed in Fig. 4. In this case, the time separation of the effects is clearly exhibited. Note that while for h E the contributions of waveform and echo are identical in shape (besides the amplitude), for h J they are fundamentally different. Note also that if the no-overlap assumption is violated, the steplike shape of the nonlinear memory receives non-negligible corrections from the mixing terms ḣ∞ ḣecho þ c:c: which are generally of oscillatory nature. The latter is a fundamental manifestation of the memory's nonlinear nature. Thus, a echo could lead to a significant alteration of the "known" memory contribution associated to BBH waveforms.

B. Computing fluxes

The balance flux laws not only provide a framework to incorporate the memory into the GW strain, but they also yield explicit fluxes of physical properties, such as energy and angular momentum, that pass through a detector placed at I þ or fall into the BH. Using Eq. ( 33), we can extract both the energy and angular momentum fluxes carried by FIG. 4. Selected modes of the memory corrected waveform. The dashed and dotted lines mark the memory contributions h E and h J from ECO and QBH directly related to the echoes. The "no-echo" counterpart is displayed in both the bottom and top plots as the solid orange line. The parameter choice and the event are copied from Fig. 3.

the GWs. Since these equations are separable into components corresponding to the initial strain and the echo-related contributions, so are the flux formulas. Concretely, for a generic strain, the (dimensionless) energy and momentum flux per unit time and angle read

The flux formula (42) are defined for asymptotic strains at I þ , hence the choice of the time coordinate, u. Strictly speaking, in Eq. ( 42) (and all memory definitions before), one should replace strain h with the asymptotic strain h ∘ defined exclusively at I þ , and which continues into physical spacetime via

The aforementioned definition becomes especially important when computing the dimensionful versions of (42).In this case, the asymptotic strain must be replaced by the physical strain. For practical examples, such as actual measurements of BBH mergers, the radius is replaced by the estimated luminosity distance D L of the event.

Naturally, in addition to the distance information, appropriate powers of c and G have to be factored in to maintain dimensional consistency. When computing (42), it may be helpful to make use of strain's decomposition into spherical harmonics outlined, for instance, in [68] (see also Appendix A). For instance, for the memory flux associated with the echo one computes

which can be decomposed into spherical harmonics of spinweight zero,

The resulting frequency-dependent expansion coefficients, up to a prefactor depending on l and m (see Appendix A for details), read

A similar result holds for the angular momentum flux J echo , expressed in terms of coefficients β echo lmω . The coefficients α echo and β echo differ solely by a mode-dependent factor, the explicit forms of which are provided in Appendix A.I ti s important to note that this decomposition applies equally to the memory formulas (36) and (41).

Considering the process of echo production in the cavity formed by the BH barrier and the ECO's( o rQ B H 's) surface, it is natural to ask whether the involved GWs obey a form of (energy) flux conservation for gravitational radiation. This can be straightforwardly deduced by considering Fig. 5. Staying consistent in our notation, we adopt the viewpoint of the ECO, but the same formalism holds for the QBH as well: The process starts with the ingoing gravitational radiation that first interacts with the reflective surface, i.e., Y in ECO;1 lmω . This radiation is partially reflected by the surface, producing the outgoing wave Z out ECO;1

Naturally, the energy flux conservation equation can be expressed in terms of the amplitudes at each stage of the GW's interaction with the barriers or horizons. For instance, combining the configuration depicted in Fig. 5 with the radial solutions close to the ECO, one could also write

which is similar to the formulations found in [108,124]. Note that in Eq. ( 48),t heup solution was complemented by a set of ingoing waves at infinity, Y ∞ in lmω . Then,

At this point, it is crucial to emphasize that the notation for the radial solutions might vary across literature, especially when the Sasaki-Nakamura (SN) formalism [125] is used instead of the Teukolsky framework. For comprehensive treatments of echo computations using the SN formalism, we refer to works such as [108,109,124,126].

For the fluxes reaching I þ , the definitions given in Eq. ( 42) in combination with the asymptotic strain yield the correct result. In case the asymptotic strain is not directly available, we can use Eq. ( 11) to extract it from Ψ °4. Defining the flux that falls down into the BH or ECO, or the initial flux directed toward the ECO's surface, on the other hand, is less trivial. Concretely, for the energy carried by Y in ECO lmω one needs to convert Ψ 4 to Ψ 0 (as it is associated with the perturbations of the shear on the horizon) using the TS identities after which the energy flux can be derived from the area change according to Hartles formula of BH area increase. After a lengthy computation, one finds [108,124]

For the events relevant to this article, the dimensionless spin a of the ECO/BH is (mostly) negligible, thus k ≈ ω and r þ ≈ 2M. A similar result can be obtained for the reflected quantity Z out ECO lmω (see Fig. 5), i.e.,

With the defined reflectivities R ECO and R BH , the above energies can be related as follows:

and thus

Finally, we can express the energy flux across H þ in terms of the same quantity at I þ , namely,

where E in ECO can be replaced using Eq. ( 52). The energy loss per cycle in the cavity depicted in Fig. 5, as a function of frequency, then corresponds to

where jT BH j 2 ≔ 1 -jR BH j 2 . The coefficients Γ, Θ represent the fraction of energy loss due to gravitational radiation crossing the ECO/QBH horizon and the potential barrier, respectively. Equation ( 54) enables the explicit computation (and comparison) of the energy crossing the horizon H þ for a given event. The latter becomes particularly relevant in the context of ECOs as there may exist instances of those collapsing to form a BH after the first echo (first cycle) [58].

If the BH is considered in a classical framework, no echoes will be received by the detector at I þ . Considering the formed BH to be represented by a QBH, the transfer function of subsequent echoes changes, leading to a distinct echo morphology. Analogous conservation laws can be established for the angular momentum flux toward the future horizons H þ and I þ . Together, the energy and angular momentum fluxes across the BH horizon influence the entropy of the dynamical BH during the ringdown phase, thereby affecting its area in accordance with the first law of BH thermodynamics. Therefore, in principle, it should be feasible to relate the energy and angular momentum flux induced by the ringdown to the entropy of dynamical BH (for recent progress see [127]). We will defer such an investigation to future work.

IV. NUMERICAL STUDY

With the analytical descriptions of the GW echo and its associated memory in hand, we now aim to quantify its impact on the waveform detected by relevant GW interferometers. Specifically, we focus on LISA, which is particularly well suited for potential memory detections due to its high sensitivity in the low-frequency band, where the memory effect is typically stronger. To evaluate the SNR of the echoes and the corresponding memory effects, we utilize the pipeline developed in [77]. Thereby, a realistic prognosis for echo memory detection necessitates consideration of various event-related factors, including the orientation and sky position of the merger relative to the LISA frame. For the sky position and all other orientationdependent features pertinent to the SNR, we adopt the conservative baseline parameters outlined in Table I of [77].

In this section, we utilize multiple SXS simulations with vanishing or nearly vanishing spin, that are SXS:BBH:0205, 0206, 0207,1424,1448,1449,1455,1936.E x c e p t i o n st o the low-spin prescription are given SXS:BBH:0334, 1155, and 2108 which have a remnant spin amplitude j⃗ χj between 0.28 and 0.68. While we acknowledge that the algorithm computing the echoes may introduce systematic errors for these events, they will only play a secondary role in this analysis as corrections for nontrivial spin amplitudes are generally small. Regarding the reflectivity models, a nonnegligible spin introduces line-broadening as discussed above and in [31]. For the relevant cases however, the effective broadening remains close to the numerically imposed frequency resolution of the applied pipeline, i.e., Oð1Þ μHz. Thus, without loss of generality, we include the corresponding events into our consideration serving as consistency checks against potential parameter biases regarding the spin of the remnant object. For each event, the echo is numerically computed solely for the harmonic strain modes h 2;AE2 due to their overall dominance. Note however that, particularly for the flux computation, all modes are included. Thus not all strain modes subject to the calculations below carry an echo. The implications of restricting the echo to the h 2;AE2 mode will be addressed below.

A. Echo-induced gravitational wave memory

To get an intuition for the significance of the echoinduced memory based on the QBH model, we simulate different reflectivity functions and compute the memory for each listed waveform. The results are exemplarily displayed for SXS:BBH:1936 in Fig. 6. The event is chosen to match the results displayed in the companion paper [81] on investigations regarding the general echo detectability prospects. The contour plots displayed in Fig. 6 show the fraction of memory attributed solely to the echo in comparison to the classical waveform's memory (top panel) and the corresponding SNR associated with the additional memory contribution (bottom panel). 11 This SNR pertains exclusively to the detectability of the memory and is independent of the waveform itself. To facilitate a clearer interpretation, we present the relative increase in memory, thereby effectively factoring out the redshift and mass dependencies. For a comprehensive exploration of the detectability of memory concerning the latter two remnant-specific parameters, we refer to [77]. As anticipated, Fig. 6 reveals a discernible trend indicating that lower suppression from the boundary FIG. 6. Memory amplitude and SNR gain for event SXS: BBH:1936 under the reflectivity model of a QBH with α ¼ 8π, β ¼ 10 -15 and δ ¼ 0.1. The top plot shows the memory amplitude purely due to the echo as a fraction of the classical waveform memory (without echo). The bottom plot shows the gain in overall memory SNR relative to the SNR of the waveform memory without echo (at redshift z ¼ 1 and mass M ¼ 10 6 M ⊙ the total value for the SNR for the chosen baseline is roughly 10). Note that for better readability, the x-axis of the bottom plot is extended up to ϵ ¼ 4.

condition and the reflectivity function R QBH , i.e., smaller ξ and larger ϵ [compare to Eq. ( 6) with T QH ¼ ϵT H ], are associated with higher echo memory contributions. Therefore, the echo's memory contribution is proportional to its amplitude. Furthermore, the contour plot indicates that over a broad range of parameter space, the echo significantly contributes to the overall memory, achieving SNR levels comparable to the initial waveform's memory in extreme cases. While the SNR for the echo memory alone is likely to be insufficient for individual detection, it is plausible that the echo memory enhances the overall memory SNR. This enhancement is heavily contingent upon the specific event and the temporal separation between the echo and the ringdown of the waveform. The sensitivity to time separation is primarily due to LISA's response with respect to the steplike increase characteristic to the memory. The closer the two memories, sourced by the waveform and the echo, reside in time, the more pronounced the features' low frequency content appears in LISA data (see [77] for more details). Further, as elaborated above, terms such as ḣ∞ ḣecho þ c:c: now appear in the nonlinear memory leading to a nonlinear increase. Conversely, if this time separation is too pronounced, the synergistic effects between the waveform's and echo's memory diminish. This phenomenon is particularly evident in the analysis of the ECO model, which will be discussed further below.

It is crucial to highlight that, in parameter regions where the echo memory constitutes more than Oð10Þ percent of the waveform memory, the oscillatory part of the echo's time series itself becomes large enough to fall within the detectable range of LISA. This detectability heavily relies on factors such as redshift and mass, which were explored in [81]. For a reference on the pure echo SNR for a parameter region comparable to that displayed in Fig. 6,we refer to Fig. 5 therein. In this context, we find that for larger ϵ and smaller ξ's, the echo's magnitude increases significantly faster than the memory contribution. One finds that once the echo amplitude reaches a threshold of approximately ≳5% of the waveform amplitude, LISA'sS N R becomes sufficiently high to enable individual echo detection (given the right redshift and mass range of merger events), rendering the echo-induced memory as a primary tool of echo detection redundant. Nevertheless, we find it surprising that an event of such short time span compared to the original waveform exerts such large contributions to the memory.

Transitioning to the ECO model, we present the parameter dependence of the echo memory contribution, its SNR, and the echo amplitude in Fig. 7. The maximal echo amplitude here serves as a gauge for comparisons between echo memory and pure echo SNR. Due to its simplified parameter dependence, instead of focusing on a single event, we include all SXS simulations listed above in the plot. It is evident that while all events exhibit a trend similar to that of SXS:BBH:1936 for the QBH, there are notable event-dependent variations in the echo and echo memory amplitude ratios. These variations arise from distinct simulation parameters such as mass ratios and initial spins, leading to a diverse range of QNM content and, consequently, echoes. Figure 7 further displays crucial difference between the reflectivity models and their associated parametrization. Specifically, while the response and transfer functions for the ECO and QBH are relatively similar (as illustrated in Fig. 2), an increase in T QH for the ECO not only affects the exponential suppression but also compresses the time interval between the waveform ringdown and the echo. This results in a nonlinear response of LISA in terms of the memory SNR in relation to the linear increase in the amplitude of the echo memory contribution when T QH is varied. A similar phenomenon can be observed for the QBH if β is adjusted accordingly.

Most importantly, when comparing the results for the QBH model and the ECO, we find no significant differences in the memory contributions, provided that the parameters are selected to ensure comparable amplitudes of the transfer functions. This observation suggests that, despite the distinct underlying physics of the two models, their impact on gravitational wave memory is nearly identical. Validating this hypothesis further, we examine the remaining parameters, i.e., γ in R ECO as well as α, β and δ in R QBH .Whileγ and β primarily influence the QNMs of the ECO and the QBH, respectively, (not to be confused with the QNMs of the ringdown, which feed the echo), α and δ modulate the position of the characteristic frequencies as well as the distinctiveness of their cusps. Our findings confirm that, as anticipated, none of these parameters significantly affect the FIG. 7. Comparison of the echoes maximal magnitude, memory, and memory SNR with respect to corresponding quantities for the classical waveform (without memory) for the ECO model. Each line corresponds to a given event of the list above, where the color intensity is arranged such that the lines fade according their value for the relevant fraction. The graphs represent the ECO analog to Fig. 6.

echo memory with the exception of δ. If large enough, δ can reduce the amplitude of the echo memory due to the increased absorption by the BH horizon. Although this feature would enable a distinction between reflectivity models solely based on the amplitude of the resulting echo memory, given the current detection estimates, it is unlikely that this amplitude can be determined with sufficient accuracy.

For the tests of different reflectivity model parameters, the following parameter ranges were considered: γ ∈ ½10 -15 ; 10 -4 , α ∈ ½4 log 2; 8π,andδ ∈ ½0; 1. The upper bound for γ was chosen such that the separability condition remains intact, the lower bound is arbitrary as it would simply delay the echoes arrival further. For α,w e chose the interval depending on relevant values in literature [31]. Finally, for δ the lower bound marks the level at which ECO and QBH transfer functions are indistinguishable at high frequencies. The upper bound establishes a critical value beyond which the characteristic frequencies' cusps are sufficiently pronounced to impact the transfer function's overall amplitude.

B. Energy and momentum flux across I +

To obtain a first characterization of the echo in terms of fluxes, we compare the time-integrated energy and momentum flux of the echo for ECO and QBH against the corresponding fluxes of the waveform at I þ . The relevant parameters are varied just as in the previous subsection. For the QBH, Fig. 8 displays the integrated energy and momentum fluxes measured morally at (u → þ∞). Numerically, this boundary is approximated by the maximum time extent of the simulated SXS waveform.

The values presented in Fig. 8 indicate that, within the parameter region for which the memory contributions is non-negligible, the fractional integrated energy flux exhibits behavior that closely resembles the fractional memory. Specifically, the energy flux associated with the echo can comprise a considerable fraction of the total energy carried by the complete waveform. A similar pattern is observed for the momentum flux; however, it is noteworthy that, for the tested events, the echo carries less momentum than energy. A similar relation between energy and momentum flux is evident for events with higher remnant spins as well, as demonstrated by Fig. 9. Throughout these events, waveforms with larger energy flux also exhibit larger angular momentum flux, although the latter is generally much smaller in magnitude when compared to the full waveform's flux. We find that this behavior is uncorrelated with respect to the remnants spin. Physically, the observation that the echo carries less angular momentum flux than the initial waveform is sensible due to the generating FIG. 9. Integrated energy and momentum flux for the listed events under the reflectivity model of an ECO. Again, each line corresponds to a given event of the list above, where the opacities are arranged such that the lines fade according to the largest value of the relevant fraction. The graphs represent the ECO analog to Fig. 8.

process of the respective strain signals. While a coalescing binary is generally expected to radiate away large amounts of angular momentum, the perturbed Schwarzschild BH is not.

In Fig. 9 the analysis of the ECO reflectivity model was extended to more extreme ratios T QH =T H to demonstrate that the energy flux due to the echo can as well dominate the one of the waveforms. It is important to highlight that in this scenario, the individual echoes computed for the ECO model [see Eq. ( 18)] transition into a continuous signal due to the decreasing time separation between individual echoes for larger T QH . The strain associated with this situation is exemplarily illustrated in Fig. 12 of [66].

As for the memory, our findings indicate that the parameters γ, α,a n dβ do not substantially impact the energy or angular momentum flux when varied within the previously defined domains. Instead, Fig. 9 demonstrates that the fluxes carried by the echo strongly vary among simulated waveforms, leading to the conclusion that the direct and indirect features of GW echoes are fundamentally determined by the classical QNM content of the remnant, given the reflectivity models presented here capture the physical reality of these systems sufficiently accurately. The latter is supported by the equations above, where Ψ 0;lm ðvÞF ðvÞ represents the information encapsulated in the ingoing QNMs. This observations implies that, generally, the echo, its memory and associated fluxes are primarily determined by the physical properties of the remnant compact object. The reflectivity primarily affects the amplitude. Further support of this statement is provided by Appendix C.T h e r e f o r e ,t h e selected models can be considered as robust in terms of potential errors due to additional phenomenology.

The above result equally holds when the echo is computed for all numerically accessible strain modes instead of only h 2;AE2 . Naturally, the fluxes increase slightly when the echo's mode content is extended. This, however, has to be anticipated given the sum over strain modes in the flux determining factors α echo (β echo ) in Eq. (A2) [Eqs. (A4) and (A6)].

Finally, we turn to the computation of an energy balance (54) that can be established for the above fluxes. We are particularly interested in how the energy measured at the GW detector compares to the energy absorbed by the compact object or the BH at H þ . To establish a flux balance, we compute Eq. ( 49) and, in combination with Eqs. ( 51)- (54), estimate the energy consumed by the QBH/ECO in comparison to the energy radiated to I þ . Our findings reveal that, even when adjusting the reflectivity of the remnant object to boost the emission of large echoes, the overwhelming majority of energy-over 99% of the total energy loss per cycle-is absorbed by the celestial body itself. In other words, for a stable cavity between the event horizon and the surrounding potential barrier, nearly all of the energy E in ECO is returned to the compact object, whether it is a QBH or an ECO.

To quantify the impact of characteristic absorption frequencies as the major distinguishing factor between the reflectivity models studied in this work, the reflectivity parameter δ is varied. Thereby, the remaining parameters are fixed such that the ECO and QBH models yield identical transfer functions up to the cusps associated with ω N . It is found that the flux is highly sensitive to δ. Already for δ ¼ 0.01, the QBHs flux decreases by Oð1Þ percent with respect to the ECOs flux. This trend continues roughly linearly and holds equally for both energy as well as momentum flux. By increasing δ further, the QBH reflectivity dramatically decreases, until, for δ ≫ 1,t h e remnant turns into a classical BH as all energy is absorbed by H þ . The sensitivity to δ suggests that the flux measurement or equally the amplitude of the echo-related memory contribution could in principle be used to distinguish the origin of echo-like features in the gravitational waveform, given the strain time series of such an echo reaches sufficient SNR. It is however unlikely that the precision on the echo memory amplitude outperforms the SNR of characteristic frequency features in LISA'st i m e delay interferometry (TDI) data, as it is investigated in [81]. A investigation of synergizing echo strain and memory in LISA data as well as detection prospects for more memory focused instruments is left to future work.

Finally, we test which of the cavity's "walls" acts as more confining regarding the energy traversing it. The latter can be quantified by the factors Γ and Θ in Eq. ( 54) denoting the energy contributions toward the BH horizon, H þ , and future null infinity, I þ , respectively. We find that the signal entering the cavity as Y in ECO lmω in fact encounters the horizon of ECO and QBH and the potential barrier with roughly equal integrated transmitivities. 12 The reflectivity, however, integrated over the frequency band is less than a Oð1Þ percent, leaving the echo with a minimal amount of energy after escaping the cavity toward I þ ,as depicted in Fig. 5. The large amount of absorbed energy implies for a large class of events, according to [58],t h e ECO may transform into a QBH within the first cycle of radiation traversing the cavity. This implies that even when features of the characteristic frequencies are detected either in form of cusps in TDI data [81] or in reduced echo memories, it is not immediately clear if the remnant object initially was a QBH or a ECO. Therefore, the detection of echos with characteristic frequencies cannot rule out the existence of ECOs.

V. DISCUSSION

The GW memory effect will likely enter the regime of detectability of future space-based detectors such as LISA [76,77] as well as of next-generation ground-based detectors such as Cosmic Explorer [128] or the Einstein Telescope [129]. Besides testing the very fundamental nonlinearity of GR, it simultaneously establishes a smoking gun for deviations beyond GR [17] and other subtle effects that can affect the permanent displacement of freely floating test masses. The primary challenge in detecting new and subtle features in gravitational waveforms lies in the need to generate highly precise GW templates for extracting signals from noise-plagued data. Each additional feature incorporated into the template significantly expands the parameter space of the simulated waveform, resulting in a steep increase in computational costs. In this work, we explored one type of feature induced by quantum phenomenology and the potential existence of ECOs that can easily be added to the waveform without drastic changes on waveform templates. The quantum phenomenological approach was rooted in basic area quantization arguments, which facilitated the development of a reflectivity model for QBHs. This model was compared against established models for ECOs. For both models, corresponding GW echoes were constructed using the methodology outlined in [66]. We subsequently computed the memory and flux contributions of these echoes. We provided semianalytical expressions for the gravitational memory effect induced by the echo as well as the corresponding fluxes. We further formulated an energy balance for the ingoing radiation toward H þ during and after the ringdown phase of the merger. Numerical investigation of both the echo-induced memory, focusing on its amplitude and SNR, and the flux balance were conducted, considering in-depth different scenarios within a reasonable parameter space and the resulting phenomenological consequences.

Our results most notably reveal that memory corrections induced by GW echoes can reach levels of significant contributions with respect to the main waveform'smemory . Although, for a large region of parameter space, the echo memory is too small to be detectable individually, it can however boost the total memories of SNR considerably. Further synergies are observed when the time separation between the echo and the merger waveform are small. Physically, this corresponds to the reflective shell being located closer to the BH potential barrier. Thus, deviation from the expected memory could imply the existence of echoes, offering an indirect detection avenue that does not require extensive waveform model adjustments. Given the results of [81], in practice a combined search of memoryand strain time series-related features of echoes is likely to provide the most robust detection prospects, given the dominance of time series SNR in comparison to the induced memory SNR for a given echo. Generally, the time separation of echoes is highly model dependent, and there is no consensus on its exact magnitude. Some studies even propose searching for rogue echoes, which cannot be clearly associated with a specific merger event [82]. Similar searches testing mainly outlier events or instrumental glitches have been proposed in other contexts as well, such as, for instance, in investigations of topological dark matter (see for instance [130] and references therein).

Regarding the flux analysis, our investigation demonstrates that the majority of energy and angular momentum stored in the ingoing Ψ 0 heading toward H þ will ultimately be absorbed by the BH or ECO. Thus, one does not expect drastic phenomenological or astrophysical consequences from the reflective shells installed in this work besides the subtle echos in the interferometer data. Note further that the frequency spectrum considered in this work applies only to GWs. For other types of radiation, the implications of, for instance, the area quantization arguments above may be dissimilar, see [31] and references therein. Yet, despite transporting only a negligible fraction of energy and momentum flux of the ingoing Ψ 0 to future null infinity I þ , the echo's contribution to the total energy and angular momentum of the waveform can still be significant. We further note that the echo consistently carries less angular momentum than energy compared to the classical waveform across all tested events. This trend is largely independent of the specific reflectivity parameters of the system and matches the expectations based on the particular dynamics that give rise to the echo.

It is important to acknowledge that the echo memory's manifestation within the gravitational waveform is roughly model independent. Solely its amplitude can differ based on the selection of the reflectivity model. In this work, the distinction between models was encapsulated in the reflectivity parameter δ. Therefore, in principle, the echo-induced memory corrections establish an easy-to-add but powerful feature in waveform templates for future GW detections. In this context, our findings suggest that the reflectivity models only play a secondary role for the exact shape of the echo time series as well as the associated memory. It is rather the dynamics of the binary configuration and the resulting QNM spectrum of the remnant body that determines the fluxes and h echo . Measuring an echo therefore, in principle, establishes another pathway toward studying QNMs in the context of quantum BHs.

The model-independence of the memory does not come without a price. We emphasize at this point that a major challenge of identifying the echo-induced memory is its distinction with respect to other memory corrections. Such include significant corrections due to deviations from GR [17] and/or quantum effects unrelated to echoes [25,26]. Further theoretical and analytical investigations of such effects would help identify distinct fingerprints of the potential features in both the memory and the oscillating part of the strain time series. Another theoretical aspect that has been only partially addressed in this work, as well as in the companion paper [81], and that warrants further exploration, is the BH reflectivity model itself. Although the arguments laid out in Sec. II (primarily based on [81]) are consistent with phenomenological and observational constraints, advancing the study of gravitational wave GW echoes demands reflectivity models rooted in quantum information theory as well as astrophysics. Specifically, it is crucial to ensure that these models do not violate the core principles of BH information theory. In particular, the boundary conditions for QBHs should be the subject of further investigation. One has to acknowledge that the prescription given in Sec. II B is oversimplified. The ingoing radiation toward the BH horizon might, in fact, introduce a feedback term in the Teukolsky equation. Another alternative could involve GW scattering, rather than reflection off the horizon. These crucial aspects of QBH models must be thoroughly explored in future research to ensure a more accurate understanding of their physical behavior and the corresponding GW signatures.

In conclusion, we find that echoes and their signatures provide a smoking gun for quantum corrections to the BH's horizon as well as the existence of ECOs, provided a sufficient measurement precision. As their features appear in the nonlinear GW memory, potentially among numerous other phenomenological imprints, the memory's detection has the potential to be the next milestone in gravitational physics and to fundamentally change our understanding of BHs and other compact stellar objects.

which is of spin-weight 1 and where we find The latter two decompositions cover all terms appearing within the angular momentum flux (42). The four terms present in second equation in (42) differ only by a time derivative of h l 1 m 1 or h l 2 m 2 .

APPENDIX B: PARAMETER DEPENDENCE OF QBH ECHOES

The QBH echo is exemplarily shown for different parameter configurations in Fig. 10. The small noiselike oscillations are due to numerical issues with handling the transfer functions poles, representing the QNMs of the QBH. For the results of Sec. IV, these oscillations are irrelevant.

APPENDIX C: QUASINORMAL MODES AND OVERTONES

As described in the main text, the radiation that initially enters the cavity, Y in lmω can be rewritten in terms of Schwarzchild QNMs ω n and the corresponding overtones. One way of expanding Y in lmω in the relevant time domain, i.e., for v>v H Σ is given by

ðA n e -iω n v þ B n e i ωn v Þ; ðC1Þ

where ω n as well as the overtones A n , B n are complex. The latter equation is used to determine v H Σ in the echo construction pipeline via a least mismatch scheme. Using Eq. (C1),t h e overtones can directly be related to the echo strain as

ðA n e -iω n v þ B n e i ωn v ÞðvÞ ¼ 1 ω 2 KðωÞ

This implies that the amplitude and phase of the overtones is nontrivially modulated by the factor ω -2 KðωÞ. Since in this particular work, KðωÞ is not analytical, as exact analytical mapping between the modes and overtones of h echo 22 and Y in 22 is postponed to future explorations. Numerically, the latter is worked out in Eq. (C1) by extracting the resulting echo modes via a similar fit analysis as for Y in lmω in the main text. In practice, this corresponds to casting the echo time series into a form such as

ðA echo n e -iω n u þ B echo n e i ωn u Þ; ðC3Þ

where u echo marks time after the merger at which the echo arrives at the detector. Note that, when applying this fit function the residual error in the mismatch is comparably large, i.e., of order Oð0.1Þ percent. This is partially due to including only overtones up to n max < 6 in the fitting for both Y in 22 an h echo 22 . Therefore, the values for the echo overtones displayed in Fig. 11 should be taken with a grain of salt.