Introduction

The LIGO, Virgo, and KAGRA collaboration [1][2][3] has detected about 100 gravitational wave (GW) events [4][5][6][7], providing us with unprecedented opportunities to explore the strong gravitational elds surrounding compact binary systems. This marks a new era in testing general relativity (GR) [4,.

To 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 rst numerical simulation of a binary black hole (BBH) merger in GR [29], evolving modi ed theories of gravity remains challenging, as reviewed in [30]. The primary dif culty lies in the potential modi cation of the principal part of the Einstein equations, rendering the continuum equations unable to admit a well-posed initial value problem. Several approaches have been proposed to address this issue. The rst approach involves directly formulating well-posed schemes for evolving the full set of evolution equations in certain speci c theories, such as weakly-coupled Horndeski gravity theories (e.g. Einsteinscalar-Gauss-Bonnet) [31][32][33][34][35] based on the modi ed generalized harmonic formulation [36,37], as well as Damour-Esposito-Farèse scalar-tensor theories [38][39][40][41][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][44][45] and scalar Gauss-Bonnet gravity [46,47]. The nal treatment introduces auxiliary variables to x high-energy (short length/timescale) degrees of freedom, while leaving low-energy parts unchanged, see [48][49][50]. Recently, Corman et al [51] compared the results of these three approaches under various physical setups.

The 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 in nitely large space, the spatial Cauchy domain is typically truncated at a nite distance from the source, which prevents direct access to GWs at future null in nity. One approach to address this limitation is to measure waveform quantities at nite radii within the Cauchy domain and then extrapolate them to in nity [52]. However, this extrapolation does not guarantee that the Einstein equations are satis ed, making it fail to capture memory effects [53]. A more rigorous solution is to use Cauchy-characteristic evolution (CCE) and Cauchycharacteristic matching (CCM) 6 [56], which can faithfully extract GWs at future null in nity up to Bondi-Metzner-Sachs freedom. Recently, a CCE [57,58] and a CCM [59] algorithm for GR have been built in a NR code SpECTRE [60]. Well-posedness analysis of CCE and CCM in Bondi gauges can be found in [61][62][63][64]. However, a robust CCE/CCM code for beyond-GR theories is still missing.

In this paper, we take a step to let the SpECTRE CCE system capture beyond-GR features. Speci cally, we choose to include a (massless) scalar eld, which serves as an important ingredient in theories such as Horndeski [65] and dynamical Chern-Simons gravity [66,67]. We will be testing the accuracy of our code and demonstrating its ability to reveal ringdown tails [68,69] in the scalar eld, as well as the memory effects sourced by the scalar eld.

This paper is organized as follows. In section 2, 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 modi ed in the presence of a scalar eld. Next in section 3, we focus on the memory effects sourced by the scalar eld and the corresponding balance laws. We then test our code by evolving scalar elds on two prescribed spacetime backgrounds (section 4) and performing the full CCE procedure to compute tails and memory effects (section 5). Finally, we summarize the results in section 6. Throughout this paper, complex conjugates are represented by overlines.

Einstein-Klein-Gordon system via CCE

The Einstein-Klein-Gordon action offers the most basic description for a system containing both scalar and tensor elds. For instance, after performing a conformal transformation [70], Damour-Esposito-Farèse scalar-tensor theories [70][71][72] in the Einstein frame are governed by:

where the real-valued scalar eld ψ is minimally coupled to the metric. The action leads to

Here, we see that the Einstein eld equations obtain a source term, while the scalar eld 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 complexvalued massive scalar eld with an arbitrary potential.

Equations of motion

By adopting the Bondi-Sachs metric of an asymptotically at spacetime [73][74][75] ds 2 = -e 2β (rW + 1) du 2

with x A = (θ, ϕ) standing for angular coordinates, the Einstein eld equations can be converted to a set of hypersurface equations [76,77]

with

The expressions of S β , S Q , S U , S W , S H , L H , and LH as functions of (J, β, Q, U, W) can be found in section IV of [57]. The spin-weighted derivative operators ð and ð are de ned to be

where D A is the angular covariant derivative compatible with the unit sphere metric. On the other hand, the KG equation for the scalar eld ψ becomes [76] 2r∂

where we have de ned an auxiliary variable Π = ∂ u ψ, and

In equation ( 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 rst u = const null slice, then radially integrate equations ( 4) and ( 6) by order to determine β, Q, U, W, H, and nally Π. Following this, we advance J and ψ to the subsequent null slice based on the values of H = ∂ u J and Π = ∂ u ψ.

SpECTRE CCE system

A spectral CCE system has been implemented for GR in SpECTRE [57,58]. To facilitate spectral implementations, the authors introduced numerically adapted coordinates (ȗ,y,x Ȃ):

where 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 ȗ = const null slice, the simulation domain spans radially from the worldtube at y = -1 (r = R wt ) to future null in nity at y = 1 (r = ∞). With this new coordinate system, the time (ȗ) derivative of J, H = ∂ ȗJ, is related to the original one H, via [57]

The second term arises as the Jacobian of the coordinate transformation. Similarly, the time derivative of the scalar eld, Π = ∂ ȗψ, transforms in the same way:

The values of other variables are not affected by the coordinate transformation, namely F = F, for all F ∈ {ψ, J, β, Q, U, W}.

Inserting equation (10) into equation ( 6) yields the hypersurface equation for Π:

with

Here we have used two identities:

Notice that in equations ( 11)-( 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 [57]:

To compute ð ψ, we follow [58] and use libsharp routines [78,79] to obtain a modal decomposition of ψ in terms of spin-weighted spherical harmonics; subsequently, we apply the differential matrix.

The transformation of the metric sector in equation ( 4) has been extensively discussed in [58]. Here we simply present additional terms contributed by ψ

(15)

Worldtube transformation

In order to integrate the hypersurface equation for Π in equation (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 [80]. 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.

It turns out that the transformation has a concise expression:

Derivation details can be found in appendix. In equation ( 16), ∂ t ′ ψ is the Cauchy time derivative of ψ; ðψ is the angular derivative of ψ (see equation ( 5)) with respect to the numerically adapted coordinates; and U (0) is de ned in equation (34b) of [57], with a spin weight of 1.

Physically speaking, the numerically adapted angular coordinates xȂ on the worldtube evolve over time with respect to the Cauchy coordinates, via (equation ( 28) of [57])

where the quantity U (0) captures the rate of the coordinate motion. Therefore, U (0) corresponds to the shift vector pulled back 7 to the 3D timelike worldtube. On the other hand, the numerically adapted time ȗ is chosen to ow at the same rate as Cauchy's time t ′ -the lapse is manually set to 1 [57]. As a result, the worldtube transformation in equation ( 16) represents the Lie derivative of ψ along ∂ ȗ, expressed in terms of Cauchy coordinates. Since ψ 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 equation ( 16), we simply need to communicate ∂ t ′ ψ from the Cauchy to the characteristic system.

After setting the boundary condition Π| wt , we integrate the hypersurface equation in equation (11) and evolve the scalar eld forward in time. On the worldtube, the scalar eld computed by the Cauchy system, ψ Cauchy , should agree with the one obtained from the characteristic system ψ Char . This yields a worltube constraint

In actual numerical simulations, the constraint is normally nonvanishing since the two elds 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 speci c examination for the worldtube transformation in equation ( 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 .

Computation of Weyl scalars at future null in nity

The SpECTRE CCE system [57,58] uses the Newman-Penrose (NP) formalism to compute Weyl scalars. The relevant equations are 8 [81]

where 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 [57]

Compared with their vacuum GR counterparts, the equations gain additional source terms Φ 01 , Λ, and Φ 21 . They are related to the scalar eld ψ through [81]

The asymptotic expansion of these quantities at future null in nity I + is of particular interest. Hereafter, we denote the terms in the expansion with the following notation:

where f represents any variable under investigation, and n is an integer. As for the massless scalar eld ψ, 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 obeysfoot_3 ψ ∼ r -1 . Substituting the asymptotic behavior of ψ into equation ( 20) leads to

where 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:

In the code, we compute the value of Π (1) , ψ (1) , and ðψ (1) via

The asymptotic expansions of Wely scalars in GR have been shown in equation ( 93) of [57].

With the presence of ψ, their expressions are modi ed to

Notice 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 equation (19a) is valid geometrically and does not rely on the detail of a gravity theory, the asymptotic expansion in [57] (see their equation (93a)) implicitly assumes the vanishing of the stress-energy tensor. Instead, here we directly expand equation (19a) and obtain:

The value of β (2) can be obtained by expanding the rr component of Einstein's equations (see equation (2.9c) in [75] or equation (2.9a) in [77]

Finally, the Weyl scalar Ψ

4 and the Bondi-Sachs News take the same forms as in GR, as shown in equations (93e) and ( 42) of [57], respectively.

Bianchi identities at future null in nity

One way to validate our CCE code is to examine Bianchi identities at future null in nity, 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 [83]:

where the source terms Φ 22 , Φ 00 , Φ 11 , Φ 02 , and Φ 20 are given by

In equation (30) we have used

to obtain the asymptotic expansions. The overhead dots denote the time derivative ∂ u .

Meanwhile, the Bianchi identities involve NP spin coef cients, whose expressions in terms of Bondi-Sachs variables can be found in equation ( 86) of [57]. Near future null in nity, their asymptotic behaviors read 10

Plugging equations ( 30)- (32) into equation (29), we obtain the following Bianchi identities at future null in nity

) 2 πψ (1) ððψ (1) π ḣ(1) ψ (1) 2 -√ 2π h(1) ψ (1) ψ( 1) .

(33e) 10 Our results are consistent with those in [84,85] after switching the convention and performing a tetrad transformation.

We will use the PYTHON package scri [86][87][88][89][90] to verify the identities in our following calculations. The package adopts a different convention (hereafter MB) [90,91] compared to SpECTRE CCE. The conversion factors read [52]

MB = √ 2σ (2) .

Then equations (33) become 1) ψ( 1) . (35e)

Scalar-induced memory effects

A massless scalar eld can propagate to future null in nity [92] and leave nontrivial imprints on memory effects. Here we extend the discussions in [53,93] and derive the memory contribution from the scalar eld. Since most of our discussions in this section adopt only the leading-order asymptotic behavior near future null in nity, we suppress the superscript de ned in equation (22) for conciseness, unless stated otherwise. To be consistent with previous discussions, we still use the convention of SpECTRE CCE, as opposed to MB (equation ( 34)).

The computation of memory effects can be achieved via the balance laws [75] (see their equations (2.11a) and (C6))

where m = -1 2 W (2) is the Bondi mass aspect; NA is related to the angular momentum aspect

Finally, C AB is the O(r -1 ) order of the angular metric h AB de ned in equation ( 3), and N AB = ∂ u C AB is the Bondi News tensor. Following [75,93], we decompose C AB into an electric (Φ) and a magnetic (Ψ) potential

where ϵ CA = i 2 q C ∧ qA is the volume form compatible with the unit sphere metric.

Electric memory

We insert equation (38) into equation (36a) and obtain [93]

Here α(θ, ϕ) 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 ux ∼ ´ψ 2 du, compared to its GR counterpart [93]. On the other hand, the Bondi mass aspect contributes to the ordinary (linear) memory, and is given by [85]

Notice that the scalar ordinary memory vanishes in a nonradiative regime ( ψ ∼ 0). The potential Φ is related to the electric memory J E via [93]

where the scalar pieces read

The tensor contributions J ET o and J ET null are the same as GR, and we provide their expressions below for completeness

Magnetic memory

The magnetic memory can be obtained by projecting both sides of equation (36b) with ϵ AE D E , and the corresponding expression for Ψ in GR is provided in equation ( 45) of [93]. With the presence of the scalar eld, a convenient way to compute Ψ is by replacing Ṅ = q A ṄA in the GR version with the following combination:

which yields

To obtain N, we rst contract equation (37) with q A and get

In GR, the angular momentum aspect N is solely determined by Ψ 1 , e.g. see appendix B of [93]. As for the Einstein-Klein-Gordon system, an extra term shows up. To see this, we use equations (B2), (B5) and (B6) in [93], namely

together with the asymptotic expansion for Ψ 1 in equation ( 26) and the one for β in equation ( 28), then we arrive at

Plugging equations ( 46) and ( 48) into equation ( 45), and by virtue of the Bianchi identity in equation (33d), we nd the scalar contribution to Ψ vanishes identically. The corresponding magnetic memory

with

The expression can be further simpli ed by replacing the term ð( hð ḣḣðh) in J B o with its negative complex conjugate 11 . After combing with J B null , we obtain

where we have used the fact that the Bondi mass aspect is real-valued, namely [84] Im

Notice that the right-hand side of equation ( 51) is linear in h, as opposed to the electric memory in equation ( 41) that depends on energy uxes, the balance law for the angular momentum 11 Notice that only the imaginary part contributes to J B . aspect in equation (36b) in fact de nes a projection operator for the magnetic component of h. In particular, equation ( 51) implies that

It corresponds to the radiative current moment, see e.g. equation ( 30) in [94], which plays an important role in building up gravitational recoils [89,94,95] and the excitation of quadratic quasinormal modes [96,97]. Under parity conjugation, J B lm transforms as follows

where we have used equation (C10d) in [89]. 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 [98].

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 [58], our tests primarily focus on the new features introduced by the scalar eld. This involves two main parts: solving the KG equation in equation (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.

Below in this section, we test the implementation of the KG equation by evolving scalar elds 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 equation (16). For (a), we compare the computed scalar modes at future null in nity to expected (analytical) behaviors. For (b), we monitor the worldtube constraint de ned in equation ( 18), as discussed in section 2.3.

Bouncing Schwarzchild BH

Following [58,99], we consider a Schwarzchild BH in an oscillating coordinate system-a time-dependent transformation is applied to KS coordinates {t KS , x KS , y KS , z KS } [58,99]:

where a and b are two constants. Since the transformation is a pure gauge effect, one expects no GWs at future null in nity, 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 in nity. This makes the bouncing BH system the most demanding test [58].

For our testing purposes, it is preferable to construct a scalar pro le that can be controlled analytically, allowing us to compare numerical results with the expected analytical expression. To do this, we rst work in the outgoing Eddington-Finkelstein coordinate system {u, r, θ, ϕ}. The retarded time u is related to the KS time t KS via where we have used equation ( 56) and

In our tests, we follow [58] and set the oscillation period b to 40M and the amplitude a to 2M. Four worldtube radii, 15M, 20M, 25M, 30M, are used for CCE. We choose the absolute tolerance of SpECTRE's time stepper residual to be 10 -foot_4 to ensure that differences between our results and the expected one ψ 0 (u) = sin u are dominated by spatial resolution. Figure 1(a) 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 gure 1(b), the worldtube constraint de ned in equation ( 18) follows a similar convergence behavior and nally levels off at ∼10 -12 .

Linearized Bondi-Sachs metric

In the second test, we choose the linearized Bondi-Sachs metric [100] as our spacetime background. The corresponding worldtube data read [99]

Following [58], we consider a l = m = 2 angular pro le, with s Z lm given by 12

The expression of J l (r), U l (r), β l (r) and W l (r) can be found in equations ( 139) and ( 140) of [99]. In particular, we have

Here B 2 , C 2a and C 2b are three free complex constants. Below we set B 2 to 0 for simplicity, which subsequently yields 3C 2a = ν 2 C 2b due to the asymptotic atness condition J l=2 ∼ O(r -1 ). Therefore, the size of the linearized Bondi-Sachs wave is controlled by a single parameter C 2a . As for the scalar sector, we choose the worldtube data to be

When the linearized Bondi-Sachs wave vanishes (C 2a = 0), the scalar eld propagates through at spacetime. The consequent wavefunction at future null in nity is simply given by rψ(u)| I + = sin u. However, for a nite C 2a , an exact analytical solution is not available. Instead, we can consider a perturbative expansion:

To the rst order in C 2a , ψ (1) (u, θ, ϕ) is generated by the coupling between metric components (characterized by the angular pro le s Z l=m=2 ) and the leading-order evolution of ψ, namely sin u (characterized by the angular pro le 0 Z l=m=0 ). The angular selection rule immediately leads to

In other words, the amplitude of the l = m = 2 harmonic of the scalar eld is expected to scale linearly with C 2a . By contrast, other harmonics of the scalar eld are generated by nonlinear couplings, making their amplitudes depend at least quadratically on C 2a .

To perform our tests, we x the frequency of the linearized Bondi-Sachs wave at ν = 0.2, while varying the wave amplitude C 2a from 10 -6 to 1. To ensure the accuracy of our analysis, we still set the absolute tolerance of the time stepper to 10 -12 . Meanwhile, we nd different choices of angular resolution l max lead to only a ∼10 -12 difference in nal products. Figure 2(a) shows the amplitude of ψ l=2,m=0 (blue dots) and ψ l=m=2 (orange dots) as a function of C 2a , measured at future null in nity. As expected, they depend quadratically and linearly on C 2a , respectively. In gure 2(b) we provide the worldtube constraint (equation ( 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 2a , consistent with the behavior at future null in nity.

Cauchy characteristic evolution

After examining the characteristic evolution of the KG equation in equation (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 rst construct initial data by solving the extended conformal thin-sandwich equations with the spectral elliptic solver in SpEC [101,102]. Next, we evolve the system nonlinearly using a rst-order generalized harmonic formulation [103]. Our Cauchy evolution for the scalar eld follows the method in [80]. Finally, we use the dumped worldtube data for CCE, as outlined in [57,58] and earlier in this paper. For conciseness, we always set the initial BH mass to unity and use it as our code unit.

Below 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 eld decays as a power law [68,69]. In particular, the power law has different exponents at timelike in nity and null in nity [68,69,[104][105][106]. The late-time tail at timelike in nity was investigated by Scheel et al [80]. Below in section 5.1, we focus on the behavior at null in nity to demonstrate the accuracy of our CCE system. (b) As discussed in section 3, the GW (tensor emission) of a BH stirred by a scalar eld consists of memory effects. They are governed by the balance laws in equations ( 41) and (49). In the second part of this section (section 5.2), we will examine the deviations from the balance laws, as well as the Bianchi identities derived in section 2.5.

Throughout this section, we initialize the scalar eld with a linear ansatz:

where the coef cient A is determined from the worldtube data obtained via the Cauchy evolution. We note that this initial data is constructed fully empirically. Developing an ab initio initial-data solver is beyond the scope of this paper and is left for future work.

Late-time tail

In the rst test, we place a Schwarzschild BH at the coordinate center. Initially, the scalar eld pro le is constructed as follows: with 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, wherefoot_5

is an auxiliary variable introduced by Scheel et al [80]. The choice of Π ′ represents an infalling scalar pulse into the BH. As shown in gure 3(a), the blue curve corresponds to the l = m = 1 spherical harmonic of the extracted scalar eld 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 eld is extracted at future null in nity, where the Bondi radius r approaches ∞ while the retarded time u remains xed. By comparison, we also obtain the scalar radiation at an approximate timelike in nity by measuring late-Cauchy time ψ at a xed Cauchy radius, as done by Scheel et al [80]. The orange curve in gure 3(a) shows the corresponding l = m = 1 scalar harmonic. To make a fair comparison with CCE, our extraction radius is still at 75. We can see that ψ 11 at both timelike and null in nity evolves similarly during the excitation and ringdown phases, whereas their late-time tails exhibit distinct decay rates. Speci cally, in gures 3(b) and (c) we t 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 t:

where t i 's are time samples. To improve the t 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 t to a one-dimensional minimization problem, and we deal with it using the Nelder-Mead method in SciPy [107]. We have checked that the results are insensitive to the initial guess. The best ts are plotted as orange dashed lines in gures 3(b) and (c). We nd the tail exponents at null in nity (µ n ) and timelike in nity (µ t ) to be 2.93 and 4.92, respectively. They are consistent with the predictions of BH perturbation theory [68,69,[104][105][106]: µ n = l + 2 and µ t = 2l + 3.

Memory effects

In the second test, we consider a Kerr BH with a dimensionless spin χ = 0.7, initially surrounded by a scalar eld with Π ′ = 0 [de ned in equation (73)]. The initial pro le of ψ is chosen to be the same as equation (72), with A = 1, r ′ 0 = 20 and w = 2. The CCE worldtube is placed at r ′ wt = 100M. Figure 4 shows the l = m = 1 harmonic of the scalar eld measured at future null in nity. The scalar eld initially contains an ingoing and an outgoing component in the radial direction. The outgoing piece reaches the CCE worldtube and appears at future null in nity at a time of r ′ wt -r ′ 0 = 80M; whereas the ingoing wave rst strikes the Kerr BH at r ′ 0 = 20M, the excited quasinormal modes are then transmitted to future null in nity by CCE at r ′ wt + r ′ 0 = 120M.

Due to the no-hair theorem [26,38,[108][109][110], a Kerr BH does not carry a scalar charge. Consequently, the scalar eld vanishes when the BH settles to stationary. During the dynamical process, the scalar eld stirs spacetime nonlinearly (equation ( 4)) and excites GWs. The blue curves in gures 5(a) and (b) 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 section 3, a strain can be decomposed into an electric and a magnetic memory. In this case, we nd that the magnetic sector vanishes. Three dominant electric components (see equation ( 41)) are plotted in orange, green, and red in gures 5(a) and (b). 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 in nity. Unlike in BBH systems, where the tensor-driven null component J ET null induces a strong memory effect [93], its contribution is negligible in this process. Two lower panels of gures 5(a) and (b) 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 gure 5(c). 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, gure 5(d) shows the (l = 3, m = 2) harmonic of the spin memory [111], as a time integration of h 32 .

As discussed in section 2.5, the Bianchi identities in equation (35) yield relations between the Weyl scalars, the strain, and the scalar eld. Meanwhile, equation ( 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 rst 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 ).

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 eld in the equations of motion (section 2.1), the computation of Weyl scalars (section 2.4), the Bianchi identities (section 2.5), as well as the memory effects (section 3). In addition, we reformulated the characteristic form of the KG equation using the numerically adapted coordinates [equation ( 8)], facilitating its implementation in our NR code SpECTRE. We also derived a concise worldtube transformation to construct the boundary condition for the characteristic KG equation from a generic Cauchy evolution (section 2.3).

To evaluate the accuracy of our code, we designed various test systems. We rst focused on the implementation of the KG equation and the worldtube transformation by evolving scalar elds 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 eld decays as a power law u -l-2 , consistent with Price's law at future null in nity. Furthermore, we veri ed that the tensor emissions extracted with CCE are consistent with the balance laws and the Bianchi identities. 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). (c) and (d): 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 equation (50).

The 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 example, 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 [112][113][114], allowing us to make comparisons with the post-Newtonian approximation [115][116][117]. These numerical studies will improve [77,[118][119][120]. lays a foundation for performing CCE simulations modi ed theories of gravity, e.g. [35,50], paving the way for studying the associated memory effects [114].

Our simulations show tails in scalar elds. A possible avenue for future work is to look for tails in breathing modes (e.g. [42]) and tensor modes [121][122][123] emitted by binary systems. The accurate simulation of these tails is sensitive to the choice of boundary conditions [124]. 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 [59].

Finally, it would also be interesting to apply our CCE algorithm to extract scalar emissions from intermediate-mass-ratio inspirals. Recently, a worldtube excision method [125,126] 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 nite radius to perturbative calculations. Our CCE algorithm could further improve this comparison. Additionally, the algorithm allows for verifying balance laws at future null in nity.