Cauchycharacteristic evolution CCE is a powerful method for accurately extracting gravitational waves at future null infinity. In this work we extend the previously implemented CCE system within the numerical relativity code SpECTRE by incorporating a scalar field. This allows the system to capture features of beyondgeneralrelativity theories. We derive scalar contributions to the equations of motion Weyl scalar computations Bianchi identities and balance laws at future null infinity. Our algorithm tested across various scenarios accurately reveals memory effects induced by both scalar and tensor fields and captures Prices powerlaw tail ul2 in scalar fields at future null infinity in contrast to the t2l3 tail at future timelike infinity.
['General Relativity and Quantum Cosmology']
EinsteinKleinGordon system via Cauchycharacteristic evolution Computation of memory and ringdown tail
{'Einstein-Klein-Gordon system via Cauchy-characteristic evolution: Computation of memory and ringdown tail': "Sizheng Ma , 1, ∗ Kyle C. Nelli , 2 Jordan Moxon , 2 Mark A. Scheel , 2 Nils Deppe , 3, 4, 5 Lawrence E. Kidder , 5 William Throwe , 5 and Nils L. Vu 2 1 Perimeter Institute for Theoretical Physics, Waterloo, ON N2L2Y5, Canada 2 TAPIR 350-17, California Institute of Technology, 1200 E California Boulevard, Pasadena, CA 91125, USA 3 Laboratory for Elementary Particle Physics, Cornell University, Ithaca, New York 14853, USA 4 Department of Physics, Cornell University, Ithaca, New York 14853, USA 5 Cornell Center for Astrophysics and Planetary Science, Cornell University, Ithaca, New York 14853, USA (Dated: September 11, 2024)  \nCauchy-characteristic evolution (CCE) is a powerful method for accurately extracting gravitational waves at future null infinity. In this work, we extend the previously implemented CCE system within the numerical relativity code SpECTRE by incorporating a scalar field. This allows the system to capture features of beyond-general-relativity theories. We derive scalar contributions to the equations of motion, Weyl scalar computations, Bianchi identities, and balance laws at future null infinity. Our algorithm, tested across various scenarios, accurately reveals memory effects induced by both scalar and tensor fields and captures Price's power-law tail ( u -l -2 ) in scalar fields at future null infinity, in contrast to the t -2 l -3 tail at future timelike infinity.", 'I. INTRODUCTION': 'The LIGO, Virgo, and KAGRA collaboration [1-3] has detected about 100 gravitational wave (GW) events [4-7], providing us with unprecedented opportunities to explore the strong gravitational fields surrounding compact binary systems. This marks a new era in testing general relativity (GR) [4, 8-28].  \nThe progress mentioned above primarily centers on the Cauchy formalism, where the equations of motion are formulated as initial boundary value problems. Since it is not feasible to simulate an infinitely large space, the spatial Cauchy domain is typically truncated at a finite distance from the source, which prevents direct access to GWs at future null infinity. One approach to address this limitation is to measure waveform quantities at finite radii within the Cauchy domain and then extrapolate them to infinity [52]. However, this extrapolation does not guarantee that the Einstein equations are satisfied, making it fail to capture memory effects [53]. A more rigorous solution is to use Cauchy-characteristic evolution (CCE) and Cauchy-characteristic matching (CCM) [54], which can faithfully extract GWs at future null infinity up to Bondi-Metzner-Sachs freedom. Recently, a CCE [55, 56] and a CCM [57] algorithm for GR have been built in a NR code SpECTRE [58]. However, a robust CCE/CCM code for beyond-GR theories is still missing.  \nTo rigorously test GR, it is crucial to have precise GW predictions for both GR and alternative theories of gravity. Numerical relativity (NR) stands as the only ab initio method capable of producing complete inspiral-merger-ringdown waveforms to achieve this goal. Although two decades have passed since the first numerical simulation of a binary black hole (BBH) merger in GR [29], evolving modified theories of gravity remains challenging, as reviewed in [30]. The primary difficulty lies in the potential modification of the principal part of the Einstein equations, rendering the original numerical methods for GR ill-posed. Several approaches have been proposed to address this issue. The first approach involves directly formulating well-posed schemes for evolving the full set of evolution equations in certain specific theories, such as weakly-coupled Horndeski gravity theories (e.g. Einstein-scalar-Gauss-Bonnet) [31-35] based on the modified generalized harmonic formulation [36, 37], as well as Damour-Esposito-Farèse scalar-tensor theories [38-42]. The second approach is to expand the equations perturbatively in the coupling constant and solve them order by order. This order reduction scheme has been applied to the evolution of dynamical Chern-Simons [43-45] and scalar Gauss-Bonnet gravity [46, 47]. The final treatment introduces auxiliary variables to fix high-energy (short length/timescale) degrees of freedom, while leaving low-energy parts unchanged, see [48-50]. Recently, Corman et al. [51] compared the results of these three approaches under various physical setups.  \nIn this paper, we take a step to let the SpECTRE CCE system capture beyond-GR features. Specifically, we choose to include a (massless) scalar field, which serves as an important ingredient in theories such as Horndeski [59] and  \ndynamical Chern-Simons gravity [60, 61]. We will be testing the accuracy of our code and demonstrating its ability to reveal ringdown tails [62, 63] in the scalar field, as well as the memory effects sourced by the scalar field.  \nThis paper is organized as follows. In Sec. II, we outline the details of CCE for an Einstein-Klein-Gordon system. In particular, we show how the equations of motion, Weyl scalar computations, and Bianchi identities are modified in the presence of a scalar field. Next in Sec. III, we focus on the memory effects sourced by the scalar field and the corresponding balance laws. We then test our code by evolving scalar fields on two prescribed spacetime backgrounds (Sec. IV) and performing the full CCE procedure to compute tails and memory effects (Sec. V). Finally, we summarize the results in Sec. VI. Throughout this paper, complex conjugates are represented by overlines.', 'II. EINSTEIN-KLEIN-GORDON SYSTEM VIA CCE': 'The Einstein-Klein-Gordon action offers the most basic description for a system containing both scalar and tensor fields. For instance, after performing a conformal transformation [64], Damour-Esposito-Farèse scalar-tensor theories [64-66] in the Einstein frame are governed by:  \nS = ∫ d 4 x √ -g [ R 16 π -1 2 ∇ µ ψ ∇ µ ψ + V ( ψ ) ] , (1)  \nwhere the real-valued scalar field ψ is minimally coupled to the metric. The action leads to  \nR µν = 8 π ∇ µ ψ ∇ ν ψ, □ ψ + ∂V ∂ψ = 0 .  \n(2)  \nHere, we see that the Einstein field equations obtain a source term, while the scalar field obeys the Klein-Gordon (KG) equation. In our following discussions, we assume V ( ψ ) = 0 for simplicity (rendering ψ massless), though these discussions can be easily extended to a complex-valued massive scalar field with an arbitrary potential.', 'A. Equations of motion': 'By adopting the Bondi-Sachs metric of an asymptotically flat spacetime [67-69]  \nds 2 = -e 2 β ( rW +1) du 2 -2 e 2 β dudr + r 2 h AB ( dx A -U A du )( dx B -U B du ) , (3)  \nwith x A = ( θ, φ ) standing for angular coordinates, the Einstein field equations can be converted to a set of hypersurface equations [70, 71]  \n∂ r β = S β ( J ) + 2 πr ( ∂ r ψ ) 2 , ∂ r ( r 2 Q ) = S Q ( J, β ) + 16 πr 2 ∂ r ψ ð ψ, ∂ r U = S U ( J, β, Q ) , ∂ r ( r 2 W ) = S W ( J, β, Q, U ) + 2 πe 2 β [ J ( ¯ ð ψ ) 2 + ¯ J ( ð ψ ) 2 -2 K ¯ ð ψ ð ψ ] , ∂ r ( rH ) + L H ( J, β, Q, U, W ) H + L ¯ H ( J, β, Q, U, W ) ¯ H = S H ( J, β, Q, U, W ) + 4 π e 2 β r ( ð ψ ) 2 , (4)  \nwith U = U A q A , Q = r 2 e -2 β q A h AB ∂ r U B , H = ∂ u J , and q A = ( -1 , -i csc θ ). The expressions of S β , S Q , S U , S W , S H , L H , and L ¯ H as functions of ( J, β, Q, U, W ) can be found in Sec. IV of [55]. The spin-weighted derivative operators ð and ¯ ð are defined to be  \nð ψ = q A D A ψ, ¯ ð ψ = ¯ q A D A ψ (5)  \nwhere D A is the angular covariant derivative compatible with the unit sphere metric. On the other hand, the KG equation for the scalar field ψ becomes [70]  \n2 r∂ r Π+2Π = 1 r ∂ r [ r 2 ( rW +1) ∂ r ψ ] + e 2 β r ( N ψ 1 -N ψ 2 + N ψ 3 ) -N ψ 5 -r 2 N ψ 4 , (6)  \nwhere we have defined an auxiliary variable Π = ∂ u ψ , and  \nN ψ 1 = 2 K Re ( ð β ¯ ð ψ ) + K ¯ ðð ψ, N ψ 2 = Re ( ¯ J ðð ψ +2 ¯ ð β ¯ ð ψJ + ¯ ð J ¯ ð ψ ) , N ψ 3 = Re ( ð K ¯ ð ψ ) , N ψ 4 = 2Re ( ¯ ð ψ∂ r U + ð ¯ U∂ r ψ +2 ¯ U ð ∂ r ψ ) , N ψ 5 = 2Re ( U ¯ ð ψ ) . (7)  \nwith  \nΛ ˘ Π = 1 2 (1 -˘ y ) ∂ ˘ y ˘ W∂ ˘ y ˘ ψ + (1 -˘ y ) 2 4 R wt ∂ 2 ˘ y ˘ ψ + ∂ ˘ u R wt R wt (1 -˘ y ) ∂ 2 ˘ y ˘ ψ + 1 2 (1 -˘ y ) ˘ W∂ 2 ˘ y ˘ ψ + 1 2 ˘ W∂ ˘ y ˘ ψ -1 2 Re ( ¯ ð ˘ ψ∂ ˘ y ˘ U + ð ¯ U∂ ˘ y ˘ ψ +2 ˘ U ¯ ð ∂ ˘ y ˘ ψ +2 ð R wt R wt ¯ ˘ U∂ ˘ y ˘ ψ ) . (12)  \nHere we have used two identities:  \nð ∂ r = (1 -˘ y ) 2 2 R wt ( ð ∂ ˘ y + ˘ ð R wt R wt ∂ ˘ y ) , ∂ r = (1 -˘ y ) 2 2 R wt ∂ ˘ y . (13)  \nNotice that in Eqs. (11), (12), and (13) the spin-weighted derivative ð is still evaluated with respect to the original angular system x A . It is related to the numerically adapted version ˘ ð by the following relation [55]:  \nð ˘ ψ = ˘ ð ˘ ψ -(1 -˘ y ) ˘ ð R wt R wt ∂ ˘ y ˘ ψ. (14)  \nTo compute ˘ ð ˘ ψ , we follow [56] and use libsharp routines [72, 73] to obtain a modal decomposition of ˘ ψ in terms of spin-weighted spherical harmonics; subsequently, we apply the differential matrix.  \nIn Eq. (4), the presence of ψ does not alter the left-hand side of the evolution equations, making their hierarchical structure unchanged: the right-hand side of the β equation solely comprises J, ψ and their radial derivatives, the right-hand side of the Q equation comprises solely J, ψ, β and their derivatives, and likewise for the other evolved variables. To evolve the system, we provide initial data for J and ψ on the first u = const null slice, then radially integrate Eqs. (4) and (6) by order to determine β, Q, U, W, H , and finally Π. Following this, we advance J and ψ to the subsequent null slice based on the values of H = ∂ u J and Π = ∂ u ψ .', 'B. SpECTRE CCE system': 'A spectral CCE system has been implemented for GR in SpECTRE [55, 56]. To facilitate spectral implementations, the authors introduced numerically adapted coordinates (˘ u, ˘ y, ˘ x ˘ A ):  \n˘ u = u, ˘ y = 1 -2 R wt r , ˘ θ = θ, ˘ φ = φ, (8)  \nwhere R wt is the Bondi radius of a worldtube for CCE, and it is not to be confused with the Ricci scalar R . On a ˘ u = const null slice, the simulation domain spans radially from the worldtube at ˘ y = -1( r = R wt ) to future null infinity at ˘ y = 1( r = ∞ ). With this new coordinate system, the time (˘ u ) derivative of J , ˘ H = ∂ ˘ u J , is related to the original one H , via [55]  \n˘ H = H +(1 -˘ y ) ∂ ˘ u R wt R wt ∂ ˘ y J. (9)  \nThe second term arises as the Jacobian of the coordinate transformation. Similarly, the time derivative of the scalar field, ˘ Π = ∂ ˘ u ψ , transforms in the same way:  \n˘ Π = Π + (1 -˘ y ) ∂ ˘ u R wt R wt ∂ ˘ y ψ. (10)  \nThe values of other variables are not affected by the coordinate transformation, namely ˘ F = F , for all F ∈ { ψ, J, β, Q, U, W } .  \nInserting Eq. (10) into Eq. (6) yields the hypersurface equation for ˘ Π:  \n(1 -˘ y ) ∂ ˘ y ˘ Π+ ˘ Π = -Re ( ˘ U ¯ ð ˘ ψ ) +(1 -˘ y ) [ Λ ˘ Π + e 2 ˘ β 4 R wt ( ˘ N ψ 1 -˘ N ψ 2 + ˘ N ψ 3 ) ] , (11)  \nR wt  \nThe transformation of the metric sector in Eq. (4) has been extensively discussed in [56]. Here we simply present additional terms contributed by ˘ ψ  \n˘ = ( ˘ ) + 2 (1 ˘ )( ˘ )  \n∂ ˘ y β S ˘ β J π -y ∂ ˘ y ψ 2 , (1 -˘ y ) ∂ ˘ y ˘ Q +2 ˘ Q = S ˘ Q ( ˘ J, ˘ β ) + 16 π (1 -˘ y ) ð ˘ ψ∂ ˘ y ˘ ψ, (1 -˘ y ) ∂ ˘ y ˘ W +2 ˘ W = S ˘ W ( ˘ J, ˘ β, ˘ Q, ˘ U ) + πe 2 ˘ β (1 -˘ y ) R wt ( ˘ J ( ¯ ð ˘ ψ ) 2 + ¯ ˘ J ( ð ˘ ψ ) 2 -2 ˘ K ð ˘ ψ ¯ ð ˘ ψ ) , (1 -˘ y ) ∂ ˘ y ˘ H + ˘ H + L ˘ H ( ˘ J, ˘ β, ˘ Q, ˘ U, ˘ W ) ˘ H + L ¯ ˘ H ( ˘ J, ˘ β, ˘ Q, ˘ U, ˘ W ) ¯ ˘ H = S ˘ H ( ˘ J, ˘ β, ˘ Q, ˘ U, ˘ W ) + 2 π (1 -˘ y ) e 2 ˘ β ( ð ˘ ψ ) 2 . (15)', 'C. Worldtube transformation': "In order to integrate the hypersurface equation for ˘ Π in Eq. (11), we need to supply a boundary condition at a worldtube R wt . Normally it comes from a separate Cauchy evolution of the KG equation, e.g., see [74]. Since the two simulation methods adopt different gauges, a worldtube transformation is needed to construct the value of ˘ Π at the worldtube out of Cauchy variables.  \nIt turns out that the transformation has a concise expression:  \n˘ Π | wt = ∂ t ' ψ +Re ( U (0) ¯ ˘ ð ψ ) . (16)  \nDerivation details can be found in Appendix A. In Eq. (16), ∂ t ' ψ is the Cauchy time derivative of ψ ; ¯ ˘ ð ψ is the angular derivative of ψ [see Eq. (5)] with respect to the numerically adapted coordinates; and U (0) is defined in Eq. (34b) of [55], with a spin weight of 1. Physically speaking, the numerically adapted angular coordinates ˘ x ˘ A on the worldtube evolve over time with respect to the Cauchy coordinates, via [Eq. (28) of [55]]  \n∂ t ' ˘ x ˘ A = -1 2 ( U (0) ¯ ˘ q ˘ A + ¯ U (0) ˘ q ˘ A ) , (17)  \nwhere the quantity U (0) captures the rate of the coordinate motion. Therefore, U (0) corresponds to the shift vector pulled back 1 to the 3D timelike worldtube. On the other hand, the numerically adapted time ˘ u is chosen to flow at the same rate as Cauchy's time t ' - the lapse is manually set to 1 [55]. As a result, the worldtube transformation in Eq. (16) is nothing but the Lie derivative of ψ .  \nAfter setting the boundary condition ˘ Π | wt , we integrate the hypersurface equation in Eq. (11) and evolve the scalar field forward in time. On the worldtube, the scalar field computed by the Cauchy system, ψ Cauchy , should agree with the one obtained from the characteristic system ψ Char . This yields a worltube constraint  \nSince ψ is an evolved variable of the characteristic system, we can directly compute the value of ¯ ˘ ð ψ on the characteristic grid using the libsharp routines, without taking information from the Cauchy side. Meanwhile, the spin-weighted variable U (0) is already available while evolving the GR system. Therefore, in order to perform the worldtube transformation in Eq. (16), we simply need to communicate ∂ t ' ψ from the Cauchy to the characteristic system.  \nC wt ≡ ψ Char -ψ Cauchy . (18)  \nIn actual numerical simulations, the constraint is normally nonvanishing since the two fields are evolved independently; and it should converge to 0 in the continuum limit. Therefore, tracking its evolution provides a diagnosis of our simulations. In fact, C wt offers a specific examination for the worldtube transformation in Eq. (16), since the characteristic evolution of ψ Char is fully controlled by the boundary condition ˘ Π | wt - a wrong transformation would lead to a wrong characteristic evolution of ψ Char , which drives it away from ψ Cauchy .", 'D. Computation of Weyl scalars at future null infinity': "The SpECTRE CCE system [55, 56] uses the Newman-Penrose (NP) formalism to compute Weyl scalars. The relevant equations are 2 [75]  \nΨ 0 = Dσ -δκ -σ ( ρ + ¯ ρ +3 /epsilon1 -¯ /epsilon1 ) + κ ( τ -¯ π NP + ¯ α +3 β NP ) , (19a)  \nΨ 1 = Dτ -∆ κ -( τ + ¯ π NP ) ρ -(¯ τ + π NP ) σ -( /epsilon1 -¯ /epsilon1 ) τ +(3 γ + ¯ γ ) κ -Φ 01 , (19b)  \nΨ 2 = Dµ -δπ NP -(¯ ρµ + σλ ) -π NP (¯ π NP -¯ α + β NP ) + µ ( /epsilon1 +¯ /epsilon1 ) + νκ -2Λ , (19c)  \nΨ 3 = ¯ δµ -δλ + ν ( ρ -¯ ρ ) + π NP ( µ -¯ µ ) + µ ( α + ¯ β NP ) + λ (¯ α -3 β NP ) + Φ 21 . (19d)  \nwhere D = l a ∇ a , ∆ = n a ∇ a , δ = m a ∇ a are the derivative operators used in the NP formalism. The tetrad ( l a , n a , m a ) is chosen to be [55]  \nl a = 1 √ 2 ∂ a r , n a = √ 2 e -2 β ( ∂ a u -V 2 r ∂ a r + 1 2 ¯ Uq a + 1 2 U ¯ q a ) , m a = -1 √ 2 r ( √ K +1 2 q a -J √ 2(1 + K ) ¯ q a ) . (20)  \nCompared with their vacuum GR counterparts, the equations gain additional source terms Φ 01 , Λ, and Φ 21 . They are related to the scalar field ψ through [75]  \nΦ 01 = 1 2 R ml = 4 πδψDψ, Φ 21 = 1 2 R n ¯ m = 4 π ¯ δψ ∆ ψ, Λ = 1 24 R = -1 12 ( R ln -R m ¯ m ) = -2 π 3 [ Dψ ∆ ψ -δψ ¯ δψ ] . (21)  \nThe asymptotic expansion of these quantities at future null infinity I + is of particular interest. Hereafter, we denote the terms in the expansion with the following notation:  \nf = f ( n ) r n + f ( n +1) r n +1 + O ( 1 r n +2 ) , (22)  \nwhere f represents any variable under investigation, and n is an integer. As for the massless scalar field ψ , its value at I + , namely ψ (0) , is a constant and can be set to 0 without changing physical aspects of the system. Therefore, its radial falloff obeys ψ ∼ r -1 . Substituting the asymptotic behavior of ψ into Eq. (20) leads to  \nDψ = -1 √ 2 ψ (1) r 2 + O ( 1 r 3 ) , ∆ ψ = √ 2 e -2 β (0) Π (1) r + O ( 1 r 2 ) , δψ = -1 √ 2 ð ψ (1) r 2 + O ( 1 r 3 ) , (23)  \nwhere we have used the fact that U (0) = J (0) = 0, K (0) = 1, and V ∼ r . The radial falloff rate for Φ 01 , Λ, and Φ 21 can be obtained accordingly:  \nΦ 01 = 2 π ψ (1) ð ψ (1) r 4 + O ( 1 r 5 ) , Λ = 2 π 3 e -2 β (0) Π (1) ψ (1) r 3 + O ( 1 r 4 ) , Φ 21 = -4 πe -2 β (0) Π (1) ¯ ð ψ (1) r 3 + O ( 1 r 4 ) . (24)  \nIn the code, we compute the value of Π (1) , ψ (1) , and ð ψ (1) via  \nΠ (1) = -2 R wt ( ∂ ˘ y ˘ Π+ ∂ ˘ u R wt R wt ∂ ˘ y ˘ ψ )∣ ∣ ∣ ∣ ˘ y =1 , ψ (1) = -2 R wt ∂ ˘ y ˘ ψ, ð ψ (1) = -2 ( ð R wt R wt ∂ ˘ y ˘ ψ + ð ∂ ˘ y ˘ ψ )∣ ∣ ∣ ∣ ˘ y =1 . (25)  \nThe asymptotic expansions of Wely scalars in GR have been shown in Eq. (93) of [55]. With the presence of ψ , their expressions are modified to  \nΨ (4) 1 = Ψ (4)GR 1 -Φ (4) 01 , Ψ (3) 2 = Ψ (3)GR 2 -2Λ (3) , Ψ (2) 3 = Ψ (2)GR 3 . (26)  \nNotice that the expression for Ψ (2) 3 is the same as its GR counterpart, since the additional source term Φ 21 ∼ O ( r -3 ) decays faster than that of Ψ 3 . The computation of Ψ 0 requires more attention. Although the formula in Eq. (19a)  \nis valid geometrically and does not rely on the detail of a gravity theory, the asymptotic expansion in [55] [see their Eq. (93a)] implicitly assumes the vanishing of the stress-energy tensor. Instead, here we directly expand Eq. (19a) and obtain:  \nΨ (5) 0 = β (2) J (1) + ¯ J (1) J (1)2 4 -3 2 J (3) . (27)  \nThe value of β (2) can be obtained by expanding the rr component of Einstein's equations [see Eq. (2.9c) in [69] or Eq. (2.9a) in [71]]  \nβ (2) = -J (1) ¯ J (1) 16 -πψ (1)2 . (28)  \nFinally, the Weyl scalar Ψ (1) 4 and the Bondi-Sachs News take the same forms as in GR, as shown in Eq. (93e) and (42) of [55], respectively.", 'E. Bianchi identities at future null infinity': 'One way to validate our CCE code is to examine Bianchi identities at future null infinity, which establish connections between Weyl scalars Ψ (5) ... (1) 0 ... 4 and the strain h (1) . Their expressions in GR can be found in, e.g. [52]. To extend the identities to our case, we consider the following NP equations [76]:  \n-∆ λ + ¯ δν -λ ( µ + ¯ µ ) -(3 γ -¯ γ ) λ +(3 α + ¯ β NP + π NP -¯ τ ) ν -Ψ 4 = 0 , (29a)  \n- -∆Ψ 3 + δ Ψ 4 +(4 β NP -τ )Ψ 4 -2(2 µ + γ )Ψ 3 +3 ν Ψ 2 -¯ δ Φ 22 +∆Φ 21 +(¯ τ -2 ¯ β NP -2 α )Φ 22 +2( γ + ¯ µ )Φ 21\n- +2 λ Φ 12 -2 ν Φ 11 -¯ ν Φ 20 = 0 ,  \n(29b)  \n- -∆Ψ 2 + δ Ψ 3 +2 ν Ψ 1 -3 µ Ψ 2 +2( β NP -τ )Ψ 3 + σ Ψ 4 -D Φ 22 + δ Φ 21 +2(¯ π NP + β NP )Φ 21 -2 µ Φ 11 -¯ λ Φ 20\n- +2 π NP Φ 12 +(¯ ρ -2 /epsilon1 -2¯ /epsilon1 )Φ 22 -2∆Λ = 0 ,  \n(29c)  \n- -∆Ψ 1 + δ Ψ 2 + ν Ψ 0 +2( γ -µ )Ψ 1 -3 τ Ψ 2 +2 σ Ψ 3 +∆Φ 01 -¯ δ Φ 02 +2(¯ µ -γ )Φ 01 -2 ρ Φ 12 -¯ ν Φ 00 +2 τ Φ 11  \n+(¯ τ -2 ¯ β NP +2 α )Φ 02 +2 δ Λ = 0 , (29d)  \n-  \n∆Ψ  \n0  \n+  \nδ  \nΨ  \n1  \n+(4  \nγ  \n-  \nµ  \n- +(¯ ρ +2 /epsilon1 -2¯ /epsilon1 )Φ 02 = 0 .  \n(29e)  \nwhere the source terms Φ 22 , Φ 00 , Φ 11 , Φ 02 , and Φ 20 are given by  \nΦ 22 = 1 2 R nn = 4 π (∆ ψ ) 2 = 8 π ˙ ψ (1) 2 r 2 + O ( 1 r 3 ) , Φ 00 = 1 2 R ll = 4 π ( Dψ ) 2 = 2 π ψ (1) 2 r 4 + O ( 1 r 5 ) , Φ 11 = 1 4 ( R ln + R m ¯ m ) = 2 π ( Dψ ∆ ψ + δψ ¯ δψ ) = -2 π ψ (1) ˙ ψ (1) r 3 + O ( 1 r 4 ) , ¯ Φ 20 = Φ 02 = 1 2 R mm = 4 π ( δψ ) 2 = 2 π ( ð ψ (1) ) 2 r 4 + O ( 1 r 5 ) . (30)  \nIn Eq. (30) we have used  \nδ = -1 √ 2 r ð + O ( 1 r 2 ) , ∆ = √ 2 ∂ u + O ( 1 r ) , (31)  \nto obtain the asymptotic expansions. The overhead dots denote the time derivative ∂ u . Meanwhile, the Bianchi identities involve NP spin coefficients, whose expressions in terms of Bondi-Sachs variables can be found in Eq. (86) of  \n)Ψ  \n0  \n-  \nτ  \n2(2  \n+  \nβ  \nNP  \n)Ψ  \n1  \n+3  \nσ  \n2  \nΨ  \n-  \nD  \n02  \nΦ  \n+  \nδ  \n01  \nΦ  \nπ  \n+2(¯  \nNP  \n-  \nβ  \nNP  \n)Φ  \n01  \n-  \n2  \nκ  \n12  \nΦ  \n-  \n¯  \nλ  \n00  \nΦ  \n+2  \nσ  \n11  \nΦ  \n[55]. Near future null infinity, their asymptotic behaviors read 3  \nJ = ¯ h (1) r + O ( 1 r 2 ) , σ = ¯ h (1) 2 √ 2 r 2 + O ( 1 r 3 ) , µ = -1 √ 2 r + O ( 1 r 2 ) , β NP = cot θ 2 √ 2 r + O ( 1 r 2 ) , λ (1) = 1 √ 2 ˙ h (1) + O ( 1 r 2 ) , ρ = -1 √ 2 r , κ = 0 , τ ∼ O ( 1 r 2 ) , /epsilon1 ∼ O ( 1 r 2 ) , π NP ∼ O ( 1 r 2 ) , ν ∼ O ( 1 r 2 ) , γ ∼ O ( 1 r 2 ) . (32)  \nPlugging Eqs. (30), (31), and (32) into Eq. (29), we obtain the following Bianchi identities at future null infinity  \nΨ (1) 4 = -h (1) , (33a)  \n˙ Ψ (3) 2 = -1 2 ð Ψ (2) 3 + 1 4 ¯ h (1) Ψ (1) 4 + 8 π 3 ˙ ψ (1) 2 -4 π 3 ψ (1) ¨ ψ (1) , (33c)  \n˙ Ψ (2) 3 = -1 2 ð Ψ (1) 4 , (33b)  \n˙ Ψ (4) 1 = -1 2 ð Ψ (3) 2 + 1 2 ¯ h (1) Ψ (2) 3 -8 π 3 ˙ ψ (1) ð ψ (1) + 4 π 3 ψ (1) ð ˙ ψ (1) , (33d)  \n˙ Ψ (5) 0 = -1 2 ð Ψ (4) 1 + 3 4 ¯ h (1) Ψ (3) 2 +2 π ( ð ψ (1) ) 2 -πψ (1) ðð ψ (1) -π ˙ ¯ h (1) ψ (1) 2 -√ 2 π ¯ h (1) ψ (1) ˙ ψ (1) , (33e)  \nWe will use the PYTHON package scri [79-83] to verify the identities in our following calculations. The package adopts a different convention (hereafter MB) [83, 84] compared to SpECTRE CCE. The conversion factors read [52]  \nΨ [MB] n = ( -√ 2) 2 -n Ψ n , ð [MB] = 1 √ 2 ð , σ (2)MB = √ 2 σ (2) . (34)  \nThen Eqs. (33) become  \nΨ (1)[MB] 4 = -¨ ¯ σ (2)[MB] , (35a)  \n˙ Ψ (2)[MB] 3 = ð [MB] Ψ (1)[MB] 4 , (35b)  \n˙ Ψ (4)[MB] 1 = ð [MB] Ψ (3)[MB] 2 +2 σ (2)[MB] Ψ (2)[MB] 3 + 16 π 3 ˙ ψ (1) ð [MB] ψ (1) -8 π 3 ψ (1) ð [MB] ˙ ψ (1) , (35d)  \n˙ Ψ (3)[MB] 2 = ð [MB] Ψ (2)[MB] 3 + σ (2)[MB] Ψ (1)[MB] 4 + 8 π 3 ˙ ψ (1) 2 -4 π 3 ψ (1) ¨ ψ (1) , (35c)  \n˙ Ψ (5)[MB] 0 = ð [MB] Ψ (4)[MB] 1 +3 σ (2)[MB] Ψ (3)[MB] 2 +8 π ( ð [MB] ψ (1) ) 2 -4 πψ (1) ð [MB] ð [MB] ψ (1) -4 π ˙ σ (2)[MB] ψ (1) 2 -4 πσ (2)[MB] ψ (1) ˙ ψ (1) . (35e)', 'III. SCALAR-INDUCED MEMORY EFFECTS': 'A massless scalar field can propagate to future null infinity [85] and leave nontrivial imprints on memory effects. Here we extend the discussions in [53, 86] and derive the memory contribution from the scalar field. Since most of our discussions in this section adopt only the leading-order asymptotic behavior near future null infinity, we suppress the superscript defined in Eq. (22) for conciseness, unless stated otherwise. To be consistent with previous discussions, we still use the convention of SpECTRE CCE, as opposed to MB [Eq. (34)].  \nThe computation of memory effects can be achieved via the balance laws [69] [see their Eq. (2.11a) and (C6)]  \n1 1  \n˙ ˆ N A = -6 π ˙ ψD A ψ +2 πψ∂ A ˙ ψ +4 πuD A ˙ ψ 2 + u 8 D A ( N BC N BC ) -3 8 N AB D C C BC + 3 8 C AB D C N BC + 1 8 N BC D B C AC -1 8 C BC D B N AC + 1 4 D B D A D C C BC -1 4 D 2 D C C AC -1 4 uD A D B D C N BC , (36b)  \n˙ m = -4 π ˙ ψ 2 -8 N AB N AB + 4 D A D B N AB , (36a)  \nwhere m = -1 2 W (2) is the Bondi mass aspect; ˆ N A is related to the angular momentum aspect N A via  \nˆ N A = N A -uD A m -1 16 D A ( C BC C BC ) -1 4 C AB D C C BC . (37)  \nFinally, C AB is the O ( r -1 ) order of the angular metric h AB defined in Eq. (3), and N AB = ∂ u C AB is the Bondi News tensor.  \nFollowing [69, 86], we decompose C AB into an electric (Φ) and a magnetic (Ψ) potential  \nC AB = ( D A D B -1 2 q AB D 2 ) Φ+ /epsilon1 C ( A D B ) D C Ψ , (38)  \nwhere /epsilon1 CA = i 2 q C ∧ ¯ q A is the volume form compatible with the unit sphere metric.', 'A. Electric memory': 'We insert Eq. (38) into Eq. (36a) and obtain [86]  \nΦ = D -1 [ m + ∫ ( 4 π ˙ ψ 2 + 1 4 ˙ h ˙ ¯ h ) du + α ( θ, φ ) ] . (39)  \nHere α ( θ, φ ) is an arbitrary function on a two-sphere and D ≡ 1 8 D 2 ( D 2 +2) is an angular differential operator. We see that the null (nonlinear) memory obtains an additional contribution from the scalar energy flux ∼ ∫ ˙ ψ 2 du , compared to its GR counterpart [86]. On the other hand, the Bondi mass aspect contributes to the ordinary (linear) memory, and is given by [78]  \nm = -Re ( Ψ 2 + 1 4 ˙ h ¯ h ) -4 π 3 ψ ˙ ψ. (40)  \nNotice that the scalar ordinary memory vanishes in a nonradiative regime ( ˙ ψ ∼ 0).  \nThe potential Φ is related to the electric memory J E via [86]  \nJ E = 1 2 ¯ ð 2 Φ = J ET o + J ET null + J ES o + J ES null , (41)  \nwhere the scalar pieces read  \nJ ES o = 1 2 ¯ ð 2 D -1 ( -4 π 3 ψ ˙ ψ ) , J ES null = 1 2 ¯ ð 2 D -1 (∫ 4 π ˙ ψ 2 du ) . (42)  \nThe tensor contributions J ET o and J ET null are the same as GR, and we provide their expressions below for completeness  \nJ ET o = 1 2 ¯ ð 2 D -1 [ -Re ( Ψ 2 + 1 4 ˙ h ¯ h )] , J ET null = 1 2 ¯ ð 2 D -1 [∫ 1 4 ˙ h ˙ ¯ hdu + α ( θ, φ ) ] . (43)', 'B. Magnetic memory': 'The magnetic memory can be obtained by projecting both sides of Eq. (36b) with /epsilon1 AE D E , and the corresponding expression for Ψ in GR is provided in Eq. (45) of [86]. With the presence of the scalar field, a convenient way to compute Ψ is by replacing ˙ ˆ N = q A ˙ ˆ N A in the GR version with the following combination:  \n˙ ˆ N → ˙ ˆ N +6 π ˙ ψ ð ψ -2 πψ ð ˙ ψ, (44)  \nwhich yields  \nD 2 D Ψ = -Im [ ¯ ð ˙ ˆ N +8 π ¯ ð ˙ ψ ð ψ + 1 8 ð ( 3 h ¯ ð ˙ ¯ h -3 ˙ h ¯ ð ¯ h + ˙ ¯ h ¯ ð h -¯ h ¯ ð ˙ h ) ] . (45)  \nTo obtain ˆ N , we first contract Eq. (37) with q A and get  \nˆ N = N -u ð m -1 8 ð ( h ¯ h ) -1 4 ¯ h ð h. (46)  \nIn GR, the angular momentum aspect N is solely determined by Ψ 1 , e.g., see Appendix B of [86]. As for the Einstein-Klein-Gordon system, an extra term shows up. To see this, we use Eqs. (B2), (B5) and (B6) in [86], namely  \nU (3) = -1 3 ( Q (2) -¯ h ¯ Q (1) ) , U (3) = -2 3 N + 1 8 ð ( h ¯ h ) + 1 2 ¯ h ð h, ¯ Q (1) = ð h, (47)  \ntogether with the asymptotic expansion for Ψ 1 in Eq. (26) and the one for β in Eq. (28), then we arrive at  \nN = 2Ψ 1 -2 πψ ð ψ. (48)  \nPlugging Eqs. (46) and (48) into Eq. (45), and by virtue of the Bianchi identity in Eq. (33d), we find the scalar contribution to Ψ vanishes identically. The corresponding magnetic memory J B = -i 2 ¯ ð 2 Ψ reads [86]  \nJ B = J B o + J B null , (49)  \nwith  \nJ B o = i 2 ¯ ð 2 D -1 D -2 Im { ¯ ð [ -ð Ψ 2 + 1 4 ( ¯ h ð ˙ h -˙ ¯ h ð h ) ]} , J B null = i 2 ¯ ð 2 D -1 D -2 Im { 1 8 ð ( 3 h ¯ ð ˙ ¯ h -3 ˙ h ¯ ð ¯ h + ˙ ¯ h ¯ ð h -¯ h ¯ ð ˙ h ) } . (50)  \nThe expression can be further simplified by replacing the term ¯ ð ( ¯ h ð ˙ h -˙ ¯ h ð h ) in J B o with its negative complex conjugate 4 . After combing with J B null , we obtain  \nJ B = -i 2 ¯ ð 2 D -1 D -2 Im ¯ ðð [ Ψ 2 + 1 4 ˙ h ¯ h ] = i 8 ¯ ð 2 D -1 Im ð 2 h, (51)  \nwhere we have used the fact that the Bondi mass aspect is real-valued, namely [77]  \nIm ( Ψ 2 + 1 4 ð 2 h + 1 4 ˙ h ¯ h ) = 0 . (52)  \nNotice that the right-hand side of Eq. (51) is linear in h , as opposed to the electric memory in Eq. (41) that depends on energy fluxes, the balance law for the angular momentum aspect in Eq. (36b) in fact defines a projection operator for the magnetic component of h . In particular, Eq. (51) implies that  \nJ B lm ∼ h lm -( -1) m ¯ h l, -m . (53)  \nIt corresponds to the radiative current moment, see e.g., Eq. (30) in [87], which plays an important role in building up gravitational recoils [82, 87, 88] and the excitation of quadratic quasinormal modes [89, 90]. Under parity conjugation, J B lm transforms as follows  \nJ B lm → ( -1) l +1 J B lm , (54)  \nwhere we have used Eq. (C10d) in [82]. Therefore, J B lm carries odd parity. In the context of Schwarzschild perturbation theory, the dynamics of J B is described by the Regge-Wheeler equation [91].', 'IV. TESTING THE CHARACTERISTIC EVOLUTION': 'Having outlined CCE for the Einstein-Klein-Gordon system, we now proceed to test our code built in SpECTRE. Given that the metric evolution in GR has been thoroughly examined in [56], our tests primarily focus on the new  \nfeatures introduced by the scalar field. This involves two main parts: solving the KG equation in Eq. (11) and conducting the full CCE procedure for the coupled (scalar+tensor) system. We address these two parts separately in this and the subsequent sections.  \nBelow in this section, we test the implementation of the KG equation by evolving scalar fields on two prescribed spacetime backgrounds (namely without scalar-induced backreaction into the metric sector). We pay particular attention to two aspects: (a) the correct volume integration of the KG equation, and (b) the accurate implementation of the worldtube transformation in Eq. (16). For (a), we compare the computed scalar modes at future null infinity to expected (analytical) behaviors. For (b), we monitor the worldtube constraint defined in Eq. (18), as discussed in Sec. II C.', 'A. Bouncing Schwarzchild BH': "Following [56, 92], we consider a Schwarzchild BH in an oscillating coordinate system - a time-dependent transformation is applied to Kerr-Schild coordinates { t KS , x KS , y KS , z KS } [56, 92]:  \nt ' = t KS , x ' = x KS -a sin 4 ( 2 πt KS b ) , y ' = y KS , z ' = z KS , (55a)  \nwhere a and b are two constants. Since the transformation is a pure gauge effect, one expects no gravitational waves at future null infinity, even though each metric component evolves with time t ' nontrivially. The whole characteristic system must be properly established to cancel out the gauge effect at future null infinity. This makes the bouncing BH system the most demanding test [56].  \nFor our testing purposes, it is preferable to construct a scalar profile that can be controlled analytically, allowing us to compare numerical results with the expected analytical expression. To do this, we first work in the outgoing Eddington-Finkelstein coordinate system { u, r, θ, φ } . The retarded time u is related to the Kerr-Schild time t KS via  \nu = t KS -r -4 M ln ( r 2 M -1 ) , (56)  \nwith M being the mass of the BH. The KG equation □ ψ = 0 in the Eddington-Finkelstein frame reads  \n2 ∂ r ∂ u ( rψ ) = ∂ r [( 1 -2 M r ) ∂ r ( rψ ) ] -2 M r 3 ( rψ ) + ð ¯ ð ( rψ ) r 2 , (57)  \nwhere ð and ¯ ð are derivatives associated with the outgoing Eddington-Finkelstein angular coordinates { θ, φ } . Considering a spherical profile, namely ð ¯ ð ( rψ ) = 0, the solution to the KG equation has the following form  \nrψ = ∞ ∑ n =0 ψ n ( u ) r n . (58)  \nThe coefficients ψ n 's are functions of the retarded time u . Plugging Eq. (58) into (57) yields a three-term recurrence relation  \ndψ 1 du = 0 , 2 dψ n +1 du = 2 M n +1 n 2 ψ n -1 -nψ n , n > 0 . (59)  \nIt has a solution of (assuming M = 1)  \nψ 0 = sin u, ψ 1 = 0 , ψ 2 = -1 2 cos u, ψ 3 = 1 2 sin u, ψ 4 = 3 4 cos u -9 8 sin u, ψ 5 = -77 20 cos u -3 2 sin u, ψ 6 = 15 16 cos u + 51 4 sin u, ψ 7 = 1287 28 cos u -1809 80 sin u, ψ 8 = -12579 80 cos u -19857 128 sin u, ψ 9 = -73557 160 cos u + 133813 140 sin u, ψ 10 = 49797063 8960 cos u + 1272267 1600 sin u, ψ 11 = -116136241 24640 cos u -57286503 1792 sin u, ψ 12 = -16472195091 89600 cos u + 419653067 5120 sin u, ψ 13 = 197057851611 232960 cos u + 517853793843 492800 sin u, ψ 14 = 23027722022071 3942400 cos u -2420206444191 313600 sin u, ψ 15 = -166944468961581 2464000 cos u -218435295225339 7321600 sin u. (60)  \nFigure 1. The convergence of numerical error with angular resolution ( l max ) in the bouncing Schwarzschild BH test. The left panel shows the difference between the computed ψ at future null infinity and the expected value ψ 0 ( u ) = sin u , while the right panel displays the worldtube constraint as defined in Eq. (18).  \n<!-- image -->  \nHere we truncate the solution at n = 15, which is sufficient for our following tests. At future null infinity, the solution is given by  \nrψ | I + = ψ 0 ( u ) = sin u. (61)  \nAs mentioned in Sec. II C, we need to provide a boundary condition for ∂ t ' ψ at a worldtube to perform CCE. This can be obtained by differentiating Eq. (58):  \n∂ t ' ψ = ∞ ∑ n =0 ˙ ψ n ( u ) r n +1 du dt ' -∞ ∑ n =0 ( n +1) ψ n ( u ) r n +2 dr dt ' , (62)  \nwith  \ndu dt ' = 1 -r +2 M r -2 M dr dt ' , dr dt ' = 8 πa br [ x ' + a sin 4 ( 2 πt ' b )] sin 3 ( 2 πt ' b ) cos ( 2 πt ' b ) , (63)  \nwhere we have used Eq. (56) and  \nr 2 = [ x ' + a sin 4 ( 2 πt ' b )] 2 + y ' 2 + z ' 2 . (64)  \nIn our tests, we follow [56] and set the oscillation period b to 40 M and the amplitude a to 2 M . Four worldtube radii, 15 M, 20 M, 25 M, 30 M , are used for CCE. We choose the absolute tolerance of SpECTRE 's time stepper residual to be 10 -12 to ensure that differences between our results and the expected one ψ 0 ( u ) = sin u are dominated by spatial resolution. Figure 1a provides the differences with angular resolution l max ranging from 8 to 22. We can see that they converge exponentially to the level of 10 -11 -10 -10 . Similarly, as shown in Fig. 1b, the worldtube constraint defined in Eq. (18) follows a similar convergence behavior and finally levels off at ∼ 10 -12 .", 'B. Linearized Bondi-Sachs metric': "In the second test, we choose the linearized Bondi-Sachs metric [93] as our spacetime background. The corresponding worldtube data read [92]  \nJ lin lm = √ ( l +2)! ( l -2)! 2 Z lm Re [ J l ( r ) e iνu ] , U lin lm = √ l ( l +1) 1 Z lm Re [ U l ( r ) e iνu ] , (65a)  \nβ lin lm = 0 Z lm Re [ β l ( r ) e iνu ] , W lin lm = 0 Z lm Re [ W l ( r ) e iνu ] . (65b)  \nFigure 2. The linearized Bondi-Sachs test. The left panel shows the amplitude of ψ l =2 ,m =0 (blue dots) and ψ l = m =2 (orange dots) as a function of C 2 a . They depend quadratically and linearly on C 2 a , as expected. The right panel provides the worldtube constraint as defined in Eq. (18).  \n<!-- image -->  \nFollowing [56], we consider a l = m = 2 angular profile, with s Z lm given by 5  \ns Z lm ( θ, φ ) = 1 √ 2 [ s Y lm ( θ, φ ) + ( -1) m s Y l -m ( θ, φ )] . (66)  \nThe expression of J l ( r ) , U l ( r ) , β l ( r ) and W l ( r ) can be found in Eqs. (139) and (140) of [92]. In particular, we have  \nJ l =2 ( r ) = 24 B 2 +3 iνC 2 a -iν 3 C 2 b 36 + C 2 a 4 r -C 2 b 12 r 3 . (67)  \nHere B 2 , C 2 a and C 2 b are three free complex constants. Below we set B 2 to 0 for simplicity, which subsequently yields 3 C 2 a = ν 2 C 2 b due to the asymptotic flatness condition J l =2 ∼ O ( r -1 ). Therefore, the size of the linearized Bondi-Sachs wave is controlled by a single parameter C 2 a . As for the scalar sector, we choose the worldtube data to be  \nψ ( u ) | worldtube = sin u r , ∂ t ' ψ | worldtube = cos u r . (68)  \nWhen the linearized Bondi-Sachs wave vanishes ( C 2 a = 0), the scalar field propagates through flat spacetime. The consequent wavefunction at future null infinity is simply given by rψ ( u ) | I + = sin u . However, for a finite C 2 a , an exact analytical solution is not available. Instead, we can consider a perturbative expansion:  \nrψ ( u, θ, φ ) | I + = sin u + C 2 a ψ (1) ( u, θ, φ ) + C 2 2 a ψ (2) ( u, θ, φ ) + O ( C 3 2 a ) . (69)  \nTo the first order in C 2 a , ψ (1) ( u, θ, φ ) is generated by the coupling between metric components (characterized by the angular profile s Z l = m =2 ) and the leading-order evolution of ψ , namely sin u (characterized by the angular profile 0 Z l = m =0 ). The angular selection rule immediately leads to  \nψ (1) ( u, θ, φ ) ∼ 0 Z 22 ( θ, φ ) . (70)  \nIn other words, the amplitude of the l = m = 2 harmonic of the scalar field is expected to scale linearly with C 2 a . By contrast, other harmonics of the scalar field are generated by nonlinear couplings, making their amplitudes depend at least quadratically on C 2 a .  \nTo perform our tests, we fix the frequency of the linearized Bondi-Sachs wave at ν = 0 . 2, while varying the wave amplitude C 2 a from 10 -6 to 1. To ensure the accuracy of our analysis, we still set the absolute tolerance of the time  \nstepper to 10 -12 . Meanwhile, we find different choices of angular resolution l max lead to only a ∼ 10 -12 difference in final products. Figure 2a shows the amplitude of ψ l =2 ,m =0 (blue dots) and ψ l = m =2 (orange dots) as a function of C 2 a , measured at future null infinity. As expected, they depend quadratically and linearly on C 2 a , respectively. In Fig. 2b we provide the worldtube constraint [Eq. (18)] associated with ( l = 2 , m = 0) and ( l = m = 2). We see that the constraint violation for ψ l =2 ,m =2 is always on the order of 10 -15 , whereas for ψ l =2 ,m =2 it scales quadratically with C 2 a , consistent with the behavior at future null infinity.", 'V. CAUCHY CHARACTERISTIC EVOLUTION': 'After examining the characteristic evolution of the KG equation in Eq. (11), we now proceed to test the full CCE procedure for the Einstein-Klein-Gordon system, which involves both Cauchy and characteristic evolutions. Here we consider a simple setup where a scalar pulse strikes a BH. To perform the simulation, we first construct initial data by solving the extended conformal thin-sandwich equations with the spectral elliptic solver in SpEC [94, 95]. Next, we evolve the system nonlinearly using a first-order generalized harmonic formulation [96]. Our Cauchy evolution for the scalar field follows the method in [74]. Finally, we use the dumped worldtube data for CCE, as outlined in [55, 56] and earlier in this paper. For conciseness, we always set the initial BH mass to unity and use it as our code unit.  \nBelow we will consider two distinct features for tests. (a) It is expected that the scalar emission from a perturbed BH could undergo a tail phase at late times, where the scalar field decays as a power law [62, 63]. In particular, the power law has different exponents at timelike infinity and null infinity [62, 63, 97-99]. The late-time tail at timelike infinity was investigated by Scheel et al. [74]. Below in Sec. V A, we focus on the behavior at null infinity to demonstrate the accuracy of our CCE system. (b) As discussed in Sec. III, the GW (tensor emission) of a BH stirred by a scalar field consists of memory effects. They are governed by the balance laws in Eqs. (41) and (49). In the second part of this section (Sec. V B), we will examine the deviations from the balance laws, as well as the Bianchi identities derived in Sec. II E.', 'A. Late-time tail': "In the first test, we place a Schwarzschild BH at the coordinate center. Initially, the scalar field profile is constructed as follows:  \nψ = A e -( r ' -r ' 0 ) 2 /w 2 r ' Y l = m =1 ( θ ' , φ ' ) , (71)  \nwith A = 0 . 1 , r ' 0 = 50 and w = 5. The angular dependence is set to be the l = m = 1 spherical harmonic. Meanwhile, we choose Π ' = ∂ r ' ψ at t ' = 0, where 6  \nΠ ' ≡ -1 α ' ( ∂ t ' ψ -β i ' ∂ i ' ψ ) , (72)  \nis an auxiliary variable introduced by Scheel et al. [74]. The choice of Π ' represents an infalling scalar pulse into the BH.  \nWe can see that ψ 11 at both timelike and null infinity evolves similarly during the excitation and ringdown phases, whereas their late-time tails exhibit distinct decay rates. Specifically, in Figs. 3b and 3c we fit the late-time portion of ψ 11 (600 < time < 1050, i.e., the gray-shaded regions) to a power law A ( t + t 0 ) -µ , with A,t 0 , and µ being three constants. We consider the following cost function for the fit:  \nAs shown in Fig. 3a, the blue curve corresponds to the l = m = 1 spherical harmonic of the extracted scalar field as a function of the retarded time u , with the CCE worltube placed at a radius of 75. The evolution consists of three stages: the excitation phase (time ≲ 125), the quasinormal ringing phase (125 ≲ time ≲ 250), and the tail phase (time ≳ 250). We note that the scalar field is extracted at future null infinity, where the Bondi radius r approaches ∞ while the retarded time u remains fixed. By comparison, we also obtain the scalar radiation at an approximate timelike infinity by measuring late-Cauchy time ψ at a fixed Cauchy radius, as done by Scheel et al. [74]. The orange curve in Fig. 3a shows the corresponding l = m = 1 scalar harmonic. To make a fair comparison with CCE, our extraction radius is still at 75.  \n∑ t i [log | ψ 11 ( t i ) | -log A + µ log( t i + t 0 )] 2 , (73)  \n(a) Null infinity versus timelike infinity  \n<!-- image -->  \nFigure 4 shows the l = m = 1 harmonic of the scalar field measured at future null infinity. Since the initial data is time-symmetric, it radially contains an ingoing and an outgoing component. The outgoing piece reaches the CCE worldtube and appears at future null infinity at a time of r ' wt -r ' 0 = 80 M ; whereas the ingoing wave first strikes the Kerr BH at r ' 0 = 20 M , the excited quasinormal modes are then transmitted to future null infinity by CCE at  \n<!-- image -->  \n(b) Null infinity  \n<!-- image -->  \n(c) Timelike infinity  \nFigure 3. Resolving ringdown tail in the scalar emission from a Schwarzschild BH. The BH is struck by a scalar field composed of a l = m = 1 spherical harmonic [Eq. (71)]. Figure 3a shows the extracted ψ 11 at null infinity (in blue, computed with CCE) and at an approximate timelike infinity (in orange, measured at a fixed Cauchy radius). In Figs. 3b and 3c, the late-time portion (in gray) of ψ 11 is fitted to a power law. The best fits are plotted as two orange dashed curves. The green curves represent the difference between ψ 11 and the best fits.  \nwhere t i 's are time samples. To improve the fit performance, we notice that for a given nonlinear parameter t 0 , two linear parameters log A and µ can be determined uniquely using ordinary least squares linear regression. This reduces the fit to a one-dimensional minimization problem, and we deal with it using the Nelder-Mead method built in SciPy [100]. We have checked that the results are insensitive to the initial guess. The best fits are plotted as orange dashed lines in Figs. 3b and 3c. We find the tail exponents at null infinity ( µ n ) and timelike infinity ( µ t ) to be 2.93 and 4.92, respectively. They are consistent with the predictions of BH perturbation theory [62, 63, 97-99]: µ n = l +2 and µ t = 2 l +3.", 'B. Memory effects': "In the second test, we consider a Kerr BH with a dimensionless spin χ = 0 . 7, initially surrounded by a time-symmetric scalar field [Π ' = 0, defined in Eq. (72)]. The initial profile of ψ is chosen to be the same as Eq. (71), with A = 1 , r ' 0 = 20 and w = 2. The CCE worldtube is placed at r ' wt = 100 M .  \nFigure 4. The l = m = 1 spherical harmonic of the scalar emission from a Kerr BH, extracted at future null infinity.  \n<!-- image -->  \nr ' wt + r ' 0 = 120 M . Due to the no-hair theorem [26, 38, 101-103], a Kerr BH does not carry a scalar charge. Consequently, the scalar field vanishes when the BH settles to stationary.  \nAs discussed in Sec. II E, the Bianchi identities in Eq. (35) yield relations between the Weyl scalars, the strain, and the scalar field. Meanwhile, Eq. (52) imposes a further constraint between Ψ 2 and h . Figure 6 shows the L 2 -norms of deviations from these constraints, computed at three numerical resolutions. In the first two rows, the convergence of the deviations with increasing resolution validates our CCE procedure. In the bottom row, the constraint violations for Ψ 4 [(35a)] and Ψ 3 [(35b)] do not converge, presumably because they are already small ( ∼ 10 -7 ).  \nDuring the dynamical process, the scalar field stirs spacetime nonlinearly [Eq. (4)] and excites GWs. The blue curves in Figs. 5a and 5b represent the corresponding ( l = 2 , m = 0) and ( l = m = 2) harmonics, respectively. Both curves exhibit a permanent jump after the passage of the scalar wave, with a quasinormal-mode ringing superimposed. As discussed in Sec. III, a strain can be decomposed into an electric and a magnetic memory. In this case, we find that the magnetic sector vanishes. Three dominant electric components [see Eq. (41)] are plotted in orange, green, and red in Figs. 5a and 5b. The scalar-driven null memory J ES null contributes to the overall jump, while the tensor-driven ordinary piece J ET o contributes to the quasinormal-mode ringing. The scalar-driven ordinary memory J ES o only results in a small pulse when the initial outgoing scalar wave reaches future null infinity. Unlike in BBH systems, where the tensor-driven null component J ET null induces a strong memory effect [86], its contribution is negligible in this process. Two lower panels of Figs. 5a and 5b present the difference between the total strain and the sum of all the individual memory contributions. We see that the difference is overall smaller than the numerical error (in yellow) except for the initial junk-radiation regime, thereby justifying our CCE code. To examine the magnetic memory, we consider the ( l = 3 , m = 2) harmonic in Fig. 5c. For the current system, it is solely contributed by the ordinary memory associated with the angular momentum aspect. The difference between J B o and the total strain is again smaller than the numerical error. Finally, Figure 5d shows the ( l = 3 , m = 2) harmonic of the spin memory [104], as a time integration of h 32 .", 'VI. CONCLUSION': "In this paper, we have implemented a CCE algorithm for the Einstein-Klein-Gordon system. Compared to GR, we presented the additional terms contributed by the scalar field in the equations of motion (Sec. II A), the computation of Weyl scalars (Sec. II D), the Bianchi identities (Sec. II E), as well as the memory effects (Sec. III). In addition, we reformulated the characteristic form of the KG equation using the numerically adapted coordinates [Eq. (8)], facilitating its implementation in our numerical relativity code SpECTRE . We also derived a concise worldtube transformation to construct the boundary condition for the characteristic KG equation from a generic Cauchy evolution (Sec. II C).  \nThe tests demonstrate that our CCE code can effectively reveal scalar-induced memory effects. In the future, a direct application of this work is to extract GWs emitted by binary mergers in alternative theories of gravity. For  \nTo evaluate the accuracy of our code, we designed various test systems. We first focused on the implementation of the KG equation and the worldtube transformation by evolving scalar fields on two prescribed spacetime backgrounds. As expected, we found that the numerical error decreases exponentially with increasing simulation resolution. We then examined the full CCE procedure by striking a BH with a scalar pulse. During the late-time evolution, we observed that the scalar field decays as a power law u -l -2 , consistent with Price's law at future null infinity. Furthermore, we verified that the tensor emissions extracted with CCE are consistent with the balance laws and the Bianchi identities.  \n1e 1  \n1e 1  \nFigure 5. Memory effects in the GWs emitted by a Kerr BH, perturbed by a scalar field. Figures 5a and 5b: the blue curves show the ( l = 2 , m = 0) and ( l = m = 2) harmonics of the full strain. They are compared to the three dominant electric memory components (in orange, green, and red), as evaluated based on Eqs. (42) and (43). Two lower panels display the difference between the total strain and the sum of all electric memory contributions (in black), which is further compared to numerical error (in yellow). Figures 5c and 5d: the blue curves (overlapped with the orange curves) illustrate the ( l = 3 , m = 2) harmonic of the full strain and its time integration (spin memory). They are solely contributed by the magnetic ordinary memory J B o (in green), as evaluated based on Eq. (50).  \n<!-- image -->  \nFigure 6. The L 2 -norms of the constraint violations given in Eq. (52) ( C Im ψ 2 , top left) and Eq. (35), computed at three numerical resolutions.  \n<!-- image -->  \nexample, a recent study simulated a black hole - neutron star merger in scalar-tensor theory [42], where the neutron star underwent spontaneous scalarization. Working within the Einstein frame, the current CCE algorithm can faithfully compute memory effects in both tensor and breathing modes [105-107], allowing us to make comparisons with the post-Newtonian approximation [108-110]. These numerical studies will improve our understanding of asymptotic symmetries [71, 111-113]. In addition, our implementation lays a foundation for performing CCE simulations in other modified theories of gravity, e.g., [35, 50], paving the way for studying the associated memory effects [107].  \nOur simulations show tails in scalar fields. A possible avenue for future work is to look for tails in breathing modes (e.g. [42]) and tensor modes [114-116] emitted by binary systems. The accurate simulation of these tails is sensitive to the choice of boundary conditions [117]. Incorrect boundary conditions can introduce nonoscillatory numerical artifacts during the ringdown phase, which might be mistaken for physical tails. This issue can be avoided by either positioning the outer boundaries far enough away to remain causally disconnected from the system or by adopting CCM [57].  \nFinally, it would also be interesting to apply our CCE algorithm to extract scalar emissions from intermediatemass-ratio inspirals. Recently, a worldtube excision method [118, 119] was developed to simulate a small scalar charge orbiting around a BH. By matching a perturbative description around the scalar charge to a Cauchy simulation, the method enables the evolution of a binary system in the intermediate-mass-ratio regime. The accuracy of this method was evaluated by comparing scalar modes extracted at a finite radius to perturbative calculations. Our CCE algorithm could further improve this comparison. Additionally, the algorithm allows for verifying balance laws at future null infinity.", 'ACKNOWLEDGMENTS': 'S.M. would like to thank Vijay Varma and Leo C. Stein for useful discussion. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Colleges and Universities. This material is based upon work supported by the National Science Foundation under Grants No. PHY-2407742, No. PHY- 2207342, and No. OAC-2209655 at Cornell. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. This work was supported by the Sherman Fairchild Foundation at Cornell.', 'Appendix A: Worldtube transformation': "As discussed in [55], in SpECTRE, five coordinate systems are involved in worldutbe transformations, including:  \n- · Cauchy coordinates r ' , x ' A ' , t ' .\n- · Null-radius coordinates λ, x A , u\n- · Bondi-like coordinates r, x A , u .\n- · Partially flat Bondi-like coordinates ˆ r, ˆ x ˆ A , ˆ u .\n- · Numerically adapted coordinates ˘ y, ˘ x ˘ A , ˘ u  \nThe numerically adapted coordinates are used for characteristic evolution, namely our target system. The Jacobians between these coordinate systems are summarized in Ref. [57]. Here we shall not go into full details. Instead, we list only necessary components in the following discussions.", '1. Obtaining ˘ Π | wt : from numerically adapted to partially flat Bondi-like coordinates': 'To obtain the boundary condition ˘ Π | wt , we need the following Jacobian [see Eq. (46b) of [55]]  \n∂ ˆ u = ∂ ˘ u -(1 -˘ y ) ∂ ˘ u ˆ R wt ˆ R wt ∂ ˘ y , (A1)  \nwhere ˆ R wt is the radius of the worldtube in the partially flat Bondi-like coordinate system. Applying the Jacobian to ψ on the worldtube ˘ y = -1, we obtain  \n˘ Π | wt = ( ∂ ˘ u ψ ) ˘ y = ( ∂ ˆ u ψ ) ˆ r + ∂ ˘ u ˆ R wt ˆ R wt ∂ ˘ y ψ. (A2)  \nThe second term ∂ ˘ u ˆ R wt ˆ R wt ∂ ˘ y ψ is already available in the current characteristic system. Therefore, we are left to obtain the value of ( ∂ ˆ u ψ ) ˆ r from the Cauchy system. Below we will transform ( ∂ ˆ u ψ ) ˆ r to the rest of coordinates.', '2. Obtaining ( ∂ ˆ u ψ ) ˆ r : from partially flat Bondi-like to Bondi-like coordinates': "To obtain the value of ( ∂ ˆ u ψ ) ˆ r from the Bondi-like coordinates, the necessary Jacobian is  \n( ∂ u ) r = ( ∂ ˆ u ) ˆ r + ( ∂ u ˆ x ˆ A ) r ∂ ˆ A + ˆ r ∂ u ˆ ω ˆ ω ∂ ˆ r = ( ∂ ˆ u ) ˆ r -Re U (0) ¯ ˆ ð +(1 -˘ y ) ∂ u ˆ ω ˆ ω ∂ ˘ y , (A3)  \nwhere the first line comes from Eq. (4.20) of [57]. The second line uses the evolution equation of ˆ x ˆ A :  \n( ∂ u ˆ x ˆ A ) r = -1 2 ( U (0) ¯ ˆ q ˆ A + ¯ U (0) ˆ q ˆ A ) , (A4)  \nsee Eq. (17) with t ' replaced by u and ˘ x ˘ A by ˆ x ˆ A ; as well as an identity [Eq. (46a) of [55]]  \nˆ r∂ ˆ r = (1 -˘ y ) ∂ ˘ y . (A5)  \nFinally, ˆ ω in Eq. (A3) is the conformal factor of the angular transformation ˆ x ˆ A ( x A ):  \nˆ q ˆ A ˆ B ∂ A ˆ x ˆ A ∂ B ˆ x ˆ B = 1 ˆ ω 2 q AB , (A6)  \nwhere q AB (ˆ q ˆ A ˆ B ) is the unit sphere metric of the (partially flat) Bondi-like coordinates. Applying the Jacobian in Eq. (A3) to ψ and restricting it to the worldtube, we obtain  \n( ∂ u ψ ) r = ( ∂ ˆ u ψ ) ˆ r -Re U (0) ¯ ˆ ð ψ +2 ∂ u ˆ ω ˆ ω ∂ ˘ y ψ. (A7)  \nIn practice, our CCE system does not directly provide the value of ∂ u ˆ ω . Instead, it gives us ∂ ˆ u ˆ ω . The two are related by applying the Jacobian in Eq. (A3) to ˆ ω , which yields  \n( ∂ u ˆ ω ) r = ( ∂ ˆ u ˆ ω ) ˆ r -Re U (0) ¯ ˆ ð ˆ ω. (A8)  \nNote that we have used the fact that the conformal factor ˆ ω does not depend on ˘ y , i.e., ∂ ˘ y ˆ ω = 0. Plugging Eq. (A8) into Eq. (A7), we get  \n( ∂ ˆ u ψ ) ˆ r = ( ∂ u ψ ) r +Re U (0) ¯ ˆ ð ψ -2 ( ∂ ˆ u ˆ ω ) ˆ r -Re U (0) ¯ ˆ ð ˆ ω ˆ ω ∂ ˘ y ψ. (A9)  \nThe only remaining unknown variable on the right-hand side is ( ∂ u ψ ) r . Therefore, our next task is to compute its value from Cauchy variables.", '3. Obtaining ( ∂ u ψ ) r : from Bondi-like to Cauchy coordinates': "To compute ( ∂ u ψ ) r , we need the following Jacobians  \n( ∂ u ) λ = ( ∂ u ) r +( ∂ u r ) ∂ r , ( ∂ t ' ) r ' = ( ∂ u ) λ , (A10)  \nsee Eqs. (4.6) and (4.3) of [57]. Restricting them to the worldtube, we obtain  \n( ∂ u ψ ) r = ( ∂ t ' ψ ) r ' -2 ∂ u R wt R wt ∂ ˘ y ψ, (A11)  \nwhere R wt is the radius of the worldtube in the Bondi-like coordinate system. It is related to ˆ R wt via  \nR wt = ˆ R wt ˆ ω . (A12)  \nThe first term ( ∂ t ' ψ ) r ' in Eq. (A11) is a dynamical variable of the Cauchy system, and the second term ∂ u R wt R wt ∂ ˘ y ψ is already available in the CCE system. Therefore, once we have collected the value of ( ∂ t ' ψ ) r ' from a Cauchy evolution, we can retrace our derivations presented in Secs. A 1, A 2 and A 3 to construct the boundary condition ˘ Π | wt .", '4. Simplication': "Combining Eqs. (A2), (A9) and (A11), we obtain a lengthy expression for ˘ Π | wt :  \n˘ Π | wt = ( ∂ t ' ψ ) r ' -2 ∂ u R wt R wt ∂ ˘ y ψ +Re U (0) ¯ ˆ ð ψ -2 ( ∂ ˆ u ˆ ω ) ˆ r -Re U (0) ¯ ˆ ð ˆ ω ˆ ω ∂ ˘ y ψ + ∂ ˘ u ˆ R wt ˆ R wt ∂ ˘ y ψ. (A13)  \nA pivotal step in simplifying it is to leverage an identity, see Eq. (15b) of [56] 7 :  \n∂ ˆ u ˆ R wt ˆ R wt = ∂ u R wt R wt + ∂ ˆ u ˆ ω ˆ ω +Re U (0) ˆ ð R wt R wt . (A14)  \nwhich leads to  \n˘ Π | wt = ( ∂ t ' ψ ) r ' +Re [ U (0) ( ¯ ˆ ð ψ +2 ∂ ˘ y ψ ¯ ˆ ð ˆ R wt ˆ R wt )] . (A15)  \nThe expression becomes even more simplified when we transform the covariant angular derivative ˆ ð from the partially flat Bondi-like coordinates to the numerically adapted ones ˘ ð , via [Eq. (46c) of [55]]:  \n¯ ˘ ð ψ = ¯ ˆ ð ψ +2 ∂ ˘ y ψ ( ¯ ˆ ð ˆ R wt ˆ R wt ) , (A16)  \nwhere we have used the fact ˘ y = -1 and ¯ ˆ ð ˆ R wt = ¯ ˘ ð ˆ R wt . Plugging Eq. (A16) into (A15) leads to our final result given in Eq. (16).", 'REFERENCES': "- [1] J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. 32 , 074001 (2015), arXiv:1411.4547 [gr-qc].\n- [2] F. Acernese et al. (Virgo), Class. Quant. Grav. 32 , 024001 (2015), arXiv:1408.3978 [gr-qc].\n- [3] K. Somiya (KAGRA), Class. Quant. Grav. 29 , 124007 (2012), arXiv:1111.7185 [gr-qc].\n- [4] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X6 , 041015 (2016), [erratum: Phys. Rev.X8,no.3,039903(2018)], arXiv:1606.04856 [gr-qc].\n- [5] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X9 , 031040 (2019), arXiv:1811.12907 [astro-ph.HE].\n- [6] R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X 11 , 021053 (2021), arXiv:2010.14527 [gr-qc].\n- [7] R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA), (2021), arXiv:2111.03606 [gr-qc].\n- [8] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116 , 221101 (2016), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc].\n- [9] C. M. Will, Living Rev. Rel. 17 , 4 (2014), arXiv:1403.7377 [gr-qc].\n- [10] N. Yunes and X. Siemens, Living Rev. Rel. 16 , 9 (2013), arXiv:1304.3473 [gr-qc].\n- [11] N. Yunes, K. Yagi, and F. Pretorius, Phys. Rev. D 94 , 084002 (2016), arXiv:1603.08955 [gr-qc].\n- [12] E. Berti et al. , Class. Quant. Grav. 32 , 243001 (2015), arXiv:1501.07274 [gr-qc].\n- [13] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 123 , 011102 (2019), arXiv:1811.00364 [gr-qc].\n- [14] A. Krolak, K. D. Kokkotas, and G. Schaefer, Phys. Rev. D 52 , 2089 (1995), arXiv:gr-qc/9503013.\n- [15] K. Yagi and T. Tanaka, Prog. Theor. Phys. 123\n- , 1069 (2010), arXiv:0908.3283 [gr-qc].\n- [16] S. Ma and N. Yunes, Phys. Rev. D 100 , 124032 (2019), arXiv:1908.07089 [gr-qc].\n- [17] Z. Carson, B. C. Seymour, and K. Yagi, Class. Quant. Grav. 37 , 065008 (2020), arXiv:1907.03897 [gr-qc].\n- [18] L. Sampson, N. Yunes, N. Cornish, M. Ponce, E. Barausse, A. Klein, C. Palenzuela, and L. Lehner, Phys. Rev. D 90 , 124091 (2014), arXiv:1407.7038 [gr-qc].\n- [19] P. D. Scharre and C. M. Will, Phys. Rev. D 65 , 042002 (2002), arXiv:gr-qc/0109044.\n- [20] C. M. Will and N. Yunes, Class. Quant. Grav. 21 , 4367 (2004), arXiv:gr-qc/0403100.\n- [21] E. Berti, A. Buonanno, and C. M. Will, Class. Quant. Grav. 22 , S943 (2005), arXiv:gr-qc/0504017.\n- [22] E. Berti, A. Buonanno, and C. M. Will, Phys. Rev. D 71 , 084025 (2005), arXiv:gr-qc/0411129.\n- [23] K. Yagi and T. Tanaka, Phys. Rev. D 81 , 064008 (2010), [Erratum: Phys.Rev.D 81, 109902 (2010)], arXiv:0906.4269 [gr-qc].  \n- [24] K. G. Arun, Class. Quant. Grav. 29 , 075011 (2012), arXiv:1202.5911 [gr-qc].\n- [25] V. Cardoso, S. Chakrabarti, P. Pani, E. Berti, and L. Gualtieri, Phys. Rev. Lett. 107 , 241101 (2011), arXiv:1109.6021 [gr-qc].\n- [26] N. Yunes, P. Pani, and V. Cardoso, Phys. Rev. D 85 , 102003 (2012), arXiv:1112.3351 [gr-qc].\n- [27] E. Berti, L. Gualtieri, M. Horbatsch, and J. Alsing, Phys. Rev. D 85 , 122005 (2012), arXiv:1204.4340 [gr-qc].\n- [28] S. Tuna, K. I. Ünlütürk, and F. M. Ramazanoğlu, Phys. Rev. D 105 , 124070 (2022), arXiv:2204.02138 [gr-qc].\n- [29] F. Pretorius, Phys. Rev. Lett. 95 , 121101 (2005), arXiv:gr-qc/0507014 [gr-qc].\n- [30] J. L. Ripley, Int. J. Mod. Phys. D 31 , 2230017 (2022), arXiv:2207.13074 [gr-qc].\n- [31] W. E. East and J. L. Ripley, Phys. Rev. D 103 , 044040 (2021), arXiv:2011.03547 [gr-qc].\n- [32] W. E. East and J. L. Ripley, Phys. Rev. Lett. 127 , 101102 (2021), arXiv:2105.08571 [gr-qc].\n- [33] W. E. East and F. Pretorius, Phys. Rev. D 106 , 104055 (2022), arXiv:2208.09488 [gr-qc].\n- [34] M. Corman, J. L. Ripley, and W. E. East, Phys. Rev. D 107 , 024014 (2023), arXiv:2210.09235 [gr-qc].\n- [35] M. Corman and W. E. East, (2024), arXiv:2405.18496 [gr-qc].\n- [36] A. D. Kovács and H. S. Reall, Phys. Rev. Lett. 124 , 221101 (2020), arXiv:2003.04327 [gr-qc].\n- [37] A. D. Kovács and H. S. Reall, Phys. Rev. D 101 , 124003 (2020), arXiv:2003.08398 [gr-qc].\n- [38] J. Healy, T. Bode, R. Haas, E. Pazos, P. Laguna, D. Shoemaker, and N. Yunes, Class. Quant. Grav. 29 , 232002 (2012), arXiv:1112.3928 [gr-qc].\n- [39] E. Barausse, C. Palenzuela, M. Ponce, and L. Lehner, Phys. Rev. D 87 , 081506 (2013), arXiv:1212.5053 [gr-qc].\n- [40] M. Shibata, K. Taniguchi, H. Okawa, and A. Buonanno, Phys. Rev. D 89 , 084005 (2014), arXiv:1310.0627 [gr-qc].\n- [41] M. Bezares, R. Aguilera-Miret, L. ter Haar, M. Crisostomi, C. Palenzuela, and E. Barausse, Phys. Rev. Lett. 128 , 091103 (2022), arXiv:2107.05648 [gr-qc].\n- [42] S. Ma, V. Varma, L. C. Stein, F. Foucart, M. D. Duez, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D 107 , 124051 (2023), arXiv:2304.11836 [gr-qc].\n- [43] M. Okounkova, M. Isi, K. Chatziioannou, and W. M. Farr, Phys. Rev. D 107 , 024046 (2023), arXiv:2208.02805 [gr-qc].\n- [44] M. Okounkova, L. C. Stein, M. A. Scheel, and D. A. Hemberger, Phys. Rev. D 96 , 044020 (2017), arXiv:1705.07924 [gr-qc].\n- [45] M. Okounkova, L. C. Stein, J. Moxon, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 101 , 104016 (2020), arXiv:1911.02588 [gr-qc].\n- [46] H. Witek, L. Gualtieri, P. Pani, and T. P. Sotiriou, Phys. Rev. D 99 , 064035 (2019), arXiv:1810.05177 [gr-qc].\n- [47] M. Okounkova, Phys. Rev. D 102 , 084046 (2020), arXiv:2001.03571 [gr-qc].\n- [48] J. Cayuso, N. Ortiz, and L. Lehner, Phys. Rev. D 96 , 084043 (2017), arXiv:1706.07421 [gr-qc].\n- [49] N. Franchini, M. Bezares, E. Barausse, and L. Lehner, Phys. Rev. D 106 , 064061 (2022), arXiv:2206.00014 [gr-qc].\n- [50] G. Lara, H. P. Pfeiffer, N. A. Wittek, N. L. Vu, K. C. Nelli, A. Carpenter, G. Lovelace, M. A. Scheel, and W. Throwe, Phys. Rev. D 110 , 024033 (2024), arXiv:2403.08705 [gr-qc].\n- [51] M. Corman, L. Lehner, W. E. East, and G. Dideron, (2024), arXiv:2405.15581 [gr-qc].\n- [52] D. A. B. Iozzo, M. Boyle, N. Deppe, J. Moxon, M. A. Scheel, L. E. Kidder, H. P. Pfeiffer, and S. A. Teukolsky, Phys. Rev. D 103 , 024039 (2021), arXiv:2010.15200 [gr-qc].\n- [53] K. Mitman et al. , (2024), arXiv:2405.08868 [gr-qc].\n- [54] J. Winicour, Living Rev. Rel. 12 , 3 (2009), arXiv:0810.1903 [gr-qc].\n- [55] J. Moxon, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 102 , 044052 (2020), arXiv:2007.01339 [gr-qc].\n- [56] J. Moxon, M. A. Scheel, S. A. Teukolsky, N. Deppe, N. Fischer, F. Hébert, L. E. Kidder, and W. Throwe, Phys. Rev. D 107 , 064013 (2023), arXiv:2110.08635 [gr-qc].\n- [57] S. Ma et al. , Phys. Rev. D 109 , 124027 (2024), arXiv:2308.10361 [gr-qc].\n- [58] N. Deppe, W. Throwe, L. E. Kidder, N. L. Vu, K. C. Nelli, et al. , 'Spectre,' 10.5281/zenodo.10967177 (2024).\n- [59] G. W. Horndeski, Int. J. Theor. Phys. 10 , 363 (1974).\n- [60] R. Jackiw and S. Y. Pi, Phys. Rev. D 68 , 104012 (2003), arXiv:gr-qc/0308071.\n- [61] S. Alexander and N. Yunes, Phys. Rept. 480 , 1 (2009), arXiv:0907.2562 [hep-th].\n- [62] R. H. Price, Phys. Rev. D 5 , 2419 (1972).\n- [63] R. H. Price, Phys. Rev. D 5 , 2439 (1972).\n- [64] T. Damour and G. Esposito-Farese, Class. Quant. Grav. 9 , 2093 (1992).\n- [65] P. G. Bergmann, International Journal of Theoretical Physics 1 , 25 (1968).\n- [66] R. V. Wagoner, Phys. Rev. D 1 , 3209 (1970).\n- [67] H. Bondi, M. G. J. van der Burg, and A. W. K. Metzner, Proceedings of the Royal Society of London Series A 269 , 21 (1962).\n- [68] G. Barnich and C. Troessaert, Phys. Rev. Lett. 105 , 111103 (2010), arXiv:0909.2617 [gr-qc].\n- [69] E. E. Flanagan and D. A. Nichols, Phys. Rev. D 95 , 044002 (2017), [Erratum: Phys.Rev.D 108, 069902 (2023)], arXiv:1510.03386 [hep-th].\n- [70] W. Barreto, A. Da Silva, R. Gomez, L. Lehner, L. Rosales, and J. Winicour, Phys. Rev. D 71 , 064028 (2005), arXiv:grqc/0412066.\n- [71] S. Tahura, D. A. Nichols, A. Saffer, L. C. Stein, and K. Yagi, Phys. Rev. D 103 , 104026 (2021), arXiv:2007.13799 [gr-qc].\n- [72] M. Reinecke and D. S. Seljebotn, Astronomy & Astrophysics 554 , A112 (2013).\n- [73] 'libsharp,' https://github.com/Libsharp/libsharp , accessed: 2021-08-25.\n- [74] M. A. Scheel, A. L. Erickcek, L. M. Burko, L. E. Kidder, H. P. Pfeiffer, and S. A. Teukolsky, Phys. Rev. D 69 , 104006 (2004), arXiv:gr-qc/0305027.  \n- [75] S. Chandrasekhar, The Mathematical Theory of Black Holes (Clarendon Press, 1998).\n- [76] H. Stephani, D. Kramer, M. MacCallum, C. Hoenselaers, and E. Herlt, Exact solutions of Einstein's field equations (Cambridge university press, 2009).\n- [77] T. M. Adamo, C. N. Kozameh, and E. T. Newman, Living Rev. Rel. 12 , 6 (2009), arXiv:0906.2155 [gr-qc].\n- [78] M. Scholtz and L. Holka, Gen. Rel. Grav. 46 , 1665 (2014), arXiv:1312.7084 [gr-qc].\n- [79] M. Boyle, 'Scri,' https://github.com/moble/scri .\n- [80] M. Boyle, D. Iozzo, and L. C. Stein, 'moble/scri: v1.2,' (2020).\n- [81] M. Boyle, Phys. Rev. D87 , 104006 (2013), arXiv:1302.2919 [gr-qc].\n- [82] M. Boyle, L. E. Kidder, S. Ossokine, and H. P. Pfeiffer, (2014), arXiv:1409.4431 [gr-qc].\n- [83] M. Boyle, Phys. Rev. D93 , 084031 (2016), arXiv:1509.00862 [gr-qc].\n- [84] L. Lehner and O. M. Moreschi, Phys. Rev. D 76 , 124040 (2007), arXiv:0706.1319 [gr-qc].\n- [85] A. D. Helfer, Journal of mathematical physics 34 , 3478 (1993).\n- [86] K. Mitman, J. Moxon, M. A. Scheel, S. A. Teukolsky, M. Boyle, N. Deppe, L. E. Kidder, and W. Throwe, Phys. Rev. D 102 , 104007 (2020), arXiv:2007.11562 [gr-qc].\n- [87] J. D. Schnittman, A. Buonanno, J. R. van Meter, J. G. Baker, W. D. Boggs, J. Centrella, B. J. Kelly, and S. T. McWilliams, Phys. Rev. D77 , 044031 (2008), arXiv:0707.0301 [gr-qc].\n- [88] S. Ma, M. Giesler, V. Varma, M. A. Scheel, and Y. Chen, Phys. Rev. D 104 , 084003 (2021), arXiv:2107.04890 [gr-qc].\n- [89] S. Ma and H. Yang, Phys. Rev. D 109 , 104070 (2024), arXiv:2401.15516 [gr-qc].\n- [90] P. Bourg, R. Panosso Macedo, A. Spiers, B. Leather, B. Bonga, and A. Pound, (2024), arXiv:2405.10270 [gr-qc].\n- [91] T. W. Baumgarte and S. L. Shapiro, Numerical Relativity: Solving Einstein's Equations on the Computer (Cambridge University Press, 2010).  \n[92] K. Barkett, J. Moxon, M. A. Scheel, and B. Szilágyi, Phys. Rev. D  \n102  \n, 024004 (2020), arXiv:1910.09677 [gr-qc].  \n- [93] N. T. Bishop, Class. Quant. Grav. 22 , 2393 (2005), arXiv:gr-qc/0412006.\n- [94] H. P. Pfeiffer, L. E. Kidder, M. A. Scheel, and S. A. Teukolsky, Comput. Phys. Commun. 152 , 253 (2003), arXiv:grqc/0202096 [gr-qc].\n- [95] 'The Spectral Einstein Code,' http://www.black-holes.org/SpEC.html .\n- [96] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, and O. Rinne, Class. Quant. Grav. 23 , S447 (2006), arXiv:gr-qc/0512093 [gr-qc].\n- [97] E. W. Leaver, Phys. Rev. D 34 , 384 (1986).\n- [98] C. Gundlach, R. H. Price, and J. Pullin, Phys. Rev. D 49 , 883 (1994).\n- [99] C. Gundlach, R. H. Price, and J. Pullin, Phys. Rev. D 49 , 890 (1994).\n- [100] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17 , 261 (2020).\n- [101] K. S. Thorne and J. J. Dykla, Astrophys. J. Lett. 166 , L35 (1971).\n- [102] T. P. Sotiriou and V. Faraoni, Phys. Rev. Lett. 108 , 081103 (2012), arXiv:1109.6324 [gr-qc].\n- [103] S. W. Hawking, Commun. Math. Phys. 25 , 167 (1972).\n- [104] S. Pasterski, A. Strominger, and A. Zhiboedov, JHEP 12 , 053 (2016), arXiv:1502.06120 [hep-th].\n- [105] S. M. Du and A. Nishizawa, Phys. Rev. D 94 , 104063 (2016), arXiv:1609.09825 [gr-qc].\n- [106] K. Koyama, Phys. Rev. D 102 , 021502 (2020), arXiv:2006.15914 [gr-qc].\n- [107] L. Heisenberg, N. Yunes, and J. Zosso, Phys. Rev. D 108 , 024010 (2023), arXiv:2303.02021 [gr-qc].\n- [108] R. N. Lang, Phys. Rev. D 89 , 084014 (2014), arXiv:1310.3320 [gr-qc].\n- [109] R. N. Lang, Phys. Rev. D 91 , 084027 (2015), arXiv:1411.3073 [gr-qc].\n- [110] S. Tahura, D. A. Nichols, and K. Yagi, Phys. Rev. D 104 , 104010 (2021), arXiv:2107.02208 [gr-qc].\n- [111] S. Hou and Z.-H. Zhu, JHEP 01 , 083 (2021), arXiv:2005.01310 [gr-qc].\n- [112] S. Hou and Z.-H. Zhu, Chin. Phys. C 45 , 023122 (2021), arXiv:2008.05154 [gr-qc].\n- [113] S. Hou, Astron. Nachr. 342 , 96 (2021), arXiv:2011.02087 [gr-qc].\n- [114] M. De Amicis, S. Albanesi, and G. Carullo, (2024), arXiv:2406.17018 [gr-qc].\n- [115] T. Islam, G. Faggioli, G. Khanna, S. E. Field, M. van de Meent, and A. Buonanno, (2024), arXiv:2407.04682 [gr-qc].\n- [116] G. Carullo and M. De Amicis, (2023), arXiv:2310.12968 [gr-qc].\n- [117] E. W. Allen, E. Buckmiller, L. M. Burko, and R. H. Price, Phys. Rev. D 70 , 044038 (2004), arXiv:gr-qc/0401092.\n- [118] N. A. Wittek et al. , Phys. Rev. D 108 , 024041 (2023), arXiv:2304.05329 [gr-qc].\n- [119] N. A. Wittek, A. Pound, H. P. Pfeiffer, and L. Barack, (2024), arXiv:2403.08864 [gr-qc]."}