Constraining gravitational wave velocities using gravitational and electromagnetic wave observations of white dwarf binaries

Although the general theory of relativity (GR) predicts that gravitational waves (GWs) have exactly the same propagation velocity as electromagnetic (EM) waves, many theories of gravity beyond GR expect otherwise. Accurate measurement of the difference in their propagation speed, or a tight constraint on it, could be crucial to validate or put limits on theories beyond GR. The proposed future space-borne GW detectors are poised to detect a substantial number of Galactic white dwarf binaries (GWDBs), which emit the GW as semi-monochromatic signals. Concurrently, these GWDBs can also be identified as optical variable sources. Here we proposed that allocating a GWDB's optical light curve and contemporaneous GW signal can be used to trace the difference between the velocity of GW and EM waves. Simulating GW and EM wave data from 14 verification binaries (VBs), our method constrains propagation-originated phase differences, limiting the discrepancy between the speed of light ($c$) and GW ($c_{GW}$). Through the utilization of LISA's design sensitivity and the current precision in optical observation on GWDB, our study reveals that a four-year observation of the 14 recognized VBs results in a joint constraint that confines $\Delta c/c$ ($\Delta c = c_{\mathrm{GW}} - c$) to the range of $-2.1\times10^{-12}$ and $4.8\times10^{-12}$. Additionally, by incorporating this constraint on $c_{\mathrm{GW}}$, we are able to establish boundaries for the mass of the graviton, limiting it to $m_{\mathrm{g}}\le3\times10^{-23}\,e\mathrm{V}/c^{2}$, and for the parameter associated with local Lorentz violation, $\bar{s}_{00}$, constrained within the range of $-3.4\times10^{-11}\le\bar{s}_{00}\le1.5\times10^{-11}$.


INTRODUCTION
Following Einstein's proposal of the general theory of relativity (GR) in 1915, he theoretically forecasted the existence of gravitational waves (GW) (Einstein 1916(Einstein , 1918)).According to his prediction, the speed of GW was expected to be precisely equal to the speed of light (Maggiore 2007).Further exploration into the realm of gravity has prompted certain researchers to propose alternative theories, a few of which involve modifications to the speed of GW.For example, if the graviton possesses a non-vanishing mass, its velocity would be dispersed in a vacuum as (Fan et al. 2017): where   is the mass of graviton;  and  GW are the speed of light and GW respectively; and  = ℎ  is the energy of the graviton, where ℎ is the Planck constant and  is the frequency of GW.The predicted dispersion pattern of gravitational waves (GW) has served as a means to investigate the mass of gravitons, causing phase differences among diverse frequency components of GW, and thereby ★ E-mail: kumarankit@ihep.ac.cn † E-mail: sxyi@ihep.ac.cn causing alterations in its waveform , formally a 1 post-Newtonian order term (Will 1998;Scientific et al. 2016): where  GW = ℎ/( g ) is the graviton's Compton wavelength,  is the cosmological redshift and  is the cosmological distance as defined explicitly in Will (1998).By aligning the observed waveform with its theoretical projection (Cutler et al. 1993), we can establish confidence intervals for the graviton's wavelength or mass.Employing a similar methodology, the LIGO and Virgo collaboration utilized gravitational wave events GW150914, GW170104, and GW170817 to establish upper limits for the mass of gravity, which were set at 1.2×10 −22 V/ 2 , 7.7×10 −23 V/ 2 , 9.51×10 −22 V/ 2 (Scientific et al. 2016(Scientific et al. , 2018;;Abbott et al. 2019), respectively.Correspondingly, at 1 kHz, the upper limits on Δ/ from the above mentioned GW events stands at < 4.2 × 10 −22 , 1.7 × 10 −22 , and 2.64 × 10 −22 , where Δ =  GW − .This constraint is notably stringent, nevertheless presuming  GW,∞ = , wherein  GW,∞ represents the velocity of gravitational waves at an infinitely high frequency.Subsequent discussion will unveil that additional methods, independent of this particular assumption, are still required.Moreover, if alternative methods can impose constraints on  GW in a much lower frequency regime, it could, in turn, restrict the possible range of the graviton's mass.(Liu et al. 2020), [7] (Larson & Hiscock 2000), [8] (Cutler et al. 2003).

Method Result
Refs.A theoretical framework known as the Standard Model extension (SME) characterizes general violations of Lorentz and CPT invariance in both general relativity and the standard model at attainable energies (Kosteleckỳ & Mewes 2016;Mewes 2019).In the context of SME, the GW velocity will be corrected to:

GW dispersion
where  GW,± represents the velocity of the two polarizations of GW.
If | ì  | is non-zero, birefringence will occur.With  0 exerting notable dominance over the | ì  | term, as highlighted by Kosteleckỳ & Mewes (2016), one can safely opt to exclude consideration of | ì  | in their calculationss (Liu et al. 2020).Consequently, our attention is solely on the deviation value  0 in this context.This approximation allows us to expand the GW velocity employing the spherical harmonic function: where ℜ and ℑ represent the real and imaginary parts of .  is a spherical harmonic function, therefore  should be an integer satisfying −  ≤  ≤ .Within the SME framework, it is posited that  ≤ 2. In order to constrain the nine coefficients present in the above equation, a minimum of nine GW sources, each with their respective  GW limits and corresponding sky positions, are needed.
There are several other ways to put constraints on the  GW .One of the stringent constraints arises from the time difference between GW170817 and GRB170817A (Goyal et al. 2021;Li et al. 2016).On 2017 August 17, the GW event GW170817 was observed by LIGO, and the associated gamma-ray burst(GRB) GRB170817A was observed independently by the Fermi Gamma-ray Burst Monitor(GBM).There is a 1.74 ± 0.05s delay between the prompt emission of GRB170817A and the chirp time of GW170817 (Abbott et al. 2017b).The observed time lag between the GW170817's peak and the first photon of GRB 170817A was used to constrain on the fractional speed difference to −3 × 10 −15 ≤ Δ/ ≤ +7 × 10 −16 (Abbott et al. 2017a).The result depends on the assumption of the intrinsic emission time difference between the GW and GRB, which is a debating presumption.(Ciolfi & Siegel 2015;Rezzolla & Kumar 2015;Tsang et al. 2012) A previous simulation study (Hendriks et al. 2023) showed that the joint detection of GW and its GRB prompt-emission counterpart will be a rare event.Therefore, we would expect obtaining a sample of  GW from at least nine different sources with this method, to be challenging in a short time scale.
There are also other model/source-independent methods, such as using the time difference between two detectors at large distances to the same GW signal.For example, the LIGO at the Hanford and Livingston sites have a time difference of about 10 milliseconds.Cornish et al. (2017) used a Bayesian approach, combining data from the first three LIGO GW observations, to constrain  GW to −0.45 ≤ Δ/ ≤ +0.42.Subsequently, Liu et al. (2020) further narrowed the interval to −0.03 ≤ Δ/ ≤ +0.02 with more data and constrained the coefficient s00 within −0.2 < s00 < 0.07.
Another model-independent method yielding more stringent constraints involves utilizing GW and EM signals emitted by Galactic white dwarf binaries (GWDB) (Larson & Hiscock 2000).It's believed there are roughly 11.5 billion white dwarfs (WD) in our Galaxy (Napiwotzki 2009), and a significant proportion of these will come together to form binary systems that steadily emit GW.In general, the frequency of these emitted GWs are below 10 −3 Hz (Barke et al. 2014;Amaro-Seoane et al. 2023) and cannot be detected effectively by ground-based GW detectors like LIGO, Virgo, and KAGRA.However, within a certain range, there are tens of millions of GWDBs that could be detected by space-borne GW detectors like LISA (Kupfer et al. 2018), TianQin (Huang et al. 2020;Ren et al. 2023) or Taĳi (Ruan et al. 2020).There are several known GWDBs, whose binary parameters are well measured with optical observation, and they can be effectively detected by space-borne GW detectors.Those GWDBs are usually referred to as verification GWDBs or verification binaries (VBs).
In 2000, Larson & Hiscock (2000) suggested comparing the potential phase difference between the GW and optical signals of GWDB to constrain the graviton's mass.In their work, they used the ratio between the sampling time and the total integration time to estimate the error in phase measurements.In addition, they focused on the errors caused by refraction due to interstellar medium refraction in the EM wave propagation processes and the initial phase difference at the source between the GW and EM signals.Later, Cutler et al. (2003) specifically discussed the phase error in GW measurements on GWDB, correcting the error to: where  is a correction term from the initial phase difference from the GW and EM signals.Moreover, they estimated the constrain to the graviton's Compton wavelength as  g > 35 × 10 12 km.They used analytical error propagation estimation and did not include the intrinsic phase discrepancy between GW and EM counterpart signals.All methods for measuring Δ/ and the corresponding results (including the limits from real observation and the estimated level of limits from proposed methods) are summarized in Table 1.
Here in the present work, we use simulation to consider the phase error of EM and GW altogether with more sources and also include the intrinsic phase difference between EM and GW due to polarization angle.We simulate the signals of GW and EM emitted by the VBs, and then fit the phase difference between the received signals.(Kupfer et al. 2018), [3] (Brown et al. 2011), [4] (Hermes et al. 2012), [5] (Kilic et al. 2014), [6] (Brown et al. 2020), [7] (Burdge et al. 2019b), [8] (Brown et al. 2010), [9] (Geier et al. 2013), [10] (Kilic et al. 2011) we exclude the initial phase differences and obtain an estimation of the phase differences due to propagation, and thus constraints on Δ/.The VBs utilized in this paper, along with their binary parameters, are presented in Table 2.The adopted method will be described in detail in the next section.In Section 3, we will focus on the intrinsic phase difference between GW and EM.In Section 4, we will discuss the method of simulating GW data.The process of optical light curve signal simulation will be outlined in Section 5. Section 6 will demonstrate how we match the phases using a joint fit and unveil our findings.Sections 7 and 8 are dedicated to exploring the boundaries of physical quantities and discussing the improvements.

METHODOLOGY OVERVIEW
The time variation of the strain of the two polarization modes of GW for a non-evolving circular GWDB can be approximated by the quadrupole moment, which is derived by Maggiore ( 2007) to vary sinusoidally.The response in the detector is a linear combination of the two modes in a certain proportion.(Huang et al. 2020) This proportion, known as the antenna pattern, will be provided with its specific expression in the next section.Therefore, the response in the GW detector can be represented with: while the variance in the EM properties can be represented with: with  GW =   EM .For the WDs utilized in this paper, the EM properties that can be readily measured and exhibit a simple sinusoidal period include apparent magnitude 2 (a.m.) and radial velocity 3 (r.v.).(Burdge et al. 2019a;Kupfer et al. 2018;Brown et al. 2011;Hermes et al. 2012;Kilic et al. 2014;Brown et al. 2020;Burdge et al. 2019b;Brown et al. 2010;Geier et al. 2013;Kilic et al. 2011) It is worth noting that both the GW and a.m.data vary at twice the orbital frequency, while the r.v.varies at an equal frequency.Thus in the above equation,  = 1 stands for a.m. and  = 2 for r.v..For the present scenario, we are considering the GW and EM signals received at the solar barycentre, assuming that the time and phase shift due to the relative motion between the detectors and the solar barycentre have been perfectly corrected 4 .The phase difference is: with Δ int representing the intrinsic phase difference between EM wave and GW.In the present work, we set the phase of the EM wave as fiducial, and the main intrinsic phase offset is due to the polarization 2 The a.m. is defined by , where   is the observed flux using spectral filter  and  ,0 is the reference flux (zero-point) for that photometric filter 3 The r.v. is velocity projection in the radial direction.It is derived from the Doppler red/blue-shift in high quality spectroscopic observations, which need more observation resources and thus there are less available data points in Fig. 3. 4 The relative motion between the detectors (optical telescopes and LISA) and the solar barycentre will introduce annual sinosoidal structures into the signals (due to geometrical delay and Doppler effect).If such noises are not corrected, any phase difference between GW and optical light curve can be dominated by such noise, and can not be used for our purpose.Luckily, such effects can be corrected as a standard procedure as the ephemeris of the Earth, the geolocation of the telescopes, and the orbital parameters of the LISA spacecrafts are known with high accuracies (Cutler 1998).We, therefore, assume we can directly work with data with such effects corrected.
Table 3.The table of polarization phases and the main parameters with their uncertainties: The data in each column of the table are, in order, the value and uncertainty of right ascension, declination, orbital inclination, and the polarization phases.The uncertainties of RA and DEC are from GAIA (Provencal et al. 1997;Vallenari et al. 2023), and the RA, DEC, orbital inclination and the uncertainties of inclination is from (Kupfer et al. 2018).The uncertainty of right ascension and declination is presented in the unit "mas" which refers to milliarcseconds. of the GW, which will be discussed in detail in the subsequent section.

Source
Δ prop denotes the phase difference induced during the propagation, which reflects the potential difference between  and  GW : where  is the orbital period of the GWDB.Consequently, we can obtain the desired quantity i.e.Δ as: The following sections will illustrate our approach to restricting the phases separately, allowing us to then establish a restriction on Δ.

INTRINSIC PHASE DIFFERENCE
The intrinsic phase difference Δ int comprises several components, with the polarization phase being the most significant one.The GW strain recorded by the detector can be described as a linear combination of the two GW polarizations modulated by the detector's response: where is the polarization phase and  +,× are the amplitudes of the two GW polarization, which can be represented as: where  refers to the orbital inclination and A 0 = 2( M ) 3/5  4  (  GW ) 2/3 (Huang et al. 2020).Here M is the chip mass and  is the distance.The  +,× are defined as the antenna pattern functions and can be mathematically represented as: where   is the polarization angle; (  ,   ) are the latitude and longitude of the source in the detector's coordinate frame, which can be transformed from Right Ascension (RA) and Declination (DEC) as in Huang et al. (2020).Further, we also estimated the uncertainties of the polarization phase  0 due to uncertainties in sky position, inclination, and polarization angle.We apply the Monte Carlo algorithm to calculate the phase uncertainty.We input randomly generated RA, DEC, inclination, and polarization angle values from 1000 normal distributions into the above equations.This process calculates the phase, allowing us to derive the polarization phase values with distribution.The main parameters with their uncertainties are shown in Table 3, where the uncertainties of RA and DEC are from GAIA (Provencal et al. 1997;Vallenari et al. 2023), and other data is borrowed from Kupfer et al. (2018).Besides, we take polarization angle   = 30 • with the uncertainty Δ  = 10 −3 , which is supported from Huang et al. (2020).The polarization phase and its associated uncertainties are also presented in Table 3.
Additional intrinsic phase difference will be introduced between the GW and the a.m.signals, when the VB are not in tidal synchronization (Larson & Hiscock 2000;Cutler et al. 2003).In this case, the elongated axes of WDs can be misaligned with the line of masses.According to Burkart et al. (2013); Fuller & Lai (2012); Yu et al. ( 2020), detached VBs with orbital period less than 20 min can be considered tidally synchronized, where this additional intrinsic phase difference disappears.Therefore, when considering the a.m.light curves, we limit the sample to those detached VBs with period less than 20 min.It results in only one VB in the sample: J0533+0209.We like to note that, although there are theoretical arguments that this system can be in tidal synchronization, there lacks support from observation.We therefore present two version of results, one is with the assumption that J0533+0209 is in tidal synchronization and its a.m.light curve are used, the other one is when we do not make such assumption, and we use its r.v.data instead.

GW SIMULATION
In the following section, we will employ simulation techniques to explore the accuracy with which we can determine the phase of gravitational waves.Here, we aim to implement several simplifications: First, we assume that the eccentricity of VBs has a negligible effect on the GW phase estimation.Typically, the trajectory of a GWDB deviates from a perfectly circular orbit and instead follows an elliptical trajectory characterized by a specific eccentricity.In practical terms, determining the eccentricity involves considerable uncertainty in measurement.However, the VBs utilized in this context are estimated to possess small eccentricity, allowing a circular orbit to serve as a reliable approximation.In the next section, we will observe that the radial velocities of all sources vary sinusoidally, indicating that their orbits are nearly circular.Moreover, the gravitational waves (GW) emitted by elliptical orbits have evolved from their initial simple, single-frequency sine waveforms into a more complex composition, now characterized by a superposition of multiple harmonics, as discussed in (Maggiore 2007), thereby intensifying the intricacy of the model.This paper, however, focuses on a more simplified scenario involving circular orbits.Our aim here is to test the fundamental principles of this novel method.Integrating the effects of eccentricity into the application of this method while fitting the data with the model remains a straightforward task.Thus, considering the aim of model simplification, we've opted to streamline these gravitational wave database sets to circular orbits.
The second simplification involves assuming that the orbital frequency of VBs remains static without any evolution.In reality, as GWDBs continuously emit GWs, their frequencies naturally escalate.Therefore, it becomes essential to explore how these fluctuations in frequency might impact subsequent research outcomes.Typically, LISA missions operate for a span of 4-5 years.Our main focus lies in evaluating whether the gravitational wave frequencies will exhibit substantial alterations within this 5-year time frame.
We use evol in LEGWORK (Wagg et al. 2022) to calculate the SNR.Take source J0533+0209 as an example, its orbital period change in 5 years is 4.76 × 10 −4 s, the cyclical rate of change is 3.02 × 10 −12 s/s, which is close to the results measured by Burdge et al. (2019b), i.e., 3.94 × 10 −12 s/s.The ratio of the period variation to the orbital period is 3.86 × 10 −7 , which can be considered as a constant period during the mission time period.In the subsequent processing, we will phase-fold the simulated data.
We employ LEGWORK to simulate the GW signals.For a source, the total length of probes in one mission is in the range of 4-5 years, and we therefore simulate 1.2 × 10 8 s time-scales.Since the orbital period scale is around 10 3 s, to reduce the amount of data, we take the sampling frequency as 10 −1 Hz, i.e. 10 s to measure one data point, in order to clearly visualize the GW signal.
Then we add simulated GW noise onto the signals.We consider only the noise within a frequency window around the GW signal.The width of the frequency window  is the frequency resolution 1/, where  is the duration of the data.Here, we assume that the distribution of the noise is normal.Based on the derivation in the Appendix A, the standard deviation of the noise can be given as: The sensitivity curve, as presented by Robson et al. (2019) and showcased in Figure 1, delineates the correlation between noise strain and its respective frequencies.
In figure 2, we plotted simulated signal and noise from 5000s to 15000s.The noise is filtered in the frequency range from  GW −  to  GW +  , where  = 1/ is the frequency resolution of the data.In figure 2,  ∼ 10 −4 Hz, while for the data we used below,  ∼ 10 −7 Hz.The corresponding SNR for this data segment is 15.2 in 4 years.

SIMULATION OF OPTICAL SIGNALS
For subsequent calculations, we require the optical data to vary sinusoidally as the GW, which excludes a fraction of the eclipsing binaries (though our method can be generalized easily to eclipsing binaries).A study conducted by Kupfer et al. (2018) in 2018 assembled a collection of 16 VBs.After eliminating specific eclipsing binaries and incorporating newly detected GWDBs based on the criteria previously discussed, the dataset was ultimately refined to consist of 14 VBs, as outlined in Table 2.
We also provided the basic information of the optical data of the VBs as shown in Table 4.To compare optical data quantitatively, we use relative root mean square error (RRMSE) (Loève 2017) to measure the quality of optical data:  where  is RRMES,  EM is the optical signal amplitude, Δ EM is the uncertainty, and  EM is the amount of data.RRMES is the total relative uncertainties of the data, so it can be used to qualify the accuracy of the optical data.
Here, we directly take the optical data as a sinusoidal signal and superimpose on it the error of a normally distributed random function whose standard deviation is the uncertainty of the data.For the temporal distribution, we set the data points uniformly scattered randomly within the GW observation duration.

JOINT ANALYSIS AND RESULT
Due to the long time series of the simulation and the large time interval of the optical data, the direct fitting of the signal will produce greater errors.As discussed previously, both the GW data and the optical data have a sinusoidal waveform.Therefore, before formal processing, we first choose to phase fold the signal, i.e., use the obtained time series to invert the corresponding orbital phase value at that moment.The obtained data can demonstrate the measured quantities in one cycle, as shown in Figure 3.
Following this, we applied a curve fitting method to generate phase-shifted data along with the corresponding 90% confidence interval.The fitting function utilized for this purpose is detailed as follows: where ,  are the phase and corresponding measurements, respectively, and  is a constant introduced mainly due to the inconsistent periods of GW, r.v., and a.m.In particular, for GW and a.m. = 2 and for r.v. = 1. and  are the values to be fitted, and the fitted phase  is all we need.
Here, we utilize fitting in order to obtain the expected value of the phase  with variance  2 .Regarding the fitting process, the data distribution is ideally modeled by a T-distribution.However, considering that our dataset sizes are substantial, we approximate the following aspects using a normal distribution, as detailed in (Loève 2017).
For every source, we end up with at least two sets of data, i.e., the GW signal phase distribution  GW ∼  ( GW ,  2 GW ) and the optical signal phase distribution  EM ∼  ( EM ,  2 EM ).Therefore the distribution of the phase difference between the two should also be normal.Using Eq.10 we can derive the expressions for Δ  as: Table 5.The limits on Δ/ from simulated data: The second column is the relative velocity difference using the optical r.v.data and the next two columns are the upper and lower limits of the 90% confidence interval; columns 5-7 are the same as 2-4, but using a.m.observation instead of r.v. as EM data.

Source
The 90% confidence interval of the normal distribution is ( Δ/ − 1.645 Δ/ ,  Δ/ + 1.645 Δ/ ).The upper and lower limits of the confidence intervals have been stated in Table 5.
Based on the assumption that the noises of different VBs in different frequencies are independent (Abbott et al. 2023), the joint probability of Δ/ considering all these 14 sources is: where  is the number of independent constraints.Since there are 14 sources with both optical r.v. and a.m.data, the value of  should be taken as 14.
For the product of multiple normal distributions, we determine the total expectation and variance using (Loève 2017): Consequently, the final outcome is obtained as −2.1 × 10 −12 ≤ Δ/ ≤ 4.8 × 10 −12 (−2.3 × 10 −12 ≤ Δ/ ≤ 4.9 × 10 −12 ) 5 . 5The results in the parentheses are those if we do not assume tidal synchronization of J0533+0209 and we use its r.v.instead of a.m..We present the corresponding results in parentheses too in the next section.

PHYSICS IMPLICATION
Upon acquiring the constraints on  GW , we can establish limitations on the physical parameters as outlined in Section 1.The mass of the graviton can be obtained from Eq 1 as: The mass is related to the graviton energy  = ℎ  .For GWDBs, the graviton energy is much lower than the GW events detected by LIGO due to their lower frequencies, which enables us to obtain more precise constraints.The distribution of graviton masses is derived by incorporating the constraints from each source into Eq.20.The best constraints involves utilizing source J1630-4233 with optical data of r.v. is  g ≤ 3 × 10 −23 V/ 2 .
Next, we delve into the examination of the coefficient s  .If only the simplest case is considered, which means only s00 is a non-zero coefficient, it holds: Substituting the results of the joint analysis yields s00 to −3.4 × 10 −11 ≤ s00 ≤ 1.5 × 10 −11 (−3.5 × 10 −11 ≤ s00 ≤ 1.6 × 10 −11 ), which is better than the result calculated in Liu et al. (2020).For the general case, employing 14 values for substitution into equation 4 provides limitations for all s  coefficients.However, this particular calculation will be conducted in a distinct study.

DISCUSSION
The methodology outlined in this paper showcases the potential to confine Δ/ within the realm of 10 −12 , surpassing outcomes attained through GW dispersion.However, it falls short by approximately three to four orders of magnitude compared to the constraint utilizing the time lag between GW and GRB.In the subsequent discussion, we discuss different aspects that would strengthen this method.
As depicted in Table 4, the data sizes for r.v.typically tend to fall within the range of approximately 10 to 100, whereas the data sizes for a.m.typically scale up to hundreds or even surpass a thousand.Employing source J0533+0209 once more, we experimented by adjusting the data size for a.m.within the range of 100 to 10,000, while varying the data size for r.v.from 10 to 1,000.Employing the identical methodology as in the preceding section, each dataset was simulated 50 times to derive the distribution reflecting the ultimate standard deviation of the outcomes, as depicted in Figure 4. Similarly, we varied the uncertainties of r.v. and a.m.within two magnitudes and simulated each data set 50 times to obtain the distribution of the standard deviation of Δ/ as shown in Figure 5.As the volume of data grows and quality enhancements are implemented, discernible improvements are evident in the simulation results, characterized by a decreased standard deviation and enhanced stability.Based on the above discussion, it is apparent that enhancing the quality of optical data may not significantly optimize the final results.
Increasing the number of sources presents another viable method.With reference to the relationship described in equation 19, we can roughly estimate the association between the number of sources, denoted by  s , and the standard deviation of Δ/ as: As LISA runs, more VBs will be discovered.In addition, the eclipsing binaries are not considered in this paper in order to simplify the model.In fact, if the theoretical function of the a.m.variation of the eclipsing binaries is known which can be accomplished using the Phoebe package (Prša & Zwitter 2005), it can also be fitted to constrain the phase difference.Based on our assessment, it seems plausible to include data gathered from at least 8 sources in the dataset.
Overall, our method is astrophysical model independent which constrains  GW with better accuracy.In addition, the use of more sources allows a more precise constraint on the parameter of local Lorentz violation.When analyzing the graviton mass, utilizing gravitational waves with lower frequencies leads to improved constraints on the mass.
However, the final accuracy of this method is limited by the distance of the source.Compared to GRBs, whose distances are on the order of 10 Mpc, and VBs, whose distances are generally on the order of kpc, lead to final results that are 3-4 orders of magnitude larger than GRB.Also, a drawback of this one method is that the distance of the VBs cannot reach the Mpc magnitude, so the results obtained by optimizing the data are expected to be improved by only 1-2 orders of magnitude and cannot reach the accuracy obtained by GRB.
It is worth noting that   = 10 −3 rad is the best case, while the average value of   uncertainties is 0.1 rad.However, these uncertainties are estimated solely with GW data, while the sky position can be determined much better with optical observations, which is highly correlated with   in GW inference.Therefore, the value reported from Huang et al. (2020) should be considered as a conservative upper limit on   .
The method we proposed here can be easily transplanted into super black hole binaries (SBHB), where the distance can be in the order of 100Mpc, and thus results in constraint in Δ/ orders of several orders of magnitude better.SBHB are believed to exist in the center of merged Galaxies, with masses from 106 to 10 10 solar masses.When the binary separation is close to pc scale, and the orbital period is ∼years, such systems are expected to emit GW in the nano-Hz band, which enters the range of Pulsar Timing Arrays (PTA) (Jaffe & Backer 2003).

Figure 1 .
Figure 1.The sensitivity curve of LISA depicted in (Robson et al. 2019).

Figure 3 .
Figure 3. Phase-folded signals with error bar: Upper panel: GW data with 1.2 × 10 8 s sampling frequency of 10 −1 Hz; middle panel: 35 r.v.data uniformly scattered randomly within the GW observation duration; lower panel: 5000 a.m.data for the same distribution as r.v..For clearer visualization, the GW panel has been down-sampled to 0.1% of its original resolution.The r.v.measurement requires high quality spectrum and thus the data points are much fewer compared to that of a.m..

Figure 4 .
Figure 4. Diagram of the uncertainty of Δ/ after the change of optical data volume.

Figure 5 .
Figure 5. Diagram of the uncertainty of Δ/ after the change of optical data uncertainty.

Table 2 . The observed parameters of VBs: The
contents in each column of the table are, in order, the name of the source, right ascension, declination, orbital period, GWDB masses, distance, orbital inclination, 4-year SNR in LISA, and references.Refs.: [1](Burdge et al. 2019a), [2] .