Gravitational wave as probe of superfluid dark matter

In recent years, superfluid dark matter (SfDM) has become a competitive model of emergent modified Newtonian dynamics (MOND) scenario: MOND phenomenons naturally emerge as a derived concept due to an extra force mediated between baryons by phonons as a result of axionlike particles condensed as superfluid at galactic scales; Beyond galactic scales, these axionlike particles behave as normal fluid without phonon-mediated MOND-like force between baryons, therefore SfDM also maintains the usual success of $\Lambda$CDM at cosmological scales. In this paper, we use gravitational waves (GWs) to probe the relevant parameter space of SfDM. GWs through Bose-Einstein condensate (BEC) could propagate with a speed slightly deviation from the speed-of-light due to the change in the effective refractive index, which depends on the SfDM parameters and GW-source properties. We find that Five hundred meter Aperture Spherical Telescope (FAST), Square Kilometre Array (SKA) and International Pulsar Timing Array (IPTA) are the most promising means as GW probe of relevant parameter space of SfDM. Future space-based GW detectors are also capable of probing SfDM if a multimessenger approach is adopted.


I. INTRODUCTION
Despite the success of ΛCDM model at linear order and cosmological scales, there are two discordances that haunt the cosmologists and astronomers for decades: the galactic discordances and cosmic discordances. We will argue below that recently proposed superfluid dark matter (SfDM) scenario is capable of shedding some light on both galactic and cosmic discordances.
The galactic discordances lie in the semiempirical laws that govern galactic dynamics as results of either cold dark matter (CDM) or modified Newtonian dynamics (MOND) [1][2][3]. On the one hand, the very existence of DM is irrebuttable: 1. At cosmic scales, the observations from the big bang nucleosynthesis (BBN), the cosmic microwave background (CMB) and the large scale structure (LSS) all point to some form of nonbaryonic gravitating mass, and this is the scale where MOND failed miserably.
2. At cluster scales, the discovery of bullet cluster with offset mass distributions between the baryonic mass from optical and x-ray observations and the nonbaryonic mass from weak lensing provides almost direct proof [4] of the very existence of DM.
3. At galactic scales, the flatness of galaxy rotation curves clearly evince a mass discrepancy between baryonic matter and dynamical matter, of which the discrepancy is usually attributed to the socalled DM. * cairg@itp.ac.cn † liutongbo@itp.ac.cn ‡ schwang@itp.ac.cn However, the direct N-body simulations with the use of CDM encounter with some small-scale anomalies (see recent review [5] and references therein) like core-cusp, missing satellite and too-big-to-fail problems. Although all of these small-scale anomalies can be alleviated in the paradigm of Bose-Einstein condensate (BEC) DM [6] with (known as self-interacting DM) or without selfinteractions (known as fuzzy DM), on the other hand, the MOND still seems well established [7] due to a critical acceleration scale written in the data that is otherwise unnatural to be seen in scale-free CDM : 1. Globally, regardless of the specific distributions of baryonic mass along radial direction of galaxies, the asymptotic circular velocity is solely correlated with the total enclosed baryonic mass, which is known as the baryonic Tully-Fisher relation (BTFR) [8,9] derived exactly from MOND.
2. Locally, the observed distribution of baryonic mass predicts a radial acceleration that is strongly correlated with that traced by rotation curves, which is known as the mass discrepancy-acceleration relation (MDAR) [10][11][12][13] realized recently from MOND [14] as well.

Besides the BTFR and MDAR along with other
Kepler-like galactic laws [7] leading to the same critical acceleration scale, the galaxy rotation curve can be made universal [15] if one properly normalizes the radial distance, regardless as to whether the galaxy is of high surface brightness (HSB) or low surface brightness (LSB). This universal rotation curve (URC) [16] once again indicates that somehow the dynamics know intimately about the distribution of light, which will be too fine-tuning if DM is told to do the same thing.

arXiv:1710.02425v3 [hep-ph] 7 Feb 2018
It seems like we are in a dilemma [17] regarding to above galactic discordances, however, the SfDM provides us a hybrid way out by mimicking MOND phenomenons from axionlike particles condensed as superfluid at galactic scales, and at the same time maintaining the DM nature at cluster and cosmic scales. There are currently two kinds of models of SfDM that can produce MOND-like interaction with [18] or without [19] the help of the excited phonons from the condensed superfluid. See also [20] for an explanation of MOND critical acceleration scale by coupling SfDM to Dirac-Born-Infeld (DBI) dark energy.
Apart from the reconciliation of galactic discordances, SfDM might also be capable of alleviating the cosmic discordances. The cosmic discordances lie in the inconsistent measurements between CMB and LSS. On the one hand, the Hubble expansion rate inferred from CMB observation [21] is usually smaller than the local measurement from type-IA supernova [22]; On the other hand, the matter fluctuations extrapolated from CMB constraints [21] into late-time regime is larger than that expected from low-redshift LSS observations [23,24]. Recently in [25] the cosmic shear/bulk viscosity is shown clearly to be a natural and economic reconciliation of this CMB-LSS inconsistency between high redshift, large scale and low redshift, small scale. In SfDM scenario, the galaxy is within the superfluid phase with vanishing viscosity, however beyond galaxy cluster scales, those selfinteracting axionlike particles are in normal fluid phase with viscosity that can be made in principle to reconcile the CMB-LSS tensions and even the cosmic acceleration [26]. We will explore this possibility elsewhere in future work.
In this paper, we will adopt the recent proposal [27] using the velocity change of GWs to probe the parameter space of BEC DM, which will be briefly reviewed in Sec. II. The difference here is that the BEC DM is SfDM with MOND emerging at galactic scales. Both of SfDM models considered in [18,19] are estimated in Sec. III with further considerations of two-fluid phases [28] and baryon effect. In Sec. IV, the observational perspectives of different SfDM models are presented along with a discussion on Shapiro time delay between GWs and photons due to the effect of gravitational lensing. We summarize our result in Sec. V.

II. GRAVITATIONAL-WAVE PROBE
In [29], Sabin et al. have shown that spacetime distortions can produce phonons in BEC. Thus, we can apply it to the detection of GWs. The effective metric of the excitations on the flat spacetime metric is given by wheren is the mean number density of the background field, and the differential of mean pressureP with respect to the mean energy density ρ gives rise to the square of the mean speed of sound,c 2 s = dP /dρ. So, the solutions of the Klein-Gordon equation with this effective metric describe massless excitations propagating with the speed of soundc s . As a result, we can obtain the following dispersion relation where ω k is the frequency of the mode and the momentum of mode is denoted by k.
In order to obtain the change in the speed of GWs, we must calculate the refractive index of GWs when scattering off scatters inside the medium. Here we follow [27], and apply the optical theorem, which links the index of refraction, n g to the forward scattering amplitude, f (0) as wheren =ρ/m is the mean number density of scatterers inside the medium and k is the wave number of the incident wave. Since the exchange of energy comes along with the scattering of incident GWs, the forward scatteringnf (0) is then expressed in terms of the energy density of the GWs as well as that of the phononic excitations in the ground state. Therefore, the effective refractive index is given by where, ∆k is the change of the wave number of the incident GWs and ω GW = 2πf stands for the angular frequency. As noted in [27], this effect could be enhanced sizably due to the huge occupation number in the ground state and long-range correlations of the condensate. Therefore, it can be used for us to probe the SfDM with GWs. In the case of ordinary CDM, this effect is dramatically small and can be neglected. In the following context, we will explore the relation between the energy density of the GWs and that of the phononic excitations in the ground state as mentioned above.
To proceed, we consider the GWs produced at a distance D from Earth, which is outside the galaxy. The typical energy density we adopt is where h is the strain of GWs, and the Planck mass is related to Newton constant as M 2 Pl = 1/8πG N by convention. The propagation of the GWs through the DM halo will result in the relative change in its wave number, which can be calculated as where ∆ρ represents the exchange of the energy density between GWs and phonons, and will be replaced with the energy density required for the phononic excitations, which we shall discuss later. Next, we follow [29] and assume that the condensate is contained in a 1-dimensional cavity trap. The energy spectrum of the massless modes after imposing hard-wall boundary conditions, is then given by where l ∈ {1, 2, · · · } and the denominator D halo = 4R/π is the cavity length analog for GWs propagating through a spherical halo with radius R. In order to excite the massless modes inside the medium, the minimum energy density that we need is then given by the product of the number density of phonons and the energy difference between the closest modes ∆ρ ≡n∆ω =n π 2c where ∆ω ≡ ω l+1 − ω l . This ∆ρ will compensate the change in energy density of the GWs encountering a DM halo as we noted above. The average fraction of distance the GWs propagating through the halo with a reduced speed c g = 1/n g is given by The effective speed of GWs can then be defined as where ∆τ = xD/c g +(1−x)D represents the proper time that GWs take to propagate from the source location to the detector. Therefore, the change in the speed of the GWs due to the change of refractive index, compared with the speed of light in vacuum, which we adopt as c = 1, is given by where is the change of refractive index due to the propagation of the GW through the BEC medium. The above expression depends not only on the parameters of the GW, like the frequency f and characteristic strain h, but on the profile of SfDM as well, like the enclosed mass of the DM halo M , the mass of the axionlike particle m and the characteristic energy scale Λ, which are encoded in the expression ofn andc s . So we could apply (11) (12) to constrain the (m, Λ) parameter space of SfDM for some fixed parameters of the target GWs detector, e.g. f and h. The results are given in the next section.

III. SUPERFLUID DARK MATTER MODELS
In this section, we will first study the SfDM model discussed in [18] in Sec. III A, where the case of superfluid phase without baryons is studied in Sec. III A 1, and the case of two-fluid phases without baryons is studied in Sec. III A 2, and the case of superfluid phase including baryons is studied in Sec. III A 3. The second SfDM model [19] will be studied in Sec. III B, where only the case of superfluid phase without baryons is studied.

A. Model A
The general picture of model A [18] contains two parts: a MOND effective action describing SfDM phonons, and a coupling term mediating MOND force between baryons, Here in (13), P (X) is the pressure from the effectivefield-theory point-of-view in nonrelativistic regime at the lowest order in derivative, and X ≡θ − mΦ − (∇θ) 2 /2m describes superfluid phonons θ(t, r) = µt + φ(t, r) expanded at constant chemical potential µ with phonon excitations φ(t, r) in external gravitational potential Φ. The interaction term (14) is a minimal coupling between phonon θ and baryon density ρ b with coupling constant α. The model parameters m and Λ are the mass of SfDM particle and the characteristic energy scale, respectively. How could the Lagrangian L = L θ + L θb reproduce MOND at galactic scales ? Consider a static sphericallysymmetric approximation, θ(t, r) = µt + φ(r), X(r) = µ(r) − φ 2 (r)/2m,μ(r) ≡ µ − mΦ(r), the equation-ofmotion (EOM) of θ can be readily integrated as where is total enclosed baryon mass at radius r. It was shown in [18] that only the X < 0 branch admits a MONDian regime where κ(r) μ(r) with solution φ (r) = k(r). To see this, in this MON-Dian regime, the acceleration experienced by baryons from the phonon-mediated force can be made to match the MONDian acceleration a MOND = √ a 0 a N as long as α 3 = a 0 M Pl /Λ 2 , where a 0 = 1.2 × 10 −10 m/s 2 is the critical acceleration scale in MOND scenario. Remarkably α is of order unity for Λ ∼ meV, which together with m ∼ eV also gives rise to a DM halo with mass M ∼ 10 12 M of realistic size R ∼ 10 2 kpc as shown in [18].
As noted in [18], the effective action of form X 3/2 is specifically chosen to reproduce the MOND law at galactic scales. Condensate of this form behaves like superfluid with equation-of-state P ∼ ρ 3 , which under viral expansion P = κ B T ρ + g 2 (T )ρ 2 + g 3 (T )ρ 3 + · · · suggests that, the SfDM particles have negligible 2-body interactions and interact primarily through 3-body processes. This can be compared to the usual BEC DM with equation-ofstate P ∼ ρ 2 , which is governed by the two-body interactions. More strange forms of equation-of-state have been studied before in condense matter physics, like the unitary Fermi gas with effective action of form X 5 2 . Therefore the nonanalytic form of effective action X 3 2 of SfDM is not that strange from the effective-field-theory point of view. In fact, SfDM model can be constructed in [19] for arbitrary n with effective action X n , please see Sec. III B for an introductory discussion.
More comments on the condensation of SfDM. First of all, there is no explicit self-interaction term in original paper [18] of SfDM model, and the total effective action consists of an nonanalytic kinetic term X 3 2 and a coupling term θρ b between phonons and baryons. Therefore, the phrase "strong self-interaction" is referred to the quantum effect of Bose-Einstein condensation of axionlike particles. Second, the self-interaction is not that strong, just enough for axionlike particles thermalized at galactic scales. As you can see from Eqs.(11)- (14) in the original paper [18], the lower bound for the interaction cross section satisfies the current constraints on the cross section of self-interacting dark matter (SIDM). However, as pointed out in [18], SfDM is considerably different than SIDM, therefore each constraint must be carefully revisited. Third, it is not the strong self-interaction but the phonon-mediated attraction force between baryons that is responsible for the MOND law at galactic scales. The phonon-baryon coupling term itself has already softly broken the global U (1) symmetry explicitly only at the 1/M Pl level and is therefore technically natural. Finally, such phonon-baryon coupling term can arise from baryons coupling to the vortex sector of superfluid, which would give rise to a cos θρ b operator, thereby breaking the continuous shift symmetry down to a discrete subgroup. When expanded around the state at finite chemical potential φ = θ − µt, such operator would give the phononbaryon coupling term to leading order, albeit with an oscillatory prefactor. As pointed out in [18], such phononbaryon coupling term is treated as an empirical term in the effective action necessary to obtain the MOND phenomenon.
In the following three subsections, we will derive the SfDM profile under three different circumstances, whose superfluid halo radius will be extracted to estimate the mass density, number density and sound speed in (12).

Superfluid phase without baryons
First, the equation-of-state (EOS) for the superfluid phase, can be easily obtained from where the mass density ρ = mn is in nonrelativistic case, and the number density is calculated under grand canonical ensemble n = P (µ) = P (|X|). Second, the hydrostatic equilibrium equation along with Poisson equation gives rise to After replacing ρ = mn = mP (µ) = mP (|X|) on the left-hand side (LHS) and ρ = mn = mΛ(2m) 3 2 |X| 1 2 on the right-hand side (RHS), the equation above leads to following profile equation, which can be made dimensionless by normalizing [30] r = bξ; namely, Choosing one arrives at the Lane-Emden equation The Lane-Emden equation can be solved numerically upon given boundary conditions Ξ(0) = 1 and Ξ (0) = 0 [31]. Therefore one can find the value ξ 1 with vanishing profile Ξ(ξ 1 ) = 0, and hence the size of the SfDM halo, which is determined as R = bξ 1 . Third, instead of fixing b with the central DM mass density ρ 0 , we want to use the total enclosed mass M (R) of SfDM halo. To do this, rewriting (23) as and evaluating at DM halo radius R, one finds In this subsection, we always fix halo mass at a fiducial value M (R) = 10 12 M denoted simply as M . Now, we are ready to evaluate the change of effective refractive index in (12) δn g = 5.04 × 10 −29 m eV 44 5 Λ meV The velocity of the GWs changes correspondingly according to (11), which also depends on source distance D, frequency f , and strain h. Throughout the paper, we use following illustrative configurations [32] of different GW detectors to present different GW sources : The results are presented in the first line of Fig. 2, which will be summarized along with other models in section Sec. IV.

Two-fluid phases without baryons
There is an unsatisfactory in the calculations presented in subsection III A 1. At galactic scales the axionlike particles are condensed as superfluid, while beyond galactic scales, the axionlike particles behave like normal fluid. In the case of superfluid phase alone, the SfDM halo is enclosed at a radius R where SfDM mass density vanishes. However, if we consider both superfluid and normalfluid phases [28], the SfDM halo should be enclosed at a smaller radius R c with nonvanishing mass density, where mass densities and pressures of both phases are continuous at that radius, where index s stands for the superfluid phase, and the normal fluid denoted by n, whose profile is chosen as isothermal profile for concreteness and simplicity, other DM profile like NFW profile can also be used but with more free parameters encountered. The goal is to solve the matching equations (47) (50) for R c and ρ c . The first matching condition (44) is just and the second one can be qualified by equating the superfluid pressure with the hydrostatic equilibrium equation for the pressure of normal fluid, where the total enclosed mass M (r > R c ) is computed by with abbreviation M c ≡ M (R c ). Therefore the second matching condition [33] is To solve the matching equations (47) (50), one still needs to specify the total enclosed mass M c , which is determined similarly as (30) by or in dimensionless form, Hence M c can also be expressed as a function of R c by It is worth noting that, b(M, m, Λ) is still computed according to (31) as function of M, m, Λ, where M should be chosen properly so that the total enclosed SfDM halo mass M c = 10 12 M . In fact, after solving Eqs (47), (50) and (53) The results for the change of GW velocity are presented in the second line of Fig. 2, which will be summarized along with other models in Sec. IV.

Superfluid phase including baryons
There is another unsatisfactory in the calculations presented in subsection III A 1. The baryons also contribute the Poisson equation (22), thus influence the hydrostatic equilibrium equation (21) or (23) for the pressure of superfluid phase, therefore the SfDM profile equation (24), will be changed accordingly with addition of baryons, even we assume that baryons are subdominated in halo. It seems that (59) is difficult to solve without prior knowledge of baryon distribution ρ b (r). Fortunately, the coupling term (14) that gives rise to MONDian solution φ (r) = κ(r) directly connects SfDM with baryons, therefore the baryon distribution can be expressed as d dr Combining (59), (61) with dimensionless normalizations r = bξ and |X(r)| = X 0 Ξ(ξ), one finally arrives at the profile equation with baryon correction, Here the normalization constant b is given as before by or in terms of total enclosed mass by Nevertheless, the value for ξ 1 and b are different with those in Sec. III A 1 due to baryon correction, and they should be determined after solving (62). Solving (62) is new to our knowledge (see [34] for more details on baryon-phonon coupling), and it is also tricky because ρ 0 cannot be specified for given M, m, Λ without input value of b, which itself depends on the solution of (62) through ξ 1 . We propose here an iteration algorithm described below: 1. solving (62) without the baryon correction term, and obtaining the 0th iteration solution Ξ (0) (ξ), then locating the values ξ  | smaller than given small number, and SfDM halo radius is then We demonstrate this iteration algorithm in Fig.1 with M = 10 12 M , m = 0.6 eV, Λ = 0.2 meV. As we can see, after four iterations the profile of Ξ(ξ) stays fixed, and hence the solution of (62) is found. Hence, the value for ξ 1 and b are obtained as the final iteration values ξ The results of velocity change of GW through superfluid phase with baryons are presented in the third line of Fig.2, which will be summarized along with other models in Sec. IV.

B. Model B
Although the model A presented in [18] enjoys serval appealing features: 1. DM and MOND share a common origin as different phases of axionlike particles; 2. The MOND emerges without the need of additional degrees of freedom; 3. The phonon of BEC DM is fully appreciated for generating MOND law among baryons; 4. The idea of SfDM naturally distinguishes between galaxies and galaxy clusters, there are also some drawbacks, 1. The finite temperature corrections are required to cure the instability of the wrong-sign kinetic term from perturbation around the MONDian solution at zero-temperature; 2. The value of a 0 , α and Λ depends on temperature in such a way that their values at cosmic scales are four orders of magnitude deviated from those at galactic scales; 3. The form of kinetic term of superfluid is of nonanalytic nature.
Those motivate to propose another model for SfDM [19], which makes several differences as follows: 1. The phonon excitations are no longer responsible for MOND, and thus the EOS of BEC can be of general type; 2. The zero-temperature action is fully analytic in all field variables, and stable all by itself without finite temperature corrections; 3. The MOND is realized universally for both DM and baryons.
In [19], Khoury proposed a next-to-leading order (NLO) term containing higher-derivative operators in addition to the leading order (LO) superfluid action plus minimal coupling term, (69) (70) where n = 2 for concreteness as standard BEC, and X = µ−mΦ for the absence of phonon excitations and at finite chemical potential. The LO term expanded with X = µ− mΦ at leading order, together with the minimal coupling term, gives rise to a contribution of form −(ρ s + ρ b )Φ, where ρ s = Λ 4 µ m n−1 . Furthermore, the "symmetron" field χ lives in an effective potential with Z 2 symmetry χ → −χ spontaneously broken when the effective mass square flips a sign in the MONDian regime |∇Φ| < 3a 0 with vacuum expectation value (VEV) as Expanding the NLO term around above VEV admits where the second term is subdominated in the deep MOND regime |∇Φ| a 0 . Therefore in the deep MON-Dian and symmetry-breaking phase with VEV (74), the effective action to LO in gradients is of which the EOM is of MONDian form, As one can see that, the new model of SfDM is different from the one presented in III A. In III A, only baryons experience extra phonon-mediated MONDian force that is larger than Newtonian force at galactic scales. However, in this subsection, both baryons and DM particles feel MONDian force steamed from the NLO term in symmetry broken phase within deep MONDian regime. The SfDM profile was derived in [19] for the static, spherically-symmetric halo without baryons. The hydrostatic equilibrium equation for the pressure of SfDM is of MONDian form, since both DM and baryon are coupled to MOND gravity. However, unlike the model A in Sec. III A 3 with a direct connection (61) due to the phonon-mediated MOND force among baryons, the baryon contributions have to be ignored since there is no prior knowledge of baryon distribution ρ b (r), and one can only work out the SfDM profile in DM-only calculations [35]. Parallel to the calculations in Sec. III A 1, the normalized variables Ξ = X/X 0 and ξ = br are defined in such a way with that the SfDM profile equation is dimensionless, which can be readily solved numerically upon given boundary conditions Ξ(0) = 1 and Ξ (0) = 0. The X 0 in b can be similarly expressed in terms of the enclosed SfDM halo mass as in Sec. III A 1 by noting that which gives rise to when evaluated at halo radius R = bξ 1 with Ξ(ξ 1 ) = 0. Once we have the halo radius eV. The results for δc g with n = 2 are presented in the last line of Fig.2, which will be summarized along with other models in section Sec. IV.

IV. OBSERVATIONAL PERSPECTIVES
In this section, we summarize the observational perspective of constraining different SfDM models presented in the last section.
First, how well have we know for the possible deviation of GW velocity from the speed-of-light ? Reference [36] put a bound δc g < 2 × 10 −19 ∼ 2 × 10 −15 obtained from the absence of gravitational Cherenkov radiation for the observation of the highest energy cosmic rays. However, the direct observations of GW can also put bound on the GW velocity from three different ways: 1. The simplest approach is to measure the arrival time of the GW from compact binary system and the electromagnetic (EM) waves from EM counterparts of that compact binary system, if we understand well enough about the intrinsic time-lag between GW emission and photons emission, which is often taken to be zero as ad hoc estimation for the GW velocity, namely Here δt is the arrival time difference between GW and gamma-ray burst (GRB), and D = d L /(1 + z) 2 is the physical distance estimated from the luminosity distance d L and redshift z, where d L can be directly obtained from standard siren and z is obtained from multimessenger observations. Although lacking unambiguous evidences for the correlation between the Fermi-GRB event [37] and the GW150914 [38] event, such distant GW events can in principle constrain the change of GW velocity down to δc g ≤ 10 −40 ∼ 10 −17 level [39][40][41]. Recent observed time delay (+1.74 ± 0.05) between the GW 170817 event [42] and GRB 170817 event [43,44] has put a stringent bound −3 × 10 −15 < c g < +7 × 10 −16 on the GW velocity [45].
2. Without the EM counterparts for GW events, the phase changes of GW waveform alone [46,47] can bound the Compton wavelength of graviton from massive gravity (MG), namely or equivalently the graviton mass [48] through the definition of Compton wavelength λ g = h/(m g c).
Matched filtering of the GW waveforms from inspiralling compact binaries can in principle constrain a frequency-dependent GW velocity, which manifests an offset in the relative arrival times at a detector, since the GW emitted at low frequency early during inspiralling stage will travel slightly slower than those emitted at high frequency later. However, transforming the constraint on the graviton mass to the constraint on the GW velocity is nontrivial due to the modified dispersion relation ω(k) from massive gravity that satisfies or more complicated forms in other modified gravity [49]. Nevertheless, this way of constraining the GW velocity can never reach the precision that can be achieved easily from the joined measurements with EM counterparts.
3. The joined measurement with EW counterparts of GW events can only make a radical estimation on the GW velocity, whose improvement relies on the well understanding of GRB emission relative to GW emission. There are other ways that do not heavily rely on the multimessenger observations and set an absolute upper limit on GW propagation speed. For example, [50] gave a very loose bound 0.55c < c g < 1.42c with 90% confidence on the GW velocity, if one notices that the GWs are arriving at the two detectors of LIGO with a time difference [51]. However, with more GW events and more large worldwide network of detectors, the bound can be improved significantly. [50] hence forecasts that just five GWs events by the LIGO-Virgo-Kagra network will constrain the GW velocity within 1% precision. A second example is the strongly lensed GW events [52] that can be used to produce robust constraints on GW velocity at the 10 −7 level without assuming vanishing emission lag in the source and without knowledge of the sky position of the inspiral event. Another example is to fully appreciate the longtime observations of eclipsing binaries with periodic signals under the so called phase lag test with eclipsing binaries [53], where the phase lag between the GW and the EM signals accumulates such amount that the dwarf binary system WDS J0651+2844 can be used to constrain the GW velocity at the level of few parts in a trillion.
Second, what can be read from Fig.2 ? In Fig.2, we present the δc g with respect to the SfDM model parameters m and Λ all in logarithmic unit. Different GW sources with strain h and frequency f are specified in the panels by the typical configurations of GW detectors. DO not interpret the numbers of contours as the detection ability of GW detectors. They are the sensitivity numbers required for the GW detectors in order to probe that part of parameter spaces of SfDM models. The first three lines present the SfDM model A in section Sec. III A, where Sec. III A 1, Sec. III A 2, Sec. III A 3 are presented in order in the first, second and third lines, respectively. The last line presents the SfDM model B in Sec. III B.
The ground-based GW detectors, like LIGO and ET, have to reach a sensitivity of 10 −60 δc g 10 −40 to explore the relevant parameter spaces of SfDM model A, which is difficult even with the help of multimessenger astrophysics. The future space-borne GW detectors, like LISA and BBO, have to reach a sensitivity of 10 −40 δc g 10 −30 to explore the relevant parameter spaces of SfDM model A, which is promising with help of electromagnetic counterpart. The GW detectors that are sensitive around nHz can probe most of parameter space of interest of SfDM model A with sensitivity of 10 −35 δc g 10 −10 . Likewise, the same comments for SfDM model A also apply to model B, but with much more promising perspective. However, unlike the model A, decreasing the mass of axionlike particle makes it easier for GW probe in SfDM model B.
Third, we will briefly discuss the Shapiro time delay [54] when GWs and photons encounter the gravitational potential of DM along the line of sight. Considering GWs with relatively high frequency, which are relevant for ground-based detectors, the geometrical approximation holds and we can apply the standard formula where γ PPN is parametrized post-Newtonian (PPN) parameter, D is the distance to source and b is impact parameter. In this case, both GWs and photons share the same time delay, which is given by (93). It is, however, another story when the wavelengths of GWs are larger than the size of the lensing object, i.e.   λ GW GM/c 2 , which can be rewritten for the lensing mass as M 10 5 M (f /Hz) −1 [55,56]. In this case, the geometrical approximation breaks down and we have to take its wave optics into account. The additional time delay between GW and EM signals can reach ∼ 0.1s(f/Hz) −1 . Since the mass of DM halo we take is M ∼ 10 12 M , for GWs with frequencies relevant for LIGO, ET, LISA and BBO, the geometrical approximation remains valid and we do not have to consider the arrival time difference. As for those frequencies down to 10 −7 Hz, the additional time delay should be considered when multimessenger analysis is involved.

V. CONCLUSIONS
In this paper, we have studied the possibility that probing the relevant parameter space of SfDM with GWs. The results we obtained indicate that ground-based GW detectors, like LIGO and ET, are difficult to put constraints on the parameters for all models even with the help of multimessenger approach. As for space-borne GW detectors, like LISA and BBO, the ability to constrain the parameter space will be improved with the help of electromagnetic counterpart. The GW detectors sensitive around nHz, including IPTA, FAST, and SKA, are shown to be the most promising tools to probe most of parameter space of interest.
Two comments follow. First, how to distinguish the velocity changes of GWs through BEC from those due to massive graviton ? The velocity changes of GWs from massive graviton are universal independent of the GW sources and DM halos during propagation, however, one should otherwise observe different patterns of velocity changes of GWs through different DM halos from different GW sources at different sky locations. Second, a recent paper [57] claims to rule out the dark matter emulators scenarios, like Bekenstein's TeVeS theory [58] and Moffat's Scalar-Tensor-Vector gravity theory [59]. However, the SfDM models as MOND emulators scenarios are not ruled out yet. We hope our work will shed light on the test of SfDM scenario with GWs in future.