Probing supermassive black hole binaries with orbital resonances of laser-ranged satellites

Coalescing supermassive black hole binaries (SMBHBs) are the primary source candidates for low frequency gravitational wave (GW) detections, which could bring us deep insights into galaxy evolutions over cosmic time and violent processes of spacetime dynamics. Promising candidates had been found based on optical and X-ray observations, which claims for new and ready-to-use GW detection approaches before the operations of space-borne antennas. We show that, satellite laser ranging (SLR) missions could serve as probes of coalescing SMBHBs through the GW-induced resonant effects. Lasting and characteristic imprints caused by such resonances in the residual distances or accelerations from SLR measurements are studied, and the detection SNR is analyzed with both the current and future improved ranging precisions. Within redshift z∼1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z \sim 1$$\end{document}, the threshold SNR = 5 requires 1–2 years of accumulated data for the current precision and months of data for improved precision, which are workable for the data processing of SLR missions. Meanwhile, joint detections with multiple SLR missions could further improve the total SNR and the confidence level. Such a detection scheme could fulfill the requirement of a tentative SMBHB probe during the preparing stage of LISA and Taiji, and it requires no further investment to any new and advanced facilities. It is also worthwhile to look back and re-process the archived data from the past decades, in where resonant signals from SMBHBs might be hidden.


Introduction
Ever since the landmark event GW150914 observed by Adv-LIGO [1,2], the new window to the Universe had been opened by GW detections.With the follow-on observations of nearly one hundred of events by the LIGO-VIRGO collaboration [3,4], the new era of GW astronomy has gradually started off.To enclose the exciting sources of much larger and heavier astrophysical systems, one needs to explore the low frequency end of the GW spectrum with detectors of much longer baselines.In the next decade, the first generation space-borne antennas, including LISA (Laser Interferometer Space Antenna) [5] and the LISA-like Taiji mission [6][7][8][9] would be launched and cover the mHz band [10].Mission concepts for decihertz band now include the DECIGO [11], AMIGO [12] and TianQin [13] projects.For the µHz band, the more advanced space missions with baselines of Solar system size like the µAres [14] and ASTROD [15] were suggested.
Among the candidate sources within the low frequency range, coalescing SMBHBs (10 6 −10 9 M ) are of the most exciting ones [19].Their collisions and mergers give rise to the most violent events in the visible Universe and produce the loudest GW signals within their frequency band.The knowledge of the population of coalescing SMBHBs and the detailed measurements of the entire wavetrains (inspiral, merger, and ringdown) will bring us deep insights into the growths and co-evolutions of the SMB-HBs and their host galaxies, the expansion history of our Universe, the most violent dynamic behavior of curved spacetime and the nature of gravitation [19][20][21].Therefore, coalescing SMBHBs are the primary candidates for many aforementioned missions.According to the roadmap of GW physics and astronomy [22], LISA will firstly reveal the detailed information about SMBHBs till the 2030s.On the other hand, SMBHB candidates could be identified by means of optical and X-ray observations.For example, Ref. [23] claimed that the rapid decaying binary system SDSSJ1430+2303 discovered via optical and X-ray observations is expected to merge within 3 years, which might provide us an excellent opportunity for multimessenger observations if any low frequency GW detector was ready [24,25].Having the plentiful scientific objectives and the potential opportunities, it would be of great significance if any tentative detection of coalescing SMBHB could be made within this decade before the operations of LISA and Taiji missions.Similar to the case for a bound system of charged particles, which could have resonant interactions with incident electromagnetic waves when the wave frequency matches with the energy difference between certain states of the system, the response of a self-gravitating binary system to GW is especially evident when the frequency of GW matches a harmonic of the binary's orbital frequency, thereby inducing a resonant effect.Such resonant evolution of the orbit can accumulate in time, and may eventually enter the scope of possible detections.The studies of resonant responses of self-gravitating binary systems to incident GWs can be traced back to the 1970s [26][27][28][29][30], and had raised more concerns in these years [31][32][33][34][35][36][37][38][39][40].Recently, Blas and Jenkins had made important progress for this detection scheme and developed powerful tools to estimate the resonant evolutions of the binary orbital elements under stochastic GW backgrounds [41].Sensitivities to such resonant effects from stochastic backgrounds in SLR, Lunar laser ranging (LLR) and PTA missions turn out to be promising [42], which approves the detectability of such scheme.
Most of the GW signals studied in the aforementioned literature are stochastic in nature, and still the resonant detection scheme for general deterministic and individual signals is not yet fully investigated.In this work, we aim to study the feasibility and prospect of detecting, in the circumstance where orbital resonance takes place, the deterministic GW signals from coalescing SMBHBs via the SLR technique of the present day and the next decade.The Earth-satellite distances can be measured precisely and continuously for SLR missions, which hence provides the long-term and faithful tracking of the orbital evolutions of the satellites [43,44].In fact, data from SLR missions, such as LAGEOS 1, 2 and LARES, had already been proved to be very useful in relativistic experiments [45][46][47][48][49][50][51][52][53].For GW detections, the orbital harmonics of all the in-orbit laser-range satellites [54] could form a "comb" in the sub-mHz range, which makes it possible to capture the strong chirping signals from coalescing SMBHBs consecutively by most of the missions in operation.Based on previous studies [29,[55][56][57][58], and especially [41], we give analytical and numerical analysis of the resonant evolutions of osculate orbital elements induced by GWs from coalescing SMBHBs in SLR missions and their dependence on relevant parameters.An important finding is that, when such resonance takes place, a characteristic signature is left in the orbital evolutions of the laser-ranged satellites.With the precision and multi-year data of orbit tracking, the resonance signal could be recovered with sophisticated data analysis methods.Given the joint detections of each individual signal by different SLR missions, this would fi-nally produce the high-confidence detection of a coalescing SMBHB.
Limited by the precision of SLR measurements in presentdays and in the near future, such detections may not give rise to detailed physical parameters of the sources, but the multi-year observations could probably provide us the first estimation of the population or event rate of coalescing SMBHBs within redshift z ∼ 0.1.Such a detection scheme could fulfill the requirement of a tentative SMBHB probe in the decade before the launch of LISA (and Taiji), and it requires no further investment to any new and advanced facilities.The only efforts demanded will be the thorough analysis of the data, especially the re-analysis of the archived data from the past decades in where resonant signatures from coalescing SMBHBs might be hidden.

Theoretical Tools
For clarity, we will refer to the SMBHB as the "source binary" and the Earth-satellite system as the "test binary".The phenomenon of orbital resonance can be described by the equations of motion (EoM) of the osculating orbital elements of the test binary, with GWs acting as small perturbations.
We introduce a cylindrical coordinate {r, θ, ˆ } whose origin is placed at the test binary's center of mass, and {r, θ} represent the bases of polar coordinates within the orbital plane, ˆ being the unit vector perpendicular to them.The unperturbed Keplerian motion of the satellite is characterized by six orbital elements X = {P, e, I, Ω, ω, ε}, including the orbital period, eccentricity, inclination, longitude of ascending node, argument of pericenter, and the compensated mean anomaly.
In the absence of perturbations, the separation r(t; X) of test binary, that is related to the most important observable of SLR missions, follows the Kepler's equations where the true anomaly ψ is defined as the angular position of the satellite measured from the pericenter, and E is the so-called eccentric anomaly.With perturbations induced by incident GWs, Eq.(1-3) are valid under the condition that X are regarded as the osculating orbital elements, which varies with time satisfying the RTN-type Gaussian perturbation equations [56,57,59].The effects of GWs, treated as small perturbations, can be written down in the form of Newtonian forces [29,60]: where R, T and N are the perturbing forces per unit mass in the radial, tangential and normal directions relative to the orbit e A ij being the polarization tensors of the GW (A = +, ×).Finally, in terms of transfer function, the EoM of X can be written in a compact form [41]: where nGW = nGW (ϑ, φ) denotes the direction of source.
The transfer functions T A define the linear responses of the orbital elements X to the incident GWs.Such linear responses are valid under the conditions that perturbations from GWs are small compared to that from Newtonian gravity and the back reactions from resonance to the incident GW field are ignorable.The explicit forms of the transfer functions T A and the coefficients of e A ij in the test binary frame can be found in the work [41].For SLR measurements, we are most concerned about the orbital period P (or the semi-major axis a), which is directly related to the total energy of the test binary.The transfer function of P can be expressed as where γ ≡ √ 1 − e 2 .Eqs.(6) constitute a system of ODEs, which are solved numerically in the following analysis using the initial conditions X| t=0 = X 0 .The orbital perturbation theory suggests that ψ(t) can be calculated by solving the Kepler equations Eq.(2) and Eq.(3) with the orbital elements regarded as the osculating ones.
To understand the resonant behavior qualitatively, we also derive an analytical solution under the simplifications that the orbit of test binary is nearly circular (e 1) and the incident GW is modeled as a monochromatic wave with redshifted frequency f GW , initial phase ϕ GW and amplitudes H A .At the "main" resonance frequency f GW = f res = 2/P , the variation of P is dominated by a linear drift term.The secular perturbation of P , defined as Ṗ averaged over one orbit revolution, reads where G 1A , G 2A are constants determined by the angles of the GW source and the test binary (see Appendix A), and δ A× = 1 if A = × or 0 if A = +.In the case where the source binary and test binary are face-on, Ṗsec can be expressed more concisely: with H ≡ H A (ι = 0).Depending on the values of ϕ GW , ω and ε, Ṗsec can be either positive or negative.Besides, two "secondary" resonances of order O(e) occur at f GW = 1/P and 3/P .The detailed derivation of this solution can be found in A. These discussions provide an intuitive demonstration of orbital resonance.

An Example of merging SMBHB
We assume a SMBHB with the parameters of SDSSJ1430+2303 [23] as an example for the SMBHBs that will merge in the near future.Although the interpretation and detectability of SDSSJ1430+2303 are still under discussion, we are interested in a group of similar SMBHB systems instead of this specific one.This example will be referred to as our representative Target Source in the rest of this paper (TS for short), and the variation of its parameters over a wide range will be discussed in B. Among all the SLR missions, we take the laser ranging mission LAGEOS 2 (L2 for short) as a typical representative.The resonant responses of other SLR missions, including LAGEOS 1, LARES 1/2, Ajisai and ETALON 1/2 to the same GW signal are also discussed.The parameters of these missions can be found in [52,61].
The parameters relevant to our example are listed in Tab. 1, where the initial values of P, e, I for L2 are taken from Ref. [52].Ref. [23] reported the properties of the TS constrained from electromagnetic observations.It is proposed that this system is an uneven mass-ratio, highly eccentric SMBHB.While, at the frequencies of our interest, its orbit would be sufficiently circularized.Moreover, the components masses of the TS are only determined with large uncertainty.Therefore, here we simply assume that it consists of two black holes with equal source-frame masses M bh = 4 × 10 7 M , and redshift z = 0.08105.The sky position and inclination of the TS, as well as Ω and ω are randomly selected, since Ω and ω could change in time, and we are interested in a family of SMBHBs with properties similar to SDSSJ1430+2303 rather than this specific one.A detailed discussion on the impacts of the mass ratio and other parameters is given in B. Regarding the modeling of the GW signal, we utilize the SEOBNRv4 time-domain waveform [62,63] provided by the open-source code PyCBC [64], and the phase at coalescence is set as ϕ c = 0.
Roughly speaking, the inspiral stage ends when f GW equals the GW frequency of the innermost stable circular orbit f ISCO ≡ 1/(6 3/2 πM ), where M is the total mass of the source.For the case under consideration, f res ≈ 1.5 × 10 −4 Hz > f ISCO , indicating that resonance happens mainly during the merger stage.As is shown in the top panel of Fig. 1, for L2, P exhibits a monotonic growth for ∼ 10 4 s, and finally reaches a steady value with ∆P fin /P 0 = 4.374 × 10 −14 .
The resonant responses of different SLR missions to the same GW signal from the TS are also shown in Fig. 1.As expected, the secular variation Ṗsec can be either positive or negative, depending on the parameters of satellite orbits and GW sources.The resonance frequencies of these satellites form a "comb" in the sub-mHz frequency range, and resonances would take place consecu- tively among these SLR missions when the chirping signal sweep across the "comb tooth".Hence, the correlations among such resonant events could give a high confidence level of the detection, and may even help to investigate the physical properties of the corresponding GW sources, like the TS.To make such detection scheme attainable, the more sophisticated and important observable, that of the residual separation δr(t), is defined and employed in the following discussions, see Fig. 2 and 3.

GW Detection with orbital resonance
Based on the above example, the GW-induced secular change in semi-major could only reach ∆a fin ≈ 2a 3P ∆P fin ∼ 10 −7 .Compared with the resolution of SLR distance measurements, which is at millimeter or sub-millimeter level [51], it seems difficult to identify such small changes in the semimajor out of uncertainties.
On the other hand, SLR is particularly superior in tracking the orbital dynamics of satellites.Collecting the round-trip times of laser pulses allows one to track the "normal point" distances over time.The orbital elements of the laser-ranged satellite can be derived based on such distance measurements with the help of the precise orbit determination programs, such as GEODYN [65].This inspired us to make use of more sophisticated observables instead of the averaged orbital elements, that of the residual normal point distance δr(t) or the residual acceleration δ a(t), to reveal the signatures from the GWs of coalescing SMBHBs.The residual distance is defined as where r data (t) is obtained from the SLR measurements, and r(t; X 0 ) is calculated from the initial elements X 0 via Eq.( 1).δr mod (t) consists of the contributions of all other modeled perturbations except for GWs.The residual acceleration is defined in the similar way δ a(t To make use of such data, one needs to accurately model and account for the possible gravitational and non-gravitational perturbations.A wide variety of perturbations have been investigated in the literature, such as Earth geopotential harmonics [66,67], atmospheric drag [68,69], thermal-thrust effects [70], Solar radiation pressure, dynamic solid tide and ocean tide [71,72], etc.
With these tools, one could track the long-term dynamical evolutions of the orbits and investigate in details Note that the resonance frequencies of ETALON 1/2 are much smaller than other satellites, therefore for these two satellites, we have labeled the times corresponding to fres, rather than showing them as vertical lines.Resonances take place consecutively among these SLR missions when the chirping signal sweep across the "frequency comb" in the sub-mHz frequency range.
the differences (residuals) between the observed and modeled orbits to search for the expected signals.Such data analysis method is slightly different from the one used for interferometric GW detectors [73], but has already been employed in SLR missions in a wide range of literature [49,51,69,72,74].Considering residual accelerations, after modelling the known perturbations, the residual mean accelerations deviated away from the geodesic motion for L2 (or LARES) are less than 1-2 × 10 −12 m/s 2 (or 0.5 × 10 −12 m/s 2 ) [49,74].While, for the optimal response of L2 to TS (see B for the determination of "optimal parameters"), the radial residual acceleration δa r will oscillate around −2.5×10 −13 m/s 2 at the start of the postresonance stage.This order-of-magnitude estimate gives a rather optimistic evaluation of the feasibility of this new detection method.
For the convenience of SNR estimations, we use the residual distance in the following analysis and assume the ideal case that the only perturbation to the satellite orbit is from the incident GWs of SMBHBs.And, considering the expected event rate of coalescing SMBHBs [19], resonances in SLR measurements can be treated as individual events.Shown in Fig. 2 is the optimal response of L2 to TS in terms of δr(t), and in Fig. 3 where M = 2πt/P +ε is the mean anomaly.That is, under the long-term condition (3πt/P 1), δr(t) would oscillate around ∆a fin with linearly varying amplitude, and the rate of variation is proportional to e.After subtraction of the known and modeled perturbations, if the similar behaviors as in Fig. 2 and 3 were observed in the residuals for different SLR missions, it would indicate with high confidence that the GW-induced resonance as the cause.Another conclusion which can be drawn from Fig. 3 is that the eccentricity of satellite orbit plays an important role in the post-resonance evolution, as is predicted by Eq. (11).Indeed, the growth rate of δr for L2 (e = 0.0135) is relatively large compared to, for example, LAGEOS 1 (e = 0.0045), ETALON 1 (e = 0.001) and LARES 1 (e = 0.001).
To extract the signal of GW, the residual data should be analyzed with methods such as matched filtering.Based on long-term data tracking, the SNR (dubbed ρ) for the optimal response of L2 can be approximated as where δr(f ) is the Fourier transform of δr(t), and S n (f ) represents the one-sided noise power spectral density of SLR.The resonance stage lasts less than 1 day and contributes a rather small fraction to the total SNR.Whereas, in the post-resonance stage, δr(t) oscillates with growing amplitude and results in the above polynomial SNR on the observation time T obs (see Appendix C). S n (f ) depends on the uncertainty σ of normal point measurement.Ref. [42] predicted that the precision of SLR coincides with that of LLR, and the latter will have an order-of-magnitude improvement in the next decade (which might require the installation of new retroreflectors [75]).Following their assumption, we consider two values of σ in this paper: a. current precision: σ = 3 mm, 50, 000 normal point measurements per year; b. improved precision: σ = 0.3 mm, 200, 000 measurements per year.
In Fig. 4 we plot the SNR against T obs for different missions under different laser ranging precision.As is shown, for L2, SNR > 1 can be achieved when T obs > 124 days (current precision a) or 17 days (improved precision b).Once we set a realistic threshold for GW detection to SNR = 5, T obs > 356 days are required for precision a and 49 days for precision b, which turns out to be practical and workable.
In the derivation of Eq.( 12), it is implicitly assumed that we have perfect knowledge of the unperturbed orbital period, thus S n (f ) does not account for the contribution of the prior uncertainty of P 0 (dubbed σ P0 hereafter).This can only be achieved with infinitely long in-orbit time used for the calibration of orbital elements.Therefore, it is necessary to examine whether the impact of σ P0 can be safely neglected compared to the effect of passing GWs, provided that P 0 is determined from a reasonable number of SLR measurements.To estimate σ P0 , we employ the Fisher information matrix (FIM) formalism presented in Ref. [41], which converts the uncertainties of laser ranging to those of orbital elements.Within a given time T (which does not have to coincide with T obs ), where i represents the parameters relevant to r(t), i.e. i = P, e, ε.The σ P0 − T relationship for L2, as well as the target signal (e.g. the optimal response of L2 to TS) are plotted in Fig. 5.It is clearly shown that, under precision b, the prerequisite (σ P0 ∆P ) for our SNR calculation is fully satisfied with T ∼ 10 3 days .This is achievable for missions like LAGEOS which have been operating for decades.While for precision a, the SLR data collected over a period of ∼ 10 3 days only yield a σ P0 comparable to the signal, necessitating a longer time of calibration, or the signal should be stronger by one or more orders (such as the examples given in the following analysis).
In addition, it should be noted that although the prior uncertainty of P 0 has been considered, this is still an ideal scenario in the sense that other gravitational or nongravitational perturbations are well modeled and subtracted from the data.Given the specific models of these perturbing factors to the SLR missions, the error analysis would become more realistic, while this is beyond the scope of our work and belongs to a separate research topic.
The magnitudes of resonant responses also depends on the properties of GW sources, especially the masses and redshifts (see Appendix B.4), we then go beyond TS and look for more promising candidates.The nearest SMBHB system reported so far is located in NGC 7277, with redshift 0.006 and component masses 1.54 × 10 8 M and 6.3 × 10 6 M [76].Although this system itself is not an imminent merging one, the existence of SMBHB at redshift around z ∼ 0.01 can not be ruled out.Suppose that L2, as the representative, is in resonant interaction with the GW from a SMBHB with redshift z = 0.01 and equal component masses M bh = 5.9 × 10 7 M , in the most optimistic case, SNR = 5 can be achieved after an observation time of 68 days (current precision a) or 9 days (improved precision b).Furthermore, for an imaginary SMBHB at z = 0.001, by only taking the data within resonance stage into consideration, the most optimistic SNR could reach 0.15 (current precision a) or 3 (improved precision b).To the far end, for mergers of TS-like SMBHBs at z = 0.1, data sets of 410 days (current precision a) or 56 days (improved precision b) are needed to achieve SNR > 5.The above analysis indicates that this new method could give tentative detections of the violent mergers of SMBHBs within the reach of z ∼ 0.1.Moreover, as expected, when considering the joint detection by all SLR missions in operation, both the total SNR and the confidence level could be further improved.

Concluding remarks
In this work, we have investigated the feasibility of the detection scheme for GWs from coalescing SMBHBs, through their resonant interactions with the laser-ranged satellites.The observable of residual distance or residual acceleration are introduced to make the detection attainable.The SNR of measuring the resonance-induced characteristic signals in the residual distances and the dependence on relevant parameters of GW sources and orbiters are analyzed.It turns out that, before the launches of the space-borne antennas, SLR may be the only ready-to-use approach of probing coalescing SMBHBs in the sub-mHz range.
Among the promising candidates, we take SMBHB SDSSJ1430+2303 (TS) as the representative example, which is expected to merge within 3 years.For LAGEOS 2, we discussed the SNR of detecting the merger of TS for two sets of ranging precision.In the optimistic case, SNR > 1 can be achieved when T obs > 124 days for the current precision or 17 days for the future improved precision.For the threshold SNR = 5, 356 days is required for the current precision, and 356 days for improved precision.These results are generated to similar sources and SLR missions.For TS like candidates at redshift z = 0.1, SNR = 5 requires less than two years data for the current precision and a few months data for the improved precision.These are workable for data processing of SLR missions.Moreover, as the chirping waves sweep across the"frequency comb" in the sub-mHz range, the possible joint detection by the multiple laser-ranged satellites could further improve the total SNR and the detection confidence.
To summarize, SLR missions with the resonant detection scheme could fulfill the requirement of a tentative SMBHB probe that within the reach of z ∼ 1.Not just future-oriented, the re-analysis of the archived data from the past decades with our method is also worthwhile.At last but not least, the prospect of our method also improves with the understandings of the total orbital perturbations for SLR missions.

Appendix A Derivation of the simplified analytical solution
An analytical solution can be obtained under the conditions that the test binary's orbit is near circular (e 1) and the incident GW is modeled as a monochromatic wave: where f GW and ϕ GW are the redshifted frequency and initial phase of GW, respectively.The amplitudes H A (A = +, ×) can be expressed in terms of the redshifted chirp mass M c , luminosity distance d L and inclination ι as It will be demonstrated that the resonance is strongest in the face-on case.In this scenario, the inclination angle ι equals zero, thus the polarization angle ψ P and the initial phase of GW ϕ GW are complete degenerate.Hence for simplicity we will omit the dependence on ψ P in the subsequent analysis.
In the case of small orbital variation, we can expand X in powers of h A as X = X 0 +X 1 (t; h A )+O(h 2 A ), insert this expansion into Eq.( 6) in the body of the paper, and keep only the linear order.Taking the orbital period as an example, it follows that In the rest of this appendix we will drop the subscript "0" for brevity, and one should keep in mind that X in the r.h.s.stand for X 0 .
The most cumbersome parts of T A P are e A ij ri rj and e A ij ri θj , which take the following forms in the frame of test binary: with θ ≡ ω + ψ, and sin 2ϑ cos ϕ sin I + sin 2 ϑ sin 2ϕ cos I , where we have defined ϕ ≡ φ − Ω for convenience.
For elliptic orbits, the explicit time-dependence of Kepler motion can be expressed in terms of the Hansen coefficients [77,78]: where M is the mean anomaly M = 2πt/P +ε.As a result, up to the linear order of e, we have In the first line, α ≡ f GW P denotes the ratio between GW frequency and the orbital frequency of the test binary, and δ A× = 1 if A = × or 0 if A = +.The factors F iA and G iA are combinations of C A r , C A θ , S A r and S A θ defined as By integrating Eq.( 20), it is straightforward to show that P has a secular and linear evolution when α = 2, which is of order O(e 0 ), and we will refer to α = 2 as the "main" resonance frequency.Besides, two "secondary" resonances of order O(e 1 ) occur at α = 1 and 3, indicating that the resonant responses of eccentric binaries are more complicated than the circular ones.
Within the context of our discussion, e ∼ 10 −3 − 10 −2 , thus for the simplified analytical solution we will mainly focus on the "main" resonance frequency f GW = f res = 2/P .The secular evolution rate of P , defined as Ṗ aver-  aged over one revolution, reads Eq.( 22) can be further simplified if we consider the case where the orbital planes of the source binary and the test binary are face-on i.e. {ϑ = I, φ = Ω − π/2, ι = 0}, which gives Ṗsec = 12πγH sin(ϕ GW − 2ω − 2ε), (23) where H ≡ 4M 5/3 c (πf GW ) 2/3 d −1 L .For circular orbits, ω and ε become ill-defined.This issue is usually solved by introducing the combination ξ = ε + ω, which, for a stable circular orbit, stands for the angle from the ascending node to the initial position.Therefore Ṗsec = 12πH sin(ϕ GW − 2ξ).
B A detailed discussion on the parameter space and the determination of optimal parameters The orbital resonance of the satellite depends on several parameters, including the position, orientation, redshift, and component masses of the GW source, and the initial orbital elements of the test binary.For monochromatic sources, the relationship between orbital resonance and aforementioned parameters are well manifested by the formulae deduced in Appendix A. As for chirp signal, to illustrate the impacts of these parameters, and look for the "optimal parameters" that maximize the effect of resonance, we vary some of their values while keeping others the same as Tab. 1, and calculate the evolution of P numerically.For brevity, we will leave out the subscript "0" of X 0 in the rest of this appendix.

B.1 the celestial coordinate (ϑ, φ) of GW source
The role of inclination angle ι in GW amplitude is quite straightforward.To seek for the celestial position of source which leads to maximum resonance, we set ι = 0, vary (ϑ, φ) and calculate ∆P fin /P 0 numerically.The result is visualized in Fig. 6, where the normal vector of L2 orbit is marked as a red star.As is shown, maximum resonance occurs when the source binary and the test binary are face-on.
B.2 ω of the test binary and ϕ c of the source binary Eq.( 23) indicates that the effects of ω, ε and ϕ GW are degenerate, and Ṗsec depends on the combination ϕ GW − 2ω (since we have set the initial condition ε = 0).This relationship also holds for chirping signals, only that ϕ GW −2ω should be replaced by 2(ϕ c − ω), for ϕ c is defined as the source's orbital phase at coalescence.To prove this, we vary (ω, ϕ c ) and keep their difference invariant.Shown in the upper panel of Fig. 7 are the results of 3 equivalent combinations, and other parameters take the values in Tab. 1. Obviously, the results of these combinations are indistinguishable.Furthermore, we investigate the dependence of ∆P fin /P 0 on 2(ϕ c − ω), which is equivalent to varying ϕ c and keeping ω = 0.It can be seen from the lower panel of Fig. 7 that ∆P fin /P 0 acts like a sinusoidal function of 2ϕ c .In the following we will denote the values of {ϑ, φ, ι, ω, ϕ c } which maximize ∆P fin /P 0 as the Optimal Parameters (OP hereafter).

B.3 the semi-major axis a of test binary
We consider four values of the semi-major axis: a = {0.782,1.21, 2, 3, 4} × 10 4 km.The 1st and 2nd of them correspond to the configurations of L2 and LARES 1. Besides, to illustrate the impact of a on the same basis, for each value of a we iterate over {ϑ, φ, ι, ω, ϕ c }, and then set them to the OP.Given the total mass of the test binary (which is approximately the mass of Earth), the orbital frequency, and hence the resonance frequency, are totally determined by a. Consequently, for different values of a, resonance takes place at difference stages of GW (see Fig. 8): (1) Inspiral (e.g. a = 4 × 10 4 km): f GW increases slowly, thus resonance can last for a relatively long time (∼ 10 5 s).The resulting ∆P fin /P 0 is the largest among all the a values in consideration; (2) Merger (e.g. a = 1.21×10 4 km, 2×10 4 km, 3×10 4 km): The duration of resonance is shorter than case (1), while h A and ḧA of this stage get much larger, thus the effect of resonance is only slightly weaker than case (1); (3) Coalescence and ring down (e.g. a = 0.782×10 4 km): In this extreme situation, resonance can only last for a very short duration, leading to a much smaller ∆P fin /P 0 .
For the variation of a, when a = {0.782,1.21, 2, 3, 4} × 10 4 km, ∆a fin = {0.23,5.67, 8.03, 13.89, 20.91} × 10 −7 m.B.4 the component masses M bh , redshift z and mass ratio q of SMBHB For simplicity, we first consider equal-mass SMBHBs with component mass M bh .Still, for each value of M bh , we set the angular parameters to OP.For sources at cosmological distances, M bh enters the expression of GW in the form of redshifted mass M bh,z ≡ (1 + z)M bh .Theoretically, the influence of M bh,z is twofold.Firstly, for different M bh,z , f res appears at different stages of GW; secondly, the amplitude of GW is directly related to M bh,z .In addition, the GW amplitude is also inversely proportional to d L (z), which, at low redshifts, follows the Hubble law d L ∝ z.
The relationships between ∆P fin /P 0 and M bh at redshifts z = 0.01, 0.08105, 1 are shown in Fig. 9.At the redshift of TS (z = 0.08105), range (M 1 , M 2 ) marked in Fig. 9 includes the M bh values which would cause resonance to start at the merger stage, while if M bh > M 2 , resonance occurs during the inspiral stage.A peak of ∆P fin /P 0 can be found at M bh ≈ 5.5 × 10 7 M .The TS is reported to be an uneven mass-ratio system, thus there is necessity to examine the role of mass ratio q ≡ m 2 /m 1 .By varying q and keeping the chirp mass M c fixed, the resonant responses of L2 are plotted in Fig. 10.Results show that equal mass (q = 1) turns out to be the most optimistic case.q = 1 q = 1/2 q = 1/3 q = 1/4 q = 1/5 Fig. 10.The optimal resonant responses of L2 for different values of q with fixed Mc.

C The signal-to-noise ratio of distance residual
The one-sided PSD S n (f ) is defined as twice the Fourier transform of auto-correlation function R(τ ) = n(t)n(t + τ ) , n(t) being the noise at time t.Following [41], the range measurements are assumed to be unbiased, with uncorre-lated Gaussian noise of variance σ 2 .Thus, for discrete SLR data, by denoting the time interval between two adjacent measurements as t s , R(τ ) can be modeled as thus S n (f ) = 2σ 2 t s sinc(πt s f ).(26) Considering that the maximum frequency f max of δr(f ) is usually much smaller than 1/t s , we can approximate sinc(πt s f ) to 1, and it follows that where δr i ≡ δr(t i ), N obs is the total number of normal point measurements, and T obs = t s N obs .Note that to derive the second line, we have used the Parseval's theorem.For long-term data tracking, SNR is mainly contributed by the post-resonance stage.By assuming the ideal case that the only perturbation to the satellite orbit is from the incident GW of SMBHB, during the post-resonance stage, δr(t) is the difference between two stable Keplerian orbits.We first expand r(t) in terms of the Hansen coefficients to the linear order of e as r(t; a, e, ε) = a 1 − e cos M (t; a, ε) + O(e 2 ) . (28)

Fig. 1 .
Fig.1.The responses of 7 laser-ranged satellites, including LAGEOS 1/2, LARES 1/2, ETALON 1/2 and Ajisai, to the same GW signal (the lowest panel) emitted by the TS.In each panel, the vertical line represents the time when fGW = fres.Note that the resonance frequencies of ETALON 1/2 are much smaller than other satellites, therefore for these two satellites, we have labeled the times corresponding to fres, rather than showing them as vertical lines.Resonances take place consecutively among these SLR missions when the chirping signal sweep across the "frequency comb" in the sub-mHz frequency range.

Fig. 4 .
Fig. 4. The SNRs of the most optimistic responses to TS for different satellites.The thin curves are obtained by numerical calculation, while the thick ones represent the polynomial function Eq.(12).

Fig. 5 .
Fig.5.The relationship between the error of theeorbital period σP 0 and the in-orbit time T used for calculating P0.The blue curve and the red curve represents the results based on precision a and precision b, respectively.For comparison, the optimal response of L2 to TS is also plotted with grey dashed line in the same figure.
This work is supported by the National Key Research and Development Program of China, No. 2020YFC2200601 and No. 2021YFC2201901.

Fig. 7 .
Fig. 7. Upper panel: ∆P/P0 under 3 equivalent sets of (ω, ϕc).Lower panel: the dependence of ∆P fin /P0 on ϕc with ω fixed to 0. The first set in the upper panel is marked with dashed lines.
the comparison of optimal responses of different SLR missions.During resonance (e.g.t ∈ (0.2, 0.4) day for L2), δr(t) grows in time, and finally reaches a steady value δr fin = ∆a fin .Afterwards, in the post-resonance stage, the behavior of δr(t) is in consistence with our theoretical prediction (see Appendix C for the derivation) δr(t) would oscillate around ∆a fin with a linearly varying amplitude, and the rate of variation is proportional to e.The results of numerical calculation allow us to make an examination on the magnitudes of ∆a fin , ∆e fin and ∆ε fin .In the case of L2's optimal response, the 3πt/P sin M term dominates over other O(e) terms, thus